134 lines
6.9 KiB
Text
134 lines
6.9 KiB
Text
with javascript_semantics
|
|
requires("1.0.2") -- (for in, tagstart)
|
|
function half(integer n) return floor((n-1)/2)+1 end function // convenience convert to idx
|
|
|
|
function count_primes(atom n)
|
|
// non-recursive Legendre prime counting function for a range `n`...
|
|
// has O(n^(3/4)/((log n)^2)) time complexity; O(n^(1/2)) space complexity.
|
|
if n<3 then return iff(n<2?0:1) end if // can't odd sieve for n less than 3!
|
|
integer sqrtn = trunc(sqrt(n)), // (actual limit)
|
|
mxndx = floor((sqrtn-1)/2) // odds-only limit
|
|
--
|
|
-- smalls is the current accumulated counts of odd primes 1 to sqrt(n), initialized
|
|
-- to odds-only sieving, ie {0,1,2,3,4...} meaning 0 odd primes to 1, 1 o.p to 3,...
|
|
--
|
|
-- roughs is the current odd k-rough numbers up to sqrt of range; k = 2
|
|
-- initialized to all odd positive numbers 1, 3, 5, 7, 9, 11, ... sqrt(n)
|
|
--
|
|
-- larges is an array of current phi counts for the above roughs... except they are
|
|
-- not strictly `phi`'s since they also include primes, to match `smalls` above!
|
|
-- initialized for current roughs after accounting for the even prime of two...
|
|
--
|
|
-- composite is a flag array representing odd numbers 1..sqrtn, for sieving.
|
|
-- initialized false, meaning all positive odd numbers are potentially prime
|
|
-- note that this array starts at (and keeps) 1 to match the algorithm even
|
|
-- though 1 is not actually a prime, as 1 is important in computation of phi...
|
|
--
|
|
sequence smalls = tagset(mxndx,0),
|
|
roughs = tagstart(1,mxndx+1,2),
|
|
larges = sq_floor_div(sq_sub(sq_div(n,roughs),1),2),
|
|
composite = repeat(false,mxndx+1)
|
|
|
|
integer bp = 3, // 'current' base prime
|
|
nbp = 0, // number of base primes found
|
|
mxri = mxndx, // current highest used rough index
|
|
i = 2, sqri = 4 // index and square (index-1) limit
|
|
// partial sieve loop, adjusting larges/smalls, compressing larges/roughs...
|
|
while sqri<=mxndx do // partial sieve to square index limit
|
|
if not composite[i] then
|
|
// cull from composite so they will never be found again
|
|
composite[i] = true // cull bp and multiples
|
|
for c=sqri+1 to mxndx+1 by bp do
|
|
composite[c] = true
|
|
end for
|
|
// partial sieving to current base prime is now completed!
|
|
|
|
// now adjust `larges` for latest partial sieve pass...
|
|
integer ori = 0 // compress input rough index(k) to output one
|
|
for k,q in roughs to mxri+1 do
|
|
// q is not necessarily prime but may be a product of primes not yet
|
|
// culled by partial sieving (saves ops cmprd to recursive Legendre)
|
|
// skip over values of `q` already culled in the last partial sieve:
|
|
integer qi = floor(q/2)+1; // index of always odd q!
|
|
if not composite[qi] then
|
|
// since `q` cannot be equal to bp due to cull of bp and above skip;
|
|
atom d = bp*q, // `d` is odd product of some combination of odd primes!
|
|
// the following computation is essential to the algorithm's speed,
|
|
// see the Nim entry for the full details of how this works
|
|
dadj = iff(d<=sqrtn ? larges[smalls[half(d)]-nbp+1]
|
|
: smalls[half(floor(n/d))])
|
|
ori += 1
|
|
larges[ori] = larges[k]-dadj+nbp // base primes count over subtracted!
|
|
// eliminate rough values that have been culled in partial sieve:
|
|
// note that `larges` and `roughs` indices relate to each other!
|
|
roughs[ori] = q
|
|
end if
|
|
end for
|
|
|
|
integer m = mxndx // and adjust `smalls` for latest partial sieve pass...
|
|
// this is faster than recounting over the `composite` array for each loop...
|
|
for k=(sqrtn/bp-1)||1 to bp by -2 do // k always odd!
|
|
// `c` is correction from current count to desired count...
|
|
// `e` is end limit index no correction is necessary for current cull...
|
|
integer c = smalls[half(k)]-nbp,
|
|
e = floor((k*bp)/2)
|
|
while m>=e do
|
|
smalls[m+1] -= c
|
|
m -= 1
|
|
end while
|
|
end for
|
|
|
|
nbp += 1 // increase number of found base primes
|
|
mxri = ori-1 // advance rough index for later
|
|
end if
|
|
bp += 2
|
|
sqri = (i+i)*(i+1)
|
|
i += 1
|
|
end while
|
|
|
|
// now `smalls` is a LUT of odd prime accumulated counts for all odd primes;
|
|
// `roughs` is exactly the "k-roughs" up to the sqrt of range with `k` (erm,
|
|
// mxri?) the index of the next prime above the quad root of the range;
|
|
// `larges` is the partial prime counts for each of the `roughs` values...
|
|
// note that `larges` values include the count of the odd base primes!!!
|
|
// - and `composite` is never used again!
|
|
|
|
// the following does the top-most "phi tree" calculation:
|
|
// the answer to here is all valid `phis`, combined here by subtraction,
|
|
// + compensate for included odd base prime counts over subracted above:
|
|
atom result = larges[1] - sum(larges[2..mxri+1])
|
|
+ trunc((mxri+1 + 2*(nbp-1))*mxri/2)
|
|
+ 1 // include the only even prime, ie 2
|
|
|
|
// This loop adds the counts due to the products of the `roughs` primes,
|
|
// of which we only use two different ones at a time, as all the
|
|
// combinations with lower primes than the cube root of the range have
|
|
// already been computed and included with the previous major loop...
|
|
// see text description in the Nim entry for how this works...
|
|
for ri,p in roughs from 2 do // for all `roughs` (now prime) bar '1':
|
|
atom m = trunc(n/p), // `m` is the `p` quotient
|
|
// so that the end limit `e` can be calculated based on `n`/(`p`^2)
|
|
e = smalls[half(floor(m/p))+1]-nbp+1
|
|
// the following test is equivalent to non-splitting optmization:
|
|
if e<=ri then exit end if // quit when no more pairs! - aka stop
|
|
// at about `p` of cube root of range!
|
|
for k=ri+1 to e do // for all `roughs` greater than `p` to limit:
|
|
result += smalls[half(floor(m/roughs[k]))];
|
|
end for
|
|
// compensate for all the extra base prime counts just added!
|
|
result -= (e-ri)*(nbp+ri-2)
|
|
end for
|
|
return result
|
|
end function
|
|
|
|
atom t = time()
|
|
constant expected = {0,4,25,168,1229,9592,78498,664579,5761455,
|
|
50847534,455052511,4118054813,37607912018,
|
|
346065536839,3204941750802}
|
|
for i=0 to iff(platform()=JS?11:14) do -- (sp: keep js under 2s)
|
|
atom c = count_primes(power(10,i))
|
|
assert(c==expected[i+1])
|
|
string e = elapsed(time()-t,0.1," (%s)")
|
|
printf(1,"10^%d = %d%s\n",{i,c,e})
|
|
end for
|
|
printf(1,"\nTook %s\n",elapsed(time()-t))
|