/src/cryptopp/nbtheory.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | // nbtheory.cpp - originally written and placed in the public domain by Wei Dai |
2 | | |
3 | | #include "pch.h" |
4 | | |
5 | | #ifndef CRYPTOPP_IMPORTS |
6 | | |
7 | | #include "nbtheory.h" |
8 | | #include "integer.h" |
9 | | #include "modarith.h" |
10 | | #include "algparam.h" |
11 | | #include "smartptr.h" |
12 | | #include "misc.h" |
13 | | #include "stdcpp.h" |
14 | | #include "trap.h" |
15 | | |
16 | | #ifdef _OPENMP |
17 | | # include <omp.h> |
18 | | #endif |
19 | | |
20 | | NAMESPACE_BEGIN(CryptoPP) |
21 | | |
22 | | // Keep sync'd with primetab.cpp |
23 | | const unsigned int maxPrimeTableSize = 3511; |
24 | | const word s_lastSmallPrime = 32719; |
25 | | |
26 | | const word16 * GetPrimeTable(unsigned int &size) |
27 | 0 | { |
28 | 0 | extern const word16 precomputedPrimeTable[maxPrimeTableSize]; |
29 | 0 | size = maxPrimeTableSize; |
30 | 0 | return precomputedPrimeTable; |
31 | 0 | } |
32 | | |
33 | | bool IsSmallPrime(const Integer &p) |
34 | 0 | { |
35 | 0 | unsigned int primeTableSize; |
36 | 0 | const word16 * primeTable = GetPrimeTable(primeTableSize); |
37 | |
|
38 | 0 | if (p.IsPositive() && p <= primeTable[primeTableSize-1]) |
39 | 0 | return std::binary_search(primeTable, primeTable+primeTableSize, (word16)p.ConvertToLong()); |
40 | 0 | else |
41 | 0 | return false; |
42 | 0 | } |
43 | | |
44 | | bool TrialDivision(const Integer &p, unsigned bound) |
45 | 0 | { |
46 | 0 | unsigned int primeTableSize; |
47 | 0 | const word16 * primeTable = GetPrimeTable(primeTableSize); |
48 | |
|
49 | 0 | CRYPTOPP_ASSERT(primeTable[primeTableSize-1] >= bound); |
50 | |
|
51 | 0 | unsigned int i; |
52 | 0 | for (i = 0; primeTable[i]<bound; i++) |
53 | 0 | if ((p % primeTable[i]) == 0) |
54 | 0 | return true; |
55 | | |
56 | 0 | if (bound == primeTable[i]) |
57 | 0 | return (p % bound == 0); |
58 | 0 | else |
59 | 0 | return false; |
60 | 0 | } |
61 | | |
62 | | bool SmallDivisorsTest(const Integer &p) |
63 | 0 | { |
64 | 0 | unsigned int primeTableSize; |
65 | 0 | const word16 * primeTable = GetPrimeTable(primeTableSize); |
66 | 0 | return !TrialDivision(p, primeTable[primeTableSize-1]); |
67 | 0 | } |
68 | | |
69 | | bool IsFermatProbablePrime(const Integer &n, const Integer &b) |
70 | 0 | { |
71 | 0 | if (n <= 3) |
72 | 0 | return n==2 || n==3; |
73 | | |
74 | 0 | CRYPTOPP_ASSERT(n>3 && b>1 && b<n-1); |
75 | 0 | return a_exp_b_mod_c(b, n-1, n)==1; |
76 | 0 | } |
77 | | |
78 | | bool IsStrongProbablePrime(const Integer &n, const Integer &b) |
79 | 0 | { |
80 | 0 | if (n <= 3) |
81 | 0 | return n==2 || n==3; |
82 | | |
83 | 0 | CRYPTOPP_ASSERT(n>3 && b>1 && b<n-1); |
84 | |
|
85 | 0 | if ((n.IsEven() && n!=2) || GCD(b, n) != 1) |
86 | 0 | return false; |
87 | | |
88 | 0 | Integer nminus1 = (n-1); |
89 | 0 | unsigned int a; |
90 | | |
91 | | // calculate a = largest power of 2 that divides (n-1) |
92 | 0 | for (a=0; ; a++) |
93 | 0 | if (nminus1.GetBit(a)) |
94 | 0 | break; |
95 | 0 | Integer m = nminus1>>a; |
96 | |
|
97 | 0 | Integer z = a_exp_b_mod_c(b, m, n); |
98 | 0 | if (z==1 || z==nminus1) |
99 | 0 | return true; |
100 | 0 | for (unsigned j=1; j<a; j++) |
101 | 0 | { |
102 | 0 | z = z.Squared()%n; |
103 | 0 | if (z==nminus1) |
104 | 0 | return true; |
105 | 0 | if (z==1) |
106 | 0 | return false; |
107 | 0 | } |
108 | 0 | return false; |
109 | 0 | } |
110 | | |
111 | | bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds) |
112 | 0 | { |
113 | 0 | if (n <= 3) |
114 | 0 | return n==2 || n==3; |
115 | | |
116 | 0 | CRYPTOPP_ASSERT(n>3); |
117 | |
|
118 | 0 | Integer b; |
119 | 0 | for (unsigned int i=0; i<rounds; i++) |
120 | 0 | { |
121 | 0 | b.Randomize(rng, 2, n-2); |
122 | 0 | if (!IsStrongProbablePrime(n, b)) |
123 | 0 | return false; |
124 | 0 | } |
125 | 0 | return true; |
126 | 0 | } |
127 | | |
128 | | bool IsLucasProbablePrime(const Integer &n) |
129 | 0 | { |
130 | 0 | if (n <= 1) |
131 | 0 | return false; |
132 | | |
133 | 0 | if (n.IsEven()) |
134 | 0 | return n==2; |
135 | | |
136 | 0 | CRYPTOPP_ASSERT(n>2); |
137 | |
|
138 | 0 | Integer b=3; |
139 | 0 | unsigned int i=0; |
140 | 0 | int j; |
141 | |
|
142 | 0 | while ((j=Jacobi(b.Squared()-4, n)) == 1) |
143 | 0 | { |
144 | 0 | if (++i==64 && n.IsSquare()) // avoid infinite loop if n is a square |
145 | 0 | return false; |
146 | 0 | ++b; ++b; |
147 | 0 | } |
148 | | |
149 | 0 | if (j==0) |
150 | 0 | return false; |
151 | 0 | else |
152 | 0 | return Lucas(n+1, b, n)==2; |
153 | 0 | } |
154 | | |
155 | | bool IsStrongLucasProbablePrime(const Integer &n) |
156 | 0 | { |
157 | 0 | if (n <= 1) |
158 | 0 | return false; |
159 | | |
160 | 0 | if (n.IsEven()) |
161 | 0 | return n==2; |
162 | | |
163 | 0 | CRYPTOPP_ASSERT(n>2); |
164 | |
|
165 | 0 | Integer b=3; |
166 | 0 | unsigned int i=0; |
167 | 0 | int j; |
168 | |
|
169 | 0 | while ((j=Jacobi(b.Squared()-4, n)) == 1) |
170 | 0 | { |
171 | 0 | if (++i==64 && n.IsSquare()) // avoid infinite loop if n is a square |
172 | 0 | return false; |
173 | 0 | ++b; ++b; |
174 | 0 | } |
175 | | |
176 | 0 | if (j==0) |
177 | 0 | return false; |
178 | | |
179 | 0 | Integer n1 = n+1; |
180 | 0 | unsigned int a; |
181 | | |
182 | | // calculate a = largest power of 2 that divides n1 |
183 | 0 | for (a=0; ; a++) |
184 | 0 | if (n1.GetBit(a)) |
185 | 0 | break; |
186 | 0 | Integer m = n1>>a; |
187 | |
|
188 | 0 | Integer z = Lucas(m, b, n); |
189 | 0 | if (z==2 || z==n-2) |
190 | 0 | return true; |
191 | 0 | for (i=1; i<a; i++) |
192 | 0 | { |
193 | 0 | z = (z.Squared()-2)%n; |
194 | 0 | if (z==n-2) |
195 | 0 | return true; |
196 | 0 | if (z==2) |
197 | 0 | return false; |
198 | 0 | } |
199 | 0 | return false; |
200 | 0 | } |
201 | | |
202 | | struct NewLastSmallPrimeSquared |
203 | | { |
204 | | Integer * operator()() const |
205 | 0 | { |
206 | 0 | return new Integer(Integer(s_lastSmallPrime).Squared()); |
207 | 0 | } |
208 | | }; |
209 | | |
210 | | bool IsPrime(const Integer &p) |
211 | 0 | { |
212 | 0 | if (p <= s_lastSmallPrime) |
213 | 0 | return IsSmallPrime(p); |
214 | 0 | else if (p <= Singleton<Integer, NewLastSmallPrimeSquared>().Ref()) |
215 | 0 | return SmallDivisorsTest(p); |
216 | 0 | else |
217 | 0 | return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p); |
218 | 0 | } |
219 | | |
220 | | bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level) |
221 | 0 | { |
222 | 0 | bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1); |
223 | 0 | if (level >= 1) |
224 | 0 | pass = pass && RabinMillerTest(rng, p, 10); |
225 | 0 | return pass; |
226 | 0 | } |
227 | | |
228 | | unsigned int PrimeSearchInterval(const Integer &max) |
229 | 0 | { |
230 | 0 | return max.BitCount(); |
231 | 0 | } |
232 | | |
233 | | static inline bool FastProbablePrimeTest(const Integer &n) |
234 | 0 | { |
235 | 0 | return IsStrongProbablePrime(n,2); |
236 | 0 | } |
237 | | |
238 | | AlgorithmParameters MakeParametersForTwoPrimesOfEqualSize(unsigned int productBitLength) |
239 | 0 | { |
240 | 0 | if (productBitLength < 16) |
241 | 0 | throw InvalidArgument("invalid bit length"); |
242 | | |
243 | 0 | Integer minP, maxP; |
244 | |
|
245 | 0 | if (productBitLength%2==0) |
246 | 0 | { |
247 | 0 | minP = Integer(182) << (productBitLength/2-8); |
248 | 0 | maxP = Integer::Power2(productBitLength/2)-1; |
249 | 0 | } |
250 | 0 | else |
251 | 0 | { |
252 | 0 | minP = Integer::Power2((productBitLength-1)/2); |
253 | 0 | maxP = Integer(181) << ((productBitLength+1)/2-8); |
254 | 0 | } |
255 | |
|
256 | 0 | return MakeParameters("RandomNumberType", Integer::PRIME)("Min", minP)("Max", maxP); |
257 | 0 | } |
258 | | |
259 | | class PrimeSieve |
260 | | { |
261 | | public: |
262 | | // delta == 1 or -1 means double sieve with p = 2*q + delta |
263 | | PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta=0); |
264 | | bool NextCandidate(Integer &c); |
265 | | |
266 | | void DoSieve(); |
267 | | static void SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv); |
268 | | |
269 | | Integer m_first, m_last, m_step; |
270 | | signed int m_delta; |
271 | | word m_next; |
272 | | std::vector<bool> m_sieve; |
273 | | }; |
274 | | |
275 | | PrimeSieve::PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta) |
276 | | : m_first(first), m_last(last), m_step(step), m_delta(delta), m_next(0) |
277 | 0 | { |
278 | 0 | DoSieve(); |
279 | 0 | } |
280 | | |
281 | | bool PrimeSieve::NextCandidate(Integer &c) |
282 | 0 | { |
283 | 0 | bool safe = SafeConvert(std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin(), m_next); |
284 | 0 | CRYPTOPP_UNUSED(safe); CRYPTOPP_ASSERT(safe); |
285 | 0 | if (m_next == m_sieve.size()) |
286 | 0 | { |
287 | 0 | m_first += long(m_sieve.size())*m_step; |
288 | 0 | if (m_first > m_last) |
289 | 0 | return false; |
290 | 0 | else |
291 | 0 | { |
292 | 0 | m_next = 0; |
293 | 0 | DoSieve(); |
294 | 0 | return NextCandidate(c); |
295 | 0 | } |
296 | 0 | } |
297 | 0 | else |
298 | 0 | { |
299 | 0 | c = m_first + long(m_next)*m_step; |
300 | 0 | ++m_next; |
301 | 0 | return true; |
302 | 0 | } |
303 | 0 | } |
304 | | |
305 | | void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv) |
306 | 0 | { |
307 | 0 | if (stepInv) |
308 | 0 | { |
309 | 0 | size_t sieveSize = sieve.size(); |
310 | 0 | size_t j = (word32(p-(first%p))*stepInv) % p; |
311 | | // if the first multiple of p is p, skip it |
312 | 0 | if (first.WordCount() <= 1 && first + step*long(j) == p) |
313 | 0 | j += p; |
314 | 0 | for (; j < sieveSize; j += p) |
315 | 0 | sieve[j] = true; |
316 | 0 | } |
317 | 0 | } |
318 | | |
319 | | void PrimeSieve::DoSieve() |
320 | 0 | { |
321 | 0 | unsigned int primeTableSize; |
322 | 0 | const word16 * primeTable = GetPrimeTable(primeTableSize); |
323 | |
|
324 | 0 | const unsigned int maxSieveSize = 32768; |
325 | 0 | unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong(); |
326 | |
|
327 | 0 | m_sieve.clear(); |
328 | 0 | m_sieve.resize(sieveSize, false); |
329 | |
|
330 | 0 | if (m_delta == 0) |
331 | 0 | { |
332 | 0 | for (unsigned int i = 0; i < primeTableSize; ++i) |
333 | 0 | SieveSingle(m_sieve, primeTable[i], m_first, m_step, (word16)m_step.InverseMod(primeTable[i])); |
334 | 0 | } |
335 | 0 | else |
336 | 0 | { |
337 | 0 | CRYPTOPP_ASSERT(m_step%2==0); |
338 | 0 | Integer qFirst = (m_first-m_delta) >> 1; |
339 | 0 | Integer halfStep = m_step >> 1; |
340 | 0 | for (unsigned int i = 0; i < primeTableSize; ++i) |
341 | 0 | { |
342 | 0 | word16 p = primeTable[i]; |
343 | 0 | word16 stepInv = (word16)m_step.InverseMod(p); |
344 | 0 | SieveSingle(m_sieve, p, m_first, m_step, stepInv); |
345 | |
|
346 | 0 | word16 halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p; |
347 | 0 | SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv); |
348 | 0 | } |
349 | 0 | } |
350 | 0 | } |
351 | | |
352 | | bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector) |
353 | 0 | { |
354 | 0 | CRYPTOPP_ASSERT(!equiv.IsNegative() && equiv < mod); |
355 | |
|
356 | 0 | Integer gcd = GCD(equiv, mod); |
357 | 0 | if (gcd != Integer::One()) |
358 | 0 | { |
359 | | // the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv) |
360 | 0 | if (p <= gcd && gcd <= max && IsPrime(gcd) && (!pSelector || pSelector->IsAcceptable(gcd))) |
361 | 0 | { |
362 | 0 | p = gcd; |
363 | 0 | return true; |
364 | 0 | } |
365 | 0 | else |
366 | 0 | return false; |
367 | 0 | } |
368 | | |
369 | 0 | unsigned int primeTableSize; |
370 | 0 | const word16 * primeTable = GetPrimeTable(primeTableSize); |
371 | |
|
372 | 0 | if (p <= primeTable[primeTableSize-1]) |
373 | 0 | { |
374 | 0 | const word16 *pItr; |
375 | |
|
376 | 0 | --p; |
377 | 0 | if (p.IsPositive()) |
378 | 0 | pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong()); |
379 | 0 | else |
380 | 0 | pItr = primeTable; |
381 | |
|
382 | 0 | while (pItr < primeTable+primeTableSize && !(*pItr%mod == equiv && (!pSelector || pSelector->IsAcceptable(*pItr)))) |
383 | 0 | ++pItr; |
384 | |
|
385 | 0 | if (pItr < primeTable+primeTableSize) |
386 | 0 | { |
387 | 0 | p = *pItr; |
388 | 0 | return p <= max; |
389 | 0 | } |
390 | | |
391 | 0 | p = primeTable[primeTableSize-1]+1; |
392 | 0 | } |
393 | | |
394 | 0 | CRYPTOPP_ASSERT(p > primeTable[primeTableSize-1]); |
395 | |
|
396 | 0 | if (mod.IsOdd()) |
397 | 0 | return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector); |
398 | | |
399 | 0 | p += (equiv-p)%mod; |
400 | |
|
401 | 0 | if (p>max) |
402 | 0 | return false; |
403 | | |
404 | 0 | PrimeSieve sieve(p, max, mod); |
405 | |
|
406 | 0 | while (sieve.NextCandidate(p)) |
407 | 0 | { |
408 | 0 | if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p)) |
409 | 0 | return true; |
410 | 0 | } |
411 | | |
412 | 0 | return false; |
413 | 0 | } |
414 | | |
415 | | // the following two functions are based on code and comments provided by Preda Mihailescu |
416 | | static bool ProvePrime(const Integer &p, const Integer &q) |
417 | 0 | { |
418 | 0 | CRYPTOPP_ASSERT(p < q*q*q); |
419 | 0 | CRYPTOPP_ASSERT(p % q == 1); |
420 | | |
421 | | // this is the Quisquater test. Numbers p having passed the Lucas - Lehmer test |
422 | | // for q and verifying p < q^3 can only be built up of two factors, both = 1 mod q, |
423 | | // or be prime. The next two lines build the discriminant of a quadratic equation |
424 | | // which holds iff p is built up of two factors (exercise ... ) |
425 | |
|
426 | 0 | Integer r = (p-1)/q; |
427 | 0 | if (((r%q).Squared()-4*(r/q)).IsSquare()) |
428 | 0 | return false; |
429 | | |
430 | 0 | unsigned int primeTableSize; |
431 | 0 | const word16 * primeTable = GetPrimeTable(primeTableSize); |
432 | |
|
433 | 0 | CRYPTOPP_ASSERT(primeTableSize >= 50); |
434 | 0 | for (int i=0; i<50; i++) |
435 | 0 | { |
436 | 0 | Integer b = a_exp_b_mod_c(primeTable[i], r, p); |
437 | 0 | if (b != 1) |
438 | 0 | return a_exp_b_mod_c(b, q, p) == 1; |
439 | 0 | } |
440 | 0 | return false; |
441 | 0 | } |
442 | | |
443 | | Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits) |
444 | 0 | { |
445 | 0 | Integer p; |
446 | 0 | Integer minP = Integer::Power2(pbits-1); |
447 | 0 | Integer maxP = Integer::Power2(pbits) - 1; |
448 | |
|
449 | 0 | if (maxP <= Integer(s_lastSmallPrime).Squared()) |
450 | 0 | { |
451 | | // Randomize() will generate a prime provable by trial division |
452 | 0 | p.Randomize(rng, minP, maxP, Integer::PRIME); |
453 | 0 | return p; |
454 | 0 | } |
455 | | |
456 | 0 | unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36); |
457 | 0 | Integer q = MihailescuProvablePrime(rng, qbits); |
458 | 0 | Integer q2 = q<<1; |
459 | |
|
460 | 0 | while (true) |
461 | 0 | { |
462 | | // this initializes the sieve to search in the arithmetic |
463 | | // progression p = p_0 + \lambda * q2 = p_0 + 2 * \lambda * q, |
464 | | // with q the recursively generated prime above. We will be able |
465 | | // to use Lucas tets for proving primality. A trick of Quisquater |
466 | | // allows taking q > cubic_root(p) rather than square_root: this |
467 | | // decreases the recursion. |
468 | |
|
469 | 0 | p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2); |
470 | 0 | PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2); |
471 | |
|
472 | 0 | while (sieve.NextCandidate(p)) |
473 | 0 | { |
474 | 0 | if (FastProbablePrimeTest(p) && ProvePrime(p, q)) |
475 | 0 | return p; |
476 | 0 | } |
477 | 0 | } |
478 | | |
479 | | // not reached |
480 | 0 | return p; |
481 | 0 | } |
482 | | |
483 | | Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits) |
484 | 0 | { |
485 | 0 | const unsigned smallPrimeBound = 29, c_opt=10; |
486 | 0 | Integer p; |
487 | |
|
488 | 0 | unsigned int primeTableSize; |
489 | 0 | const word16 * primeTable = GetPrimeTable(primeTableSize); |
490 | |
|
491 | 0 | if (bits < smallPrimeBound) |
492 | 0 | { |
493 | 0 | do |
494 | 0 | p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2); |
495 | 0 | while (TrialDivision(p, 1 << ((bits+1)/2))); |
496 | 0 | } |
497 | 0 | else |
498 | 0 | { |
499 | 0 | const unsigned margin = bits > 50 ? 20 : (bits-10)/2; |
500 | 0 | double relativeSize; |
501 | 0 | do |
502 | 0 | relativeSize = std::pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1); |
503 | 0 | while (bits * relativeSize >= bits - margin); |
504 | |
|
505 | 0 | Integer a,b; |
506 | 0 | Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize)); |
507 | 0 | Integer I = Integer::Power2(bits-2)/q; |
508 | 0 | Integer I2 = I << 1; |
509 | 0 | unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt); |
510 | 0 | bool success = false; |
511 | 0 | while (!success) |
512 | 0 | { |
513 | 0 | p.Randomize(rng, I, I2, Integer::ANY); |
514 | 0 | p *= q; p <<= 1; ++p; |
515 | 0 | if (!TrialDivision(p, trialDivisorBound)) |
516 | 0 | { |
517 | 0 | a.Randomize(rng, 2, p-1, Integer::ANY); |
518 | 0 | b = a_exp_b_mod_c(a, (p-1)/q, p); |
519 | 0 | success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1); |
520 | 0 | } |
521 | 0 | } |
522 | 0 | } |
523 | 0 | return p; |
524 | 0 | } |
525 | | |
526 | | Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u) |
527 | 0 | { |
528 | | // isn't operator overloading great? |
529 | 0 | return p * (u * (xq-xp) % q) + xp; |
530 | | /* |
531 | | Integer t1 = xq-xp; |
532 | | cout << hex << t1 << endl; |
533 | | Integer t2 = u * t1; |
534 | | cout << hex << t2 << endl; |
535 | | Integer t3 = t2 % q; |
536 | | cout << hex << t3 << endl; |
537 | | Integer t4 = p * t3; |
538 | | cout << hex << t4 << endl; |
539 | | Integer t5 = t4 + xp; |
540 | | cout << hex << t5 << endl; |
541 | | return t5; |
542 | | */ |
543 | 0 | } |
544 | | |
545 | | Integer ModularSquareRoot(const Integer &a, const Integer &p) |
546 | 0 | { |
547 | | // Callers must ensure p is prime, GH #1249 |
548 | 0 | CRYPTOPP_ASSERT(IsPrime(p)); |
549 | |
|
550 | 0 | if (p%4 == 3) |
551 | 0 | return a_exp_b_mod_c(a, (p+1)/4, p); |
552 | | |
553 | 0 | Integer q=p-1; |
554 | 0 | unsigned int r=0; |
555 | 0 | while (q.IsEven()) |
556 | 0 | { |
557 | 0 | r++; |
558 | 0 | q >>= 1; |
559 | 0 | } |
560 | |
|
561 | 0 | Integer n=2; |
562 | 0 | while (Jacobi(n, p) != -1) |
563 | 0 | ++n; |
564 | |
|
565 | 0 | Integer y = a_exp_b_mod_c(n, q, p); |
566 | 0 | Integer x = a_exp_b_mod_c(a, (q-1)/2, p); |
567 | 0 | Integer b = (x.Squared()%p)*a%p; |
568 | 0 | x = a*x%p; |
569 | 0 | Integer tempb, t; |
570 | |
|
571 | 0 | while (b != 1) |
572 | 0 | { |
573 | 0 | unsigned m=0; |
574 | 0 | tempb = b; |
575 | 0 | do |
576 | 0 | { |
577 | 0 | m++; |
578 | 0 | b = b.Squared()%p; |
579 | 0 | if (m==r) |
580 | 0 | return Integer::Zero(); |
581 | 0 | } |
582 | 0 | while (b != 1); |
583 | | |
584 | 0 | t = y; |
585 | 0 | for (unsigned i=0; i<r-m-1; i++) |
586 | 0 | t = t.Squared()%p; |
587 | 0 | y = t.Squared()%p; |
588 | 0 | r = m; |
589 | 0 | x = x*t%p; |
590 | 0 | b = tempb*y%p; |
591 | 0 | } |
592 | | |
593 | 0 | CRYPTOPP_ASSERT(x.Squared()%p == a); |
594 | 0 | return x; |
595 | 0 | } |
596 | | |
597 | | bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p) |
598 | 0 | { |
599 | | // Callers must ensure p is prime, GH #1249 |
600 | 0 | CRYPTOPP_ASSERT(IsPrime(p)); |
601 | |
|
602 | 0 | Integer D = (b.Squared() - 4*a*c) % p; |
603 | 0 | switch (Jacobi(D, p)) |
604 | 0 | { |
605 | 0 | default: |
606 | 0 | CRYPTOPP_ASSERT(false); // not reached |
607 | 0 | return false; |
608 | 0 | case -1: |
609 | 0 | return false; |
610 | 0 | case 0: |
611 | 0 | r1 = r2 = (-b*(a+a).InverseMod(p)) % p; |
612 | 0 | CRYPTOPP_ASSERT(((r1.Squared()*a + r1*b + c) % p).IsZero()); |
613 | 0 | return true; |
614 | 0 | case 1: |
615 | 0 | Integer s = ModularSquareRoot(D, p); |
616 | 0 | Integer t = (a+a).InverseMod(p); |
617 | 0 | r1 = (s-b)*t % p; |
618 | 0 | r2 = (-s-b)*t % p; |
619 | 0 | CRYPTOPP_ASSERT(((r1.Squared()*a + r1*b + c) % p).IsZero()); |
620 | 0 | CRYPTOPP_ASSERT(((r2.Squared()*a + r2*b + c) % p).IsZero()); |
621 | 0 | return true; |
622 | 0 | } |
623 | 0 | } |
624 | | |
625 | | Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq, |
626 | | const Integer &p, const Integer &q, const Integer &u) |
627 | 0 | { |
628 | | // Callers must ensure p and q are prime, GH #1249 |
629 | 0 | CRYPTOPP_ASSERT(IsPrime(p) && IsPrime(q)); |
630 | | |
631 | | // GCC warning bug, https://stackoverflow.com/q/12842306/608639 |
632 | | #ifdef _OPENMP |
633 | | Integer p2, q2; |
634 | | #pragma omp parallel |
635 | | #pragma omp sections |
636 | | { |
637 | | #pragma omp section |
638 | | p2 = ModularExponentiation((a % p), dp, p); |
639 | | #pragma omp section |
640 | | q2 = ModularExponentiation((a % q), dq, q); |
641 | | } |
642 | | #else |
643 | 0 | const Integer p2 = ModularExponentiation((a % p), dp, p); |
644 | 0 | const Integer q2 = ModularExponentiation((a % q), dq, q); |
645 | 0 | #endif |
646 | |
|
647 | 0 | return CRT(p2, p, q2, q, u); |
648 | 0 | } |
649 | | |
650 | | Integer ModularRoot(const Integer &a, const Integer &e, |
651 | | const Integer &p, const Integer &q) |
652 | 0 | { |
653 | | // Callers must ensure p and q are prime, GH #1249 |
654 | 0 | CRYPTOPP_ASSERT(IsPrime(p) && IsPrime(q)); |
655 | |
|
656 | 0 | Integer dp = EuclideanMultiplicativeInverse(e, p-1); |
657 | 0 | Integer dq = EuclideanMultiplicativeInverse(e, q-1); |
658 | 0 | Integer u = EuclideanMultiplicativeInverse(p, q); |
659 | 0 | CRYPTOPP_ASSERT(!!dp && !!dq && !!u); |
660 | 0 | return ModularRoot(a, dp, dq, p, q, u); |
661 | 0 | } |
662 | | |
663 | | /* |
664 | | Integer GCDI(const Integer &x, const Integer &y) |
665 | | { |
666 | | Integer a=x, b=y; |
667 | | unsigned k=0; |
668 | | |
669 | | CRYPTOPP_ASSERT(!!a && !!b); |
670 | | |
671 | | while (a[0]==0 && b[0]==0) |
672 | | { |
673 | | a >>= 1; |
674 | | b >>= 1; |
675 | | k++; |
676 | | } |
677 | | |
678 | | while (a[0]==0) |
679 | | a >>= 1; |
680 | | |
681 | | while (b[0]==0) |
682 | | b >>= 1; |
683 | | |
684 | | while (1) |
685 | | { |
686 | | switch (a.Compare(b)) |
687 | | { |
688 | | case -1: |
689 | | b -= a; |
690 | | while (b[0]==0) |
691 | | b >>= 1; |
692 | | break; |
693 | | |
694 | | case 0: |
695 | | return (a <<= k); |
696 | | |
697 | | case 1: |
698 | | a -= b; |
699 | | while (a[0]==0) |
700 | | a >>= 1; |
701 | | break; |
702 | | |
703 | | default: |
704 | | CRYPTOPP_ASSERT(false); |
705 | | } |
706 | | } |
707 | | } |
708 | | |
709 | | Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b) |
710 | | { |
711 | | CRYPTOPP_ASSERT(b.Positive()); |
712 | | |
713 | | if (a.Negative()) |
714 | | return EuclideanMultiplicativeInverse(a%b, b); |
715 | | |
716 | | if (b[0]==0) |
717 | | { |
718 | | if (!b || a[0]==0) |
719 | | return Integer::Zero(); // no inverse |
720 | | if (a==1) |
721 | | return 1; |
722 | | Integer u = EuclideanMultiplicativeInverse(b, a); |
723 | | if (!u) |
724 | | return Integer::Zero(); // no inverse |
725 | | else |
726 | | return (b*(a-u)+1)/a; |
727 | | } |
728 | | |
729 | | Integer u=1, d=a, v1=b, v3=b, t1, t3, b2=(b+1)>>1; |
730 | | |
731 | | if (a[0]) |
732 | | { |
733 | | t1 = Integer::Zero(); |
734 | | t3 = -b; |
735 | | } |
736 | | else |
737 | | { |
738 | | t1 = b2; |
739 | | t3 = a>>1; |
740 | | } |
741 | | |
742 | | while (!!t3) |
743 | | { |
744 | | while (t3[0]==0) |
745 | | { |
746 | | t3 >>= 1; |
747 | | if (t1[0]==0) |
748 | | t1 >>= 1; |
749 | | else |
750 | | { |
751 | | t1 >>= 1; |
752 | | t1 += b2; |
753 | | } |
754 | | } |
755 | | if (t3.Positive()) |
756 | | { |
757 | | u = t1; |
758 | | d = t3; |
759 | | } |
760 | | else |
761 | | { |
762 | | v1 = b-t1; |
763 | | v3 = -t3; |
764 | | } |
765 | | t1 = u-v1; |
766 | | t3 = d-v3; |
767 | | if (t1.Negative()) |
768 | | t1 += b; |
769 | | } |
770 | | if (d==1) |
771 | | return u; |
772 | | else |
773 | | return Integer::Zero(); // no inverse |
774 | | } |
775 | | */ |
776 | | |
777 | | int Jacobi(const Integer &aIn, const Integer &bIn) |
778 | 639 | { |
779 | 639 | CRYPTOPP_ASSERT(bIn.IsOdd()); |
780 | | |
781 | 639 | Integer b = bIn, a = aIn%bIn; |
782 | 639 | int result = 1; |
783 | | |
784 | 519k | while (!!a) |
785 | 518k | { |
786 | 518k | unsigned i=0; |
787 | 1.08M | while (a.GetBit(i)==0) |
788 | 566k | i++; |
789 | 518k | a>>=i; |
790 | | |
791 | 518k | if (i%2==1 && (b%8==3 || b%8==5)) |
792 | 94.6k | result = -result; |
793 | | |
794 | 518k | if (a%4==3 && b%4==3) |
795 | 128k | result = -result; |
796 | | |
797 | 518k | std::swap(a, b); |
798 | 518k | a %= b; |
799 | 518k | } |
800 | | |
801 | 639 | return (b==1) ? result : 0; |
802 | 639 | } |
803 | | |
804 | | Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n) |
805 | 0 | { |
806 | 0 | unsigned i = e.BitCount(); |
807 | 0 | if (i==0) |
808 | 0 | return Integer::Two(); |
809 | | |
810 | 0 | MontgomeryRepresentation m(n); |
811 | 0 | Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two()); |
812 | 0 | Integer v=p, v1=m.Subtract(m.Square(p), two); |
813 | |
|
814 | 0 | i--; |
815 | 0 | while (i--) |
816 | 0 | { |
817 | 0 | if (e.GetBit(i)) |
818 | 0 | { |
819 | | // v = (v*v1 - p) % m; |
820 | 0 | v = m.Subtract(m.Multiply(v,v1), p); |
821 | | // v1 = (v1*v1 - 2) % m; |
822 | 0 | v1 = m.Subtract(m.Square(v1), two); |
823 | 0 | } |
824 | 0 | else |
825 | 0 | { |
826 | | // v1 = (v*v1 - p) % m; |
827 | 0 | v1 = m.Subtract(m.Multiply(v,v1), p); |
828 | | // v = (v*v - 2) % m; |
829 | 0 | v = m.Subtract(m.Square(v), two); |
830 | 0 | } |
831 | 0 | } |
832 | 0 | return m.ConvertOut(v); |
833 | 0 | } |
834 | | |
835 | | // This is Peter Montgomery's unpublished Lucas sequence evaluation algorithm. |
836 | | // The total number of multiplies and squares used is less than the binary |
837 | | // algorithm (see above). Unfortunately I can't get it to run as fast as |
838 | | // the binary algorithm because of the extra overhead. |
839 | | /* |
840 | | Integer Lucas(const Integer &n, const Integer &P, const Integer &modulus) |
841 | | { |
842 | | if (!n) |
843 | | return 2; |
844 | | |
845 | | #define f(A, B, C) m.Subtract(m.Multiply(A, B), C) |
846 | | #define X2(A) m.Subtract(m.Square(A), two) |
847 | | #define X3(A) m.Multiply(A, m.Subtract(m.Square(A), three)) |
848 | | |
849 | | MontgomeryRepresentation m(modulus); |
850 | | Integer two=m.ConvertIn(2), three=m.ConvertIn(3); |
851 | | Integer A=m.ConvertIn(P), B, C, p, d=n, e, r, t, T, U; |
852 | | |
853 | | while (d!=1) |
854 | | { |
855 | | p = d; |
856 | | unsigned int b = WORD_BITS * p.WordCount(); |
857 | | Integer alpha = (Integer(5)<<(2*b-2)).SquareRoot() - Integer::Power2(b-1); |
858 | | r = (p*alpha)>>b; |
859 | | e = d-r; |
860 | | B = A; |
861 | | C = two; |
862 | | d = r; |
863 | | |
864 | | while (d!=e) |
865 | | { |
866 | | if (d<e) |
867 | | { |
868 | | swap(d, e); |
869 | | swap(A, B); |
870 | | } |
871 | | |
872 | | unsigned int dm2 = d[0], em2 = e[0]; |
873 | | unsigned int dm3 = d%3, em3 = e%3; |
874 | | |
875 | | // if ((dm6+em6)%3 == 0 && d <= e + (e>>2)) |
876 | | if ((dm3+em3==0 || dm3+em3==3) && (t = e, t >>= 2, t += e, d <= t)) |
877 | | { |
878 | | // #1 |
879 | | // t = (d+d-e)/3; |
880 | | // t = d; t += d; t -= e; t /= 3; |
881 | | // e = (e+e-d)/3; |
882 | | // e += e; e -= d; e /= 3; |
883 | | // d = t; |
884 | | |
885 | | // t = (d+e)/3 |
886 | | t = d; t += e; t /= 3; |
887 | | e -= t; |
888 | | d -= t; |
889 | | |
890 | | T = f(A, B, C); |
891 | | U = f(T, A, B); |
892 | | B = f(T, B, A); |
893 | | A = U; |
894 | | continue; |
895 | | } |
896 | | |
897 | | // if (dm6 == em6 && d <= e + (e>>2)) |
898 | | if (dm3 == em3 && dm2 == em2 && (t = e, t >>= 2, t += e, d <= t)) |
899 | | { |
900 | | // #2 |
901 | | // d = (d-e)>>1; |
902 | | d -= e; d >>= 1; |
903 | | B = f(A, B, C); |
904 | | A = X2(A); |
905 | | continue; |
906 | | } |
907 | | |
908 | | // if (d <= (e<<2)) |
909 | | if (d <= (t = e, t <<= 2)) |
910 | | { |
911 | | // #3 |
912 | | d -= e; |
913 | | C = f(A, B, C); |
914 | | swap(B, C); |
915 | | continue; |
916 | | } |
917 | | |
918 | | if (dm2 == em2) |
919 | | { |
920 | | // #4 |
921 | | // d = (d-e)>>1; |
922 | | d -= e; d >>= 1; |
923 | | B = f(A, B, C); |
924 | | A = X2(A); |
925 | | continue; |
926 | | } |
927 | | |
928 | | if (dm2 == 0) |
929 | | { |
930 | | // #5 |
931 | | d >>= 1; |
932 | | C = f(A, C, B); |
933 | | A = X2(A); |
934 | | continue; |
935 | | } |
936 | | |
937 | | if (dm3 == 0) |
938 | | { |
939 | | // #6 |
940 | | // d = d/3 - e; |
941 | | d /= 3; d -= e; |
942 | | T = X2(A); |
943 | | C = f(T, f(A, B, C), C); |
944 | | swap(B, C); |
945 | | A = f(T, A, A); |
946 | | continue; |
947 | | } |
948 | | |
949 | | if (dm3+em3==0 || dm3+em3==3) |
950 | | { |
951 | | // #7 |
952 | | // d = (d-e-e)/3; |
953 | | d -= e; d -= e; d /= 3; |
954 | | T = f(A, B, C); |
955 | | B = f(T, A, B); |
956 | | A = X3(A); |
957 | | continue; |
958 | | } |
959 | | |
960 | | if (dm3 == em3) |
961 | | { |
962 | | // #8 |
963 | | // d = (d-e)/3; |
964 | | d -= e; d /= 3; |
965 | | T = f(A, B, C); |
966 | | C = f(A, C, B); |
967 | | B = T; |
968 | | A = X3(A); |
969 | | continue; |
970 | | } |
971 | | |
972 | | CRYPTOPP_ASSERT(em2 == 0); |
973 | | // #9 |
974 | | e >>= 1; |
975 | | C = f(C, B, A); |
976 | | B = X2(B); |
977 | | } |
978 | | |
979 | | A = f(A, B, C); |
980 | | } |
981 | | |
982 | | #undef f |
983 | | #undef X2 |
984 | | #undef X3 |
985 | | |
986 | | return m.ConvertOut(A); |
987 | | } |
988 | | */ |
989 | | |
990 | | Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u) |
991 | 0 | { |
992 | | // Callers must ensure p and q are prime, GH #1249 |
993 | 0 | CRYPTOPP_ASSERT(IsPrime(p) && IsPrime(q)); |
994 | | |
995 | | // GCC warning bug, https://stackoverflow.com/q/12842306/608639 |
996 | | #ifdef _OPENMP |
997 | | Integer d = (m*m-4), p2, q2; |
998 | | #pragma omp parallel |
999 | | #pragma omp sections |
1000 | | { |
1001 | | #pragma omp section |
1002 | | { |
1003 | | p2 = p-Jacobi(d,p); |
1004 | | p2 = Lucas(EuclideanMultiplicativeInverse(e,p2), m, p); |
1005 | | } |
1006 | | #pragma omp section |
1007 | | { |
1008 | | q2 = q-Jacobi(d,q); |
1009 | | q2 = Lucas(EuclideanMultiplicativeInverse(e,q2), m, q); |
1010 | | } |
1011 | | } |
1012 | | #else |
1013 | 0 | const Integer d = (m*m-4); |
1014 | 0 | const Integer t1 = p-Jacobi(d,p); |
1015 | 0 | const Integer p2 = Lucas(EuclideanMultiplicativeInverse(e,t1), m, p); |
1016 | |
|
1017 | 0 | const Integer t2 = q-Jacobi(d,q); |
1018 | 0 | const Integer q2 = Lucas(EuclideanMultiplicativeInverse(e,t2), m, q); |
1019 | 0 | #endif |
1020 | |
|
1021 | 0 | return CRT(p2, p, q2, q, u); |
1022 | 0 | } |
1023 | | |
1024 | | unsigned int FactoringWorkFactor(unsigned int n) |
1025 | 0 | { |
1026 | | // extrapolated from the table in Odlyzko's "The Future of Integer Factorization" |
1027 | | // updated to reflect the factoring of RSA-130 |
1028 | 0 | if (n<5) return 0; |
1029 | 0 | else return (unsigned int)(2.4 * std::pow((double)n, 1.0/3.0) * std::pow(log(double(n)), 2.0/3.0) - 5); |
1030 | 0 | } |
1031 | | |
1032 | | unsigned int DiscreteLogWorkFactor(unsigned int n) |
1033 | 0 | { |
1034 | | // assuming discrete log takes about the same time as factoring |
1035 | 0 | if (n<5) return 0; |
1036 | 0 | else return (unsigned int)(2.4 * std::pow((double)n, 1.0/3.0) * std::pow(log(double(n)), 2.0/3.0) - 5); |
1037 | 0 | } |
1038 | | |
1039 | | // ******************************************************** |
1040 | | |
1041 | | void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits) |
1042 | 0 | { |
1043 | | // no prime exists for delta = -1, qbits = 4, and pbits = 5 |
1044 | 0 | CRYPTOPP_ASSERT(qbits > 4); |
1045 | 0 | CRYPTOPP_ASSERT(pbits > qbits); |
1046 | |
|
1047 | 0 | if (qbits+1 == pbits) |
1048 | 0 | { |
1049 | 0 | Integer minP = Integer::Power2(pbits-1); |
1050 | 0 | Integer maxP = Integer::Power2(pbits) - 1; |
1051 | 0 | bool success = false; |
1052 | |
|
1053 | 0 | while (!success) |
1054 | 0 | { |
1055 | 0 | p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12); |
1056 | 0 | PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta); |
1057 | |
|
1058 | 0 | while (sieve.NextCandidate(p)) |
1059 | 0 | { |
1060 | 0 | CRYPTOPP_ASSERT(IsSmallPrime(p) || SmallDivisorsTest(p)); |
1061 | 0 | q = (p-delta) >> 1; |
1062 | 0 | CRYPTOPP_ASSERT(IsSmallPrime(q) || SmallDivisorsTest(q)); |
1063 | 0 | if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p)) |
1064 | 0 | { |
1065 | 0 | success = true; |
1066 | 0 | break; |
1067 | 0 | } |
1068 | 0 | } |
1069 | 0 | } |
1070 | |
|
1071 | 0 | if (delta == 1) |
1072 | 0 | { |
1073 | | // find g such that g is a quadratic residue mod p, then g has order q |
1074 | | // g=4 always works, but this way we get the smallest quadratic residue (other than 1) |
1075 | 0 | for (g=2; Jacobi(g, p) != 1; ++g) {} |
1076 | | // contributed by Walt Tuvell: g should be the following according to the Law of Quadratic Reciprocity |
1077 | 0 | CRYPTOPP_ASSERT((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4); |
1078 | 0 | } |
1079 | 0 | else |
1080 | 0 | { |
1081 | 0 | CRYPTOPP_ASSERT(delta == -1); |
1082 | | // find g such that g*g-4 is a quadratic non-residue, |
1083 | | // and such that g has order q |
1084 | 0 | for (g=3; ; ++g) |
1085 | 0 | if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2) |
1086 | 0 | break; |
1087 | 0 | } |
1088 | 0 | } |
1089 | 0 | else |
1090 | 0 | { |
1091 | 0 | Integer minQ = Integer::Power2(qbits-1); |
1092 | 0 | Integer maxQ = Integer::Power2(qbits) - 1; |
1093 | 0 | Integer minP = Integer::Power2(pbits-1); |
1094 | 0 | Integer maxP = Integer::Power2(pbits) - 1; |
1095 | |
|
1096 | 0 | do |
1097 | 0 | { |
1098 | 0 | q.Randomize(rng, minQ, maxQ, Integer::PRIME); |
1099 | 0 | } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q)); |
1100 | | |
1101 | | // find a random g of order q |
1102 | 0 | if (delta==1) |
1103 | 0 | { |
1104 | 0 | do |
1105 | 0 | { |
1106 | 0 | Integer h(rng, 2, p-2, Integer::ANY); |
1107 | 0 | g = a_exp_b_mod_c(h, (p-1)/q, p); |
1108 | 0 | } while (g <= 1); |
1109 | 0 | CRYPTOPP_ASSERT(a_exp_b_mod_c(g, q, p)==1); |
1110 | 0 | } |
1111 | 0 | else |
1112 | 0 | { |
1113 | 0 | CRYPTOPP_ASSERT(delta==-1); |
1114 | 0 | do |
1115 | 0 | { |
1116 | 0 | Integer h(rng, 3, p-1, Integer::ANY); |
1117 | 0 | if (Jacobi(h*h-4, p)==1) |
1118 | 0 | continue; |
1119 | 0 | g = Lucas((p+1)/q, h, p); |
1120 | 0 | } while (g <= 2); |
1121 | 0 | CRYPTOPP_ASSERT(Lucas(q, g, p) == 2); |
1122 | 0 | } |
1123 | 0 | } |
1124 | 0 | } |
1125 | | |
1126 | | NAMESPACE_END |
1127 | | |
1128 | | #endif |