\ mostly from ( 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 ( 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 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