\ Slightly Better Prime Counting Program. It is derived from YAPCP
\ posted by Albert van der Horst
\ Usage (after loading):
\ 1000 pi . \ counts and prints the number of primes <=1000.
\ 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 ;
: sqrt ( u1 -- u2 )
0 d>f fsqrt f>d drop ;
\ 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 sqrt 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 ;
0 value extra-size \ additional max-table size
-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) ( n -- primes )
\ PRIMES is the number of primes <= n; works for 2R \ Use R@ as a local, representing IP.
R@ PRIME[] /
DUP R@ PRIME[] DUP * < IF
PI% R@ PRIME[] PI% - 2 + \ 2 represents P and P^2
ELSE
DUP
R@ 1 ?DO
OVER I DISMISS% -
LOOP NIP
THEN
RDROP ;
: PI ( n -- primes )
\ PRIMES is the number of primes <= n;
DUP MAX-TABLE < IF
PI(small)
ELSE
DUP 1- \ N PItobe
OVER PI(sqrt) 1+ 1 DO
OVER I DISMISS - 1+ \ Dismiss multiples, add 1 for prime itself.
LOOP NIP
THEN
;
' PI IS PI%
' DISMISS IS DISMISS%
: PI ( n -- primes )
dup sqrt over 0 d>f fln 0.65e f* fexp f>d drop ( min-elems opt-elems )
table-limit min max 9 max make-tables
dup 2 > if
PI
else
1- \ give same result as yapcp.fs for pi(0..2)
then
primes here - allot
pi-table free throw ;