首页  编辑  

determine if a number is prime, quickly (2)

Tags: /超级猛料/Alogrith.算法和数据结构/数值运算/   Date Created:

...determine if a number is prime, quickly (2)?  

Autor: John A. Stockton  

[Write new comment]

[ Print tip ]    

Tip Rating (3):  

     

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

*** }