Prime Number Generation ======================= In issue 50 R.Harker (D6E) wrote a short and simple prime number generator (program Primes1). It went up all integers in sequence testing if any other number less than it would divide it evenly, ie if it had any factors. REM > Primes1 REM D6E (R. Harker)'s original version MODE0 place=1 divisable=0 number=0 T%=TIME:REPEAT FOR divisable=1 TO number IF (number MOD divisable)=0 THEN a=a+1 IF a>2 THEN divisable=number IF (number MOD 2=0 AND number>2) THEN divisable=number NEXT divisable IF a=2 THEN PRINT number;" Is prime number ";place;" after ";(TIME-T%)/100;" secs" IF a=2 THEN place=place+1 a=0 number=number+1 UNTIL FALSE n 10 100 1000 Prime(n) 29 541 7919 Run time 1.79s 221.61s 27500s In issue 51, Peter Davy wrote an improved version (program Primes2). The major improvement was that, obviously, all prime number over 2 are odd, so only odd numbers were tested. Also, all the factors of a number must be less than its square root. For example, 36 has factors 3, 2, 3, 4, 6, 9, 12 and 18. The unique factors are 2 and 3 which are both less than or equal to the square root, 6. So the next logical improvement is to only test up the the square root. The timing table shown below shows that this improves the speed about 18 fold. REM > Primes2 REM K2K (Peter Davy)'s modification of REM D6E (R. Harker)'s Prime Number program MODE0 PRINT 2;" is prime number ";1:PRINT 3;" is prime number ";2 number%=5 try%=1 order%=2 T%=TIME:REPEAT REPEAT try%=try%+2 IF (number% MOD try%)=0 THEN factorfound=TRUE ELSE factorfound=FALSE UNTIL factorfound OR try%>SQR(number%) IF NOT factorfound AND try%>SQR(number%) THEN order%=order%+1:PRINT number%; " is prime number ";order%;" after ";(TIME-T%)/100;" secs" number%=number%+2:try%=1 UNTIL FALSE n 10 100 1000 Prime(n) 29 541 7919 Run time 0.46s 11.96s 331.74s All factors of a number can be reduced themselves to prime numbers. 36 is 9 * 4 which is 3 * 3 * 2 * 2. The next logical step is to test only with prime numbers up to the square root. This process is called the sieve of Erastores. This would involve storing a list of all the prime numbers encountered so far. However, by examining primes it can be seen that all prime number over 3 are members of the sequence 6n-1, 6n+1. Note that this does not mean that all members of the sequence 6n-1, 6n+1 are primes. All cats have four legs, that dog has four legs, therefore it is a cat. Program Primes3 uses this knowlege to extend the seive process to test number in the sequence 6n-1, 6n+1; most of which are primes; with numbers in the sequence 6n-1, 6n+1. The slight inefficiency cause by the sequence generating a number of non-primes is countered by the simplicity of the formula. The timing table shows that Primes3 is about 20 times faster than R.Harker's original. REM > Primes3 MODE0 REM W79 (J.G.Harston)'s modification to REM K2K (Peter Davy)'s modification to REM D6E (R. Harker)'s Primes program : PRINT 2;" is prime number ";1:PRINT 3;" is prime number ";2 number%=5:number_inc%=2 order%=2 T%=TIME:REPEAT try%=1:try_inc%=4 factorfound=(number% MOD 2)=0 OR (number% MOD 3)=0 IF NOT factorfound THEN REPEAT:try%=try%+try_inc%:try_inc%=6-try_inc%:factorfound= (number% MOD try%)=0:UNTIL factorfound OR try%>SQR(number%) IF NOT (factorfound AND try%<=SQR(number%)) order%=order%+1:PRINT number%; " is prime number ";order%;" after ";(TIME-T%)/100;" secs" number%=number%+number_inc%:number_inc%=6-number_inc% UNTIL FALSE n 10 100 1000 Prime(n) 29 541 7919 Run time 0.44s 10.31s 268.75s Rewriting the program using resident integer variable speeds thing up by another third. The reduction in the length of the names and the fact that they do not need to be searched for in BASIC's variable area but are in fixed locations gives the extra bust of speed. The table shows that Primes4 is 30 times as fast as the original. REM > Primes4 REM 21-09-96 JGH: Fast 6n-1, 6n+1 prime lister : MODE128:PRINT2;" is prime number 1":PRINT3;" is prime number 2" N%=5:I%=2:O%=2:T%=TIME:REPEAT D%=1:A%=4:F%=FALSE:IF(N%MOD3):REPEATD%=D%+A%:A%=6-A%:F%=(N%MODD%)=0:UNTILF% OR D%>=SQRN% IFNOT(F% AND D%<=SQRN%):O%=O%+1:PRINTN%;" is prime number ";O%;" after ";(TIME-T%)/100;" secs" N%=N%+I%:I%=6-I%:UNTILFALSE n 10 100 1000 Prime(n) 29 541 7919 Run time 0.32s 7.43s 175.77s In the eighteenth century the Swiss mathematician Leonhard Eular discovered that the formula: f(n)=n^2+n+41 returns a large number of prime numbers. For all values of n from 0 to 39 the result is a prime. It fails at n=40 as f(n)=40^2+40+41 which is 41^2. However, this formula generates a larger number of primes than any other quadratic (Mathematics: The New Gold Age; Delvin, Keith; p53). REM > Primes5 REM 19-09-96 JGH: Fast quadratic combination : MODE128:O%=3 PRINT2;" is prime 1":PRINT3;" is prime 2" N%=5:I%=2:O%=2:T%=TIME:REPEAT F%=TRUE:IFN%>42 AND N%<1602:N=(SQR(4*N%-163)-1)/2:F%=N<>INTN F%=F%OR((N%MOD3)=0):D%=1:J%=4:IFF%:REPEAT D%=D%+J%:J%=6-J%:F%=(N% MOD D%)=0:UNTILF% OR D%>=SQR(N%) IFNOT(F% AND D%<=SQR(N%)):O%=O%+1:PRINTN%;" is prime ";O%;" after ";(TIME-T%)/100;" secs" N%=N%+I%:I%=6-I%:UNTILFALSE n 10 100 1000 Prime(n) 29 541 7919 Run time 0.32s 7.98s 185.05s Unfortunately, this cannot be used on its own to generate primes as it skips many primes. For example, f(3) is 53, but f(4) is 61, missing out 59. However, re-arranging the equation gives a method of testing for primeness over a limited range. REM > Primes6 REM 21-09-96 JGH: Function to test for primeness : MODE0:N%=5:I%=2:O%=2:T%=TIME:REPEAT IFFNNum_IsPrime(N%):O%=O%+1:PRINTN%;" is prime ";O%;" after ";(TIME-T%)/100;" secs" N%=N%+I%:I%=6-I%:UNTILFALSE : DEFFNNum_IsPrime(N%):LOCAL N,D%,K%,F%:IF(N%MOD2)=0 OR (N%MOD3)=0:=FALSE IFN%>42 AND N%<1602:N=(SQR(4*N%-163)-1)/2:IFN=INTN:=TRUE D%=1:J%=4:REPEAT D%=D%+J%:J%=6-J%:F%=(N% MOD D%)=0:UNTILF% OR D%>=SQR(N%):=NOT(F% AND D%<=SQR(N%)) n 10 100 1000 Prime(n) 29 541 7919 Run time 0.36s 8.63s 195.74s 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 Comparative speed: 1 18 20 30 29 26