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