PrimeInt.Pas Program Listing

{ Start of file PrimeInt.Pas ***********************************************}

{$A+,B-,D+,E+,F-,I+,L+,N+,O-,R-,S+,V+}
{$M 16384,0,655360}

unit PrimeInt; { Converted from Apple II+ Pascal to Turbo Pascal 5.0 }
               { and changed from SextPrime to PrimeInt }
               { 90/07/08 by Harry Smith, Saratoga CA }

{ Algorithm 357 - Collected Algorithms from ACM }
{ An Efficient Prime Number Generator using the sieve of Eratosthenes }

{ Driver can generates all the primes <= 2147107031 = 2^31 - 376617 }
{
  Name    = 'PrimeInt - Generates Primes';
  Version = 'Version 1.00, last revised: 90/09/15 0600 hours';
  Author  = 'Copyright (c) 1981-1990 by author: Harry J. Smith,';
  Address = '19628 Via Monte Dr., Saratoga CA 95070.  All rights reserved.';
}
{***************************************************************************}

{ Developed in TURBO Pascal 5.0 }

interface

const
  JLim = 6500;       { Size of working array }
  MaxT = 2147107031; { Approximately the max prime that can be generated }

type
  PrimeType  = LongInt;
  PrimeStore = array[1..JLim] of PrimeType; { Working storage for primes }

var
  IP : PrimeStore; { Working storage for primes }

{ The following PROCEDURE(s) can be called externally: }

  function NPrime(M : Integer) : Integer;
  { M = Number of new primes desired this call, must be >= 2 on 1st call }
  {     keep M <= JLim / 2 for speed, must be < JLim. }
  { The current length of the tables in array IQ and JQ is returned }
  {   This is approx. the number of primes < Sqrt(current prime) }
  {   It must not reach QS = 4793 because the 4793rd prime = 46349 and }
  {   this squared is greater than 2^32 - 1, the max LongInt. }
  {   If no more primes can be generated, -4793 is returned. }

{************************************************************}

implementation

const
  QS = 4793; { Size of IQ and JQ arrays, must be large enough to hold }
             {   all primes <= Sqrt(MaxT) }

type
  WordStore  = array[1..QS] of Word;      { Secondary table of primes }
  LongStore  = array[1..QS] of PrimeType; { Secondary table of primes squared
                                           etc. = next composite to cross off }
var
  IQ : LongStore; { Secondary table of primes squared etc. }
  JQ : WordStore; { Secondary table of primes < Sqrt(current prime) }
                  {   2, 3, 5, 7, ... 46349 (the 4793rd prime) }

{--------------------------------------}
function NPrime(M : Integer) : Integer;
label 10, 20, 30, 40, 50, 60;
const                  { Typed constants are Static variables }
  IJ  : Word      = 0; { Number of Primes in secondary prime table }
  IK  : Word      = 0; { Number of Primes for testing secondary primes }
  INC : Word      = 0; { = 2, 4, 2, 4, ... }
  J   : Word      = 0; { Array index, sieve index on exit }
  NJ  : Word      = 0; { Test number for secondary table of primes }
  NI  : PrimeType = 0; { Test number - 2 * J for primes }

{ note:  IQ[1] is used as a Static variable to save highest number in sieve }

var
  I   : Word; { Array index }
  JQI : Word; { JQ[I] to store in IP[J] }
  K   : Word; { Number of Primes in primary prime table }

begin
  K:= 0;  if (IP[1] >= 0) and (IJ < QS) then  GoTo 60;
  if IP[1] < 0 then begin
    { Set initial conditions }
    IP[1]:=  2;  JQ[1]:= 2;  IK:=    2;  INC:= 2;
    IQ[2]:=  9;  JQ[2]:= 3;  IQ[1]:= 3;  IJ:=  3;
    IQ[3]:= 25;  JQ[3]:= 5;  NJ:=    5;  K:=   1;
  end;
  if IJ >= QS then begin
    for I:= 1 to M do  IP[I]:= 0; { Clear primes that cannot be generated }
    NPrime:= -QS;  Exit;          { Indicate no more primes can be generated }
  end;
  { Prepare to delete a sequence of composite numbers }
10:
  J:= K + 1;  NI:= IQ[1] - J - J;
  IQ[1]:= NI + JLim + JLim;
  for I:= J to JLim do  IP[I]:= 0;
20:
  I:= IJ;  if IQ[IJ] >= IQ[1] then  GoTo 50;
  { Extend the list of primes in array JQ counting so as to omit
    multiples of 2 and 3 }
30:
  NJ:= NJ + INC;  INC:= 6 - INC;
  if Sqr( JQ[IK+1]) <= NJ then  IK:= IK + 1;
  for J:= 3 to IK do
    if (NJ DIV JQ[J]) * JQ[J] = NJ then  GoTo 30;
  IJ:= IJ + 1;  JQ[IJ]:= NJ;  IQ[IJ]:= PrimeType(NJ) * NJ;
  if IJ >= QS then begin
    { Clear primes that cannot be generated }
    for I:= (K + 1) to M do  IP[I]:= 0;
    NPrime:= -QS;  Exit;     { Indicate no more primes can be generated }
  end;
  GoTo 20;
  { If J+J+NI is composite, enter its smallest prime factor in IP[J].
    If J+J+NI is prime, then set IP[J] to 0 }
40:
  IP[J]:= JQI;  J:= J + JQI;
  if J < JLim then  GoTo 40;
  IQ[I]:= NI + J + J;
50:
  I:= I - 1;  JQI:= JQ[I];  J:= (IQ[I] - NI) DIV 2;
  if J < JLim then  GoTo 40;
  if I <> 1 then  GoTo 50;  J:= K;
  { Pack the next M primes in IP[1], ..., IP[M] }
60:
  J:= J + 1;  if IP[J] <> 0 then  GoTo 60; { IP[JLim] is always 0 }
  if J = JLim then  GoTo 10;  
  K:= K + 1;  IP[K]:= NI + J + J;
  if K <> M then  GoTo 60;
  NPrime:= IJ
end; { NPrime }

{ PrimeInt Unit initialization }
begin
  IP[1]:= -1;
end.

{ End of file PrimeInt.Pas *************************************************}