\ 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