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