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