117 lines
3.5 KiB
Text
117 lines
3.5 KiB
Text
-- Naieve stuff, mostly, but coded with enthuiasm!
|
|
-- Following the idea behind (but not the code from!) the J submission:
|
|
-- Generate 10 primes (kept in p10) -- (print K=1)
|
|
-- Multiply each of them by the first ten primes
|
|
-- Sort and find unique values, take the first ten of those -- (print K=2)
|
|
-- Multiply each of them by the first ten primes
|
|
-- Sort and find unique values, take the first ten of those -- (print K=3)
|
|
-- ...
|
|
-- However I just keep a "top 10", using a bubble insertion, and stop
|
|
-- multiplying as soon as everything else for p10[i] will be too big.
|
|
|
|
-- (as calculated earlier from this routine,
|
|
-- or that "return 1" in pi() works just fine.)
|
|
--constant f17={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59}
|
|
constant f17={2,3,5,7,11,13,17}
|
|
|
|
function pi(integer n)
|
|
-- approximates the number of primes less than or equal to n
|
|
-- if n<=10 then return 4 end if
|
|
-- -- best estimate
|
|
-- return floor(n/(log(n)-1))
|
|
-- if n<=20 then return 1 end if -- (or use a table:)
|
|
if n<17 then
|
|
for i=1 to length(f17) do
|
|
if n<=f17[i] then return i end if
|
|
end for
|
|
end if
|
|
-- -- upper bound for n>=17 (Rosser and Schoenfeld 1962):
|
|
-- return floor(1.25506*n/log(n))
|
|
-- lower bound for n>=17 (Rosser and Schoenfeld 1962):
|
|
return floor(n/log(n))
|
|
end function
|
|
|
|
function primes(integer n)
|
|
-- return the first n prime numbers (tested 0 to 20,000, which took ~86s)
|
|
sequence prime
|
|
integer count = 0
|
|
integer lowN, highN, midN
|
|
|
|
-- First, iteratively estimate the sieve size required
|
|
lowN = 2*n
|
|
highN = n*n+1
|
|
while lowN<highN do
|
|
midN = floor((lowN+highN)/2)
|
|
if pi(midN)>n then
|
|
highN = midN
|
|
else
|
|
lowN = midN+1
|
|
end if
|
|
end while
|
|
-- Then apply standard sieve and store primes as we find
|
|
-- them towards the (no longer used) start of the sieve.
|
|
prime = repeat(1,highN)
|
|
for i=2 to highN do
|
|
if prime[i] then
|
|
count += 1
|
|
prime[count] = i
|
|
if count>=n then exit end if
|
|
for k=i+i to highN by i do
|
|
prime[k] = 0
|
|
end for
|
|
end if
|
|
end for
|
|
return prime[1..n]
|
|
end function
|
|
|
|
procedure display(integer k, sequence kprimes)
|
|
printf(1,"%d: ",k)
|
|
for i=1 to length(kprimes) do
|
|
printf(1,"%5d",kprimes[i])
|
|
end for
|
|
puts(1,"\n")
|
|
end procedure
|
|
|
|
function bubble(sequence next, integer v)
|
|
-- insert v into next (discarding next[$]), keeping next in ascending order
|
|
-- (relies on next[1] /always/ being smaller that anything that we insert.)
|
|
for i=length(next)-1 to 1 by -1 do
|
|
if v>next[i] then
|
|
next[i+1] = v
|
|
exit
|
|
end if
|
|
next[i+1] = next[i]
|
|
end for
|
|
return next
|
|
end function
|
|
|
|
procedure almost_prime()
|
|
sequence p10 = primes(10)
|
|
sequence apk = p10 -- (almostprime[k])
|
|
sequence next = repeat(0,length(p10))
|
|
integer high, test
|
|
for k=1 to 5 do
|
|
display(k,apk)
|
|
if k=5 then exit end if
|
|
next = apk
|
|
for i=1 to length(p10) do
|
|
-- next[i] = apk[i]*p10[1]
|
|
next[i] = apk[i]*2
|
|
end for
|
|
high = next[$]
|
|
for i=2 to length(p10) do
|
|
for j=1 to length(next) do
|
|
test = apk[j]*p10[i]
|
|
if not find(test,next) then
|
|
if test>high then exit end if
|
|
next = bubble(next,test)
|
|
high = next[$]
|
|
end if
|
|
end for
|
|
end for
|
|
apk = next
|
|
end for
|
|
if getc(0) then end if
|
|
end procedure
|
|
|
|
almost_prime()
|