5.9 An Extended Example (Finding Prime Numbers)

In previous discussion on the greatest common divisor, it was mentioned that the GCD could be found by factoring the two numbers into products of primes, that is numbers with exactly two different positive divisors. The sequence of prime numbers

2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, and so on

has long fascinated mathematicians because of their usefulness and their unique properties.

Although this is a sequence, because it can be numbered off in order with the natural numbers, it is possible to prove that there does not exist a formula that can be used to compute the ith term of the sequence. In order to obtain a list of primes, it is therefore necessary to resort to computational methods that produce them one at a time. A possible algorithm for determining if a number is prime is:

let a trial divisor have the value 2
set a boolean Prime to true
while Prime or the trial divisor  the square root of the test number
  divide the number being tested by the trial divisor
  if there is no remainder, set Prime to false
  increment the trial divisor

Once the square root of the test number has been reached, there is no point in checking for more divisors, as there can be a divisor larger than this only if there was also one smaller. Here is some Modula-2 code to implement and test this algorithm:

MODULE CheckPrime;

(* Written by R.J. Sutcliffe *)
(* using ISO Modula-2 *)
(* last revision 1991 03 01 *)

FROM STextIO IMPORT
  WriteString, WriteLn, ReadChar, SkipLine;
FROM SWholeIO IMPORT
  ReadCard, WriteCard;
FROM RealMath IMPORT
  sqrt;
  
VAR
  theNum : CARDINAL;
  ans : CHAR;

PROCEDURE IsPrime (testNum : CARDINAL) : BOOLEAN;

VAR
  prime : BOOLEAN;
  trialDiv, stop : CARDINAL;

BEGIN
  IF (testNum = 0) OR (testNum = 1)
    THEN
      prime := FALSE
    ELSE
      prime := TRUE;
      trialDiv := 2;
      stop := TRUNC (sqrt (FLOAT (testNum)));
      WHILE (trialDiv <= stop) AND prime
        DO
        IF testNum MOD trialDiv = 0
          THEN
            prime := FALSE
          END;
        INC (trialDiv);
      END; (* while *)
    END;  (* if *)
  RETURN prime;
END IsPrime;

BEGIN (* main *)
  REPEAT
    WriteString ("The number you type will be checked ");
    WriteString ("to see if it is prime ==>");
    ReadCard (theNum);
    SkipLine;
    WriteLn;
    WriteString ("the number ");
    WriteCard (theNum, 0);
    WriteString (" is ");
    IF NOT IsPrime (theNum)
      THEN
        WriteString ("not");
      END;
    WriteString (" prime.");
    WriteLn;
    WriteString ("Do another? Y/N  ");
    ReadChar (ans);
    SkipLine;
  UNTIL CAP (ans) = "N";
END CheckPrime.

This code will do, if all that is required is to test a single number. However, it has two drawbacks. First, it is inefficient. It is not necessary to use every number from 2 through the square root of the number being tested as a trial divisor; only the primes in that range. If a way can be found to do this, the number of computations is greatly reduced. Second, no record or list of primes is kept by this program. A number that is given to it twice must be tested twice. Although this algorithm cannot be replaced by a simple call to a closed formula, it can be improved upon considerably.

A much more efficient method of finding primes is often given in elementary mathematics textbooks. It is called the Sieve of Eratosthenes (after the Greek mathematician) and when done by hand it works like this:

1. Write down the numbers to be tested in convenient rows

2. Cross out 1 and circle the number 2--it is prime.

3. Cross out all multiples of 2; they cannot be prime.

4. The next number not crossed out is a prime; circle it and cross out all its multiples that are not already crossed out. When looking for multiples, start with the square of the number just circled; any lesser multiples have already been crossed out.

5. Repeat steps 3 and 4 until all the numbers less than the square root of the largest number that was written down are either circled (prime) or crossed out (composite).

6. All remaining numbers of those written down are prime, and may now be circled.

A count of the circled numbers reveals that there are 25 primes less than 100. This hand method works very well for checking up to a few hundred numbers, but becomes rather tedious beyond that point. The calculations are, however, very simple to computerize. Using the numbers above for the steps:

1. Create an array of booleans and set them all to true at first. (true = prime)
2. Set array element 1 to false. Now 2 is prime.
3. Set the values whose index in the array is a multiple of the last prime found to false.
4. The next index where the array holds the value true is the next prime.
5. Repeat steps 3 and 4 until the last prime found is greater than the square root of the largest number in the array.

The program that follows tested the numbers to 1000 and printed out the primes it found, ten to a line.

Code:

MODULE Sieve;

(* Written by R.J. Sutcliffe *)
(* using ISO Modula-2 *)
(* last revision 1993 03 01 *)

FROM STextIO IMPORT
  WriteString, WriteLn, ReadChar;
FROM SWholeIO IMPORT
  WriteCard;
FROM RealMath IMPORT
  sqrt;

CONST
  size = 1000;

VAR
  primeFlags : ARRAY [1 .. size] OF BOOLEAN;
  numFound,
  numToTest, (* the candidate for primality *)
  count, lineCount,        (* Loop counters *)
  mul,   (* for multiples *)
  max    (* largest to test *) : CARDINAL;
  key : CHAR;

BEGIN
  numFound := 0;
  max := TRUNC (sqrt (FLOAT(size))); 
  primeFlags [1] := FALSE;
  FOR count := 2 TO size
    DO    (* Initialize variables *)
      primeFlags [count] := TRUE
    END;

  (* Find primes *)
  FOR numToTest := 2 TO max
    DO
      IF primeFlags [numToTest]
        THEN  (* got the next prime here *)
          INC (numFound);  (* count them as we go *)
          mul := numToTest * numToTest;  (* start at its square *)
          WHILE mul <= size
            DO
              primeFlags [mul] := FALSE;   (* Cross out *)
              INC (mul, numToTest)    (* Do all multiples *)
            END  (* while *)
        END  (* if *)
    END;  (* for *)
  
  (* now, count the rest *)
  FOR numToTest := max + 1 TO size
    DO
      IF primeFlags [numToTest]
        THEN  (* got the next prime here *)
          INC (numFound);  (* count them as we go *)
        END  (* if *)
    END;  (* for *)

  (* Now print them out 10 to a line *)

  lineCount := 0;
  FOR count := 2 TO size
    DO
      IF primeFlags [count]
        THEN
          WriteCard (count, 6);
          INC (lineCount)
        END;
      IF lineCount = 10 
        THEN
          lineCount := 0;  (* reset line count *)
          WriteLn;   (* and start new line *)
        END
    END;   (* for *)
  WriteLn;
  WriteString (" ... a total of ");
  WriteCard (numFound, 6);
  WriteString (" primes ");
  WriteString ("Press a key to continue ==> ");
  ReadChar (key);
END Sieve.

Results:

     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   101   103   107   109   113
   127   131   137   139   149   151   157   163   167   173
   179   181   191   193   197   199   211   223   227   229
   233   239   241   251   257   263   269   271   277   281
   283   293   307   311   313   317   331   337   347   349
   353   359   367   373   379   383   389   397   401   409
   419   421   431   433   439   443   449   457   461   463
   467   479   487   491   499   503   509   521   523   541
   547   557   563   569   571   577   587   593   599   601
   607   613   617   619   631   641   643   647   653   659
   661   673   677   683   691   701   709   719   727   733
   739   743   751   757   761   769   773   787   797   809
   811   821   823   827   829   839   853   857   859   863
   877   881   883   887   907   911   919   929   937   941
   947   953   967   971   977   983   991   997
 ... a total of    168 primes 

NOTES: 1. The value of size could have been set much higher, but there was no need to print any more primes here.

2. In some implementations of Modula-2, the maximum array index range is limited to a number much less than MAX (CARDINAL), and if size were to be greater than this, the program would not compile. This is due to limitations on the total size of a variable.

3. A file of primes could be kept on hand for use in other problems.


Contents