6502.org Forum  Projects  Code  Documents  Tools  Forum
It is currently Sun Nov 24, 2024 10:14 am

All times are UTC




Post new topic Reply to topic  [ 42 posts ]  Go to page Previous  1, 2, 3  Next
Author Message
PostPosted: Wed Sep 30, 2015 6:46 pm 
Offline
User avatar

Joined: Thu Dec 11, 2008 1:28 pm
Posts: 10986
Location: England
(continued the VTL porting discussion over here.)


Top
 Profile  
Reply with quote  
PostPosted: Thu Oct 01, 2015 3:29 pm 
Offline
User avatar

Joined: Thu Dec 11, 2008 1:28 pm
Posts: 10986
Location: England
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:
   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:
   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


Top
 Profile  
Reply with quote  
PostPosted: Thu Oct 01, 2015 4:59 pm 
Offline

Joined: Wed Jan 08, 2014 3:31 pm
Posts: 578
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:
\ 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 ;


Top
 Profile  
Reply with quote  
PostPosted: Thu Oct 01, 2015 5:52 pm 
Offline
User avatar

Joined: Sun Jun 30, 2013 10:26 pm
Posts: 1950
Location: Sacramento, CA, USA
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.


Top
 Profile  
Reply with quote  
PostPosted: Thu Oct 01, 2015 7:08 pm 
Offline
User avatar

Joined: Thu Dec 11, 2008 1:28 pm
Posts: 10986
Location: England
For interest and for comparison I converted Andrew's code into Basic, with a couple of necessary tweaks:
Code:
   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.)


Top
 Profile  
Reply with quote  
PostPosted: Fri Oct 02, 2015 12:07 am 
Offline
User avatar

Joined: Tue Mar 02, 2004 8:55 am
Posts: 996
Location: Berkshire, UK
BigEd wrote:
For interest and for comparison I converted Andrew's code into Basic, with a couple of necessary tweaks:
Code:
   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.

_________________
Andrew Jacobs
6502 & PIC Stuff - http://www.obelisk.me.uk/
Cross-Platform 6502/65C02/65816 Macro Assembler - http://www.obelisk.me.uk/dev65/
Open Source Projects - https://github.com/andrew-jacobs


Top
 Profile  
Reply with quote  
PostPosted: Fri Oct 02, 2015 3:50 am 
Offline
User avatar

Joined: Sun Jun 30, 2013 10:26 pm
Posts: 1950
Location: Sacramento, CA, USA
BigEd wrote:
For interest and for comparison I converted Andrew's code into Basic, with a couple of necessary tweaks:
Code:
   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:
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.


Top
 Profile  
Reply with quote  
PostPosted: Sat Oct 03, 2015 11:33 am 
Offline
User avatar

Joined: Thu Dec 11, 2008 1:28 pm
Posts: 10986
Location: England
(I've ported Andrew's code to easy6502, it runs much quicker than visual6502:
Code:
; 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)

)


Top
 Profile  
Reply with quote  
PostPosted: Sat Oct 03, 2015 12:33 pm 
Offline

Joined: Sat Jul 28, 2012 11:41 am
Posts: 442
Location: Wiesbaden, Germany
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:
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

_________________
6502 sources on GitHub: https://github.com/Klaus2m5


Last edited by Klaus2m5 on Mon Oct 05, 2015 2:24 pm, edited 2 times in total.

Top
 Profile  
Reply with quote  
PostPosted: Sat Oct 03, 2015 2:39 pm 
Offline
User avatar

Joined: Thu Dec 11, 2008 1:28 pm
Posts: 10986
Location: England
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.


Top
 Profile  
Reply with quote  
PostPosted: Sun Oct 04, 2015 5:35 am 
Offline

Joined: Sat Jul 28, 2012 11:41 am
Posts: 442
Location: Wiesbaden, Germany
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.

_________________
6502 sources on GitHub: https://github.com/Klaus2m5


Top
 Profile  
Reply with quote  
PostPosted: Sun Oct 04, 2015 3:35 pm 
Offline
User avatar

Joined: Thu Dec 11, 2008 1:28 pm
Posts: 10986
Location: England
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:
Attachment:
Lehmer-Sieve-1928.png
Lehmer-Sieve-1928.png [ 532.7 KiB | Viewed 1204 times ]


Last edited by BigEd on Tue Oct 25, 2022 9:06 pm, edited 1 time in total.

Top
 Profile  
Reply with quote  
PostPosted: Fri Oct 23, 2015 8:32 pm 
Offline

Joined: Sun Feb 22, 2004 9:01 pm
Posts: 111
In the 8BS disk magazine for the BBC computer there was a discussion on comparative prime number generating algorithms: here.

_________________
--
JGH - http://mdfs.net


Top
 Profile  
Reply with quote  
PostPosted: Fri Oct 23, 2015 8:44 pm 
Offline
User avatar

Joined: Thu Dec 11, 2008 1:28 pm
Posts: 10986
Location: England
That's interesting... but ouch, all that repeated calculating of the same square root!


Top
 Profile  
Reply with quote  
PostPosted: Fri Oct 23, 2015 9:33 pm 
Offline
User avatar

Joined: Sun Jun 30, 2013 10:26 pm
Posts: 1950
Location: Sacramento, CA, USA
The Beeb is a 2MHz machine, right?
Code:
                     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.


Top
 Profile  
Reply with quote  
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 42 posts ]  Go to page Previous  1, 2, 3  Next

All times are UTC


Who is online

Users browsing this forum: No registered users and 6 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to: