\ 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 : 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 ; \ With this sqrt(MAX-INT) we can handle all positive numbers. -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 ; \ : ?PRIME ( n -- f ?) \ DUP 2 DO \ DUP I /MOD \ I < IF 2DROP -1 LEAVE THEN \ 0= IF DROP 0 LEAVE THEN \ LOOP ; \ : make-tables ( 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 ; \ --------------------------------------------------- : BIN-SEARCH ( nmin nmax xt -- nres ) \ !! undocumented >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 I return the ith prime NUMBER. : PRIME[] ( n -- nth-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< ( n -- f ) PRIME[] LIMIT @ > 0= ; \ Auxiliary for BIN-SEARCH : PI(small)x ( n -- primes ) \ PRIMES is the number of primes <= n; works for 0 ) LIMIT ! 0 PIhalf ['] p< BIN-SEARCH 1- ; : pi(small) ( n -- primes ) \ PRIMES is the number of primes <= n; works for 2 0= ; \ Auxiliary for BIN-SEARCH : PI(sqrt) ( n -- primes ) \ PRIMES is the number of primes <= sqrt(n); works for 0R \ 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 ; : PI ( n -- primes ) \ PRIMES is the number of primes <= n; 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% : PI ( n -- primes ) dup sqrt make-tables 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