Monday, November 23, 2015

Small Prime Sieve in FORTH

An explanation of the algorithm here can be found in a blog entry on the Sieve of Eratosthenes in my math and English blog, and in a previous entry in this blog on primes in C, along with my reasons for doing this, such as I have.

FORTH requires a lot of explanation, so I'm going to show you a C-like version first.

Comments in FORTH are between parentheses. Or you can start a comment with a backslash and it will end at the end of the line.

So, just read the code naturally and read the explanation in the comments. It's surprisingly readable (because I have written it to be readable). It's also full of hidden surprises for those used to algol-like in-fix languages including C and Ada.

Here's the C-like source:



( Archetypical implementation of the sieve of eratosthenes in FORTH -- gforth --
\ Copyright Joel Rees, 2015
\ By Joel Rees, Amagasaki, Japan, 2015
\ All rights reserved.
\ Permission granted by the author to use this code
\ for any purpose, 
\ on condition that substantial use 
\ shall retain this copyright and permission notice.
\
\ Timing results with MAXSIEVE set to 4194304 and output commented out:
\ me@fun:~/work/mathgames/sieve$ time gforth forthsieve.fs -e bye

\ real 0m1.061s
\ user 0m0.988s
\ sys 0m0.012s
\
\ Comparing to the C:
\
\ me@fun:~/work/mathgames/sieve$ time ./sieve_c
\ 
\ real    0m0.457s
\ user    0m0.436s
\ sys    0m0.020s
\
\ A bit more than double the run-times, 
\ so, really, in the same ballpark,
\ if not quite in the same league for speed.
)


256 constant MAXSIEVE
MAXSIEVE 1- 2 / constant FINALPASS

5 constant DISPWIDTH ( enough digits to display MAXSIEVE )

create sieve MAXSIEVE allot


: sieveMain ( -- )
0 sieve c!      ( 0 is not prime. )
0 sieve 1+ c!   ( 1 is not prime. )
sieve MAXSIEVE 2 ?do    ( set flags to true for 2 to FINALPASS. )
   -1 over i + c! loop  ( sieve pointer remains on stack. ) 
FINALPASS 2 ?do ( clear flags at all multiples. )
   dup i + c@ if      ( don't bother if not prime. )
      MAXSIEVE i dup + ?do    ( start at first multiple -- double. )
         0 over i + c!      ( clear at this multiple. )
      j +loop           ( sieve still on stack. )
   then
loop            ( sieve still on stack. )
MAXSIEVE 0 ?do
   i DISPWIDTH .r ." : is "
   dup i + c@ 0= if 
      ." not " 
   then
   ." prime." cr
loop drop ;

sieveMain




Okay, it looks like it ought to make sense, right?

But, if you are new to FORTH, and used to algol-like languages like C and Ada, it doesn't. So, let's look at it again, with comments that explain the FORTH itself.

But before we do that I should explain the stack.

FORTH's native interpreter operates in a post-fix grammar. In other words,
MAXSIEVE - 1
is written
MAXSIEVE 1 -
When the interpreter sees the constant MAXSIEVE, it pushes it on the stack. Then it sees 1 and pushes it on the stack. Then it sees the subtraction operator and subtracts the top item on the stack from the second:

inputstack
MAXSIEVEMAXSIEVE
1MAXSIEVE 1
-(result of MAXSIEVE - 1)

Or, borrowing from the code:
 MAXSIEVE 1- 2 / constant FINALPASS

inputstackwhat happens
MAXSIEVEMAXSIEVEConstant pushed on stack.
1-(result of decrement)Top item of stack gets decremented.
2(result1) 2 Numeric constant pushed on stack.
/(result of division)2nd on stack gets divided by top.
constant FINALPASS(empty)Constant FINALPASS generated in symbol table.

Now we can look at the code again:



256 constant MAXSIEVE
\ Defines MAXSIEVE as a true constant.
\ Incidentally, gforth is case insensitive, so we could call it "maxsieve", too.
\ But we won't.

MAXSIEVE 1- 2 / constant FINALPASS
\ See the explanation above.

5 constant DISPWIDTH ( enough columns to display MAXSIEVE )
\ We want this for integer output.

create sieve MAXSIEVE allot
\ This names and allocates a byte array named "sieve".
\ "sieve" in the input will now leave the array address on the stack:
\ ( -- adr )


: sieveMain ( --)
\ The colon operator defines a new word in the FORTH dictionary.
\ Or you can think of it as a new function or operator.
\ It is a symbol table entry with code that can be executed.
\ The ( -- ) is just a comment that shows tha there are no stack effects.
\ Some dialects of FORTH actually parse the stack effects, by the way.

\ "c!" stores the 2nd on stack at the address on top:
\ The stack effect is ( n adr -- )

0 sieve c!      ( 0 is not prime. )\ So, this stores a 0 in the first entry in the array,
\ for the integer zero.

0 sieve 1+ c!   ( 1 is not prime. ) \ And this stores a 0 in the second entry in the array,
\ for the integer one.


\ ?do 

\ Now we put the address on the stack again,
\ and set up a counted loop
\ with the limit and the initial count on stack
\ and start the loop, testing the limit against the count:
sieve MAXSIEVE 2 ?do    ( set flags to true for 2 to FINALPASS. )

\ "-1" is read as the numeric constant negative one.

\ "over" copies the 2nd on the stack and pushes it on top,
\ the stack effect is ( adr n -- adr n adr )

\ "i" pushes the current value of the loop counter on the stack.
\ The stack effect is ( -- counter )

\ "+" adds the top two numbers on stack, of course.
\ The stack effect is ( n1 n2 -- sum ).

\ "loop" marks the end of the loop.
\ It has no parameter stack effect: ( -- ),
\ but it removes the limit and current loop counter
\ from where they are stored.
\ It also branches back to the beginning of the loop,
\ if the limit has not been exceeded.

\ So, in the following code,
\ the loop counter is added to the array address,
\ and -1 is stored into the array as a byte value
\ at each byte in the array.

   -1 over i + c! loop  ( sieve pointer -- still on stack. )

\ This is the nested loop that clears the non-primes
\ by synthetic multiplication.
\ Note that "i" is always the innermost loop counter.
\ When we need to access the outer loop counter in the inner loop,
\ we do so with "j".

\ "?do" tests the limit against the initial value
\ and does not enter the loop if the limit is already exceeded.
\ "dup" duplicates the top of stack, effect is ( n -- n n ).
\ "c@" gets the byte value at the address on stack
\ and replaces the address with the value obtained.
\ Stack effect here is ( adr -- value )

\ "if" sets up a conditional evaluation.
\ Stack effect is ( flag -- ).
\ If the flag is zero, code up to "then" is skipped.
\ If the flag is non-zero, execution continues,
\ up to an optional "else".
\ The postfix-influenced interpretation of "then"
\ may be a little hard to get used to.
\ It is often aliased to "endif".

\ "+loop" is like loop, but it counts by the increment on top of stack,
\ stack effect is ( count -- )

FINALPASS 2 ?do ( clear flags at all multiples. )
\ Since the flag clearing pass starts at double each prime,
\ any prime over half of MAXSIEVE will never need to be looked at.

   dup i + c@ if      ( don't bother if not prime. )
      MAXSIEVE i dup + ?do    ( start at first multiple -- double. )
         0 over i + c!      ( clear at this multiple. )
      j +loop           ( sieve still on stack. )
\ Since j is the outer loop counter,
\ it is the current prime being scanned.
   then
loop            ( sieve still on stack. )

\ ".r" prints the 2nd integer on the stack
\ in a field width given by the top of stack.
\ "."" (sorry about the quotes) prints the following string,
\ up to the next quote.

\ "0=" inverts the flag on top of stack.
\ Non-zero become zero, zero becomes -1 .

\ "cr" writes a carriage return out.

\ "drop" drops a value from the stack,
\ stack effect is ( n -- ).

\ ";" terminates the colon definition.

MAXSIEVE 0 ?do
   i DISPWIDTH .r ." : is "
   dup i + c@ 0= if
      ." not "
   then
   ." prime." cr
loop drop ;

\ Now we can execute the newly defined function:

sieveMain
\ which generates the sieve and prints it.



As I said, this FORTH source looks like it was written by a C programmer. (It was. I am far more familiar with C than FORTH.)

You notice that there don't seem to be any local variables in FORTH. Modern FORTHs do have named local variables, but most local variables in FORTH are just parameters of expressions, anonymous on the stack.

Instead of local variables, FORTH traditionally cut functions up in smaller pieces. Some exercise of care allowed heavy re-use of those pieces, yielding more compact, highly readable code. (And then you could go to far and the code could become too dense to read.)

FORTH allows one to make new vocabularies, and that would allow us to hide the small pieces.

Finally, trying to use this on my bif-c revealed some of the bugs I haven't squashed yet. Trying to use it on the fig-FORTH for 6800 I transcribed some time back shows the wisdom of using the small pieces. The small pieces are far easier to debug. This works with the fig-FORTH, and, with the vocabulary stuff commented out, with my bif-c. (And now I have some debugging work to do on the bif-c.)


( Archetypical implementation of the sieve of eratosthenes in FORTH -- fig, bif-c --  )
( using more of the FORTH idiom. )
( Copyright Joel Rees, 2015 )
( By Joel Rees, Amagasaki, Japan, 2015 )
( All rights reserved. )
( Permission granted by the author to use this code )
( for any purpose,  )
( on condition that substantial use  )
( shall retain this copyright and permission notice. )



VOCABULARY sieve-local ( Make a local symbol table. )
sieve-local DEFINITIONS ( Switch to the local vocabulary. )


256 CONSTANT MAXSIEVE
MAXSIEVE 1 - 2 / CONSTANT FINALPASS

5 CONSTANT DISPWIDTH ( enough digits to display MAXSIEVE )


0 VARIABLE sieve ( Old FORTHs don't provide a default behavior for CREATE )
( gforth will leave the zero there. Old FORTHs need an initial value. )

   HERE sieve - DUP ( Old FORTHs don't provide a CELL width. )
   MAXSIEVE SWAP - ALLOT ( Allocate the rest of the byte array. )

   CONSTANT CELLWIDTH ( To show how this can be done. )

: sieveInit ( -- adr )
0 sieve C!      ( 0 is not prime. )
0 sieve 1+ C!   ( 1 is not prime. )
sieve MAXSIEVE 2 DO    ( set flags to true for 2 to FINALPASS. )
   -1 OVER I + C! LOOP  ( sieve pointer -- still on stack. )
;

: primePass ( adr prime -- adr )
MAXSIEVE OVER DUP + DO    ( start at first multiple -- double. )
   OVER I + 0 SWAP C!     ( clear at this multiple. )
   DUP +LOOP              ( next multiple )
DROP ;      ( sieve address still on stack. )

: findPrimes ( adr -- adr )
FINALPASS 2 DO   ( clear flags at all multiples. )
   DUP I + C@ IF ( don't bother if not prime. )
      I primePass
   ENDIF
LOOP ;           ( sieve still on stack. )

: printPrimes ( adr -- )
MAXSIEVE 0 DO
   I DISPWIDTH .R ." : is "
   DUP I + C@ 0= IF
      ." not "
   ENDIF
   ." prime." CR
LOOP DROP ;


FORTH DEFINITIONS

: sieveMain ( -- )
[ sieve-local ] sieveInit
findPrimes
printPrimes ;


sieveMain



If you are wondering about the speed of this version compared to the long version, it's pretty close. With the print commented out and the size set to 4194304, here's what we get in gforth:
me@fun:~/work/mathgames/sieve$ time gforth forthsieve.fs -e bye

real    0m1.096s
user    0m1.080s
sys    0m0.008s
which looks about right for 7 instructions vs. 6 in the inner loop.

(In bif-c, with the buggy parts commented out, too, it takes about 5 seconds. I guess, while I'm debugging it, I should add an option to parse an input file and go, so I can get accurate timings. But timings really aren't meaningful yet, there.)

You might also be interested in the FORTH language sieve entry in Rosetta Code.
Comparing the times for the current example using the dictionary buffer without allocating it:
me@fun:~/work/mathgames/sieve$ time gforth rs_sieve.fs -e bye

real    0m1.463s
user    0m1.444s
sys    0m0.012s
in gforth.

[JMR 2 Jan 2020: Forth makes it fairly easy to query the sieve. See how to do that, and a comparable approach in C, here.]

No comments:

Post a Comment