\ Improved Slightly Better Prime Counting Program. It is derived from YAPCP
\ posted by Albert van der Horst
\ Last change(AH): added TABLE6 to short cut the dismissal of primes <=13.
\ Known bugs: this program doesn't like small numbers,
\ it gives a memory error.
\ Before TABLE6
\ 10 PI went wrong on dec alpha gforth 0.5.0, 100 PI and beyond was okay.
\ Now TABLE6
\ 10^4 PI goes wrong on dec alpha 10^5 ...10^10 works out right.
\ Usage (after loading):
\ 100000 pi . \ 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.
\ 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 + @ IF 2* dup BEGIN < WHILE
\ swap ! over REPEAT 2drop THEN LOOP , rot - +LOOP here 2dup cell+ /
\ >r r@ * ELSE 1- 1+ DO ' allot
\ from CORE-EXT :
\ ?DO Value nip
\ from BLOCK-EXT :
\ \
\ from EXCEPTION :
\ throw
\ from FILE :
\ S" (
\ from FLOAT :
\ d>f f>d
\ from FLOAT-EXT :
\ fsqrt
\ from LOCAL :
\ TO
\ from MEMORY :
\ allocate free
\ from non-ANS :
\ { -- } bounds Defer rdrop IS
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
25000000 value table-limit \ about 109MB with 64-bit cells
\ table-limit 2/ cells -> pi-table, + some for primes
: 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 )
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/ cells allocate throw { flags }
flags u 2/ cells -1 fill
u 0 dsqrt drop 2/ cells 0 ?do
flags i cells + @ if
i 2* 3 + dup i + begin ( prime index )
dup u 2/ < while
dup cells flags + 0 swap ! over +
repeat
2drop
then
loop
flags ;
-1 value max-table \ limit for table lookup
0 value primes \ lookup table
0 value pi-table \ table for small pi lookups
-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 cells + @ if
i 2* 3 + ,
then
loop ;
: sieve-to-pi ( addr u -- )
\ take a sieve with u elements and convert it to a pi-table
1 rot rot cells bounds ?do
i @ - dup i !
1 cells +loop
drop ;
: make-tables ( n -- )
\ set up primes, pi-table etc.
to max-table
here to primes
max-table sieve to pi-table
pi-table max-table 2/ 2dup sieve-to-primes \ a little too big
here primes - 0 cell+ / to PIhalf
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-table free throw ;
: test
100000 1 do i 0 pi i 7 .r 7 d.r cr loop ;