5 #ifndef CRYPTOPP_IMPORTS
19 NAMESPACE_BEGIN(CryptoPP)
21 const word s_lastSmallPrime = 32719;
25 std::vector<word16> * operator()()
const
27 const unsigned int maxPrimeTableSize = 3511;
29 std::auto_ptr<std::vector<word16> > pPrimeTable(
new std::vector<word16>);
30 std::vector<word16> &primeTable = *pPrimeTable;
31 primeTable.reserve(maxPrimeTableSize);
33 primeTable.push_back(2);
34 unsigned int testEntriesEnd = 1;
36 for (
unsigned int p=3; p<=s_lastSmallPrime; p+=2)
39 for (j=1; j<testEntriesEnd; j++)
40 if (p%primeTable[j] == 0)
42 if (j == testEntriesEnd)
44 primeTable.push_back(p);
45 testEntriesEnd = UnsignedMin(54U, primeTable.size());
49 return pPrimeTable.release();
53 const word16 * GetPrimeTable(
unsigned int &size)
56 size = (
unsigned int)primeTable.size();
57 return &primeTable[0];
60 bool IsSmallPrime(
const Integer &p)
62 unsigned int primeTableSize;
63 const word16 * primeTable = GetPrimeTable(primeTableSize);
65 if (p.IsPositive() && p <= primeTable[primeTableSize-1])
66 return std::binary_search(primeTable, primeTable+primeTableSize, (word16)p.
ConvertToLong());
71 bool TrialDivision(
const Integer &p,
unsigned bound)
73 unsigned int primeTableSize;
74 const word16 * primeTable = GetPrimeTable(primeTableSize);
76 assert(primeTable[primeTableSize-1] >= bound);
79 for (i = 0; primeTable[i]<bound; i++)
80 if ((p % primeTable[i]) == 0)
83 if (bound == primeTable[i])
84 return (p % bound == 0);
89 bool SmallDivisorsTest(
const Integer &p)
91 unsigned int primeTableSize;
92 const word16 * primeTable = GetPrimeTable(primeTableSize);
93 return !TrialDivision(p, primeTable[primeTableSize-1]);
101 assert(n>3 && b>1 && b<n-1);
102 return a_exp_b_mod_c(b, n-1, n)==1;
110 assert(n>3 && b>1 && b<n-1);
112 if ((n.IsEven() && n!=2) || GCD(b, n) != 1)
124 Integer z = a_exp_b_mod_c(b, m, n);
125 if (z==1 || z==nminus1)
127 for (
unsigned j=1; j<a; j++)
146 for (
unsigned int i=0; i<rounds; i++)
148 b.Randomize(rng, 2, n-2);
149 if (!IsStrongProbablePrime(n, b))
155 bool IsLucasProbablePrime(
const Integer &n)
169 while ((j=Jacobi(b.Squared()-4, n)) == 1)
179 return Lucas(n+1, b, n)==2;
182 bool IsStrongLucasProbablePrime(
const Integer &n)
196 while ((j=Jacobi(b.Squared()-4, n)) == 1)
220 z = (z.Squared()-2)%n;
239 if (p <= s_lastSmallPrime)
240 return IsSmallPrime(p);
242 return SmallDivisorsTest(p);
244 return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p);
249 bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1);
251 pass = pass && RabinMillerTest(rng, p, 10);
255 unsigned int PrimeSearchInterval(
const Integer &max)
260 static inline bool FastProbablePrimeTest(
const Integer &n)
262 return IsStrongProbablePrime(n,2);
267 if (productBitLength < 16)
272 if (productBitLength%2==0)
274 minP =
Integer(182) << (productBitLength/2-8);
280 maxP =
Integer(181) << ((productBitLength+1)/2-8);
283 return MakeParameters(
"RandomNumberType", Integer::PRIME)(
"Min", minP)(
"Max", maxP);
291 bool NextCandidate(
Integer &c);
294 static void SieveSingle(std::vector<bool> &sieve, word16 p,
const Integer &first,
const Integer &step, word16 stepInv);
296 Integer m_first, m_last, m_step;
299 std::vector<bool> m_sieve;
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)
308 bool PrimeSieve::NextCandidate(
Integer &c)
310 bool safe = SafeConvert(std::find(m_sieve.begin()+m_next, m_sieve.end(),
false) - m_sieve.begin(), m_next);
312 if (m_next == m_sieve.size())
314 m_first += long(m_sieve.size())*m_step;
315 if (m_first > m_last)
321 return NextCandidate(c);
326 c = m_first + long(m_next)*m_step;
332 void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word16 p,
const Integer &first,
const Integer &step, word16 stepInv)
336 size_t sieveSize = sieve.size();
337 size_t j = (word32(p-(first%p))*stepInv) % p;
339 if (first.
WordCount() <= 1 && first + step*long(j) == p)
341 for (; j < sieveSize; j += p)
346 void PrimeSieve::DoSieve()
348 unsigned int primeTableSize;
349 const word16 * primeTable = GetPrimeTable(primeTableSize);
351 const unsigned int maxSieveSize = 32768;
352 unsigned int sieveSize = STDMIN(
Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong();
355 m_sieve.resize(sieveSize,
false);
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]));
365 Integer qFirst = (m_first-m_delta) >> 1;
366 Integer halfStep = m_step >> 1;
367 for (
unsigned int i = 0; i < primeTableSize; ++i)
369 word16 p = primeTable[i];
370 word16 stepInv = (word16)m_step.
InverseMod(p);
371 SieveSingle(m_sieve, p, m_first, m_step, stepInv);
373 word16 halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p;
374 SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
381 assert(!equiv.IsNegative() && equiv < mod);
387 if (p <= gcd && gcd <= max && IsPrime(gcd) && (!pSelector || pSelector->IsAcceptable(gcd)))
396 unsigned int primeTableSize;
397 const word16 * primeTable = GetPrimeTable(primeTableSize);
399 if (p <= primeTable[primeTableSize-1])
405 pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.
ConvertToLong());
409 while (pItr < primeTable+primeTableSize && !(*pItr%mod == equiv && (!pSelector || pSelector->IsAcceptable(*pItr))))
412 if (pItr < primeTable+primeTableSize)
418 p = primeTable[primeTableSize-1]+1;
421 assert(p > primeTable[primeTableSize-1]);
424 return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector);
433 while (sieve.NextCandidate(p))
435 if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p))
454 if (((r%q).Squared()-4*(r/q)).IsSquare())
457 unsigned int primeTableSize;
458 const word16 * primeTable = GetPrimeTable(primeTableSize);
460 assert(primeTableSize >= 50);
461 for (
int i=0; i<50; i++)
463 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
465 return a_exp_b_mod_c(b, q, p) == 1;
476 if (maxP <=
Integer(s_lastSmallPrime).Squared())
479 p.Randomize(rng, minP, maxP, Integer::PRIME);
483 unsigned int qbits = (pbits+2)/3 + 1 + rng.
GenerateWord32(0, pbits/36);
484 Integer q = MihailescuProvablePrime(rng, qbits);
496 p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2);
497 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
499 while (sieve.NextCandidate(p))
501 if (FastProbablePrimeTest(p) && ProvePrime(p, q))
512 const unsigned smallPrimeBound = 29, c_opt=10;
515 unsigned int primeTableSize;
516 const word16 * primeTable = GetPrimeTable(primeTableSize);
518 if (bits < smallPrimeBound)
522 while (TrialDivision(p, 1 << ((bits+1)/2)));
526 const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
529 relativeSize = pow(2.0,
double(rng.
GenerateWord32())/0xffffffff - 1);
530 while (bits * relativeSize >= bits - margin);
533 Integer q = MaurerProvablePrime(rng,
unsigned(bits*relativeSize));
536 unsigned int trialDivisorBound = (
unsigned int)STDMIN((
unsigned long)primeTable[primeTableSize-1], (
unsigned long)bits*bits/c_opt);
537 bool success =
false;
540 p.Randomize(rng, I, I2, Integer::ANY);
541 p *= q; p <<= 1; ++p;
542 if (!TrialDivision(p, trialDivisorBound))
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);
556 return p * (u * (xq-xp) % q) + xp;
575 return a_exp_b_mod_c(a, (p+1)/4, p);
586 while (Jacobi(n, p) != -1)
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;
609 for (
unsigned i=0; i<r-m-1; i++)
617 assert(x.Squared()%p == a);
623 Integer D = (b.Squared() - 4*a*c) % p;
624 switch (Jacobi(D, p))
632 r1 = r2 = (-b*(a+a).InverseMod(p)) % p;
633 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
636 Integer s = ModularSquareRoot(D, p);
637 Integer t = (a+a).InverseMod(p);
640 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
641 assert(((r2.Squared()*a + r2*b + c) % p).IsZero());
654 p2 = ModularExponentiation((a % p), dp, p);
656 q2 = ModularExponentiation((a % q), dq, q);
658 return CRT(p2, p, q2, q, u);
664 Integer dp = EuclideanMultiplicativeInverse(e, p-1);
665 Integer dq = EuclideanMultiplicativeInverse(e, q-1);
666 Integer u = EuclideanMultiplicativeInverse(p, q);
667 assert(!!dp && !!dq && !!u);
668 return ModularRoot(a, dp, dq, p, q, u);
795 while (a.GetBit(i)==0)
799 if (i%2==1 && (b%8==3 || b%8==5))
802 if (a%4==3 && b%4==3)
809 return (b==1) ? result : 0;
820 Integer v=p, v1=m.Subtract(m.Square(p), two);
828 v = m.Subtract(m.Multiply(v,v1), p);
830 v1 = m.Subtract(m.Square(v1), two);
835 v1 = m.Subtract(m.Multiply(v,v1), p);
837 v = m.Subtract(m.Square(v), two);
840 return m.ConvertOut(v);
1002 #pragma omp parallel
1003 #pragma omp sections
1008 p2 = Lucas(EuclideanMultiplicativeInverse(e,p2), m, p);
1013 q2 = Lucas(EuclideanMultiplicativeInverse(e,q2), m, q);
1016 return CRT(p2, p, q2, q, u);
1019 unsigned int FactoringWorkFactor(
unsigned int n)
1024 else return (
unsigned int)(2.4 * pow((
double)n, 1.0/3.0) * pow(log(
double(n)), 2.0/3.0) - 5);
1027 unsigned int DiscreteLogWorkFactor(
unsigned int n)
1031 else return (
unsigned int)(2.4 * pow((
double)n, 1.0/3.0) * pow(log(
double(n)), 2.0/3.0) - 5);
1036 void PrimeAndGenerator::Generate(
signed int delta,
RandomNumberGenerator &rng,
unsigned int pbits,
unsigned int qbits)
1040 assert(pbits > qbits);
1042 if (qbits+1 == pbits)
1046 bool success =
false;
1050 p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
1051 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
1053 while (sieve.NextCandidate(p))
1055 assert(IsSmallPrime(p) || SmallDivisorsTest(p));
1057 assert(IsSmallPrime(q) || SmallDivisorsTest(q));
1058 if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
1070 for (g=2; Jacobi(g, p) != 1; ++g) {}
1072 assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
1076 assert(delta == -1);
1080 if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
1093 q.Randomize(rng, minQ, maxQ, Integer::PRIME);
1094 }
while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
1101 Integer h(rng, 2, p-2, Integer::ANY);
1102 g = a_exp_b_mod_c(h, (p-1)/q, p);
1104 assert(a_exp_b_mod_c(g, q, p)==1);
1111 Integer h(rng, 3, p-1, Integer::ANY);
1112 if (Jacobi(h*h-4, p)==1)
1114 g = Lucas((p+1)/q, h, p);
1116 assert(Lucas(q, g, p) == 2);