Оформил: DeeCo
Автор: http://www.swissdelphicenter.ch

{ *** HTPrimes 
************************************************************** 

      version:  1.1, 13 June 2002 
                (Version history at end of document.) 

      author:   john.a.stockton@btopenworld.com 
                Comments, suggestions, bug reports -- all welcome. 

      purpose:  Suite of functions to assist exploring the prime integers to 
                High(Cardinal), ie. up to 4,294,967,295. Includes an 
efficient 
                primality test combining factorisation and SPRP approaches, 
                implemented (almost) entirely in Pascal. 

      licence:  Free for non-commercial use. 
                A postcard would be very much appreciated from commercial 
                users (royalty cheques also gratefully received - if your 
                conscience gets the better of you). 


*************************************************************************** 

  NOTE: The StrongProbablePrime and WeakProbablePrime functions depend on an 
        ancillary function, MulMod. MulMod is a kludge that uses an 
assembler 
        shortcut. It might not work on all compilers. It has been tested on 
        Delphi 5 (dcc32 v13.0) and Delphi 6 (dcc32 v14.0). 

        The current implementation of MulMod is also prone to trigger #DE 
        exceptions (expressed as EIntOverflow by Delphi), for large values 
of 
        B in calls to StrongProbablePrime and WeakProbablePrime. A slower, 
but 
        safe version of the function is included in a comment. 


*************************************************************************** 
}

 unit HTPrimes;

 interface { *************************************************** Interface 
*** }

 const
   MaxPrime = High(Cardinal) - 4;
     { The highest prime in the range of the Cardinal type. 
      P(203,280,221) = 4,294,967,291. }

 type
   TPrimeFactor = record
     p, q: Cardinal;
   end;

   TPrimeFactorArray = array of TPrimeFactor;
     { Record and dynamic array types used by the GetPrimeFactors function 
      (see below). }

   TPrimeArray = array of Cardinal;
     { Dynamic array type used by the GetGoldbachPairs function (see 
below). }

 { ************************************************* Function Declarations 
*** }

 function IsPrime(Nmbr: Cardinal): Boolean;
   { True if Nmbr is prime, False otherwise. }

 function StrongProbablePrime(N, B: Cardinal): Boolean;
   { True if N is a strong probable prime (SPRP) base B. False if N is 
    composite, or if either N or B < 2. 

    SPRP base B - Write: 
      N - 1 = 2^s * d, with s >= 0 and d odd, 
    then N is probably prime if either: 
      B^d = 1 (mod N), or 
      (B^d)^(2^r) = -1 (mod N), for 0 <= r < s. }

 function WeakProbablePrime(N, B: Cardinal): Boolean;
   { True if N is a probable prime (PRP) base N. False if N is composite, 
    or if either N or B < 2. 

    PRP base B - If B^(N-1) = 1 (mod N), then N is probably prime. }

 function Pseudoprime(N, B: Cardinal): Boolean;
   { True if N is a PRP base B, but composite. }

 function StrongPseudoprime(N, B: Cardinal): Boolean;
   { True if N is a SPRP base B, but composite. }

 function IsComposite(Nmbr: Cardinal): Boolean;
   { True if Nmbr is composite. This isn''t equivalent to not IsPrime(Nmbr), 
    because 0 and 1 are neither prime nor composite. }

 function NextPrime(Nmbr: Cardinal): Cardinal;
   { Returns the first prime greater than Nmbr, or 0 if the search exceeds 
    MaxPrime. }

 function PreviousPrime(Nmbr: Cardinal): Cardinal;
   { Returns the greatest prime less than Nmbr, or 0 if no such prime 
exists. }

 function RandomPrime(Range: Cardinal): Cardinal;
   { Returns a random prime, X, such that 2 <= X < Range for Range > 2. 
    As with the system function Random, initialise the random number 
generator 
    by calling Randomize, or assigning a value to RandSeed. 

    Returns 0 for Range <= 2. }

 function IsJuniorTwin(Nmbr: Cardinal): Boolean;
   { True if Nmbr and Nmbr + 2 are both prime, ie. if they form a 'twin pair'. 
    False if either are composite, or if Nmbr >= MaxPrime. }

 function IsSeniorTwin(Nmbr: Cardinal): Boolean;
   { True if Nmbr and Nmbr - 2 form a twin pair. }

 function GapToNextPrime(Nmbr: Cardinal): Cardinal;
   { Returns the gap to the next prime, where Gap = NextPrime - Nmbr. 
    Returns 0 if the next prime exceeds MaxPrime. }

 function CountPrimesInRange(P, Q: Cardinal): Cardinal;
   { Returns the count of primes, X, that satisfy P <= X <= Q. }

 function RangeContainsPrime(P, Q: Cardinal): Boolean;
   { True if there is at least one prime, X, that satisfies P <= X <= Q. 
    Equivalent to CountPrimesInRange(P, Q) > 0, but faster. }

 function IsPrimeFactor(Nmbr, Fctr: Cardinal): Boolean;
   { True if Fctr is prime, and a factor of Nmbr. }

 function GetPrimeFactors(Nmbr: Cardinal; var Fctrs: TPrimeFactorArray): Boolean;
   { True if Nmbr is greater than one. Fctrs will contain Nmbr''s prime 
factors 
    in TPrime records such that p = the factor, q = exponent that p is 
raised 
    to in factorization of Nmbr. It is not necessary to initialise the array 
    referenced by Fctrs, it will be dynamically sized as required. 

    False if Nmbr is less than two. }

 function CountPrimeFactors(Nmbr: Cardinal): Word;
   { Equivalent to summing exponents after a call to GetPrimeFactors. 
However, 
    this function does not require memory for storing the factors and is 
faster 
    if you only need to know how many factors there are. }

 function CountUniqueFactors(Nmbr: Cardinal): Word;
   { Equivalent to Length(Fctrs) after a call to GetPrimeFactors, but 
faster. }

 function LeastPrimeFactor(Nmbr: Cardinal): Cardinal;
   { Returns the smallest prime factor of Nmbr. 0 if Nmbr is less than 2. 
    If you''re only interested in the least prime factor of Nmbr, this 
function 
    is faster than calling GetPrimeFactors and then examining Fctrs[0]. }

 function GreatestPrimeFactor(Nmbr: Cardinal): Cardinal;
   { Returns the greatest prime factor of Nmbr. 0 if Nmbr is less than 2. 
    As with LeastPrimeFactors, this function is faster than calling 
    GetPrimeFactors and then examining Fctrs[High(Fctrs)]. }

 function GreatestCommonDivisor(P, Q: Cardinal): Cardinal; overload;
   { Returns the greatest common divisor of P and Q, ie. the largest integer 
    that divides them both. }

 function GreatestCommonDivisor(Nmbrs: array of Cardinal): Cardinal;
   overload;
   { Returns the greatest common divisor of the elements of Nmbrs, ie. the 
    largest integer that divides them all. 

    Returns 0 if Nmbrs has less than two elements. }

 function LeastCommonMultiple(P, Q: Cardinal): Int64; overload;
   { Returns the lowest positive integer that both P and Q divide. Returns 0 
if 
    either P or Q are 0. }

 function LeastCommonMultiple(Nmbrs: array of Cardinal): Int64; overload;
   { Returns the lowest positive integer that can be divided by all elements 
of 
    Nmbrs. Returns 0 if any element of Nmbrs is 0, or if Nmbrs has less than 
    2 elements. 

    Note: The LCM of even a small collection of integers can be surprisingly 
          large, and neither version of this function makes any effort to 
deal 
          with overflows beyond High(Int64). Use with caution. }

 function AreRelativelyPrime(P, Q: Cardinal): Boolean;
   { True if P and Q are relatively prime, ie. their greatest common divisor 
is 
    one. }

 function AreMutuallyPrime(Nmbrs: array of Cardinal): Boolean;
   { True if the elements of Nmbrs are mutually relatively prime, ie. if 
there 
    is no integer that divides them all (other than one). 

    False if Nmbrs are not mutually relatively prime, or if Nmbrs has less 
than 
    two elements. }

 function ArePairwisePrime(Nmbrs: array of Cardinal): Boolean;
   { True if the elements of Nmbrs are pairwise relatively prime, ie. if 
every 
    unique pairing of elements in Nmbrs is relatively prime. 

    False if Nmbrs are not pairwise relatively prime, or if Nmbrs has less
 than
     two elements. }

 function GetGoldbachPairs(Nmbr: Cardinal; var Primes: TPrimeArray): Boolean;
   { True if Nmbr is even and greater than 2. The Primes array will be 
populated 
    with the low members of all prime pairs that satisfy P1 + P2 = Nmbr. 

    False if Nmbr is odd or 2. The length of Pairs will be 0. }

 function CountGoldbachPairs(Nmbr: Cardinal): Cardinal;
   { If Nmbr is even and greater than 2, returns the number of prime pairs 
that 
    satisfy P1 + P2 = Nmbr. Equivalent to Length(Primes) after a call to 
    GetGoldbachPairs, but faster and without the memory overhead. 

    0 if Nmbr is odd or 2 (or if Nmbr refutes Goldbach''s Conjecture). }

 function IsGoldbachPair(P, Q: Cardinal): Boolean;
   { True if both P and Q are prime, and P + Q is even and greater than 2. }

 function HasGoldbachPair(Nmbr: Cardinal): Boolean;
   { True if Nmbr is even, greater than two and is the sum of two primes, in 
at 
    least one way. Equivalent to CountGoldbachPairs(Nmbr) > 0, but much 
    faster. }

 implementation { ***************************************** implementation 
*** }

 uses Math;

 {$B-}

 { ****************************************************** function IsPrime 
*** }
 function IsPrime(Nmbr: Cardinal): Boolean;
 const
   PF: array[0..8] of Cardinal = (2, 3, 5, 7, 11, 13, 17, 19, 23);
 var
   i: Cardinal;
 begin
   case Nmbr of
     0..1: IsPrime := False;
     { 0 and 1 are neither prime nor composite. }
     2, 7, 61: IsPrime := True;
     { These three rogues break the test, so deal with them here. }
     else
       begin
     { Do a quick factorisation test using the primes in PF. The list is so 
      short that restricting the factor search to values <= Sqrt(Nmbr) is 
      counter-productive, because of the overhead of calculating the square 
      root. If you extend the factor list, try changing the last inequality 
      in the repeat...until test to (PF[i] > Trunc(Sqrt(Nmbr))). }
         i := 0;
         repeat
           Result := ((Nmbr mod PF[i]) > 0);
           Inc(i);
         until (not Result) or (i > High(PF)) or (PF[i] >= Nmbr);

     { If Result is still True, hammer Nmbr with some SPRP tests. 
      If Nmbr < 4,759,123,141 and is a 2, 7 and 61 SPRP, then Nmbr is 
prime. }
         Result := Result and
           StrongProbablePrime(Nmbr, 2) and
           StrongProbablePrime(Nmbr, 7) and
           StrongProbablePrime(Nmbr, 61);
       end;
   end;
 end; { function IsPrime }

 { ****************************** function MulMod (fast but risky version) 
*** }
 function MulMod(x, y, m: Cardinal): Cardinal;
   { Assumes that on entry: eax = x, edx = y and ecx = m. Self destructs with 
a 
    #DE (divide error) exception when x and y are large. }
 asm
   mul edx
   div ecx
   mov eax, edx
 end; { function MulMod }

 { ******************************* function MulMod (safe but slow version) 
*** 
function MulMod(x, y, m: Cardinal): Cardinal; 
var 
  i, j: Cardinal; 
begin 
  i := x mod m; 
  j := y mod m; 
  asm 
    mov eax, i 
    mul j 
    div m 
    mov i, edx 
  end; 
  MulMod := i; 
end; }

 { ****************************************** function StrongProbablePrime 
*** }
 function StrongProbablePrime(N, B: Cardinal): Boolean;
 var
   d, i, r, s: Cardinal;
 begin
   if (Min(N, B) < 2) then
     StrongProbablePrime := False
   else
   begin
     { Find d and s that satisfy N - 1 = 2^s * d, where d is odd and s >= 0. }
     d := N - 1;
     s := 0;
     while ((d and 1) = 0) do
     begin
       d := d shr 1;
       Inc(s);
     end;

     { Use right-to-left binary exponentiation to Calculate B^d mod N, 
      result in i. }
     i := 1;
     r := B;
     while (d > 1) do
     begin
       if ((d and 1) = 1) then
         i := MulMod(i, r, N);
       d := d shr 1;
       r := MulMod(r, r, N);
     end;
     i := MulMod(i, r, N);

     { If i = 1 or N - 1, then N is a SPRP base B. }
     Result := (i = 1) or (i = N - 1);

     { Otherwise, see if i^(2^r) mod N = N - 1, where 0 <= r < s. }
     r := 0;
     while ((not Result) and (r + 1 <= s)) do
     begin
       i := MulMod(i, i, N);
       Result := (i = N - 1);
       Inc(r);
     end;
   end;
 end; { function StrongProbablePrime }

 { ******************************************** function WeakProbablePrime 
*** }

 function WeakProbablePrime(N, B: Cardinal): Boolean;
 var
   d, i, r: Cardinal;
 begin
   if (Min(N, B) < 2) then
     WeakProbablePrime := False
   else
   begin
     { Use right-to-left binary exponentiation to calculate B^(N-1) mod N, 
      result in i. }
     d := N - 1;
     i := 1;
     r := B;
     while (d > 1) do
     begin
       if ((d and 1) = 1) then
         i := MulMod(i, r, N);
       d := d shr 1;
       r := MulMod(r, r, N);
     end;
     i := MulMod(i, r, N);

     WeakProbablePrime := (i = 1);
   end;
 end; { function WeakProbablePrime }

 { ************************************************** function PseudoPrime 
*** }
 function Pseudoprime(N, B: Cardinal): Boolean;
 begin
   Pseudoprime := WeakProbablePrime(N, B) and IsComposite(N);
 end; { function Pseudoprime }

 { ******************************************** function StrongPseudoPrime 
*** }
 function StrongPseudoprime(N, B: Cardinal): Boolean;
 begin
   StrongPseudoPrime := StrongProbablePrime(N, B) and IsComposite(N);
 end; { function StrongPseudoprime }

 { ************************************************** function IsComposite 
*** }
 function IsComposite(Nmbr: Cardinal): Boolean;
 begin
   { Not quite as trivial as negating IsPrime(Nmbr), but almost. }
   IsComposite := (Nmbr > 3) and (not IsPrime(Nmbr));
 end; { function IsComposite }

 { **************************************************** function NextPrime 
*** }
 function NextPrime(Nmbr: Cardinal): Cardinal;
 begin
   if (Nmbr >= MaxPrime) then
     { Return 0 as an error result if Nmbr >= MaxPrime. }
     NextPrime := 0
   else
   begin
     { Deal with Nmbr < 2 separately, because the main search only considers 
      odd candidates and will therefore miss 2. }
     if (Nmbr < 2) then
       NextPrime := 2
     else
     begin
       { Set Result to the next odd number. }
       Result := Nmbr + 1;
       if ((Result and 1) = 0) then
         Inc(Result);

       { Then test consecutive odd numbers until we find a prime. }
       while (not IsPrime(Result)) do
         Inc(Result, 2);
     end;
   end
 end; { function NextPrime }

 { ************************************************ function PreviousPrime 
*** }
 function PreviousPrime(Nmbr: Cardinal): Cardinal;
 begin
   case Nmbr of
     { 0, 1 and 2 have no previous prime, so return 0 as an error result. }
     0..2: PreviousPrime := 0;
     { As in NextPrime, the main search only considers odd candidates so 
      Start = 3 is a special case and has to be dealt with separately. }
     3: PreviousPrime := 2;
     else
       begin
         { Set Result to the greatest preceding odd number. }
         Result := Nmbr - 1;
         if ((Result and 1) = 0) then
           Dec(Result);

     { Then test consecutive odd numbers (counting down) until we discover a 
      prime. }
         while (not IsPrime(Result)) do
           Dec(Result, 2);
       end;
   end;
 end; { function PreviousPrime }

 { ************************************************** function RandomPrime 
*** }
 function RandomPrime(Range: Cardinal): Cardinal;
 begin
   if (Range <= 2) then
     RandomPrime := 0
   else
     repeat
       Result := NextPrime(Random(Range - 1));
     until (Result > 0) and (Result < Range);
 end; { function RandomPrime }

 { ************************************************* function IsJuniorTwin 
*** }
 function IsJuniorTwin(Nmbr: Cardinal): Boolean;
 begin
   Result := (Nmbr < MaxPrime);
   if Result then
     Result := IsPrime(Nmbr) and IsPrime(Nmbr + 2);
 end; { function IsJuniorTwin }

 { ************************************************* function IsSeniorTwin 
*** }
 function IsSeniorTwin(Nmbr: Cardinal): Boolean;
 begin
   Result := (Nmbr > 2);
   if Result then
     Result := IsPrime(Nmbr) and IsPrime(Nmbr - 2);
 end; { function IsSeniorTwin }

 { *********************************************** function GapToNextPrime 
*** }
 function GapToNextPrime(Nmbr: Cardinal): Cardinal;
 begin
   Result := NextPrime(Nmbr);
   if (Result > 0) then
     Result := Result - Nmbr;
 end; { function GapToNextPrime }

 { ******************************************* function CountPrimesInRange 
*** }
 function CountPrimesInRange(P, Q: Cardinal): Cardinal;
 begin
   Result := 0;
   if (P <= Q) then
   begin
     if (P > 0) then
       Dec(P);
     while (P <= Q) do
     begin
       P := NextPrime(P);
       Inc(Result);
     end;
     Dec(Result);
   end;
 end; { function CountPrimesInRange }

 { ******************************************* function RangeContainsPrime 
*** }
 function RangeContainsPrime(P, Q: Cardinal): Boolean;
 begin
   RangeContainsPrime := False;
   if ((P <= Q) and (P < MaxPrime)) then
   begin
     if (P > 0) then
       Dec(P);
     Result := (NextPrime(P) <= Q);
   end;
 end; { function RangeContainsPrime }

 { ************************************************ function IsPrimeFactor 
*** }
 function IsPrimeFactor(Nmbr, Fctr: Cardinal): Boolean;
 begin
   IsPrimeFactor := ((Nmbr mod Fctr) = 0) and IsPrime(Fctr);
 end; { function IsPrimeFactor }

 { ********************************************** function GetPrimeFactors 
*** }
 function GetPrimeFactors(Nmbr: Cardinal; var Fctrs: TPrimeFactorArray): Boolean;
 var
   i, k: Cardinal;
 begin
   SetLength(Fctrs, 0);
   Result := (Nmbr > 1);
   if Result then
   begin
     k := 0;
     repeat
       { If Nmbr is prime at this stage, it must be the original Nmbr''s greatest 
        prime factor. If Nmbr isn''t prime, search for a prime that divides it, 
        starting from two. }
       if IsPrime(Nmbr) then
         i := Nmbr
       else
       begin
         i := 2;
         while ((Nmbr mod i) > 0) do
           i := NextPrime(i);
       end;

       { Store the found factor. Append it to Fctrs if it''s not already known, 
        otherwise increment the exponent. }
       if (i > k) then
       begin
         k := i;
         SetLength(Fctrs, Length(Fctrs) + 1);
         Fctrs[High(Fctrs)].p := i;
         Fctrs[High(Fctrs)].q := 1;
       end
       else
         Inc(Fctrs[High(Fctrs)].q);

       { Divide Nmbr by the last factor found. If this doesn''t reduce Nmbr to 
        one then go back to find the lowest prime factor of the new value. }
       Nmbr := Nmbr div i;
     until (Nmbr = 1);
   end;
 end; { function GetPrimeFactors }

 { ******************************************** function CountPrimeFactors 
*** }
 function CountPrimeFactors(Nmbr: Cardinal): Word;
 var
   i: Cardinal;
 begin
   Result := 0;
   while (Nmbr > 1) do
   begin
     if IsPrime(Nmbr) then
       i := Nmbr
     else
     begin
       i := 2;
       while ((Nmbr mod i) > 0) do
         i := NextPrime(i);
     end;
     Inc(Result);
     Nmbr := Nmbr div i;
   end;
 end; { function CountPrimeFactors }

 { ******************************************* function CountUniqueFactors 
*** }
 function CountUniqueFactors(Nmbr: Cardinal): Word;
 var
   i, j: Cardinal;
 begin
   Result := 0;
   j := 0;
   while (Nmbr > 1) do
   begin
     if IsPrime(Nmbr) then
     begin
       if (Nmbr <> j) then
         Inc(Result);
       i := Nmbr;
     end
     else
     begin
       i := 2;
       while ((Nmbr mod i) > 0) do
         i := NextPrime(i);
       if (i > j) then
       begin
         j := i;
         Inc(Result);
       end;
     end;
     Nmbr := Nmbr div i;
   end;
 end; { function CountUniqueFactors }

 { ********************************************* function LeastPrimeFactor 
*** }
 function LeastPrimeFactor(Nmbr: Cardinal): Cardinal;
 begin
   if (Nmbr < 2) then
     LeastPrimeFactor := 0
   else if IsPrime(Nmbr) then
     Result := Nmbr
   else
   begin
     Result := 2;
     while ((Nmbr mod Result) > 0) do
       Result := NextPrime(Result);
   end;
 end; { function LeastPrimeFactor }

 { ****************************************** function GreatestPrimeFactor 
*** }
 function GreatestPrimeFactor(Nmbr: Cardinal): Cardinal;
 begin
   if (Nmbr < 2) then
     GreatestPrimeFactor := 0
   else
   begin
     Result := Nmbr;
     while (not IsPrime(Result)) do
       Result := Result div LeastPrimeFactor(Result);
   end;
 end; { function GreatestPrimeFactor }

 { **************************************** function GreatestCommonDivisor 
*** }
 function GreatestCommonDivisor(P, Q: Cardinal): Cardinal;
 var
   i: Cardinal;
 begin
   Result := 0;
   if (Min(P, Q) > 0) then
   begin
     { Use the Euclidean Algorithm to find the GCD of P and Q. Much faster 
than 
      calculating the product of the prime factors that they have in common, 
      but much less interesting. }
     while (Q > 0) do
     begin
       i := P mod Q;
       P := Q;
       Q := i;
     end;
     Result := P;
   end;
 end; { function GreatestCommonDivisor(P, Q: Cardinal) }

 function GreatestCommonDivisor(Nmbrs: array of Cardinal): Cardinal;
 var
   i, j: Cardinal;
 begin
   if (Length(Nmbrs) < 2) then
     Result := 0
   else
   begin
     j      := Low(Nmbrs);
     Result := GreatestCommonDivisor(Nmbrs[j], Nmbrs[j + 1]);
     for i := j + 2 to High(Nmbrs) do
       Result := GreatestCommonDivisor(Result, Nmbrs[i]);
   end;
 end; { function GreatestCommonDivisor(Nmbrs: array of Cardinal) }

 { ****************************************** function LeastCommonMultiple 
*** }
 function LeastCommonMultiple(P, Q: Cardinal): Int64;
 begin
   Result := GreatestCommonDivisor(P, Q);
   if (Result > 0) then
   begin
     { LCM(a, b) = ab/GCD(a, b) }
     P := P div Result;
     LeastCommonMultiple := Int64(P) * Int64(Q);
   end;
 end; { function LeastCommonMultiple(P, Q: Cardinal) }

 function LeastCommonMultiple(Nmbrs: array of Cardinal): Int64;
 var
   i, j: Cardinal;
 begin
   Result := 0;
   if (Length(Nmbrs) > 1) then
   begin
     j := Low(Nmbrs);
     Result := LeastCommonMultiple(Nmbrs[j], Nmbrs[j + 1]);
     for i := j + 2 to High(Nmbrs) do
       Result := LeastCommonMultiple(Result, Nmbrs[i]);
   end;
 end; { LeastCommonMultiple(Nmbrs: array of Cardinal) }

 { ******************************************* function AreRelativelyPrime 
*** }
 function AreRelativelyPrime(P, Q: Cardinal): Boolean;
 begin
   AreRelativelyPrime := (GreatestCommonDivisor(P, Q) = 1);
 end; { function AreRelativelyPrime }

 { ********************************************* function AreMutuallyPrime 
*** }
 function AreMutuallyPrime(Nmbrs: array of Cardinal): Boolean;
 begin
   AreMutuallyPrime := (GreatestCommonDivisor(Nmbrs) = 1)
 end; { function AreMutuallyPrime }

 { ********************************************* function ArePairwisePrime 
*** }
 function ArePairwisePrime(Nmbrs: array of Cardinal): Boolean;
 var
   i, j: Cardinal;
 begin
   Result := (Length(Nmbrs) > 1);
   if Result then
   begin
     i := Low(Nmbrs);
     while (Result and (i <= High(Nmbrs))) do
     begin
       j := i + 1;
       while (Result and (j <= High(Nmbrs))) do
       begin
         Result := (Result and AreRelativelyPrime(Nmbrs[i], Nmbrs[j]));
         Inc(j);
       end;
       Inc(i);
     end;
   end;
 end; { function ArePairwisePrime }

 { ********************************************* function GetGoldbachPairs 
*** }
 function GetGoldbachPairs(Nmbr: Cardinal; var Primes: TPrimeArray): Boolean;
 var
   i: Cardinal;
 begin
   SetLength(Primes, 0);
   Result := (Nmbr > 2) and ((Nmbr and 1) = 0);
   if Result then
   begin
     { Starting from 2, subtract consecutive primes from Nmbr and check if the 
      result is prime. Stop when Nmbr div 2 is exceeded, because then we''re
       just discovering the reversed versions of the pairs we have already. }
     i := 2;
     while (i <= (Nmbr div 2)) do
     begin
       if IsPrime(Nmbr - i) then
       begin
         SetLength(Primes, Length(Primes) + 1);
         Primes[High(Primes)] := i;
       end;
       i := NextPrime(i);
     end;
   end;
 end; { function GetGoldbachPairs }

 { ******************************************* function CountGoldbachPairs 
*** }
 function CountGoldbachPairs(Nmbr: Cardinal): Cardinal;
 var
   i: Cardinal;
 begin
   Result := 0;
   if ((Nmbr > 2) and ((Nmbr and 1) = 0)) then
   begin
     i := 2;
     while (i <= (Nmbr div 2)) do
     begin
       if IsPrime(Nmbr - i) then
         Inc(Result);
       i := NextPrime(i);
     end;
   end;
 end; { function CountGoldbachPairs }

 { ********************************************** function HasGoldbachPair 
*** }
 function HasGoldbachPair(Nmbr: Cardinal): Boolean;
 var
   i: Cardinal;
 begin
   Result := ((Nmbr > 2) and ((Nmbr and 1) = 0));
   if Result then
   begin
     i := 2;
     Result := False;
     while (not Result and (i <= (Nmbr div 2))) do
       if IsPrime(Nmbr - i) then
         Result := True
     else
       i := NextPrime(i);
   end;
 end; { function HasGoldbachPair }

 { *********************************************** function IsGoldbachPair 
*** }
 function IsGoldbachPair(P, Q: Cardinal): Boolean;
 begin
   IsGoldbachPair := (((P + Q) and 1) = 0) and (P + Q > 2) and
     IsPrime(P) and IsPrime(Q);
 end; { function IsGoldbachPair }

 { ******************************************************* Version History 
*** 

  To 
do --------------------------------------------------------------------- 

  Any suggestions? 

  v1.1, 13 June 
2002 -------------------------------------------------------- 

  Additions: 
    function WeakProbablePrime 
    function PseudoPrime 
    function StrongPseudoprime 
    function RandomPrime 

  Changes: 
    01  Moved IsPrime''s declaration to make it more obvious. 
    02  MulMod: Promoted so WeakProbablePrime can use it. 
    03  Mulmod: 'Safe' version re-written to permit SPRP and PRP with large 
        bases. 
    04  IsPrime: Shortened the prime factor list. The shorter list seems to 
        be more efficient in combination with the SPRP tests. 
    05  IsPrime: Scrapped the restriction that (trial factor) <= Sqrt(Nmbr). 
        With the short prime factor list, calculating Trunc(Sqrt(Nmbr)) 
        is an unwarranted overhead. Also replaced the while loop with a 
        repeat until loop, because it looks neater and doesn''t damage 
        performance. 

  v1.0, 12 June 
2002 -------------------------------------------------------- 

  Original release. }

 end. { ********************************************************** THE END 
*** }