Crypto++  5.6.5
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 
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 
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
244  return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p);
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 
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  Integer p2, q2;
650  #pragma omp parallel
651  #pragma omp sections
652  {
653  #pragma omp section
654  p2 = ModularExponentiation((a % p), dp, p);
655  #pragma omp section
656  q2 = ModularExponentiation((a % q), dq, q);
657  }
658  return CRT(p2, p, q2, q, u);
659 }
660 
661 Integer ModularRoot(const Integer &a, const Integer &e,
662  const Integer &p, const Integer &q)
663 {
664  Integer dp = EuclideanMultiplicativeInverse(e, p-1);
665  Integer dq = EuclideanMultiplicativeInverse(e, q-1);
666  Integer u = EuclideanMultiplicativeInverse(p, q);
667  CRYPTOPP_ASSERT(!!dp && !!dq && !!u);
668  return ModularRoot(a, dp, dq, p, q, u);
669 }
670 
671 /*
672 Integer GCDI(const Integer &x, const Integer &y)
673 {
674  Integer a=x, b=y;
675  unsigned k=0;
676 
677  CRYPTOPP_ASSERT(!!a && !!b);
678 
679  while (a[0]==0 && b[0]==0)
680  {
681  a >>= 1;
682  b >>= 1;
683  k++;
684  }
685 
686  while (a[0]==0)
687  a >>= 1;
688 
689  while (b[0]==0)
690  b >>= 1;
691 
692  while (1)
693  {
694  switch (a.Compare(b))
695  {
696  case -1:
697  b -= a;
698  while (b[0]==0)
699  b >>= 1;
700  break;
701 
702  case 0:
703  return (a <<= k);
704 
705  case 1:
706  a -= b;
707  while (a[0]==0)
708  a >>= 1;
709  break;
710 
711  default:
712  CRYPTOPP_ASSERT(false);
713  }
714  }
715 }
716 
717 Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b)
718 {
719  CRYPTOPP_ASSERT(b.Positive());
720 
721  if (a.Negative())
722  return EuclideanMultiplicativeInverse(a%b, b);
723 
724  if (b[0]==0)
725  {
726  if (!b || a[0]==0)
727  return Integer::Zero(); // no inverse
728  if (a==1)
729  return 1;
730  Integer u = EuclideanMultiplicativeInverse(b, a);
731  if (!u)
732  return Integer::Zero(); // no inverse
733  else
734  return (b*(a-u)+1)/a;
735  }
736 
737  Integer u=1, d=a, v1=b, v3=b, t1, t3, b2=(b+1)>>1;
738 
739  if (a[0])
740  {
741  t1 = Integer::Zero();
742  t3 = -b;
743  }
744  else
745  {
746  t1 = b2;
747  t3 = a>>1;
748  }
749 
750  while (!!t3)
751  {
752  while (t3[0]==0)
753  {
754  t3 >>= 1;
755  if (t1[0]==0)
756  t1 >>= 1;
757  else
758  {
759  t1 >>= 1;
760  t1 += b2;
761  }
762  }
763  if (t3.Positive())
764  {
765  u = t1;
766  d = t3;
767  }
768  else
769  {
770  v1 = b-t1;
771  v3 = -t3;
772  }
773  t1 = u-v1;
774  t3 = d-v3;
775  if (t1.Negative())
776  t1 += b;
777  }
778  if (d==1)
779  return u;
780  else
781  return Integer::Zero(); // no inverse
782 }
783 */
784 
785 int Jacobi(const Integer &aIn, const Integer &bIn)
786 {
787  CRYPTOPP_ASSERT(bIn.IsOdd());
788 
789  Integer b = bIn, a = aIn%bIn;
790  int result = 1;
791 
792  while (!!a)
793  {
794  unsigned i=0;
795  while (a.GetBit(i)==0)
796  i++;
797  a>>=i;
798 
799  if (i%2==1 && (b%8==3 || b%8==5))
800  result = -result;
801 
802  if (a%4==3 && b%4==3)
803  result = -result;
804 
805  std::swap(a, b);
806  a %= b;
807  }
808 
809  return (b==1) ? result : 0;
810 }
811 
812 Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n)
813 {
814  unsigned i = e.BitCount();
815  if (i==0)
816  return Integer::Two();
817 
819  Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two());
820  Integer v=p, v1=m.Subtract(m.Square(p), two);
821 
822  i--;
823  while (i--)
824  {
825  if (e.GetBit(i))
826  {
827  // v = (v*v1 - p) % m;
828  v = m.Subtract(m.Multiply(v,v1), p);
829  // v1 = (v1*v1 - 2) % m;
830  v1 = m.Subtract(m.Square(v1), two);
831  }
832  else
833  {
834  // v1 = (v*v1 - p) % m;
835  v1 = m.Subtract(m.Multiply(v,v1), p);
836  // v = (v*v - 2) % m;
837  v = m.Subtract(m.Square(v), two);
838  }
839  }
840  return m.ConvertOut(v);
841 }
842 
843 // This is Peter Montgomery's unpublished Lucas sequence evalutation algorithm.
844 // The total number of multiplies and squares used is less than the binary
845 // algorithm (see above). Unfortunately I can't get it to run as fast as
846 // the binary algorithm because of the extra overhead.
847 /*
848 Integer Lucas(const Integer &n, const Integer &P, const Integer &modulus)
849 {
850  if (!n)
851  return 2;
852 
853 #define f(A, B, C) m.Subtract(m.Multiply(A, B), C)
854 #define X2(A) m.Subtract(m.Square(A), two)
855 #define X3(A) m.Multiply(A, m.Subtract(m.Square(A), three))
856 
857  MontgomeryRepresentation m(modulus);
858  Integer two=m.ConvertIn(2), three=m.ConvertIn(3);
859  Integer A=m.ConvertIn(P), B, C, p, d=n, e, r, t, T, U;
860 
861  while (d!=1)
862  {
863  p = d;
864  unsigned int b = WORD_BITS * p.WordCount();
865  Integer alpha = (Integer(5)<<(2*b-2)).SquareRoot() - Integer::Power2(b-1);
866  r = (p*alpha)>>b;
867  e = d-r;
868  B = A;
869  C = two;
870  d = r;
871 
872  while (d!=e)
873  {
874  if (d<e)
875  {
876  swap(d, e);
877  swap(A, B);
878  }
879 
880  unsigned int dm2 = d[0], em2 = e[0];
881  unsigned int dm3 = d%3, em3 = e%3;
882 
883 // if ((dm6+em6)%3 == 0 && d <= e + (e>>2))
884  if ((dm3+em3==0 || dm3+em3==3) && (t = e, t >>= 2, t += e, d <= t))
885  {
886  // #1
887 // t = (d+d-e)/3;
888 // t = d; t += d; t -= e; t /= 3;
889 // e = (e+e-d)/3;
890 // e += e; e -= d; e /= 3;
891 // d = t;
892 
893 // t = (d+e)/3
894  t = d; t += e; t /= 3;
895  e -= t;
896  d -= t;
897 
898  T = f(A, B, C);
899  U = f(T, A, B);
900  B = f(T, B, A);
901  A = U;
902  continue;
903  }
904 
905 // if (dm6 == em6 && d <= e + (e>>2))
906  if (dm3 == em3 && dm2 == em2 && (t = e, t >>= 2, t += e, d <= t))
907  {
908  // #2
909 // d = (d-e)>>1;
910  d -= e; d >>= 1;
911  B = f(A, B, C);
912  A = X2(A);
913  continue;
914  }
915 
916 // if (d <= (e<<2))
917  if (d <= (t = e, t <<= 2))
918  {
919  // #3
920  d -= e;
921  C = f(A, B, C);
922  swap(B, C);
923  continue;
924  }
925 
926  if (dm2 == em2)
927  {
928  // #4
929 // d = (d-e)>>1;
930  d -= e; d >>= 1;
931  B = f(A, B, C);
932  A = X2(A);
933  continue;
934  }
935 
936  if (dm2 == 0)
937  {
938  // #5
939  d >>= 1;
940  C = f(A, C, B);
941  A = X2(A);
942  continue;
943  }
944 
945  if (dm3 == 0)
946  {
947  // #6
948 // d = d/3 - e;
949  d /= 3; d -= e;
950  T = X2(A);
951  C = f(T, f(A, B, C), C);
952  swap(B, C);
953  A = f(T, A, A);
954  continue;
955  }
956 
957  if (dm3+em3==0 || dm3+em3==3)
958  {
959  // #7
960 // d = (d-e-e)/3;
961  d -= e; d -= e; d /= 3;
962  T = f(A, B, C);
963  B = f(T, A, B);
964  A = X3(A);
965  continue;
966  }
967 
968  if (dm3 == em3)
969  {
970  // #8
971 // d = (d-e)/3;
972  d -= e; d /= 3;
973  T = f(A, B, C);
974  C = f(A, C, B);
975  B = T;
976  A = X3(A);
977  continue;
978  }
979 
980  CRYPTOPP_ASSERT(em2 == 0);
981  // #9
982  e >>= 1;
983  C = f(C, B, A);
984  B = X2(B);
985  }
986 
987  A = f(A, B, C);
988  }
989 
990 #undef f
991 #undef X2
992 #undef X3
993 
994  return m.ConvertOut(A);
995 }
996 */
997 
998 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
999 {
1000  Integer d = (m*m-4);
1001  Integer p2, q2;
1002  #pragma omp parallel
1003  #pragma omp sections
1004  {
1005  #pragma omp section
1006  {
1007  p2 = p-Jacobi(d,p);
1008  p2 = Lucas(EuclideanMultiplicativeInverse(e,p2), m, p);
1009  }
1010  #pragma omp section
1011  {
1012  q2 = q-Jacobi(d,q);
1013  q2 = Lucas(EuclideanMultiplicativeInverse(e,q2), m, q);
1014  }
1015  }
1016  return CRT(p2, p, q2, q, u);
1017 }
1018 
1019 unsigned int FactoringWorkFactor(unsigned int n)
1020 {
1021  // extrapolated from the table in Odlyzko's "The Future of Integer Factorization"
1022  // updated to reflect the factoring of RSA-130
1023  if (n<5) return 0;
1024  else return (unsigned int)(2.4 * std::pow((double)n, 1.0/3.0) * std::pow(log(double(n)), 2.0/3.0) - 5);
1025 }
1026 
1027 unsigned int DiscreteLogWorkFactor(unsigned int n)
1028 {
1029  // assuming discrete log takes about the same time as factoring
1030  if (n<5) return 0;
1031  else return (unsigned int)(2.4 * std::pow((double)n, 1.0/3.0) * std::pow(log(double(n)), 2.0/3.0) - 5);
1032 }
1033 
1034 // ********************************************************
1035 
1036 void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits)
1037 {
1038  // no prime exists for delta = -1, qbits = 4, and pbits = 5
1039  CRYPTOPP_ASSERT(qbits > 4);
1040  CRYPTOPP_ASSERT(pbits > qbits);
1041 
1042  if (qbits+1 == pbits)
1043  {
1044  Integer minP = Integer::Power2(pbits-1);
1045  Integer maxP = Integer::Power2(pbits) - 1;
1046  bool success = false;
1047 
1048  while (!success)
1049  {
1050  p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
1051  PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
1052 
1053  while (sieve.NextCandidate(p))
1054  {
1055  CRYPTOPP_ASSERT(IsSmallPrime(p) || SmallDivisorsTest(p));
1056  q = (p-delta) >> 1;
1057  CRYPTOPP_ASSERT(IsSmallPrime(q) || SmallDivisorsTest(q));
1058  if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
1059  {
1060  success = true;
1061  break;
1062  }
1063  }
1064  }
1065 
1066  if (delta == 1)
1067  {
1068  // find g such that g is a quadratic residue mod p, then g has order q
1069  // g=4 always works, but this way we get the smallest quadratic residue (other than 1)
1070  for (g=2; Jacobi(g, p) != 1; ++g) {}
1071  // contributed by Walt Tuvell: g should be the following according to the Law of Quadratic Reciprocity
1072  CRYPTOPP_ASSERT((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
1073  }
1074  else
1075  {
1076  CRYPTOPP_ASSERT(delta == -1);
1077  // find g such that g*g-4 is a quadratic non-residue,
1078  // and such that g has order q
1079  for (g=3; ; ++g)
1080  if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
1081  break;
1082  }
1083  }
1084  else
1085  {
1086  Integer minQ = Integer::Power2(qbits-1);
1087  Integer maxQ = Integer::Power2(qbits) - 1;
1088  Integer minP = Integer::Power2(pbits-1);
1089  Integer maxP = Integer::Power2(pbits) - 1;
1090 
1091  do
1092  {
1093  q.Randomize(rng, minQ, maxQ, Integer::PRIME);
1094  } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
1095 
1096  // find a random g of order q
1097  if (delta==1)
1098  {
1099  do
1100  {
1101  Integer h(rng, 2, p-2, Integer::ANY);
1102  g = a_exp_b_mod_c(h, (p-1)/q, p);
1103  } while (g <= 1);
1104  CRYPTOPP_ASSERT(a_exp_b_mod_c(g, q, p)==1);
1105  }
1106  else
1107  {
1108  CRYPTOPP_ASSERT(delta==-1);
1109  do
1110  {
1111  Integer h(rng, 3, p-1, Integer::ANY);
1112  if (Jacobi(h*h-4, p)==1)
1113  continue;
1114  g = Lucas((p+1)/q, h, p);
1115  } while (g <= 2);
1116  CRYPTOPP_ASSERT(Lucas(q, g, p) == 2);
1117  }
1118  }
1119 }
1120 
1121 NAMESPACE_END
1122 
1123 #endif
An invalid argument was detected.
Definition: cryptlib.h:200
unsigned int WordCount() const
Determines the number of words required to represent the Integer.
Definition: integer.cpp:3301
Classes for working with NameValuePairs.
bool SafeConvert(T1 from, T2 &to)
Tests whether a conversion from -> to is safe to perform.
Definition: misc.h:557
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:291
bool IsSquare() const
Determine whether this integer is a perfect square.
Definition: integer.cpp:4352
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
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:281
signed long ConvertToLong() const
Convert the Integer to Long.
Definition: integer.cpp:2984
const Integer & Subtract(const Integer &a, const Integer &b) const
Subtracts elements in the ring.
Definition: integer.cpp:4516
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:1327
void Randomize(RandomNumberGenerator &rng, size_t bitCount)
Set this Integer to random integer.
Definition: integer.cpp:3472
Integer InverseMod(const Integer &n) const
calculate multiplicative inverse of *this mod n
Definition: integer.cpp:4384
bool IsPositive() const
Determines if the Integer is positive.
Definition: integer.h:342
static const Integer & One()
Integer representing 1.
Definition: integer.cpp:4806
Integer ConvertIn(const Integer &a) const
Reduces an element in the congruence class.
Definition: modarith.h:292
Pointer that overloads operator ->
Definition: smartptr.h:37
Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
Generates a provable prime.
Definition: nbtheory.cpp:470
a number with no special properties
Definition: integer.h:93
AlgorithmParameters MakeParameters(const char *name, const T &value, bool throwIfNotUsed=true)
Create an object that implements NameValuePairs.
Definition: algparam.h:509
bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level=1)
Verifies a prime number.
Definition: nbtheory.cpp:247
Application callback to signal suitability of a cabdidate prime.
Definition: nbtheory.h:79
static Integer Power2(size_t e)
Exponentiates to a power of 2.
Definition: integer.cpp:3051
Multiple precision integer with arithmetic operations.
Definition: integer.h:49
const T1 UnsignedMin(const T1 &a, const T2 &b)
Safe comparison of values that could be neagtive and incorrectly promoted.
Definition: misc.h:541
const Integer & Multiply(const Integer &a, const Integer &b) const
Multiplies elements in the ring.
Definition: integer.cpp:4597
bool IsPrime(const Integer &p)
Verifies a prime number.
Definition: nbtheory.cpp:237
static const Integer & Two()
Integer representing 2.
Definition: integer.cpp:4818
const Integer & Square(const Integer &a) const
Square an element in the ring.
Definition: integer.cpp:4610
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:503
#define CRYPTOPP_ASSERT(exp)
Debugging and diagnostic assertion.
Definition: trap.h:60
bool TrialDivision(const Integer &p, unsigned bound)
Definition: nbtheory.cpp:71
unsigned int BitCount() const
Determines the number of bits required to represent the Integer.
Definition: integer.cpp:3315
Classes and functions for number theoretic operations.
Integer Squared() const
Multiply this integer by itself.
Definition: integer.h:576
Performs modular arithmetic in Montgomery representation for increased speed.
Definition: modarith.h:274
An object that implements NameValuePairs.
Definition: algparam.h:426
void Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned qbits)
Generate a Prime and Generator.
Definition: nbtheory.cpp:1036
Multiple precision integer with arithmetic operations.
static const Integer & Zero()
Integer representing 0.
Definition: integer.cpp:4794
Class file for performing modular arithmetic.
Crypto++ library namespace.
bool GetBit(size_t i) const
Provides the i-th bit of the Integer.
Definition: integer.cpp:3075
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:4623