Tips (1541)

Tip Rating (11):

{ *** 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
*** }