53 lines
1.5 KiB
Text
53 lines
1.5 KiB
Text
local int = require "int"
|
|
local fmt = require "fmt"
|
|
|
|
local inc = {4, 2, 4, 2, 4, 6, 2, 6}
|
|
|
|
-- Assumes n is odd.
|
|
local function first_prime_factor(n)
|
|
if n == 1 then return 1 end
|
|
if n % 3 == 0 then return 3 end
|
|
if n % 5 == 0 then return 5 end
|
|
local k = 7
|
|
local i = 1
|
|
while k * k <= n do
|
|
if n % k == 0 then
|
|
return k
|
|
else
|
|
k = k + inc[i]
|
|
i = i % 8 + 1
|
|
end
|
|
end
|
|
return n
|
|
end
|
|
|
|
local blum = {}
|
|
local bc = 0
|
|
local counts = { [1] = 0, [3] = 0, [7] = 0, [9] = 0 }
|
|
local i = 1
|
|
while true do
|
|
local p = first_prime_factor(i)
|
|
if p % 4 == 3 then
|
|
local q = i // p
|
|
if q != p and q % 4 == 3 and int.isprime(q) then
|
|
if bc < 50 then blum:insert(i) end
|
|
counts[i % 10] += 1
|
|
bc += 1
|
|
if bc == 50 then
|
|
print("First 50 Blum integers:")
|
|
fmt.tprint("%3d ", blum, 10)
|
|
print()
|
|
elseif bc == 26828 or bc % 100_000 == 0 then
|
|
fmt.print("The %9s Blum integer is: %9s", fmt.ord(bc, true), fmt.int(i))
|
|
if bc == 400_000 then
|
|
print("\n% distribution of the first 400,000 Blum integers:")
|
|
for {1, 3, 7, 9} as j do
|
|
fmt.print(" %6.3f%% end in %d", counts[j] / 4000, j)
|
|
end
|
|
return
|
|
end
|
|
end
|
|
end
|
|
end
|
|
i = (i % 5 == 3) ? i + 4 : i + 2
|
|
end
|