\ Improved Slightly Better Prime Counting Program.  It is derived from YAPCP
\ posted by Albert van der Horst <albert@spenarnc.xs4all.nl>

\ 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 ;

: 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 ;

-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 2<n<=max-table
    3 - 2/ cells pi-table + @ ;

: pi(sqrt) ( n -- primes )
    \ PRIMES is the number of primes <= sqrt(n); works for 3<=sqrt(n)<=max-table
    sqrt pi(small) ;


: GCD ( m n -- gcd )
    \ For M and N return their GREATEST COMMON DIVIDER
    BEGIN 2DUP MOD DUP WHILE ROT DROP REPEAT DROP NIP ;

\ Tabulate the result of sieving the first 6 primes
2 3 * 5 * 7 * 11 * 13 * CONSTANT P6
:NONAME ( -- n )
    0 P6 1 DO P6 I GCD 1 = IF 1+ THEN DUP , LOOP ;
CREATE TABLE6   0 ,   EXECUTE
CONSTANT MULT6

\ A 6-prime is not divisible by any of the first 6 primes:
\  so 1 is counted as a 6-prime.
\  and the 6 primes themselves are not counted.
\ This results in a correction of 5 + for use as a prime counter.
: phi(.,6) ( n -- number )
    \ For N return the NUMBER of 6-primes less or equal.
    P6 /MOD   MULT6 *   SWAP   CELLS TABLE6 + @   + ;

\ Forward declarations
DEFER PI%      ( n -- primes )
DEFER DISMISS% ( N1 IP -- N2 )

: DISMISS ( N1 IP -- N2 )
    \ N2 are the amount of numbers <= N1 that are dismissed by the prime IP,
    \ i.e. it is divisible by PRIMES[IP] but not by a smaller prime.
    \ Compared to Deleglise : DISMISS(N,Psubn)= phi(N,Psubn) - phi(N,Psubn-1)
    \ Requires P<=N1
    >R          \ 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
        \ Initialise loop over primes, but cut short the smallest 6.
        \ Remember: highest index will be `'R@ 1-''
        R@ 7 < IF
            DUP   R@   1
        ELSE
           DUP phi(.,6) R@   7
        THEN ?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 phi(.,6) 5 +   \ N PItobe
        OVER PI(sqrt) 1+   7   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 0.6e f** f>d drop ( min-elems opt-elems )
    table-limit min max 12000 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 ;

: test
    100000 1 do i pi i 7 .r 7 .r cr loop ;