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