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