A coding challenge: prime numbers

[BigEd]

(continued the VTL porting discussion over here.)

[BigEd]

Oddly, at least one of my ideas for running faster doesn't actually run any faster, at least when I try it in Basic. And I notice Andrew's code just walks through the whole sieve without trying to be clever about stopping early, and is faster than my code which does try to stop early. There's a lesson there about premature optimisation.

For the record, here are the two approaches I wrote up in Basic before tackling the challenge in assembly language.

Here's the simpler one, expected to be slower:

Code: Select all

   10 N=100
   20 P=0
   30 DIM S(N)
   40 FOR I=2 TO N
   50 IF S(I)<>0 THEN NEXT:STOP
   60 PRINT ;I;" ";
   70 P=P+1
   80 IF I+I>N THEN GOTO 130
   90 FOR J=I+I TO N STEP I
  100 S(J)=-1
  110 NEXT
  120 NEXT I
  130 FOR K=I+1 TO N
  140 IF S(K)=0 THEN PRINT ;K;" ";:P=P+1:IF P=20 STOP
  150 NEXT K

Here's the more "clever" one, expected to be faster, but not:

Code: Select all

   10 N=100
   20 P=0
   30 DIM S(N)
   40 FOR I=2 TO N
   50 IF S(I)<>0 THEN NEXT:STOP
   60 PRINT ;I;" ";
   70 P=P+1
   80 C=1
   90 FOR J=I+I TO N STEP I
  100 C=C+1
  110 S(J)=-1
  120 NEXT J
  130 IF C>I THEN NEXT I
  140 FOR K=I+1 TO N
  150 IF S(K)=0 THEN PRINT ;K;" ";:P=P+1:IF P=20 STOP
  160 NEXT K

[Martin_H]

I wrote an ANS Forth sieve of Eratosthenes a few years back to teach myself the language. One of these days I need to get Forth working on my 6502 SBC.

Code: Select all

\ Implements the sieve of Eratosthenes for the first 1024 numbers

\ Begin borrowed code from http://www.forth.org/svfig/Len/bits.htm
\ Note: that code had a bug I fixed as indicated below.

\ Bit_array is a defining word that will create a bit array.
\ l is the number of bytes
create masks 128 c, 64 c, 32 c, 16 c, 8 c, 4 c, 2 c, 1 c,
: bit_array \ ( len -- ) ( i -- m a)
   create allot does>
      swap   \ a i
      8      \ a i 8
      /mod   \ a remainder whole
      swap   \ a whole remainder !MCH bug fix!
      masks  \ a whole remainder a
      +      \ a whole mask_a
      c@     \ a whole m
      -rot   \ m a whole
      + ;    \ m a

\ We also need words to store and fetch bits. The words .@ and .!
\ Will fetch and store, respectively, a Boolean flag on the stack to
\ a bit array. (Do not confuse .! with cset or .@ with creset.) 

: .! ( f m a -- ) rot if cset else creset then ;
: .@ ( m a -- f ) c@ and 0<> ;

\ Examples
\ 3 bit-array thflag  \ Create a bit array for 24 bit-flags
\ 11 thflag ctoggle    \ Toggle the 11th bit
\ 10 thFLAG .@ ( -- f) \ Extract the 10th bit and leave as well-formed Boolean flag

\ End of borrowed code.

\ Create a bit vector to hold the sieve.
\ 128 bytes allows for all primes less than 1024.
128 bit_array sieve

\ locates the index of the first non-zero bit in the sieve
\ with an index greater than the input.
: find_one ( n -- n )
   1+ dup 1024 swap do
      i sieve .@ \ get the bit corresponding to the integer.
      if
         leave
      else
         1+
      then
   loop ;

\ sets bits starting at n.
: set_bits ( n -- )
   1024 swap
   do
      i sieve cset
   loop ;

\ clears bits starting at n in steps of n. So n=2 starts at 2
\ and clears every other bit. While n=3 starts at three and
\ clears every third bit. However, it then resets the first
\ bit as that number is prime.
: clear_bits ( n -- )
   dup 1024 swap
   do
      i sieve creset
      dup
   +loop
   sieve cset ;

\ We'll waste two bits for integers 0 and 1 which can't be prime
: do_sieve
  2 set_bits   \ Assume all numbers 2 or greater are prime.
  2            \ Initially, let p equal 2, the first prime number.

  begin dup 1024 < while
     dup          \ n n
     clear_bits   \ n
     find_one     \ index of first one > n
  repeat ;

\ iterates through the bit vector printing the index of a prime number
: .sieve
   ." List of primes "
   1024 0 do
      i sieve .@ \ get the bit corresponding to the integer.
      if
         i . ." , "
      then
   loop ;

[barrym95838]

Both of those Beeb versions complete in less than one second at my end, Ed. I don't know what the effective clock speed is on my setup. I would go with the clearest and/or shortest, unless I was competing with someone solely for speed.

Klaus' VTL-2 code runs to completion in about 25 seconds on a 1 MHz Apple 2, but (as he mentioned), it is interpreted un-optimized brute-force ... the perfect recipe for a bit of thumb-twiddling. The code density and run-time footprint are both quite decent, though.

Mike B.

[BigEd]

For interest and for comparison I converted Andrew's code into Basic, with a couple of necessary tweaks:

Code: Select all

   10 P=20: N=100
   20 I=2: DIM S(N)
   30 REPEAT
   40 IF S(I)=0 PRINT ;I;" ";: FOR J=I TO N STEP I: S(J)=J: NEXT: P=P-1: IF P=0 STOP
   50 I=I+1
   60 UNTIL I=N

(This is a tiny bit slower than my Basic versions but still under a second.)
(The original BBC Basics don't allow for nested IF THEN ELSE so I had to improvise.)

[BitWise]

BigEd wrote:

For interest and for comparison I converted Andrew's code into Basic, with a couple of necessary tweaks:

Code: Select all

   10 P=20: N=100
   20 I=2: DIM S(N)
   30 REPEAT
   40 IF S(I)=0 PRINT ;I;" ";: FOR J=I TO N STEP I: S(J)=J: NEXT: P=P-1: IF P=0 STOP
   50 I=I+1
   60 UNTIL I=N

(This is a tiny bit slower than my Basic versions but still under a second.)
(The original BBC Basics don't allow for nested IF THEN ELSE so I had to improvise.)

It should go faster if you use integer variables like I%, P% etc. throughout -- BEEBEM makes it 0.59 secs with I, J and 0.56 sec with integers.

[barrym95838]

BigEd wrote:

For interest and for comparison I converted Andrew's code into Basic, with a couple of necessary tweaks:

Code: Select all

   10 P=20: N=100
   20 I=2: DIM S(N)
   30 REPEAT
   40 IF S(I)=0 PRINT ;I;" ";: FOR J=I TO N STEP I: S(J)=J: NEXT: P=P-1: IF P=0 STOP
   50 I=I+1
   60 UNTIL I=N

I translated this to VTL-2, with a couple of necessary tweaks:

Code: Select all

10 P=20
15 N=100
20 I=2
25 J=I
30 :J)=0
35 J=J+1
40 #=J<N*30
45 #=:I)*90
50 ?=I
55 $=32
60 J=I
65 :J)=1
70 J=J+I
75 #=N>J*65
80 P=P-1
85 #=P=0*99 
90 I=I+1
95 #=I<N*45

It runs to completion in about 2.5 seconds on a 1 MHz Apple 2.

Mike B.

[BigEd]

(I've ported Andrew's code to easy6502, it runs much quicker than visual6502:

Code: Select all

; some constants
define ZERO $30
define SPACE $20
define COUNT   20 ; Primes left to find

; some zero page locations
define CHROUT $0f ; output device (visual6502)
define INDEX  $21 ; Index in decimal
define FLAGS  $30 ; Prime flags (sieve)

  LDY   #COUNT              ; Set target count
  LDA   #2                  ; And starting index
  STA   INDEX
  SED                      ; Use decimal
repeat:
  TAX                      ; Is this number prime?
  LDA   FLAGS,X
  BNE checknext
  LDA   #SPACE                  ; Yes, print it
  wdm   0 ; or STA CHROUT
  TXA
  LSR   A
  LSR   A
  LSR   A
  LSR   A
  BEQ skip  ; tens digit zero?
  ORA   #ZERO ; to ascii
  wdm   0 ; or STA CHROUT
skip:
  TXA
  AND   #$0f  ; units digit
  ORA   #ZERO ; to ascii
  wdm   0 ; or STA CHROUT
  DEY                          ; Found twenty?
  BEQ done                     ; Yes

  CLC
  TXA
marking:
  ADC   INDEX
  BCS clearcarry ; BREAK CS at 100 (decimal)
  TAX
  STA   FLAGS,X
  BNE marking ; UNTIL EQ ; until 100

clearcarry:
  CLC
checknext:       ; Check the next index
  LDA   INDEX
  ADC   #1 ; carry is clear
  STA   INDEX
  BNE repeat ; UNTIL EQ (to 100)

done:
  NOP            ; Done (previously an RTS)

)

[Klaus2m5]

I removed the feature "thumb-twiddling" from my VTL-program. It is completely unreadable and ugly now, but still uses trial divisons (line 310 & 350) to find prime numbers.

Code: Select all

10 #=1000
100 N=N+F
110 #=S
120 N=N+G
150 !=T
200 C=!
210 #=N<U*Q
220 Q=M
230 V=V+F
250 U=V*V
300 D=H
310 A=N/D
320 #=%=O*C
330 D=D+F
340 #=D>V*R
350 A=N/D
360 #=%=O*C
370 D=D+G
380 #=D<V*P
400 ?=N
410 ?=""
420 X=X-E
430 #=X>E*C
440 #=9999
500 #=S
510 N=N+E
520 #=S
530 N=N+F
540 #=150
1000 E=1
1010 F=2
1020 G=4
1030 H=5
1050 O=0
1100 T=100
1110 S=200
1120 R=400
1130 Q=R
1140 P=310
1150 M=300
1200 V=5
1210 U=25
1220 X=20
1230 N=2
1240 #=500

Feel free to increase the count of prime numbers in line 1220.

edit: added a view improvements

[BigEd]

Hee hee. I see you are using the "come from" construction.

Back in the day, with my TI-57, I used the constant 42424626 to skip factors of 2, 3 and 5 - I extracted digits from that constant to use as deltas probably by shifting it in a circular fashion. (I'll admit I did need to reconstruct and check that number just now - I only remembered the general form of it.)

I think I wrote out the primes up to about 10,000 before my enthusiasm waned.

Funnily enough I don't think I've ever written a sieve program until this challenge.

[Klaus2m5]

Mine is a bit simpler as it does reduce the prime number canidates only to 2 of 6 dividends (no even numbers and no 3rd odd number) and I am not doing this for the divisor (room for improvment). The main performance boost resulted from the perception, that you only need to divide through odd numbers and when the dividend is >= the square of the divisor. So I only start to divide by 5 at 25, by 7 at 49 and so on.

A sieve for 10000 primes would require too much real estate for a small computer, so you probably had not other choice at the time as to do it with optimized trial division. However, I will never reach 10000 with VTL because of the 16-bit integer limit.

[BigEd]

I wasn't quite so industrious, Klaus - only primes up to 10,000, which would mean up to 9973. So VTL can do much better than I did.

BTW, you might enjoy this prime sieve in physical form:

[jgharston]

In the 8BS disk magazine for the BBC computer there was a discussion on comparative prime number generating algorithms: here.

[BigEd]

That's interesting... but ouch, all that repeated calculating of the same square root!

[barrym95838]

The Beeb is a 2MHz machine, right?

Code: Select all

                     Table of Comparative Program Times
                     ----------------------------------
                      Running time in seconds on BBC B

         n     P(n)   Prime1   Prime2   Prime3   Prime4   Prime5   Prime6
         10      29      1.79     0.46     0.44     0.32     0.32     0.36
        100     541    221.61    11.96    10.31     7.43     7.98     8.63
       1000    7919  27500      331.74   268.75   175.77   185.05   195.74

From another thread:

Quote:

cheers, Klaus

P.S.: I just ran my prime number program for 1000 primes on my emulator (~2MHz) and compared the performance. VTL02A = 72 seconds, VTL02B = 85 seconds

I didn't know that VTL02 could easily beat BASIC in execution speed. Of course, it has fewer features, but still ... are we looking at the same task here?

Mike B.