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