Crypto++ 7.0 Free C++ class library of cryptographic schemes

nbtheory.cpp
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 {
25 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
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
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
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
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
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
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
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
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
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 {
231 Integer * operator()() const
232 {
233 return new Integer(Integer(s_lastSmallPrime).Squared());
234 }
235 };
236
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
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
255 unsigned int PrimeSearchInterval(const Integer &max)
256 {
257 return max.BitCount();
258 }
259
260 static inline bool FastProbablePrimeTest(const Integer &n)
261 {
262 return IsStrongProbablePrime(n,2);
263 }
264
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
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
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
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
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
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
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
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 then 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
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
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
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
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
646 Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq,
647 const Integer &p, const Integer &q, const Integer &u)
648 {
649 Integer p2, q2;
650 #pragma omp parallel
651 #pragma omp sections
652 {
653 #pragma omp section
654 p2 = ModularExponentiation((a % p), dp, p);
655 #pragma omp section
656 q2 = ModularExponentiation((a % q), dq, q);
657 }
658 return CRT(p2, p, q2, q, u);
659 }
660
661 Integer ModularRoot(const Integer &a, const Integer &e,
662 const Integer &p, const Integer &q)
663 {
664 Integer dp = EuclideanMultiplicativeInverse(e, p-1);
665 Integer dq = EuclideanMultiplicativeInverse(e, q-1);
666 Integer u = EuclideanMultiplicativeInverse(p, q);
667 CRYPTOPP_ASSERT(!!dp && !!dq && !!u);
668 return ModularRoot(a, dp, dq, p, q, u);
669 }
670
671 /*
672 Integer GCDI(const Integer &x, const Integer &y)
673 {
674 Integer a=x, b=y;
675 unsigned k=0;
676
677 CRYPTOPP_ASSERT(!!a && !!b);
678
679 while (a[0]==0 && b[0]==0)
680 {
681 a >>= 1;
682 b >>= 1;
683 k++;
684 }
685
686 while (a[0]==0)
687 a >>= 1;
688
689 while (b[0]==0)
690 b >>= 1;
691
692 while (1)
693 {
694 switch (a.Compare(b))
695 {
696 case -1:
697 b -= a;
698 while (b[0]==0)
699 b >>= 1;
700 break;
701
702 case 0:
703 return (a <<= k);
704
705 case 1:
706 a -= b;
707 while (a[0]==0)
708 a >>= 1;
709 break;
710
711 default:
712 CRYPTOPP_ASSERT(false);
713 }
714 }
715 }
716
717 Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b)
718 {
719 CRYPTOPP_ASSERT(b.Positive());
720
721 if (a.Negative())
722 return EuclideanMultiplicativeInverse(a%b, b);
723
724 if (b[0]==0)
725 {
726 if (!b || a[0]==0)
727 return Integer::Zero(); // no inverse
728 if (a==1)
729 return 1;
730 Integer u = EuclideanMultiplicativeInverse(b, a);
731 if (!u)
732 return Integer::Zero(); // no inverse
733 else
734 return (b*(a-u)+1)/a;
735 }
736
737 Integer u=1, d=a, v1=b, v3=b, t1, t3, b2=(b+1)>>1;
738
739 if (a[0])
740 {
741 t1 = Integer::Zero();
742 t3 = -b;
743 }
744 else
745 {
746 t1 = b2;
747 t3 = a>>1;
748 }
749
750 while (!!t3)
751 {
752 while (t3[0]==0)
753 {
754 t3 >>= 1;
755 if (t1[0]==0)
756 t1 >>= 1;
757 else
758 {
759 t1 >>= 1;
760 t1 += b2;
761 }
762 }
763 if (t3.Positive())
764 {
765 u = t1;
766 d = t3;
767 }
768 else
769 {
770 v1 = b-t1;
771 v3 = -t3;
772 }
773 t1 = u-v1;
774 t3 = d-v3;
775 if (t1.Negative())
776 t1 += b;
777 }
778 if (d==1)
779 return u;
780 else
781 return Integer::Zero(); // no inverse
782 }
783 */
784
785 int Jacobi(const Integer &aIn, const Integer &bIn)
786 {
787 CRYPTOPP_ASSERT(bIn.IsOdd());
788
789 Integer b = bIn, a = aIn%bIn;
790 int result = 1;
791
792 while (!!a)
793 {
794 unsigned i=0;
795 while (a.GetBit(i)==0)
796 i++;
797 a>>=i;
798
799 if (i%2==1 && (b%8==3 || b%8==5))
800 result = -result;
801
802 if (a%4==3 && b%4==3)
803 result = -result;
804
805 std::swap(a, b);
806 a %= b;
807 }
808
809 return (b==1) ? result : 0;
810 }
811
812 Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n)
813 {
814 unsigned i = e.BitCount();
815 if (i==0)
816 return Integer::Two();
817
818 MontgomeryRepresentation m(n);
819 Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two());
820 Integer v=p, v1=m.Subtract(m.Square(p), two);
821
822 i--;
823 while (i--)
824 {
825 if (e.GetBit(i))
826 {
827 // v = (v*v1 - p) % m;
828 v = m.Subtract(m.Multiply(v,v1), p);
829 // v1 = (v1*v1 - 2) % m;
830 v1 = m.Subtract(m.Square(v1), two);
831 }
832 else
833 {
834 // v1 = (v*v1 - p) % m;
835 v1 = m.Subtract(m.Multiply(v,v1), p);
836 // v = (v*v - 2) % m;
837 v = m.Subtract(m.Square(v), two);
838 }
839 }
840 return m.ConvertOut(v);
841 }
842
843 // This is Peter Montgomery's unpublished Lucas sequence evalutation algorithm.
844 // The total number of multiplies and squares used is less than the binary
845 // algorithm (see above). Unfortunately I can't get it to run as fast as
846 // the binary algorithm because of the extra overhead.
847 /*
848 Integer Lucas(const Integer &n, const Integer &P, const Integer &modulus)
849 {
850 if (!n)
851 return 2;
852
853 #define f(A, B, C) m.Subtract(m.Multiply(A, B), C)
854 #define X2(A) m.Subtract(m.Square(A), two)
855 #define X3(A) m.Multiply(A, m.Subtract(m.Square(A), three))
856
857 MontgomeryRepresentation m(modulus);
858 Integer two=m.ConvertIn(2), three=m.ConvertIn(3);
859 Integer A=m.ConvertIn(P), B, C, p, d=n, e, r, t, T, U;
860
861 while (d!=1)
862 {
863 p = d;
864 unsigned int b = WORD_BITS * p.WordCount();
865 Integer alpha = (Integer(5)<<(2*b-2)).SquareRoot() - Integer::Power2(b-1);
866 r = (p*alpha)>>b;
867 e = d-r;
868 B = A;
869 C = two;
870 d = r;
871
872 while (d!=e)
873 {
874 if (d<e)
875 {
876 swap(d, e);
877 swap(A, B);
878 }
879
880 unsigned int dm2 = d[0], em2 = e[0];
881 unsigned int dm3 = d%3, em3 = e%3;
882
883 // if ((dm6+em6)%3 == 0 && d <= e + (e>>2))
884 if ((dm3+em3==0 || dm3+em3==3) && (t = e, t >>= 2, t += e, d <= t))
885 {
886 // #1
887 // t = (d+d-e)/3;
888 // t = d; t += d; t -= e; t /= 3;
889 // e = (e+e-d)/3;
890 // e += e; e -= d; e /= 3;
891 // d = t;
892
893 // t = (d+e)/3
894 t = d; t += e; t /= 3;
895 e -= t;
896 d -= t;
897
898 T = f(A, B, C);
899 U = f(T, A, B);
900 B = f(T, B, A);
901 A = U;
902 continue;
903 }
904
905 // if (dm6 == em6 && d <= e + (e>>2))
906 if (dm3 == em3 && dm2 == em2 && (t = e, t >>= 2, t += e, d <= t))
907 {
908 // #2
909 // d = (d-e)>>1;
910 d -= e; d >>= 1;
911 B = f(A, B, C);
912 A = X2(A);
913 continue;
914 }
915
916 // if (d <= (e<<2))
917 if (d <= (t = e, t <<= 2))
918 {
919 // #3
920 d -= e;
921 C = f(A, B, C);
922 swap(B, C);
923 continue;
924 }
925
926 if (dm2 == em2)
927 {
928 // #4
929 // d = (d-e)>>1;
930 d -= e; d >>= 1;
931 B = f(A, B, C);
932 A = X2(A);
933 continue;
934 }
935
936 if (dm2 == 0)
937 {
938 // #5
939 d >>= 1;
940 C = f(A, C, B);
941 A = X2(A);
942 continue;
943 }
944
945 if (dm3 == 0)
946 {
947 // #6
948 // d = d/3 - e;
949 d /= 3; d -= e;
950 T = X2(A);
951 C = f(T, f(A, B, C), C);
952 swap(B, C);
953 A = f(T, A, A);
954 continue;
955 }
956
957 if (dm3+em3==0 || dm3+em3==3)
958 {
959 // #7
960 // d = (d-e-e)/3;
961 d -= e; d -= e; d /= 3;
962 T = f(A, B, C);
963 B = f(T, A, B);
964 A = X3(A);
965 continue;
966 }
967
968 if (dm3 == em3)
969 {
970 // #8
971 // d = (d-e)/3;
972 d -= e; d /= 3;
973 T = f(A, B, C);
974 C = f(A, C, B);
975 B = T;
976 A = X3(A);
977 continue;
978 }
979
980 CRYPTOPP_ASSERT(em2 == 0);
981 // #9
982 e >>= 1;
983 C = f(C, B, A);
984 B = X2(B);
985 }
986
987 A = f(A, B, C);
988 }
989
990 #undef f
991 #undef X2
992 #undef X3
993
994 return m.ConvertOut(A);
995 }
996 */
997
998 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
999 {
1000 Integer d = (m*m-4);
1001 Integer p2, q2;
1002 #pragma omp parallel
1003 #pragma omp sections
1004 {
1005 #pragma omp section
1006 {
1007 p2 = p-Jacobi(d,p);
1008 p2 = Lucas(EuclideanMultiplicativeInverse(e,p2), m, p);
1009 }
1010 #pragma omp section
1011 {
1012 q2 = q-Jacobi(d,q);
1013 q2 = Lucas(EuclideanMultiplicativeInverse(e,q2), m, q);
1014 }
1015 }
1016 return CRT(p2, p, q2, q, u);
1017 }
1018
1019 unsigned int FactoringWorkFactor(unsigned int n)
1020 {
1021 // extrapolated from the table in Odlyzko's "The Future of Integer Factorization"
1022 // updated to reflect the factoring of RSA-130
1023 if (n<5) return 0;
1024 else return (unsigned int)(2.4 * std::pow((double)n, 1.0/3.0) * std::pow(log(double(n)), 2.0/3.0) - 5);
1025 }
1026
1027 unsigned int DiscreteLogWorkFactor(unsigned int n)
1028 {
1029 // assuming discrete log takes about the same time as factoring
1030 if (n<5) return 0;
1031 else return (unsigned int)(2.4 * std::pow((double)n, 1.0/3.0) * std::pow(log(double(n)), 2.0/3.0) - 5);
1032 }
1033
1034 // ********************************************************
1035
1036 void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits)
1037 {
1038 // no prime exists for delta = -1, qbits = 4, and pbits = 5
1039 CRYPTOPP_ASSERT(qbits > 4);
1040 CRYPTOPP_ASSERT(pbits > qbits);
1041
1042 if (qbits+1 == pbits)
1043 {
1044 Integer minP = Integer::Power2(pbits-1);
1045 Integer maxP = Integer::Power2(pbits) - 1;
1046 bool success = false;
1047
1048 while (!success)
1049 {
1050 p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
1051 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
1052
1053 while (sieve.NextCandidate(p))
1054 {
1055 CRYPTOPP_ASSERT(IsSmallPrime(p) || SmallDivisorsTest(p));
1056 q = (p-delta) >> 1;
1057 CRYPTOPP_ASSERT(IsSmallPrime(q) || SmallDivisorsTest(q));
1058 if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
1059 {
1060 success = true;
1061 break;
1062 }
1063 }
1064 }
1065
1066 if (delta == 1)
1067 {
1068 // find g such that g is a quadratic residue mod p, then g has order q
1069 // g=4 always works, but this way we get the smallest quadratic residue (other than 1)
1070 for (g=2; Jacobi(g, p) != 1; ++g) {}
1071 // contributed by Walt Tuvell: g should be the following according to the Law of Quadratic Reciprocity
1072 CRYPTOPP_ASSERT((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
1073 }
1074 else
1075 {
1076 CRYPTOPP_ASSERT(delta == -1);
1077 // find g such that g*g-4 is a quadratic non-residue,
1078 // and such that g has order q
1079 for (g=3; ; ++g)
1080 if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
1081 break;
1082 }
1083 }
1084 else
1085 {
1086 Integer minQ = Integer::Power2(qbits-1);
1087 Integer maxQ = Integer::Power2(qbits) - 1;
1088 Integer minP = Integer::Power2(pbits-1);
1089 Integer maxP = Integer::Power2(pbits) - 1;
1090
1091 do
1092 {
1093 q.Randomize(rng, minQ, maxQ, Integer::PRIME);
1094 } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
1095
1096 // find a random g of order q
1097 if (delta==1)
1098 {
1099 do
1100 {
1101 Integer h(rng, 2, p-2, Integer::ANY);
1102 g = a_exp_b_mod_c(h, (p-1)/q, p);
1103 } while (g <= 1);
1104 CRYPTOPP_ASSERT(a_exp_b_mod_c(g, q, p)==1);
1105 }
1106 else
1107 {
1108 CRYPTOPP_ASSERT(delta==-1);
1109 do
1110 {
1111 Integer h(rng, 3, p-1, Integer::ANY);
1112 if (Jacobi(h*h-4, p)==1)
1113 continue;
1114 g = Lucas((p+1)/q, h, p);
1115 } while (g <= 2);
1116 CRYPTOPP_ASSERT(Lucas(q, g, p) == 2);
1117 }
1118 }
1119 }
1120
1121 NAMESPACE_END
1122
1123 #endif
CRT
Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u)
Chinese Remainder Theorem.
Definition: nbtheory.cpp:553
InvalidArgument
An invalid argument was detected.
Definition: cryptlib.h:202
Integer::WordCount
unsigned int WordCount() const
Determines the number of words required to represent the Integer.
Definition: integer.cpp:3326
GetPrimeTable
const word16 * GetPrimeTable(unsigned int &size)
The Small Prime table.
Definition: nbtheory.cpp:53
algparam.h
Classes for working with NameValuePairs.
IsLucasProbablePrime
bool IsLucasProbablePrime(const Integer &n)
Determine if a number is probably prime.
Definition: nbtheory.cpp:155
SafeConvert
bool SafeConvert(T1 from, T2 &to)
Tests whether a conversion from -> to is safe to perform.
Definition: misc.h:618
NewLastSmallPrimeSquared
Definition: nbtheory.cpp:229
PrimeSieve
Definition: nbtheory.cpp:286
Integer::PRIME
a number which is probabilistically prime
Definition: integer.h:95
misc.h
Utility functions for the Crypto++ library.
Singleton
Restricts the instantiation of a class to one static object without locks.
Definition: misc.h:291
Integer::IsSquare
bool IsSquare() const
Determine whether this integer is a perfect square.
Definition: integer.cpp:4367
FirstPrime
bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector)
Finds a random prime of special form.
Definition: nbtheory.cpp:379
SmallDivisorsTest
bool SmallDivisorsTest(const Integer &p)
Tests whether a number is divisible by a small prime.
Definition: nbtheory.cpp:89
RandomNumberGenerator::GenerateWord32
virtual word32 GenerateWord32(word32 min=0, word32 max=0xffffffffUL)
Generate a random 32 bit word in the range min to max, inclusive.
Definition: cryptlib.cpp:283
IsStrongLucasProbablePrime
bool IsStrongLucasProbablePrime(const Integer &n)
Determine if a number is probably prime.
Definition: nbtheory.cpp:182
Integer::ConvertToLong
signed long ConvertToLong() const
Convert the Integer to Long.
Definition: integer.cpp:3009
ModularArithmetic::Subtract
const Integer & Subtract(const Integer &a, const Integer &b) const
Subtracts elements in the ring.
Definition: integer.cpp:4555
MaurerProvablePrime
Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
Generates a provable prime.
Definition: nbtheory.cpp:510
smartptr.h
Classes for automatic resource management.
IsSmallPrime
bool IsSmallPrime(const Integer &p)
Tests whether a number is a small prime.
Definition: nbtheory.cpp:60
Integer::IsNegative
bool IsNegative() const
Determines if the Integer is negative.
Definition: integer.h:336
RandomNumberGenerator
Interface for random number generators.
Definition: cryptlib.h:1349
stdcpp.h
Common C++ header files.
IsFermatProbablePrime
bool IsFermatProbablePrime(const Integer &n, const Integer &b)
Determine if a number is probably prime.
Definition: nbtheory.cpp:96
Integer::Randomize
void Randomize(RandomNumberGenerator &rng, size_t bitCount)
Set this Integer to random integer.
Definition: integer.cpp:3497
Jacobi
int Jacobi(const Integer &a, const Integer &b)
Calculate the Jacobi symbol.
Definition: nbtheory.cpp:785
Integer::InverseMod
Integer InverseMod(const Integer &n) const
Calculate multiplicative inverse.
Definition: integer.cpp:4407
ModularSquareRoot
Integer ModularSquareRoot(const Integer &a, const Integer &p)
Extract a modular square root.
Definition: nbtheory.cpp:572
Integer::IsPositive
bool IsPositive() const
Determines if the Integer is positive.
Definition: integer.h:342
Integer::One
static const Integer & One()
Integer representing 1.
Definition: integer.cpp:4854
InverseLucas
Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
Calculate the inverse Lucas value.
Definition: nbtheory.cpp:998
IsStrongProbablePrime
bool IsStrongProbablePrime(const Integer &n, const Integer &b)
Determine if a number is probably prime.
Definition: nbtheory.cpp:105
MontgomeryRepresentation::ConvertIn
Integer ConvertIn(const Integer &a) const
Reduces an element in the congruence class.
Definition: modarith.h:292
member_ptr
Pointer that overloads operator ->
Definition: smartptr.h:36
MihailescuProvablePrime
Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
Generates a provable prime.
Definition: nbtheory.cpp:470
GCD
Integer GCD(const Integer &a, const Integer &b)
Calculate the greatest common divisor.
Definition: nbtheory.h:142
Integer::ANY
a number with no special properties
Definition: integer.h:93
Lucas
Integer Lucas(const Integer &e, const Integer &p, const Integer &n)
Calculate the Lucas value.
Definition: nbtheory.cpp:812
MakeParameters
AlgorithmParameters MakeParameters(const char *name, const T &value, bool throwIfNotUsed=true)
Create an object that implements NameValuePairs.
Definition: algparam.h:502
VerifyPrime
bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level=1)
Verifies a number is probably prime.
Definition: nbtheory.cpp:247
PrimeSelector
Application callback to signal suitability of a cabdidate prime.
Definition: nbtheory.h:113
Integer::Power2
static Integer Power2(size_t e)
Exponentiates to a power of 2.
Definition: integer.cpp:3076
Integer
Multiple precision integer with arithmetic operations.
Definition: integer.h:49
pch.h
Precompiled header file.
UnsignedMin
const T1 UnsignedMin(const T1 &a, const T2 &b)
Safe comparison of values that could be neagtive and incorrectly promoted.
Definition: misc.h:602
MontgomeryRepresentation::Multiply
const Integer & Multiply(const Integer &a, const Integer &b) const
Multiplies elements in the ring.
Definition: integer.cpp:4636
IsPrime
bool IsPrime(const Integer &p)
Verifies a number is probably prime.
Definition: nbtheory.cpp:237
Integer::Two
static const Integer & Two()
Integer representing 2.
Definition: integer.cpp:4866
MontgomeryRepresentation::Square
const Integer & Square(const Integer &a) const
Square an element in the ring.
Definition: integer.cpp:4649
Integer::IsEven
bool IsEven() const
Determines if the Integer is even parity.
Definition: integer.h:348
STDMIN
const T & STDMIN(const T &a, const T &b)
Replacement function for std::min.
Definition: misc.h:563
CRYPTOPP_ASSERT
#define CRYPTOPP_ASSERT(exp)
Debugging and diagnostic assertion.
Definition: trap.h:60
TrialDivision
bool TrialDivision(const Integer &p, unsigned bound)
Tests whether a number is divisible by a small prime.
Definition: nbtheory.cpp:71
Integer::BitCount
unsigned int BitCount() const
Determines the number of bits required to represent the Integer.
Definition: integer.cpp:3340
nbtheory.h
Classes and functions for number theoretic operations.
DiscreteLogWorkFactor
unsigned int DiscreteLogWorkFactor(unsigned int bitlength)
Estimate work factor.
Definition: nbtheory.cpp:1027
ModularRoot
Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq, const Integer &p, const Integer &q, const Integer &u)
Extract a modular root.
Definition: nbtheory.cpp:646
EuclideanMultiplicativeInverse
Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b)
Calculate multiplicative inverse.
Definition: nbtheory.h:165
Integer::Squared
Integer Squared() const
Multiply this integer by itself.
Definition: integer.h:609
NewPrimeTable
Definition: nbtheory.cpp:23
MontgomeryRepresentation
Performs modular arithmetic in Montgomery representation for increased speed.
Definition: modarith.h:274
AlgorithmParameters
An object that implements NameValuePairs.
Definition: algparam.h:419
PrimeAndGenerator::Generate
void Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned qbits)
Generate a Prime and Generator.
Definition: nbtheory.cpp:1036
RabinMillerTest
bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds)
Determine if a number is probably prime.
Definition: nbtheory.cpp:138
integer.h
Multiple precision integer with arithmetic operations.
Integer::Zero
static const Integer & Zero()
Integer representing 0.
Definition: integer.cpp:4842
modarith.h
Class file for performing modular arithmetic.
CryptoPP
Crypto++ library namespace.
SolveModularQuadraticEquation
bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p)
Solve a Modular Quadratic Equation.
Definition: nbtheory.cpp:621
ModularExponentiation
Integer ModularExponentiation(const Integer &x, const Integer &e, const Integer &m)
Modular exponentiation.
Definition: nbtheory.h:215
Integer::GetBit
bool GetBit(size_t i) const
Provides the i-th bit of the Integer.
Definition: integer.cpp:3100
FactoringWorkFactor
unsigned int FactoringWorkFactor(unsigned int bitlength)
Estimate work factor.
Definition: nbtheory.cpp:1019
Integer::IsOdd
bool IsOdd() const
Determines if the Integer is odd parity.
Definition: integer.h:351
MontgomeryRepresentation::ConvertOut
Integer ConvertOut(const Integer &a) const
Reduces an element in the congruence class.
Definition: integer.cpp:4662