MathgeomGLS

1
{---------------------------------------------------------------------------}
2
{                                                                           }
3
{ File:       Velthuis.BigIntegers.Primes.pas                               }
4
{ Function:   Prime functions for BigIntegers                               }
5
{ Language:   Delphi version XE2 or later                                   }
6
{ Author:     Rudy Velthuis                                                 }
7
{ Copyright:  (c) 2017 Rudy Velthuis                                        }
8
{ Notes:      See http://rvelthuis.de/programs/bigintegers.html             }
9
{             See https://github.com/rvelthuis/BigNumbers                   }
10
{                                                                           }
11
{ License:    Redistribution and use in source and binary forms, with or    }
12
{             without modification, are permitted provided that the         }
13
{             following conditions are met:                                 }
14
{                                                                           }
15
{             * Redistributions of source code must retain the above        }
16
{               copyright notice, this list of conditions and the following }
17
{               disclaimer.                                                 }
18
{             * Redistributions in binary form must reproduce the above     }
19
{               copyright notice, this list of conditions and the following }
20
{               disclaimer in the documentation and/or other materials      }
21
{               provided with the distribution.                             }
22
{                                                                           }
23
{ Disclaimer: THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER "AS IS"     }
24
{             AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT     }
25
{             LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND     }
26
{             FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO        }
27
{             EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE     }
28
{             FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,     }
29
{             OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,      }
30
{             PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,     }
31
{             DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED    }
32
{             AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT   }
33
{             LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)        }
34
{             ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF   }
35
{             ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                    }
36
{                                                                           }
37
{---------------------------------------------------------------------------}
38

39
unit Velthuis.BigIntegers.Primes;
40

41
interface
42

43
uses
44
  Velthuis.BigIntegers;
45

46
type
47
  /// <summary>Type indicating primality of a number.</summary>
48
  TPrimality = (
49
    /// <summary>Number is definitely composite.</summary>
50
    primComposite,
51
    /// <summary>Number is probably prime.</summary>
52
    primProbablyPrime,
53
    /// <summary>Number is definitely prime.</summary>
54
    primPrime
55
  );
56

57
/// <summary>Detects if N is probably prime according to a Miller-Rabin test.</summary>
58
/// <param name="N">The number to be tested</param>
59
/// <param name="Precision">Determines the probability of the test, which is 1.0 - 0.25^Precision</param>
60
/// <returns>True if N is prime with the given precision, False if N is definitely composite.</returns>
61
function IsProbablePrime(const N: BigInteger; Precision: Integer): Boolean;
62

63
/// <summary>Detects if N is prime, probably prime or composite. Deterministically correct for N &lt; 341,550,071,728,321, otherwise
64
///  with a probability determined by the Precision parameter</summary>
65
/// <param name="N">The number to be tested</param>
66
/// <param name="Precision">For values >= 341,550,071,728,321, IsProbablyPrime(N, Precision) is called</param>
67
/// <returns>Returns primComposite if N is definitely composite, primProbablyPrime if N is probably prime and
68
///  primPrime if N is definitely prime.</returns>
69
function IsPrime(const N: BigInteger; Precision: Integer): TPrimality;
70

71
/// <summary>Returns a random probably prime number.</summary>
72
/// <param name="NumBits">Maximum number of bits of th random number</param>
73
/// <param name="Precision">Precision to be used for IsProbablePrime</param>
74
/// <returns> a random prime number that is probably prime with the given precision.</returns>
75
function RandomProbablePrime(NumBits: Integer; Precision: Integer): BigInteger;
76

77
/// <summary>Returns a probable prime >= N.</summary>
78
/// <param name="N">Number to start with</param>
79
/// <param name="Precision">Precision for primality test</param>
80
/// <returns>A number >= N that is probably prime with the given precision.</returns>
81
function NextProbablePrime(const N: BigInteger; Precision: Integer): BigInteger;
82

83
/// <summary>Checks if the (usually randomly chosen) number A witnesses N's compositeness.</summary>
84
/// <param name="A">Value to test with</param>
85
/// <param name="N">Number to be tested as composite</param>
86
///
87
/// <remarks>See https://en.wikipedia.org/wiki/Miller-Rabin_primality_test#Algorithm_and_running_time</remarks>
88
function IsWitness(const A, N: BigInteger): Boolean;
89

90
function IsComposite(const A, D, N: BigInteger; S: Integer): Boolean;
91

92
implementation
93

94
// See https://en.wikipedia.org/wiki/Miller-Rabin_primality_test.
95

96
uses
97
  Velthuis.RandomNumbers, System.SysUtils;
98

99
var
100
  Two: BigInteger;
101
  DeterminicityThreshold: BigInteger;
102
  CertainlyComposite: BigInteger;
103
  Random: IRandom;
104

105
// Rabin-Miller test, deterministically correct for N < 341,550,071,728,321.
106
// http://rosettacode.org/wiki/Miller-Rabin_primality_test#Python:_Proved_correct_up_to_large_N
107
// If you want to improve upon this:
108
// https://en.wikipedia.org/wiki/Miller-Rabin_primality_test#Deterministic_variants_of_the_test
109
function IsPrime(const N: BigInteger; Precision: Integer): TPrimality;
110
var
111
  R: Integer;
112
  I: Integer;
113
  PrimesToTest: Integer;
114
  D: BigInteger;
115
  N64: UInt64;
116
const
117
  CPrimesToTest: array[0..6] of Integer = (2, 3, 5, 7, 11, 13, 17);
118
  CProbabilityResults: array[Boolean] of TPrimality = (primComposite, primProbablyPrime);
119
begin
120
  if BigInteger.Compare(N, DeterminicityThreshold) > 0 then
121
    Exit(CProbabilityResults[IsProbablePrime(N, Precision)]);
122

123
  if BigInteger.Compare(N, CertainlyComposite) = 0 then
124
    Exit(primComposite);
125

126
  N64 := UInt64(N);
127
  if N64 > Int64(3474749660383) then
128
    PrimesToTest := 7
129
  else if N64 > Int64(2152302898747) then
130
    PrimesToTest := 6
131
  else if N64 > Int64(118670087467) then
132
    PrimesToTest := 5
133
  else if N64 > Int64(25326001) then
134
    PrimesToTest := 4
135
  else if N64 > Int64(1373653) then
136
    PrimesToTest := 3
137
  else
138
    PrimesToTest := 2;
139

140
  D := N - BigInteger.One;
141
  R := 0;
142

143
  while D.IsEven do
144
  begin
145
    D := D shr 1;
146
    Inc(R);
147
  end;
148

149
  for I := 0 to PrimesToTest - 1 do
150
    if IsComposite(CPrimesToTest[I], D, N, R) then
151
      Exit(primComposite);
152

153
  Result := primPrime;
154
end;
155

156
// Check if N is probably prime with a probability of at least 1.0 - 0.25^Precision
157
function IsProbablePrime(const N: BigInteger; Precision: Integer): Boolean;
158
var
159
  I: Integer;
160
  A, NLessOne: BigInteger;
161
begin
162
  if (N = 2) or (N = 3) then
163
    Exit(True)
164
  else if (N = BigInteger.Zero) or (N = BigInteger.One) or (N.IsEven) then
165
    Exit(False);
166

167
  NLessOne := N;
168
  Dec(NLessOne);
169

170
  for I := 1 to Precision do
171
  begin
172
    repeat
173
      A := BigInteger.Create(N.BitLength, Random);
174
    until (A > BigInteger.One) and (A < NLessOne);
175

176
    if IsWitness(A, N) then
177
      Exit(False);
178
  end;
179
  Result := True;
180
end;
181

182
function RandomProbablePrime(NumBits: Integer; Precision: Integer): BigInteger;
183
begin
184
  repeat
185
    Result := BigInteger.Create(NumBits, Random);
186
  until IsProbablePrime(Result, Precision);
187
end;
188

189
// Next probable prime >= N.
190
function NextProbablePrime(const N: BigInteger; Precision: Integer): BigInteger;
191
var
192
  Two: BigInteger;
193
begin
194
  Two := 2;
195
  Result := N;
196
  if Result = Two then
197
    Exit;
198
  if Result.IsEven then
199
    Result := Result.FlipBit(0);
200
  while not IsProbablePrime(Result, Precision) do
201
  begin
202
    Result := Result + Two;
203
  end;
204
  Result := N;
205
end;
206

207
function IsWitness(const A, N: BigInteger): Boolean;
208
var
209
  R: Integer;
210
  D: BigInteger;
211
begin
212
  //  Write N - 1 as (2 ^ Power) * Factor, where Factor is odd.
213
  //  Repeatedly try to divide N - 1 by 2
214
  R := 1;
215
  D := (N - BigInteger.One) shr 1;
216
  while D.IsEven do
217
  begin
218
    D := D shr 1;
219
    Inc(R);
220
  end;
221

222
  // Now check if A is a witness to N's compositeness
223
  Result := IsComposite(A, D, N, R);
224
end;
225

226
function IsComposite(const A, D, N: BigInteger; S: Integer): Boolean;
227
var
228
  I: Integer;
229
  NLessOne: BigInteger;
230
begin
231
  NLessOne := N - BigInteger.One;
232

233
  // if A^(2^0 * D) ≡ 1 (mod N) then prime.
234
  if BigInteger.ModPow(A, D, N) = BigInteger.One then
235
    Exit(False);
236

237
  for I := 0 to S - 1 do
238
    // if A^(2^I * D) ≡ N - 1 (mod N) then prime.
239
    if BigInteger.ModPow(A, (BigInteger.One shl I) * D, N) = NLessOne then
240
      Exit(False);
241

242
  Result := True;
243
end;
244

245
initialization
246
  Two := 2;
247
  DeterminicityThreshold := UInt64(341550071728321);
248
  CertainlyComposite := UInt64(3215031751);
249
  Random := TRandom.Create(Round(Now * SecsPerDay * MSecsPerSec));
250

251
end.
252

253

254