10 REM Sieve of Erastosthenes 20 MODE 6:PRINT ' TAB(9);"Sieve of Er astosthenes"':@%=10 30 PRINT "This program will calculate the prime numbers up to some given nu mber."' 40 REPEAT 50 INPUT '"Up to what number: ";N 60 IF N<4 OR INT(N)<>N OR N>5000 THE N PRINT '"Be reasonable!" 70 UNTIL N>3 AND INT(N)=N AND N<5001 80 DIM A%(N):CLS:PRINT '"Primes from 2 to ";N 90 FOR I=2 TO N 100 IF A%(I)=0 THEN PROCPrime 110 NEXT 120 PRINT CHR$(7) ''' TAB(10);"Another go? Y or N "; 130 REPEAT:G$=GET$:UNTIL G$="Y" OR G$= "N" 140 IF G$="Y" THEN RUN 150 CLS:PRINT '"Bye for now.":END 160 DEF PROCPrime 170 IF 39-POS<LEN(STR$(I)) THEN PRINT 180 PRINT ;I;" "; 190 FOR J=I TO N STEP I:A%(J)=I:NEXT 200 ENDPROC
Sieve of Erastosthenes This program will calculate the prime numbers up to some given number. Up to what number: ?100 Primes from 2 to 100 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 Another go? Y or N
Note that for large numbers the program starts off slowly; but soon starts printing primes very fast. Printing the primes between 2 and 5000 takes about 43 seconds on the BBC micro and about 59 seconds on the Electron.
The Sieve of Erastosthenes is comceptually easy. It is useful if you want a list of prime numbers. But it is not a very practical method of testing whether a number is prime.
A simple way to check if a number N is prime is to check whether it is divisible by the numbers smaller than N, step by step. The following illustrates this method. A clock is included to indicate how long it takes to test the number for primality.
Listing 9.2
LIST
10 REM Primes Version 1 20 MODE 1:COLOUR 3:PRINT '" Ineffic ient Prime tester Version 1"':@%=10:COLO UR 1 30 REPEAT 40 INPUT '"Number to be tested: ";N% 50 IF N%<4 THEN COLOUR 3:PRINT '"Be reasonable!":COLOUR 1 60 UNTIL N%>3 70 COLOUR 3:PRINT '"The number being tested is ";N% 80 TIME=0:A$="":X%=N%-1 90 FOR I%=2 TO X% 100 IF N% MOD I% = 0 THEN A$="not ":I %=X% 110 NEXT 120 COLOUR 2:PRINT '"The number is ";A $;"prime." 130 PRINT '"Time taken to test number ";INT(TIME/100+0.5);" seconds." 140 COLOUR 3:PRINT CHR$(7) '' TAB(10); "Another go? Y or N "; 150 REPEAT:G$=GET$:UNTIL G$="Y" OR G$= "N" 160 IF G$="Y" THEN RUN 170 CLS:PRINT '"Bye for now.":END
RUN
Inefficient Prime tester Version 1 Number to be tested: ?9001 The number being tested is 9001 The number is prime. Time taken to test number 20 seconds. Another go? Y or N
The program works but it really is inefficient and slow. For instance to test the primality of 9001 takes 20 seconds on the BBC micro and about 56 seconds on the Electron. To test the primality of a larger number such as 90001 takes about 205 seconds on the BBC and about 590 on the Electron. A little thought will produce enormous benefits.
For a start we needn't bother using all the numbers between 2 and N. We need only use the numbers between 2 and INT(SQR(N)). Because if M is an integer which divides N and is greater than INT(SQR(N)) then N/M is an integer that divides N and is smaller than INT(SQR(N)). This simple addition is included in version two below.
Listing 9.3
LIST
10 REM Primes Version 2 20 MODE 1:COLOUR 3:PRINT '" Ineffic ient Prime tester Version 2"':@%=10:COLO UR 1 30 REPEAT 40 INPUT '"Number to be tested: ";N% 50 IF N%<4 THEN COLOUR 3:PRINT '"Be reasonable!":COLOUR 1 60 UNTIL N%>3 70 COLOUR 3:PRINT '"The number being tested is ";N% 80 TIME=0:A$="":X%=SQR(N%) 90 FOR I%=2 TO X% 100 IF N% MOD I% = 0 THEN A$="not ":I %=X% 110 NEXT 120 COLOUR 2:PRINT '"The number is ";A $;"prime." 130 PRINT '"Time taken to test number ";INT(TIME/100+0.5);" seconds." 140 COLOUR 3:PRINT CHR$(7) '' TAB(10); "Another go? Y or N "; 150 REPEAT:G$=GET$:UNTIL G$="Y" OR G$= "N" 160 IF G$="Y" THEN RUN 170 CLS:PRINT '"Bye for now.":END
RUN
Inefficient Prime tester Version 2 Number to be tested: ?90001 The number being tested is 90001 The number is prime. Time taken to test number 1 seconds. Another go? Y or N Inefficient Prime tester Version 2 Number to be tested: ?9001 The number being tested is 9001 The number is prime. Time taken to test number 0 seconds. Another go? Y or N
This version works considerably faster. Testing the number 9001 now takes just over 1 second. A larger number, such as 987654323 takes about 87 seconds on the BBC and about 252 seconds on the Electron. Don't attempt to test such a large number with version 1 - unless you enjoy looking mindlessly at your television screen.
Remark. The addition incorporated in Prime Version 2 could also be incorporated in the program Erastosthenes to produce Erastosthenes 2. The result would improve the speed.
We can further improve the program Primes Version 2 by taking a tip from the Sieve of Erastosthenes. We can miss out all even numbers bigger than 2 and miss out every third number after 3. These suggestions have been added to the third version of the program.
Listing 9.4
LIST
10 REM Primes Version 3 20 MODE 1:COLOUR 3:PRINT ' TAB(9);"Pr ime tester Version 3"':@%=10:COLOUR 1 30 REPEAT 40 INPUT '"Number to be tested: ";N% 50 IF N%<8 THEN COLOUR 3:PRINT '"Be reasonable!":COLOUR 1 60 UNTIL N%>7 70 COLOUR 3:PRINT '"The number being tested is ";N% 80 TIME=0:A$="":X%=SQR(N%) 90 IF N% MOD 2 = 0 OR N% MOD 3 = 0 TH EN A$="not " ELSE PROCOther 100 COLOUR 2:PRINT '"The number is ";A $;"prime." 110 PRINT '"Time taken to test number ";INT(TIME/100+0.5);" seconds." 120 COLOUR 3:PRINT CHR$(7) '' TAB(10); "Another go? Y or N "; 130 REPEAT:G$=GET$:UNTIL G$="Y" OR G$= "N" 140 IF G$="Y" THEN RUN 150 CLS:PRINT '"Bye for now.":END 160 DEF PROCOther 170 FOR I%=5 TO X% STEP 6 180 IF N% MOD I% = 0 THEN A$="not ":I% =X%190 IF N% MOD (I%+2) = 0 THEN A$="not ":I%=X% 200 NEXT 210 ENDPROC
RUN
Prime tester Version 3 Number to be tested: ?90001 The number being tested is 90001 The number is prime. Time taken to test number 0 seconds. Another go? Y or N
This speeds up the program a little. Testing the prime 987654323 now takes about 30 seconds on the BBC and about 86 seconds on the Electron. The largest integer that your computer can store is 2147483647, is it a prime number?
A composite number may be written as a product of prime numbers. For instance
6 = 3 * 2
24 = 3 * 2 * 2 * 2
81018001 = 9001 * 9001
and so on. The various primes that occur are called factors of the number. By adding a few extra lines to version 3 we can have all the factors of a number displayed.
The next program prints out the factors of a number.
Listing 9.5
LIST
10 REM Prime factors 20 MODE 1:COLOUR 3:PRINT ' TAB(13);"P rime factors"':@%=10:COLOUR 1 30 REPEAT 40 INPUT '"Number to be tested: ";N% 50 IF N%<8 THEN COLOUR 3:PRINT '"Be reasonable!":COLOUR 1 60 UNTIL N%>7 70 COLOUR 3:PRINT '"The factors of "; N%;" are:"':COLOUR 2 80 TIME=0:A$="":X%=SQR(N%):S%=X%+1 90 IF N% MOD 2 = 0 THEN PROCFactor(2) 100 IF N% MOD 3 = 0 THEN PROCFactor(3) 110 FOR I%=5 TO X% STEP 6 120 IF N% MOD I% = 0 THEN PROCFactor( I%) 130 IF N% MOD (I%+2) = 0 THEN PROCFac tor(I%+2) 140 IF I%>S% THEN I%=X% 150 NEXT 160 IF 39-POS<LEN(STR$(N%)) THEN PRINT 170 IF N%>1 THEN PRINT ;N% 180 COLOUR 1:PRINT ''"Time to find fac tors: ";INT(TIME/100+0.5);" seconds." 190 COLOUR 3:PRINT CHR$(7) '' TAB(10); "Another go? Y or N "; 200 REPEAT:G$=GET$:UNTIL G$="Y" OR G$= "N" 210 IF G$="Y" THEN RUN 220 CLS:PRINT '"Bye for now.":END 230 DEF PROCFactor(J%) 240 REPEAT 250 IF 39-POS<LEN(STR$(J%)) THEN PRIN T 260 PRINT ;J%;" ";:N%=N% DIV J% 270 UNTIL N% MOD J% > 0 280 S%=SQR(N%)+1 290 ENDPROC
RUN
Prime factors Number to be tested: ?162036002 The factors of 162036002 are: 2 9001 9001 Time to find factors: 9 seconds. Another go? Y or N Prime factors Number to be tested: ?1022048 The factors of 2048 are: 2 2 2 2 2 2 2 2 2 2 2 Time to find factors: 0 seconds. Another go? Y or N
Large primes
Finding large prime numbers is a pastime for some. The largest known primes are usually Mersenne primes. Numbers of the form
2P - 1
are called Mersenne numbers because of the French monk, Father Marin Mersennes, who made some suggestions concerning the primality of such numbers.
Many early writers believed that the Mersenne numbers are prime if the exponent M is prime. For the first few cases this is indeed true.
22 - 1 = 3
23 - 1 = 7
25 - 1 = 31
27 - 1 = 127
But for P = 11 the number is not prime.
211 - 1 = 2047 = 23 * 89
For at least the last 100 years the world's largest known prime has always been a Mersenne prime (except for a short period in 1951). Prior to January 1983, the largest known prime was
244497 - 1
a Mersenne prime discovered by David Slowinski in April 1979. This was the 27th Mersenne prime known. In about 1983 Slowinski found the much larger prime
286243 - 1
which is a Mersenne prime containing 25962 decimal digits. (It would fill several pages of this book to write out the number completely.) The computer that Slowinski used was a Cray-1. This is an amazingly fast machine, even so it took 1 hour 36 minutes and 22 seconds of computer time to perform the test.
An arbitrary number of the size 2286243 - 1 would be impossible to test for primality using any of the programs in the previous section. But there are special techniques for Mersenne numbers. These were developed during the last 100 years. To test a number of the form N = 2P - 1 for primality one defines a sequence as follows:
U1 = 4
U2 = (U1*U1 - 2) MOD N
· ·
· ·
· ·
Up-1 = (Up-2*Up-2 - 2) MOD N
Then N is prime if and only if UP-1 = 0. This means that in order to test for N being prime we need to perform approximately P simple calculations. Performing 86243 operations would be quite simple for your computer if the numbers involved were small. But here one is dealing with numbers that have 86243 binary digits and your computer stores 32 binary digits. There are methods available to get around this problem but we won't go into the details here. See the next chapter.
Listing 9.6
LIST
10 REM Probabilisic primality test 20 MODE 1:COLOUR 3:PRINT ' TAB(6);"Pr obabilistic primality test"':@%=10:COLOU R 1 30 REPEAT 40 INPUT '"Number to be tested: ";N% 50 IF N%<8 THEN COLOUR 3:PRINT '"Be reasonable!":COLOUR 1 60 UNTIL N%>7 70 REPEAT 80 REM Factorise N%-1 = (2^T%)*X% 90 T%=0:X%=N%-1 100 REPEAT 110 IF X% MOD 2 = 0 THEN T%=T%+1:X%= X% DIV 2 120 UNTIL X% MOD 2 = 1 130 REM Select a base B 140 I=RND(-TIME):B%=RND(N%-3) + 1 150 COLOUR 1:PRINT '"Test using base ";B%:COLOUR 2 160 REM Raise B% to power X% 170 P%=1 180 REPEAT 190 IF X% MOD 2 = 1 THEN P%=P%*B%:P% =P% MOD N% 200 B%=B%*B%:B%=B% MOD N%:X%=X% DIV 2 210 UNTIL X%=0 220 REM Check B%^X% 230 REPEAT 240 TEST=-1 250 IF P%=1 OR P%=N%-1 THEN TEST=0 260 IF TEST AND T%<2 THEN PRINT '"Nu mber is NOT prime.":TEST=1:G$="N" 270 IF TEST THEN P%=P%*P%:P%=P% MOD N%:T%=T%-1 280 UNTIL NOT TEST 290 IF TEST=0 PRINT '"Number is PROBA BLY prime."''"Do you want to test the nu mber again with another base? Y or N" ' 300 IF TEST=0:REPEAT:G$=GET$:UNTIL G$ ="Y" OR G$="N" 310 UNTIL G$="N" 320 COLOUR 3:PRINT CHR$(7) '' TAB(10); "Another go? Y or N "; 330 REPEAT:G$=GET$:UNTIL G$="Y" OR G$= "N" 340 IF G$="Y" THEN RUN 350 CLS:PRINT '"Bye for now.":END
RUN
Probabilistic primality test Number to be tested: ?341 Test using base 29 Number is PROBABLY prime. Do you want to test the number again with another base? Y or N Test using base 31 Number is NOT prime. Another go? Y or N
Try numbers such as 341 and 561, neither of which is prime. If you find that they are a pseudoprime to some base then try another base.