\ Article: 109427 of comp.lang.forth
\ Newsgroups: comp.lang.forth
\ Path: tunews.univie.ac.at!aconews-feed.univie.ac.at!newscore.univie.ac.at!newsgate.cistron.nl!transit.news.xs4all.nl!post.news.xs4all.nl!uucp.xs4all.nl!xs4all!spenarnc.xs4all.nl!not-for-mail
\ From: Albert van der Horst
\ Subject: YAPCP (iso version)
\ Date: 27 Mar 2005 19:04:01 GMT
\ Message-ID:
\ Lines: 175
\ Organization: Dutch Forth Workshop
\ Xref: tunews.univie.ac.at comp.lang.forth:109427
\ I present here Yet Another Prime Counting Program.
\ It is supposedly ISO, except for some words every Forth has,
\ see the REQUIRE 's.
\ This algorithm is fast, or I wouldn't bother.
\ It counts to its limit 2G in a reasonable time on my POS
\ Pentium 100 Mhz on a slow Forth (ciforth, 5 minutes).
\ I'm interested in
\ - confirmation that it runs on 16 bit Forths
\ - reports of timing with fast Forth's
\ - a double precision version
\ - any tweaks Marcel Hendrix may come up with.
\ Note that this is an example of mutually recursive routines,
\ hence the DEFER (where I would prefer FORWARD).
\ It adapts to 16,32,64 Forths, but you may not like the
\ 140 M CELLS table generated during loading on 64 bit systems.
\ (And only tested on a 32 bit Forth)
\ This is not just generating a table and looking up!
\ The table is to the square root of what can be handled.
\ In ciforth the table is burned into the executable, which is
\ reasonable for 16/32 bits. See main. Then you can say
\ count 1000000000
\ [ ONLY prevents that something interesting happens for
\ count '"rm -rf /" SYSTEM' ]
\ - - - - - - - 8 < -- - - - - - - 8 < -- - - - - - - 8 < -
( 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 ;
\ ---------------------------------------------------
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 ( n -- f ?)
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 value max-table \ limit for table lookup
0 value primes \ lookup table
-1 value PIhalf \ number of primes in lookup table
: prime-table ( n -- )
to max-table
here to primes
1 BEGIN DUP ?PRIME IF DUP , THEN 1+ DUP MAX-TABLE = UNTIL DROP
here primes - 0 cell+ / to PIhalf ;
\ For I return the ith prime NUMBER.
: PRIME[] ( n -- nprime )
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%
: fsqrt ( r1 -- r2 )
\ workaround for gcc-2.95.x bug for gforth-fast (at least 0.6.2)
0 fsqrt drop ;
: PI ( n -- primes )
dup 0 d>f fsqrt f>d drop prime-table
PI
primes here - allot ;
\ 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