61 lines
1.5 KiB
Nim
61 lines
1.5 KiB
Nim
import std/[strformat, tables]
|
|
|
|
func isPrime(n: Natural): bool =
|
|
## Return "true" is "n" is prime.
|
|
if n < 2: return false
|
|
if (n and 1) == 0: return n == 2
|
|
if n mod 3 == 0: return n == 3
|
|
var d = 5
|
|
var step = 2
|
|
while d * d <= n:
|
|
if n mod d == 0:
|
|
return false
|
|
inc d, step
|
|
step = 6 - step
|
|
return true
|
|
|
|
const Inc = [4, 2, 4, 2, 4, 6, 2, 6]
|
|
|
|
func firstPrimeFactor(n: Positive): int =
|
|
## Return the first prime factor.
|
|
## Assuming "n" is odd.
|
|
if n == 1: return 1
|
|
if n mod 3 == 0: return 3
|
|
if n mod 5 == 0: return 5
|
|
var k = 7
|
|
var i = 0
|
|
while k * k <= n:
|
|
if n mod k == 0:
|
|
return k
|
|
k += Inc[i]
|
|
i = (i + 1) and 7
|
|
return n
|
|
|
|
|
|
var blum: array[50, int]
|
|
var bc = 0
|
|
var counts: CountTable[int]
|
|
var n = 1
|
|
while true:
|
|
var p = n.firstPrimeFactor
|
|
if (p and 3) == 3:
|
|
let q = n div p
|
|
if q != p and (q and 3) == 3 and q.isPrime:
|
|
if bc < 50: blum[bc] = n
|
|
counts.inc(n mod 10)
|
|
inc bc
|
|
if bc == 50:
|
|
echo "First 50 Blum integers:"
|
|
for i, val in blum:
|
|
stdout.write &"{val:3}"
|
|
stdout.write if i mod 10 == 9: '\n' else: ' '
|
|
echo()
|
|
elif bc == 26828 or bc mod 100000 == 0:
|
|
echo &"The {bc:>6}th Blum integer is: {n:>7}"
|
|
if bc == 400000:
|
|
echo "\n% distribution of the first 400_000 Blum integers:"
|
|
for i in [1, 3, 7, 9]:
|
|
echo &" {counts[i]/4000:6.3f} % end in {i}"
|
|
break
|
|
n = if n mod 5 == 3: n + 4 else: n + 2
|