79 lines
2.2 KiB
Plaintext
79 lines
2.2 KiB
Plaintext
#include "isprime.bas"
|
|
uses console
|
|
|
|
dim inc(7) as integer
|
|
inc(0) = 4: inc(1) = 2: inc(2) = 4
|
|
inc(3) = 2: inc(4) = 4: inc(5) = 6
|
|
inc(6) = 2: inc(7) = 6
|
|
|
|
function firstPrimeFactor(n as long) as long
|
|
if n = 1 then return 1
|
|
if n mod 3 = 0 then return 3
|
|
if n mod 5 = 0 then return 5
|
|
|
|
long k = 7
|
|
int idx = 0
|
|
|
|
while k * k <= n
|
|
if mod(n, k) = 0 then return k
|
|
k += inc(idx)
|
|
idx = mod((idx + 1), 8)
|
|
wend
|
|
return n
|
|
end function
|
|
|
|
sub main()
|
|
dim as long blum(49), counts(9)
|
|
long bc = 0, i = 1, pct, p, q
|
|
int j
|
|
string s, t
|
|
|
|
do
|
|
p = firstPrimeFactor(i)
|
|
if p mod 4 = 3 then
|
|
q = i \ p
|
|
if q <> p and q mod 4 = 3 and isPrime(q) then
|
|
if bc < 50 then blum(bc) = i
|
|
counts(i mod 10) = counts(i mod 10) + 1
|
|
bc += 1
|
|
|
|
if bc = 50 then
|
|
printl "First 50 Blum integers:"
|
|
for j = 0 to 49
|
|
s = str(blum(j))
|
|
while len(s) < 4: s = " " + s: wend
|
|
print s;
|
|
if mod((j + 1), 10) = 0 then printl
|
|
next
|
|
printl
|
|
elseif bc = 26828 or bc mod 100000 = 0 then
|
|
s = str(bc)
|
|
while len(s) < 6: s = " " + s: wend
|
|
t = str(i)
|
|
while len(t) < 7: t = " " + t: wend
|
|
printl "The " + s + "th Blum integer is: " + t
|
|
if bc = 400000 then
|
|
printl cr "% distribution of the first 400,000 Blum integers:"
|
|
|
|
for j = 1 to 9 step 2
|
|
if j <> 5 then
|
|
pct = counts(j)/4000
|
|
s = str(pct)
|
|
while len(s) < 5: s = " " + s: wend
|
|
printl s + "% end in " + str(j)
|
|
end if
|
|
next
|
|
end
|
|
end if
|
|
end if
|
|
end if
|
|
end if
|
|
if mod(i, 5) = 3 then i += 4 else i += 2
|
|
end do
|
|
end sub
|
|
|
|
main()
|
|
|
|
printl cr "Enter ..."
|
|
waitkey
|