A coding challenge: prime numbers

[BigEd]

This might be fun, as we have several languages to choose from on our 6502, 65816 and maybe even 65Org16 systems: the challenge from The National Museum of Computing is to write a program that prints out the first 20 prime numbers.

I had a go this afternoon, and cooked up a working version in Basic, and then rewrote it with some difficulty in 6502 assembly language. (That was 90 bytes, if anyone is counting, written for the visual6502 platform.)

(Jac's approach of printing out the right 20 numbers without computing them will get no points from me!)

Anyone for Forth?

[jac_goudsmit]

For reference, this was my solution...

Code: Select all

 
10 data 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71
20 for i=1 to 20:read p:print p:next
30 end

It meets the requirements of printing the first 20 prime numbers. You didn't say anything about calculating them

Oh well, as much as I would love to go to Bletchley Park to pick up my prize, I wouldn't be able to do that anyway.

===Jac

[BigEd]

Let's call that a reference solution!

[Rob Finch]

Quote:

This might be fun, as we have several languages to choose from on our 6502, 65816 and maybe even 65Org16 systems: the challenge from The National Museum of Computing is to write a program that prints out the first 20 prime numbers.

I found a standard sieve program written in 'C' on the NET which I modified for testing the RTF65003 core. I've left out some of the function calls like printf() and getch() that also had to be coded. The constant PRIMES should be set to 20 to find the first 20 primes. This program also tracks the number of clock ticks to complete the operation.

Code: Select all

#define LIMIT 1500000
#define PRIMES 100000

int numbers[1500000];
int primes[100000];

int main(){
    int i,j;
	int limit;
	int start_tick,end_tick;
	int ch;

	// Request I/O focus
	asm {
		jsr	(0xFFFF8014>>2)
		jsr (0xFFFF801C>>2)		; clear screen
		jsr	(0xFFFF8020>>2)		; home cursor
	}
j1:
	//disable_ints();
	firstcall {
		printf("This appears the first time only.\r\n");
	}
	start_tick = get_tick();
	printf("\r\nStart tick %d\r\n", start_tick);

	limit=LIMIT;

 	for (i=0;i<limit;i++){
		numbers[i]=i+2;
	}

    /*sieve the non-primes*/
    for (i=0;i<limit;i++){
        if (numbers[i]!=-1){
            for (j=2*numbers[i]-2;j<limit;j+=numbers[i])
                numbers[j]=-1;
        }
    }

    /*transfer the primes to their own array*/
    j = 0;
    for (i=0;i<limit&&j<PRIMES;i++)
        if (numbers[i]!=-1)
            primes[j++] = numbers[i];

	end_tick = get_tick();
	printf("Clock ticks %d\r\n", end_tick-start_tick);
	printf("Press a key to list primes.");
	getchar();

    /*print*/
	for (i=0;i<PRIMES;i++) {
		if ((i&0xf)==0) {
			ch = getcharNoWait() & 0xff;
			if (ch==3)	// CTRL-C
				break;
		}
        printf("%d\r\n",primes[i]);
	}

	//enable_ints();

return 0;
}

[BigEd]

Thanks Rob - infinitely more readable than mine!

It's a curious thing that you need to estimate the size of a sieve before you know how big it needs to be - I had to do the same, as my efforts today in Basic and 6502 were both sieves. It seems the right approach for a machine without multiply or divide.

You might lose points for printing more than 20! But it's readily fixed, of course. That very large LIMIT would stress a machine with just 16 bits of address, but of course you are targeting a bigger machine than that.

My 6502 code uses a little over 5000 clock ticks to find and print 20 primes. Printing in decimal probably took a fair number of those ticks.

[BitWise]

Slightly smaller and faster but the same type of sieving algorithm.

http://visual6502.org/JSSim/expert.html ... 8521D0CC60

Code: Select all

                            ;===============================================================================
                            ; Data Areas
                            ;-------------------------------------------------------------------------------
                            
                                            .PAGE0
                                            .ORG    $20
                                            
000020  00                : COUNT           .SPACE  1                       ; Primes left to find
000021  00                : INDEX           .SPACE  1                       ; Index in decimal
                            
                                            .ORG    $30
                            FLAGS           .SPACE  0                       ; Prime flags
                            
                            ;===============================================================================
                            ; Prime Finder
                            ;-------------------------------------------------------------------------------
                            
                                            .CODE
                                            .ORG    $1000
                            RESET:
001000  A914              :                 LDA     #20                     ; Set target count
001002  8520              :                 STA     COUNT
001004  A902              :                 LDA     #2                      ; And starting index
001006  8521              :                 STA     INDEX
                                            
001008  F8                :                 SED                             ; Use decimal
                                            REPEAT

Portable 65xx Assembler [14.10]

001009  A621              :                  LDX    INDEX                   ; Is this number prime?
00100B  B530              :                  LDA    FLAGS,X
00100D  D025              :                  IF     EQ
00100F  A920              :                   LDA   #' '                    ; Yes, print it
001011  850F              :                   STA   $0F
                                              
001013  8A                :                   TXA
001014  4A                :                   LSR   A
001015  4A                :                   LSR   A
001016  4A                :                   LSR   A
001017  4A                :                   LSR   A
001018  F004              :                   IF    NE
00101A  0930              :                    ORA  #'0'
00101C  850F              :                    STA  $0F
                                              ENDIF
                                              
00101E  8A                :                   TXA
00101F  290F              :                   AND   #$0F
001021  0930              :                   ORA   #'0'
001023  850F              :                   STA   $0F
                                             
001025  C620              :                   DEC   COUNT                   ; Found twenty?
001027  F014              :                   BREAK EQ                      ; Yes
                                             
001029  18                :                   CLC
                                              REPEAT                        ; No, mark multiples
00102A  8A                :                    TXA
00102B  6521              :                    ADC  INDEX
00102D  B005              :                    BREAK CS
00102F  AA                :                    TAX
001030  F630              :                    INC  FLAGS,X
001032  D0F6              :                   UNTIL EQ
                                             ENDIF
                                             
001034  18                :                  CLC                            ; Check the next index
001035  A521              :                  LDA    INDEX
001037  6901              :                  ADC    #1
001039  8521              :                  STA    INDEX
00103B  D0CC              :                 UNTIL EQ
                                            
00103D  60                :                 RTS                             ; Done

[BitWise]

Changing the INC FLAGS,X to STA FLAGS,X makes the key loop a few cycles shorter on each iteration.

http://visual6502.org/JSSim/expert.html ... 8521D0CC60

.. And the initial LDX INDEX could be replaced with a TAX

http://visual6502.org/JSSim/expert.html ... 8521D0CD60

[MichaelM]

For the FORTH afficionados, I posted the source for two different Sieve programs. (Note: the second program is substantially faster than the first one given at the link.)

Although not exactly in compliance with the challenge, they do appear to find all of the primes from 0 to the limit specified. I imagine that either could be modified to comply with the challenge, but the modifications required are presently beyond my skill level in Forth. I would love to see the code for a modified version that meets the guidelines for the challenge.

[theGSman]

For completeness, here is a BASIC listing and run session as entered on the Commodore 64.
(Remember that the C64 screen is 40 characters wide).

Code: Select all

LIST

100 MA=100:MP=20:T1=TI:DIM NP(MA)
110 FOR I=2 TO MA:NP(I)=0:NEXT
120 D=2
130 IF D*D>MA THEN 170
140 IF NP(D) THEN 160
150 FOR I=2*D TO MA STEP D:NP(I)=1:NEXT
160 D=D+1:GOTO 130
170 PRINT "CALCULATION TOOK";TI-T1;" JIFFY 
 TICKS":PRINT
180 PRINT"PRIMES:";:I=2:C=0 
190 IF NP(I) THEN I=I+1:GOTO 190
200 PRINT I;:I=I+1
210 C=C+1:IF C<MP THEN PRINT",";:GOTO 190

220 PRINT
READY.
RUN
CALCULATION TOOK 76 JIFFY TICKS

PRIMES: 2 , 3 , 5 , 7 , 11 , 13 , 17 , 1
9 , 23 , 29 , 31 , 37 , 41 , 43 , 47 , 5
3 , 59 , 61 , 67 , 71

READY.

[BigEd]

Very nice theGSman - here's a conversion to Beeb: http://goo.gl/1OcLtu
Edit: original link

Thanks for your two offerings Michael - I notice one prints out 1900 which is I think the correct number of primes less than 16384, but the other prints out 1899. One thing I like about the 20 primes challenge is that it could easily show up an off-by-one error in the counting. Points deducted for off-by-one errors!

(I could just about run your SIEVE and sieve words at https://repl.it/languages/forth)

Andrew, that's a very short program indeed! Well done! I see SED was put to good use. And the structured control flow in your portable assembler looks very handy.

[theGSman]

MichaelM wrote:

For the FORTH afficionados, I posted the source for two different Sieve programs.

Your second "sieve" routine seems to be similar to the way I wrote mine but I can't for the life of me figure out how the printing ("primes") routine does the job.

For comparison, this is how I coded the routines in FORTH. Note that I like to intersperse my code with stack diagrams and try to find meaningful names for the variables in the hope that it makes the code easier to follow.

Code: Select all

255 CONSTANT TRUE 0 CONSTANT FALSE
100 CONSTANT MAX_NUM
20 CONSTANT MAX_PRIMES
CREATE isPrime MAX_NUM ALLOT

: sieve ( -- )
   MAX_NUM 2 DO
      TRUE isPrime I + !
   LOOP
   2 ( start of divisor )
   BEGIN
      DUP isPrime + C@ ( d isPrime[d] )
      IF
         MAX_NUM OVER 2* ( d MAX_NUM 2d ) 
         DO
            FALSE I isPrime + C!
            DUP
         +LOOP
      THEN
      1+ DUP DUP * MAX_NUM >=
   UNTIL DROP
;

: printSieve ( -- )
   2 ( first prime number )
   0 ( count )
   CR ." PRIMES: "
   BEGIN
      OVER isPrime + C@ ( n c isPrime[n] )
      IF
         OVER .
         1+ DUP MAX_PRIMES >= ( n c+1 c+1>=MAX_PRIMES? )
         IF
             CR 2DROP EXIT
         ELSE
             ." , "
         THEN
      THEN
      SWAP 1+ SWAP
   AGAIN
;

: erasthosenes
   sieve
   printSieve
;

[Klaus2m5]

In VTL (brute force and not optimized):

Code: Select all

10 X=0
20 N=1
100 N=N+1
110 A=N-1
120 A=A-1
130 #=A<2*200
140 B=N/A
150 #=%=0*100
160 #=120
200 ?=N
210 ?=""
220 X=X+1
230 #=X<20*100

[BigEd]

That's great! (Recipe to run VTL in Kowalski here - I tried and so far failed to port VTL to Beeb.)

[BigEd]

Andrew, I think I found some more small optimisations of your code:
- use Y for the COUNT.
- hoist the TXA out of the "mark multiples" repeat loop
- hoist the CLC for "check next index" one line, into the "found a prime" branch

http://visual6502.org/JSSim/expert.html ... 8521d0ce60

There are a couple of algorithmic optimisations which might or might not cost more than they save, for such a small sieve size, but would surely be worth it for a bigger problem:
- only store odd numbers in the sieve, printing 2 as a special case. Jac would approve.
- quit scanning for new primes when the index is beyond the square root. No need for anything stronger than counting: count the steps in the marking loop and compare it to the index. (And account for the half-size sieve...)

If keeping the sieve at full size, so it contains even and odd numbers, once 2 is dealt with the marking loop need only visit odd multiples of the index, by doubling the stride. Again it means special-casing 2.

With your cunning use of decimal mode, I'm not sure how easy those might be to implement.

[barrym95838]

BigEd wrote:

... I tried and so far failed to port VTL to Beeb ...

Is it the I/O or page-zero conflicts that's throwing a wrench in the works? The Apple 1 and Apple 2 monitors only need about a dozen (or less) zp locations, and they're all in the lower half, so it was safe for me to usurp $80 to $ff for the interpreter. That may not be a universally satisfactory arrangement.

Mike B.

P.S. It seems that Andrew's the man for finding cunning uses of decimal mode.