\ Improved Slightly Better Prime Counting Program, double-cell version \ It is derived from YAPCP posted by Albert van der Horst \ \ Improvements by Albert van der Horst and Anton Ertl \ Usage (after loading): \ 100000. pi d. \ counts and prints the number of primes <=100000. \ Improvements: \ - Compiles on Gforth (but some Gforthisms were introduced) \ If you don't like the Gforthisms, try \ http://www.complang.tuwien.ac.at/forth/programs/sbpcp-1.0.fs \ (slower, but as non-standard as yapcp) \ - Is usable on 64-bit systems \ - It is faster, because it uses direct lookup tables instead of \ binary search. Speedup factor 2.5-3 for the ranges I tested. \ - added stack effect comments and other code cleanup. \ - TABLE6 to short cut the dismissal of primes <=13 (speedup factor 5) \ - table size for high speed up to TABLE-LIMIT (for limiting RAM consumption) \ - table compression (for factor 4-5 bigger tables with the same RAM) \ This program uses techniques from the world record holding program \ for counting primes by Deleglise e.a., but it is O(N) not O(N^2/3). \ See also : \ http://home.hccnet.nl/a.w.m.van.der.horst/benchpin.frt \ http://home.hccnet.nl/a.w.m.van.der.horst/pipi.pas (More explanation). \ The program uses the following words \ from CORE : \ environment? drop : ; 2/ cells fill i + c@ IF 2* dup BEGIN < WHILE \ swap c! over REPEAT 2drop THEN LOOP , and 0= rshift ! - here 2dup \ cell+ / 1+ @ mod rot * Constant DO = Create execute um/mod m* >r r@ \ um* ELSE 2over ' min r> max allot cr \ from CORE-EXT : \ Value ?DO nip :noname .r \ from BLOCK-EXT : \ \ \ from DOUBLE : \ d+ m*/ d< d- d.r \ from EXCEPTION : \ throw \ from FILE : \ S" ( \ from FLOAT : \ d>f f>d \ from FLOAT-EXT : \ fsqrt f** \ from LOCAL : \ TO \ from MEMORY : \ allocate free \ from non-ANS : \ { W: -- } Defer 2nip rdrop IS d> s" X:deferred" environment? drop \ ask for deferred words \ the table size will not grow beyond the following limit unless necessary \ set it as required for your RAM size 100000000 value table-limit \ about 93MB with 64-bit cells \ table-limit 2/ -> pi-lo-table, + some for prime : fsqrt ( r1 -- r2 ) \ workaround for gcc-2.95.x bug for gforth-fast (at least 0.6.2) 0 fsqrt drop ; : f** ( r1 r2 -- r ) \ workaround for gcc-2.95.x bug for gforth-fast (at least 0.6.2) 0 f** drop ; : dsqrt ( d1 -- d2 ) d>f fsqrt f>d ; \ this is adapted from the Byte Sieve : sieve { u -- addr } \ compute the primes up to U, in the form of an array of U/2 cells \ at ADDR, each containing -1 if the number is prime and 0 if it is not u 2/ allocate throw { flags } flags u 2/ 1 fill u 0 dsqrt drop 2/ 0 ?do flags i + c@ if i 2* 3 + dup i + begin ( prime index ) dup u 2/ < while dup flags + 0 swap c! over + repeat 2drop then loop flags ; -1 value max-table \ limit for table lookup 0 value primes \ lookup table 0 value pi-lo-table \ table for small pi lookups, offset since last hi entry 0 value pi-hi-table \ table for small pi lookups, entries for 1024|(N-3) -1 value PIhalf \ number of primes in lookup table : sieve-to-primes { addr u -- } \ take a sieve at addr with u elements, and generate a primes \ table HERE from it 1 , 2 , u 0 ?do addr i + c@ if i 2* 3 + , then loop ; \ our pi-table has a funny format (for more compact storage): \ for every odd number N there's a byte, containing pi(N)-pi((N-3)&~$3ff) \ for every N (where N mod 1024 = 3), theres a cell containing pi(N) \ so you get pi(N)=pi_hi((N-3)&~$3ff)+pi_lo(N) (see pi(small)). : sieve-to-pi { addr u -- } \ take a sieve with u elements and convert it to a pi-table 1 { lastbase } 1 u 0 ?do i addr + c@ + i $1ff and 0= if \ i is already divided by 2 dup to lastbase dup i 9 rshift cells pi-hi-table + ! then dup lastbase - i addr + c! loop drop ; : make-tables ( n -- ) \ set up primes, pi-table etc. to max-table here to primes max-table sieve to pi-lo-table pi-lo-table max-table 2/ 2dup sieve-to-primes \ a little too big here primes - 0 cell+ / to PIhalf max-table 10 rshift 1+ cells allocate throw to pi-hi-table sieve-to-pi ; \ For I return the ith prime NUMBER. : PRIME[] ( n -- nth-prime ) CELLS PRIMES + @ ; : pi(small) ( d -- dprimes ) \ PRIMES is the number of primes <= n; works for 2R \ Use R@ as a local, representing IP. 1 R@ PRIME[] m*/ 2DUP R@ PRIME[] DUP um* d< IF ( d ) PI% R@ PRIME[] 0 PI% d- 2. d+ \ 2 represents P and P^2 ELSE \ Initialise loop over primes, but cut short the smallest 6. \ Remember: highest index will be `'R@ 1-'' R@ 7 < IF 2DUP R@ 1 ELSE 2DUP phi(.,6) R@ 7 THEN ( d1 d2 n1 n2 ) ?DO 2OVER I DISMISS% d- LOOP 2NIP THEN RDROP ; : PI ( d -- dprimes ) \ PRIMES is the number of primes <= n; 2dup max-table 0 d< if pi(small) else 2DUP phi(.,6) 5. d+ \ N PItobe 2OVER PI(sqrt) drop 1+ 7 DO 2OVER I DISMISS d- 1. d+ \ Dismiss multiples, add 1 for prime itself. LOOP 2NIP THEN ; ' PI IS PI% ' DISMISS IS DISMISS% : PI ( d -- dprimes ) 2dup dsqrt drop >r 2dup d>f 0.6e f** f>d drop ( d dopt-elems r: dmin-elems) table-limit min r> max 12000 max make-tables 2dup 2. d> if PI else 1. - \ give same result as yapcp.fs for pi(0..2) then primes here - allot pi-lo-table free throw pi-hi-table free throw ; : test 100000 1 do i 0 pi i 7 .r 7 d.r cr loop ;