\ mostly from <ie0yap.dfr@spenarnc.xs4all.nl>
( Copyright{2005}: Albert van der Horst, HCC FIG Holland by GNU Public License)
( $Id: count.frt,v 1.14.1.1 2005/03/27 17:41:11 albert Exp $ )

\ This program is an improved version of my prime counting benchmark,
\ using 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 number of primes less or equal than n is a function denoted
\ by the greek letter pi, hence the names.

: \D
        POSTPONE \                       \ Comment out if debugging
; IMMEDIATE

\ You may comment this out, consider it a shopping list.
\  REQUIRE DEFER
\  REQUIRE RDROP
\  REQUIRE :NONAME
\  REQUIRE NIP

\ 74 LOAD ( Just kidding).
\ ---------------------------------------------------
( BIN-SEARCH binary_search_by stack ) \ AvdH
( nmin nmax xt -- nres )
\ See description in next screen.
\

: BIN-SEARCH    >R
    BEGIN       \ Loop variant IMAX - IMIN
        2DUP  <> WHILE
        2DUP + 2/  ( -- ihalf )
        DUP R@ EXECUTE IF
           1+  SWAP ROT DROP \ Replace IMIN
        ELSE
           SWAP DROP \ Replace IMAX
        THEN
    REPEAT
DROP RDROP ;

\ from  <ieafyt.nnw@spenarnc.xs4all.nl>
( binary_search_description )
\  ( BIN-SEARCH    : n IMIN, n IMAX, xt COMP -- n IRES )
\  Uses a comparison routine with execution token `COMP'
\  `COMP' must have the stack diagram ( IT -- flag) , where flag
\  typically means that IT compares lower or equal to some fixed
\  value. It may be  TRUE , FALSE or undefined for `IMIN' , but
\  it must be monotonic down in the range [IMIN,IMAX), i.e.
\  if IMIN<=IX<=IY<IMAX then if IX COMP gives false, IY COMP
\  cannot give true.

\  BIN-SEARCH finds the first index `IT' between `IMIN'
\  and `IMAX' (exclusive) for which `COMP' returns false.
\  or else ``IMAX''.
\  An empty range is possible, (`IMIN' and `IMAX' are equal.)
\  See also  binary_search_test in the examples section.

\ ---------------------------------------------------

8 CELLS CONSTANT MAX-BITS            \ # bits length of primes we have in table.

\ For N: "it IS prime".
\ Cases 0 1 are in fact illegal but return true.
: ?PRIME    DUP 2 DO
                DUP I /MOD
                I <  IF 2DROP -1 LEAVE THEN
                0= IF DROP  0 LEAVE THEN
            LOOP ;

\ With this sqrt(MAX-INT) we can handle all positive numbers.
1 MAX-BITS 2/ 1- LSHIFT 142 100 */ CONSTANT MAX-TABLE

\ You may tweek MAX-TABLE by hand, no harm done.
\D 1000 CONSTANT MAX-TABLE       \ Trim table during testing, or for 64 bit systems.

\ Makes 2 the first prime, and 1 the zero-th ...
\ For 32 bit Forth's we need (7/5)*0x8000/15.5*.7=4K400 cells (or 16 bit words).
\ For 64 bits we need (7/5)*0x8000,0000/31.5*.7 140M cells (or 32 bit words).
:NONAME 1 BEGIN DUP ?PRIME IF DUP , THEN 1+ DUP MAX-TABLE = UNTIL DROP ;
CREATE PRIMES EXECUTE

HERE PRIMES - 0 CELL+ / CONSTANT PIhalf   \ Number of primes in table.

\ For I return the ith prime NUMBER.
: PRIME[]   CELLS PRIMES + @ ;

\ For N return the NUMBER of primes less than or equal to it.
\ Use this only for small numbers < ``MAX-TABLE''
    VARIABLE LIMIT              \ Auxiliary for BIN-SEARCH
    : p< PRIME[] LIMIT @ > 0= ;    \ Auxiliary for BIN-SEARCH
: PI(small)   LIMIT !   0 PIhalf ['] p< BIN-SEARCH 1- ;

\ For N return the NUMBER of primes less than or equal to its square root.
    VARIABLE LIMIT              \ Auxiliary for BIN-SEARCH
    : p< PRIME[] DUP * LIMIT @ > 0= ;    \ Auxiliary for BIN-SEARCH
: PI(sqrt)   LIMIT !   0 PIhalf ['] p< BIN-SEARCH 1- ;

DEFER PI%              \ Forward declaration
DEFER 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
: DISMISS   >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
       DUP
       R@ 1 ?DO
             OVER I DISMISS% -
       LOOP NIP
   THEN
RDROP ;

\ For N return the NUMBER of primes less or equal.
: PI
    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%

\ Redefine with error detection (outside of recursion).
\ Only needed if you trimmed max-table!
\ : PI   DUP MAX-TABLE DUP * > IF ." Too large!  " 13 ERROR THEN   PI ;

\ ---  TEST ---

VARIABLE OLD
\ Find any errors in PI for arguments < N .
: FINDPROBLEMS        ( N -- )
   1 OLD !
   3 DO
      I PI DUP OLD @ - 0= 0=
      I ?PRIME 0= 0= <> IF
         DROP ." Wrong : " I . LEAVE
\D    ELSE
\D       ." Okay  : " I .        \ Comment out as desired.
      THEN
      OLD !
\D    .S
   LOOP
;

\ : MAIN   POSTPONE ONLY    1 ARG[] EVALUATE   PI . CR ;
\ - - - - - - -  8 < -- - - - - - -  8 < -- - - - - - -  8 < -

\  Groetjes Albert

\  --
\  -- 
\  Albert van der Horst,Oranjestr 8,3511 RA UTRECHT,THE NETHERLANDS
\  Economic growth -- like all pyramid schemes -- ultimately falters.
\  albert@spenarnc.xs4all.nl http://home.hccnet.nl/a.w.m.van.der.horst