Crypto++  8.2
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 
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 
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 
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 
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
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 
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 
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 
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 
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 
624 bool 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 
649 Integer 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 
664 Integer 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 /*
675 Integer 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 
720 Integer 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 
788 int 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 
815 Integer 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 /*
851 Integer 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 
1001 Integer 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 
1022 unsigned 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 
1030 unsigned 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 
1039 void 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 
1124 NAMESPACE_END
1125 
1126 #endif
Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u)
Chinese Remainder Theorem.
Definition: nbtheory.cpp:553
An invalid argument was detected.
Definition: cryptlib.h:202
unsigned int WordCount() const
Determines the number of words required to represent the Integer.
Definition: integer.cpp:3331
const word16 * GetPrimeTable(unsigned int &size)
The Small Prime table.
Definition: nbtheory.cpp:53
Classes for working with NameValuePairs.
bool IsLucasProbablePrime(const Integer &n)
Determine if a number is probably prime.
Definition: nbtheory.cpp:155
bool SafeConvert(T1 from, T2 &to)
Tests whether a conversion from -> to is safe to perform.
Definition: misc.h:622
a number which is probabilistically prime
Definition: integer.h:95
Utility functions for the Crypto++ library.
Restricts the instantiation of a class to one static object without locks.
Definition: misc.h:263
bool IsSquare() const
Determine whether this integer is a perfect square.
Definition: integer.cpp:4390
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
bool SmallDivisorsTest(const Integer &p)
Tests whether a number is divisible by a small prime.
Definition: nbtheory.cpp:89
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
bool IsStrongLucasProbablePrime(const Integer &n)
Determine if a number is probably prime.
Definition: nbtheory.cpp:182
signed long ConvertToLong() const
Convert the Integer to Long.
Definition: integer.cpp:3012
const Integer & Subtract(const Integer &a, const Integer &b) const
Subtracts elements in the ring.
Definition: integer.cpp:4578
Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
Generates a provable prime.
Definition: nbtheory.cpp:510
Classes for automatic resource management.
bool IsSmallPrime(const Integer &p)
Tests whether a number is a small prime.
Definition: nbtheory.cpp:60
bool IsNegative() const
Determines if the Integer is negative.
Definition: integer.h:336
Interface for random number generators.
Definition: cryptlib.h:1383
Common C++ header files.
bool IsFermatProbablePrime(const Integer &n, const Integer &b)
Determine if a number is probably prime.
Definition: nbtheory.cpp:96
void Randomize(RandomNumberGenerator &rng, size_t bitCount)
Set this Integer to random integer.
Definition: integer.cpp:3503
int Jacobi(const Integer &a, const Integer &b)
Calculate the Jacobi symbol.
Definition: nbtheory.cpp:788
Integer InverseMod(const Integer &n) const
Calculate multiplicative inverse.
Definition: integer.cpp:4430
Integer ModularSquareRoot(const Integer &a, const Integer &p)
Extract a modular square root.
Definition: nbtheory.cpp:572
bool IsPositive() const
Determines if the Integer is positive.
Definition: integer.h:342
static const Integer & One()
Integer representing 1.
Definition: integer.cpp:4877
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:1001
bool IsStrongProbablePrime(const Integer &n, const Integer &b)
Determine if a number is probably prime.
Definition: nbtheory.cpp:105
Integer ConvertIn(const Integer &a) const
Reduces an element in the congruence class.
Definition: modarith.h:292
Pointer that overloads operator ->
Definition: smartptr.h:36
Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
Generates a provable prime.
Definition: nbtheory.cpp:470
Integer GCD(const Integer &a, const Integer &b)
Calculate the greatest common divisor.
Definition: nbtheory.h:142
a number with no special properties
Definition: integer.h:93
Integer Lucas(const Integer &e, const Integer &p, const Integer &n)
Calculate the Lucas value.
Definition: nbtheory.cpp:815
AlgorithmParameters MakeParameters(const char *name, const T &value, bool throwIfNotUsed=true)
Create an object that implements NameValuePairs.
Definition: algparam.h:502
bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level=1)
Verifies a number is probably prime.
Definition: nbtheory.cpp:247
Application callback to signal suitability of a cabdidate prime.
Definition: nbtheory.h:113
static Integer Power2(size_t e)
Exponentiates to a power of 2.
Definition: integer.cpp:3079
Multiple precision integer with arithmetic operations.
Definition: integer.h:49
Precompiled header file.
const T1 UnsignedMin(const T1 &a, const T2 &b)
Safe comparison of values that could be neagtive and incorrectly promoted.
Definition: misc.h:606
const Integer & Multiply(const Integer &a, const Integer &b) const
Multiplies elements in the ring.
Definition: integer.cpp:4659
bool IsPrime(const Integer &p)
Verifies a number is probably prime.
Definition: nbtheory.cpp:237
static const Integer & Two()
Integer representing 2.
Definition: integer.cpp:4889
const Integer & Square(const Integer &a) const
Square an element in the ring.
Definition: integer.cpp:4672
bool IsEven() const
Determines if the Integer is even parity.
Definition: integer.h:348
const T & STDMIN(const T &a, const T &b)
Replacement function for std::min.
Definition: misc.h:567
#define CRYPTOPP_ASSERT(exp)
Debugging and diagnostic assertion.
Definition: trap.h:69
bool TrialDivision(const Integer &p, unsigned bound)
Tests whether a number is divisible by a small prime.
Definition: nbtheory.cpp:71
unsigned int BitCount() const
Determines the number of bits required to represent the Integer.
Definition: integer.cpp:3345
Classes and functions for number theoretic operations.
unsigned int DiscreteLogWorkFactor(unsigned int bitlength)
Estimate work factor.
Definition: nbtheory.cpp:1030
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:649
Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b)
Calculate multiplicative inverse.
Definition: nbtheory.h:165
Integer Squared() const
Multiply this integer by itself.
Definition: integer.h:609
Performs modular arithmetic in Montgomery representation for increased speed.
Definition: modarith.h:274
An object that implements NameValuePairs.
Definition: algparam.h:419
void Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned qbits)
Generate a Prime and Generator.
Definition: nbtheory.cpp:1039
bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds)
Determine if a number is probably prime.
Definition: nbtheory.cpp:138
Multiple precision integer with arithmetic operations.
static const Integer & Zero()
Integer representing 0.
Definition: integer.cpp:4865
Class file for performing modular arithmetic.
Crypto++ library namespace.
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:624
Integer ModularExponentiation(const Integer &x, const Integer &e, const Integer &m)
Modular exponentiation.
Definition: nbtheory.h:215
bool GetBit(size_t i) const
Provides the i-th bit of the Integer.
Definition: integer.cpp:3103
unsigned int FactoringWorkFactor(unsigned int bitlength)
Estimate work factor.
Definition: nbtheory.cpp:1022
bool IsOdd() const
Determines if the Integer is odd parity.
Definition: integer.h:351
Integer ConvertOut(const Integer &a) const
Reduces an element in the congruence class.
Definition: integer.cpp:4685