00001
00002
00003 #include "pch.h"
00004
00005 #ifndef CRYPTOPP_IMPORTS
00006
00007 #include "nbtheory.h"
00008 #include "modarith.h"
00009 #include "algparam.h"
00010
00011 #include <math.h>
00012 #include <vector>
00013
00014 NAMESPACE_BEGIN(CryptoPP)
00015
00016 const word s_lastSmallPrime = 32719;
00017
00018 struct NewPrimeTable
00019 {
00020 std::vector<word16> * operator()() const
00021 {
00022 const unsigned int maxPrimeTableSize = 3511;
00023
00024 std::auto_ptr<std::vector<word16> > pPrimeTable(new std::vector<word16>);
00025 std::vector<word16> &primeTable = *pPrimeTable;
00026 primeTable.reserve(maxPrimeTableSize);
00027
00028 primeTable.push_back(2);
00029 unsigned int testEntriesEnd = 1;
00030
00031 for (unsigned int p=3; p<=s_lastSmallPrime; p+=2)
00032 {
00033 unsigned int j;
00034 for (j=1; j<testEntriesEnd; j++)
00035 if (p%primeTable[j] == 0)
00036 break;
00037 if (j == testEntriesEnd)
00038 {
00039 primeTable.push_back(p);
00040 testEntriesEnd = STDMIN((size_t)54U, primeTable.size());
00041 }
00042 }
00043
00044 return pPrimeTable.release();
00045 }
00046 };
00047
00048 const word16 * GetPrimeTable(unsigned int &size)
00049 {
00050 const std::vector<word16> &primeTable = Singleton<std::vector<word16>, NewPrimeTable>().Ref();
00051 size = primeTable.size();
00052 return &primeTable[0];
00053 }
00054
00055 bool IsSmallPrime(const Integer &p)
00056 {
00057 unsigned int primeTableSize;
00058 const word16 * primeTable = GetPrimeTable(primeTableSize);
00059
00060 if (p.IsPositive() && p <= primeTable[primeTableSize-1])
00061 return std::binary_search(primeTable, primeTable+primeTableSize, (word16)p.ConvertToLong());
00062 else
00063 return false;
00064 }
00065
00066 bool TrialDivision(const Integer &p, unsigned bound)
00067 {
00068 unsigned int primeTableSize;
00069 const word16 * primeTable = GetPrimeTable(primeTableSize);
00070
00071 assert(primeTable[primeTableSize-1] >= bound);
00072
00073 unsigned int i;
00074 for (i = 0; primeTable[i]<bound; i++)
00075 if ((p % primeTable[i]) == 0)
00076 return true;
00077
00078 if (bound == primeTable[i])
00079 return (p % bound == 0);
00080 else
00081 return false;
00082 }
00083
00084 bool SmallDivisorsTest(const Integer &p)
00085 {
00086 unsigned int primeTableSize;
00087 const word16 * primeTable = GetPrimeTable(primeTableSize);
00088 return !TrialDivision(p, primeTable[primeTableSize-1]);
00089 }
00090
00091 bool IsFermatProbablePrime(const Integer &n, const Integer &b)
00092 {
00093 if (n <= 3)
00094 return n==2 || n==3;
00095
00096 assert(n>3 && b>1 && b<n-1);
00097 return a_exp_b_mod_c(b, n-1, n)==1;
00098 }
00099
00100 bool IsStrongProbablePrime(const Integer &n, const Integer &b)
00101 {
00102 if (n <= 3)
00103 return n==2 || n==3;
00104
00105 assert(n>3 && b>1 && b<n-1);
00106
00107 if ((n.IsEven() && n!=2) || GCD(b, n) != 1)
00108 return false;
00109
00110 Integer nminus1 = (n-1);
00111 unsigned int a;
00112
00113
00114 for (a=0; ; a++)
00115 if (nminus1.GetBit(a))
00116 break;
00117 Integer m = nminus1>>a;
00118
00119 Integer z = a_exp_b_mod_c(b, m, n);
00120 if (z==1 || z==nminus1)
00121 return true;
00122 for (unsigned j=1; j<a; j++)
00123 {
00124 z = z.Squared()%n;
00125 if (z==nminus1)
00126 return true;
00127 if (z==1)
00128 return false;
00129 }
00130 return false;
00131 }
00132
00133 bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds)
00134 {
00135 if (n <= 3)
00136 return n==2 || n==3;
00137
00138 assert(n>3);
00139
00140 Integer b;
00141 for (unsigned int i=0; i<rounds; i++)
00142 {
00143 b.Randomize(rng, 2, n-2);
00144 if (!IsStrongProbablePrime(n, b))
00145 return false;
00146 }
00147 return true;
00148 }
00149
00150 bool IsLucasProbablePrime(const Integer &n)
00151 {
00152 if (n <= 1)
00153 return false;
00154
00155 if (n.IsEven())
00156 return n==2;
00157
00158 assert(n>2);
00159
00160 Integer b=3;
00161 unsigned int i=0;
00162 int j;
00163
00164 while ((j=Jacobi(b.Squared()-4, n)) == 1)
00165 {
00166 if (++i==64 && n.IsSquare())
00167 return false;
00168 ++b; ++b;
00169 }
00170
00171 if (j==0)
00172 return false;
00173 else
00174 return Lucas(n+1, b, n)==2;
00175 }
00176
00177 bool IsStrongLucasProbablePrime(const Integer &n)
00178 {
00179 if (n <= 1)
00180 return false;
00181
00182 if (n.IsEven())
00183 return n==2;
00184
00185 assert(n>2);
00186
00187 Integer b=3;
00188 unsigned int i=0;
00189 int j;
00190
00191 while ((j=Jacobi(b.Squared()-4, n)) == 1)
00192 {
00193 if (++i==64 && n.IsSquare())
00194 return false;
00195 ++b; ++b;
00196 }
00197
00198 if (j==0)
00199 return false;
00200
00201 Integer n1 = n+1;
00202 unsigned int a;
00203
00204
00205 for (a=0; ; a++)
00206 if (n1.GetBit(a))
00207 break;
00208 Integer m = n1>>a;
00209
00210 Integer z = Lucas(m, b, n);
00211 if (z==2 || z==n-2)
00212 return true;
00213 for (i=1; i<a; i++)
00214 {
00215 z = (z.Squared()-2)%n;
00216 if (z==n-2)
00217 return true;
00218 if (z==2)
00219 return false;
00220 }
00221 return false;
00222 }
00223
00224 struct NewLastSmallPrimeSquared
00225 {
00226 Integer * operator()() const
00227 {
00228 return new Integer(Integer(s_lastSmallPrime).Squared());
00229 }
00230 };
00231
00232 bool IsPrime(const Integer &p)
00233 {
00234 if (p <= s_lastSmallPrime)
00235 return IsSmallPrime(p);
00236 else if (p <= Singleton<Integer, NewLastSmallPrimeSquared>().Ref())
00237 return SmallDivisorsTest(p);
00238 else
00239 return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p);
00240 }
00241
00242 bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level)
00243 {
00244 bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1);
00245 if (level >= 1)
00246 pass = pass && RabinMillerTest(rng, p, 10);
00247 return pass;
00248 }
00249
00250 unsigned int PrimeSearchInterval(const Integer &max)
00251 {
00252 return max.BitCount();
00253 }
00254
00255 static inline bool FastProbablePrimeTest(const Integer &n)
00256 {
00257 return IsStrongProbablePrime(n,2);
00258 }
00259
00260 AlgorithmParameters<AlgorithmParameters<AlgorithmParameters<NullNameValuePairs, Integer::RandomNumberType>, Integer>, Integer>
00261 MakeParametersForTwoPrimesOfEqualSize(unsigned int productBitLength)
00262 {
00263 if (productBitLength < 16)
00264 throw InvalidArgument("invalid bit length");
00265
00266 Integer minP, maxP;
00267
00268 if (productBitLength%2==0)
00269 {
00270 minP = Integer(182) << (productBitLength/2-8);
00271 maxP = Integer::Power2(productBitLength/2)-1;
00272 }
00273 else
00274 {
00275 minP = Integer::Power2((productBitLength-1)/2);
00276 maxP = Integer(181) << ((productBitLength+1)/2-8);
00277 }
00278
00279 return MakeParameters("RandomNumberType", Integer::PRIME)("Min", minP)("Max", maxP);
00280 }
00281
00282 class PrimeSieve
00283 {
00284 public:
00285
00286 PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta=0);
00287 bool NextCandidate(Integer &c);
00288
00289 void DoSieve();
00290 static void SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv);
00291
00292 Integer m_first, m_last, m_step;
00293 signed int m_delta;
00294 word m_next;
00295 std::vector<bool> m_sieve;
00296 };
00297
00298 PrimeSieve::PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta)
00299 : m_first(first), m_last(last), m_step(step), m_delta(delta), m_next(0)
00300 {
00301 DoSieve();
00302 }
00303
00304 bool PrimeSieve::NextCandidate(Integer &c)
00305 {
00306 m_next = std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin();
00307 if (m_next == m_sieve.size())
00308 {
00309 m_first += m_sieve.size()*m_step;
00310 if (m_first > m_last)
00311 return false;
00312 else
00313 {
00314 m_next = 0;
00315 DoSieve();
00316 return NextCandidate(c);
00317 }
00318 }
00319 else
00320 {
00321 c = m_first + m_next*m_step;
00322 ++m_next;
00323 return true;
00324 }
00325 }
00326
00327 void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv)
00328 {
00329 if (stepInv)
00330 {
00331 unsigned int sieveSize = sieve.size();
00332 word j = word((word32(p-(first%p))*stepInv) % p);
00333
00334 if (first.WordCount() <= 1 && first + step*j == p)
00335 j += p;
00336 for (; j < sieveSize; j += p)
00337 sieve[j] = true;
00338 }
00339 }
00340
00341 void PrimeSieve::DoSieve()
00342 {
00343 unsigned int primeTableSize;
00344 const word16 * primeTable = GetPrimeTable(primeTableSize);
00345
00346 const unsigned int maxSieveSize = 32768;
00347 unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong();
00348
00349 m_sieve.clear();
00350 m_sieve.resize(sieveSize, false);
00351
00352 if (m_delta == 0)
00353 {
00354 for (unsigned int i = 0; i < primeTableSize; ++i)
00355 SieveSingle(m_sieve, primeTable[i], m_first, m_step, m_step.InverseMod(primeTable[i]));
00356 }
00357 else
00358 {
00359 assert(m_step%2==0);
00360 Integer qFirst = (m_first-m_delta) >> 1;
00361 Integer halfStep = m_step >> 1;
00362 for (unsigned int i = 0; i < primeTableSize; ++i)
00363 {
00364 word16 p = primeTable[i];
00365 word16 stepInv = m_step.InverseMod(p);
00366 SieveSingle(m_sieve, p, m_first, m_step, stepInv);
00367
00368 word16 halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p;
00369 SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
00370 }
00371 }
00372 }
00373
00374 bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector)
00375 {
00376 assert(!equiv.IsNegative() && equiv < mod);
00377
00378 Integer gcd = GCD(equiv, mod);
00379 if (gcd != Integer::One())
00380 {
00381
00382 if (p <= gcd && gcd <= max && IsPrime(gcd) && (!pSelector || pSelector->IsAcceptable(gcd)))
00383 {
00384 p = gcd;
00385 return true;
00386 }
00387 else
00388 return false;
00389 }
00390
00391 unsigned int primeTableSize;
00392 const word16 * primeTable = GetPrimeTable(primeTableSize);
00393
00394 if (p <= primeTable[primeTableSize-1])
00395 {
00396 const word16 *pItr;
00397
00398 --p;
00399 if (p.IsPositive())
00400 pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
00401 else
00402 pItr = primeTable;
00403
00404 while (pItr < primeTable+primeTableSize && !(*pItr%mod == equiv && (!pSelector || pSelector->IsAcceptable(*pItr))))
00405 ++pItr;
00406
00407 if (pItr < primeTable+primeTableSize)
00408 {
00409 p = *pItr;
00410 return p <= max;
00411 }
00412
00413 p = primeTable[primeTableSize-1]+1;
00414 }
00415
00416 assert(p > primeTable[primeTableSize-1]);
00417
00418 if (mod.IsOdd())
00419 return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector);
00420
00421 p += (equiv-p)%mod;
00422
00423 if (p>max)
00424 return false;
00425
00426 PrimeSieve sieve(p, max, mod);
00427
00428 while (sieve.NextCandidate(p))
00429 {
00430 if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p))
00431 return true;
00432 }
00433
00434 return false;
00435 }
00436
00437
00438 static bool ProvePrime(const Integer &p, const Integer &q)
00439 {
00440 assert(p < q*q*q);
00441 assert(p % q == 1);
00442
00443
00444
00445
00446
00447
00448 Integer r = (p-1)/q;
00449 if (((r%q).Squared()-4*(r/q)).IsSquare())
00450 return false;
00451
00452 unsigned int primeTableSize;
00453 const word16 * primeTable = GetPrimeTable(primeTableSize);
00454
00455 assert(primeTableSize >= 50);
00456 for (int i=0; i<50; i++)
00457 {
00458 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
00459 if (b != 1)
00460 return a_exp_b_mod_c(b, q, p) == 1;
00461 }
00462 return false;
00463 }
00464
00465 Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits)
00466 {
00467 Integer p;
00468 Integer minP = Integer::Power2(pbits-1);
00469 Integer maxP = Integer::Power2(pbits) - 1;
00470
00471 if (maxP <= Integer(s_lastSmallPrime).Squared())
00472 {
00473
00474 p.Randomize(rng, minP, maxP, Integer::PRIME);
00475 return p;
00476 }
00477
00478 unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36);
00479 Integer q = MihailescuProvablePrime(rng, qbits);
00480 Integer q2 = q<<1;
00481
00482 while (true)
00483 {
00484
00485
00486
00487
00488
00489
00490
00491 p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2);
00492 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
00493
00494 while (sieve.NextCandidate(p))
00495 {
00496 if (FastProbablePrimeTest(p) && ProvePrime(p, q))
00497 return p;
00498 }
00499 }
00500
00501
00502 return p;
00503 }
00504
00505 Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
00506 {
00507 const unsigned smallPrimeBound = 29, c_opt=10;
00508 Integer p;
00509
00510 unsigned int primeTableSize;
00511 const word16 * primeTable = GetPrimeTable(primeTableSize);
00512
00513 if (bits < smallPrimeBound)
00514 {
00515 do
00516 p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2);
00517 while (TrialDivision(p, 1 << ((bits+1)/2)));
00518 }
00519 else
00520 {
00521 const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
00522 double relativeSize;
00523 do
00524 relativeSize = pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1);
00525 while (bits * relativeSize >= bits - margin);
00526
00527 Integer a,b;
00528 Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize));
00529 Integer I = Integer::Power2(bits-2)/q;
00530 Integer I2 = I << 1;
00531 unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt);
00532 bool success = false;
00533 while (!success)
00534 {
00535 p.Randomize(rng, I, I2, Integer::ANY);
00536 p *= q; p <<= 1; ++p;
00537 if (!TrialDivision(p, trialDivisorBound))
00538 {
00539 a.Randomize(rng, 2, p-1, Integer::ANY);
00540 b = a_exp_b_mod_c(a, (p-1)/q, p);
00541 success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1);
00542 }
00543 }
00544 }
00545 return p;
00546 }
00547
00548 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u)
00549 {
00550
00551 return p * (u * (xq-xp) % q) + xp;
00552 }
00553
00554 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q)
00555 {
00556 return CRT(xp, p, xq, q, EuclideanMultiplicativeInverse(p, q));
00557 }
00558
00559 Integer ModularSquareRoot(const Integer &a, const Integer &p)
00560 {
00561 if (p%4 == 3)
00562 return a_exp_b_mod_c(a, (p+1)/4, p);
00563
00564 Integer q=p-1;
00565 unsigned int r=0;
00566 while (q.IsEven())
00567 {
00568 r++;
00569 q >>= 1;
00570 }
00571
00572 Integer n=2;
00573 while (Jacobi(n, p) != -1)
00574 ++n;
00575
00576 Integer y = a_exp_b_mod_c(n, q, p);
00577 Integer x = a_exp_b_mod_c(a, (q-1)/2, p);
00578 Integer b = (x.Squared()%p)*a%p;
00579 x = a*x%p;
00580 Integer tempb, t;
00581
00582 while (b != 1)
00583 {
00584 unsigned m=0;
00585 tempb = b;
00586 do
00587 {
00588 m++;
00589 b = b.Squared()%p;
00590 if (m==r)
00591 return Integer::Zero();
00592 }
00593 while (b != 1);
00594
00595 t = y;
00596 for (unsigned i=0; i<r-m-1; i++)
00597 t = t.Squared()%p;
00598 y = t.Squared()%p;
00599 r = m;
00600 x = x*t%p;
00601 b = tempb*y%p;
00602 }
00603
00604 assert(x.Squared()%p == a);
00605 return x;
00606 }
00607
00608 bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p)
00609 {
00610 Integer D = (b.Squared() - 4*a*c) % p;
00611 switch (Jacobi(D, p))
00612 {
00613 default:
00614 assert(false);
00615 return false;
00616 case -1:
00617 return false;
00618 case 0:
00619 r1 = r2 = (-b*(a+a).InverseMod(p)) % p;
00620 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
00621 return true;
00622 case 1:
00623 Integer s = ModularSquareRoot(D, p);
00624 Integer t = (a+a).InverseMod(p);
00625 r1 = (s-b)*t % p;
00626 r2 = (-s-b)*t % p;
00627 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
00628 assert(((r2.Squared()*a + r2*b + c) % p).IsZero());
00629 return true;
00630 }
00631 }
00632
00633 Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq,
00634 const Integer &p, const Integer &q, const Integer &u)
00635 {
00636 Integer p2 = ModularExponentiation((a % p), dp, p);
00637 Integer q2 = ModularExponentiation((a % q), dq, q);
00638 return CRT(p2, p, q2, q, u);
00639 }
00640
00641 Integer ModularRoot(const Integer &a, const Integer &e,
00642 const Integer &p, const Integer &q)
00643 {
00644 Integer dp = EuclideanMultiplicativeInverse(e, p-1);
00645 Integer dq = EuclideanMultiplicativeInverse(e, q-1);
00646 Integer u = EuclideanMultiplicativeInverse(p, q);
00647 assert(!!dp && !!dq && !!u);
00648 return ModularRoot(a, dp, dq, p, q, u);
00649 }
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765 int Jacobi(const Integer &aIn, const Integer &bIn)
00766 {
00767 assert(bIn.IsOdd());
00768
00769 Integer b = bIn, a = aIn%bIn;
00770 int result = 1;
00771
00772 while (!!a)
00773 {
00774 unsigned i=0;
00775 while (a.GetBit(i)==0)
00776 i++;
00777 a>>=i;
00778
00779 if (i%2==1 && (b%8==3 || b%8==5))
00780 result = -result;
00781
00782 if (a%4==3 && b%4==3)
00783 result = -result;
00784
00785 std::swap(a, b);
00786 a %= b;
00787 }
00788
00789 return (b==1) ? result : 0;
00790 }
00791
00792 Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n)
00793 {
00794 unsigned i = e.BitCount();
00795 if (i==0)
00796 return Integer::Two();
00797
00798 MontgomeryRepresentation m(n);
00799 Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two());
00800 Integer v=p, v1=m.Subtract(m.Square(p), two);
00801
00802 i--;
00803 while (i--)
00804 {
00805 if (e.GetBit(i))
00806 {
00807
00808 v = m.Subtract(m.Multiply(v,v1), p);
00809
00810 v1 = m.Subtract(m.Square(v1), two);
00811 }
00812 else
00813 {
00814
00815 v1 = m.Subtract(m.Multiply(v,v1), p);
00816
00817 v = m.Subtract(m.Square(v), two);
00818 }
00819 }
00820 return m.ConvertOut(v);
00821 }
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
00979 {
00980 Integer d = (m*m-4);
00981 Integer p2 = p-Jacobi(d,p);
00982 Integer q2 = q-Jacobi(d,q);
00983 return CRT(Lucas(EuclideanMultiplicativeInverse(e,p2), m, p), p, Lucas(EuclideanMultiplicativeInverse(e,q2), m, q), q, u);
00984 }
00985
00986 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q)
00987 {
00988 return InverseLucas(e, m, p, q, EuclideanMultiplicativeInverse(p, q));
00989 }
00990
00991 unsigned int FactoringWorkFactor(unsigned int n)
00992 {
00993
00994
00995 if (n<5) return 0;
00996 else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
00997 }
00998
00999 unsigned int DiscreteLogWorkFactor(unsigned int n)
01000 {
01001
01002 if (n<5) return 0;
01003 else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
01004 }
01005
01006
01007
01008 void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits)
01009 {
01010
01011 assert(qbits > 4);
01012 assert(pbits > qbits);
01013
01014 if (qbits+1 == pbits)
01015 {
01016 Integer minP = Integer::Power2(pbits-1);
01017 Integer maxP = Integer::Power2(pbits) - 1;
01018 bool success = false;
01019
01020 while (!success)
01021 {
01022 p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
01023 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
01024
01025 while (sieve.NextCandidate(p))
01026 {
01027 assert(IsSmallPrime(p) || SmallDivisorsTest(p));
01028 q = (p-delta) >> 1;
01029 assert(IsSmallPrime(q) || SmallDivisorsTest(q));
01030 if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
01031 {
01032 success = true;
01033 break;
01034 }
01035 }
01036 }
01037
01038 if (delta == 1)
01039 {
01040
01041
01042 for (g=2; Jacobi(g, p) != 1; ++g) {}
01043
01044 assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
01045 }
01046 else
01047 {
01048 assert(delta == -1);
01049
01050
01051 for (g=3; ; ++g)
01052 if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
01053 break;
01054 }
01055 }
01056 else
01057 {
01058 Integer minQ = Integer::Power2(qbits-1);
01059 Integer maxQ = Integer::Power2(qbits) - 1;
01060 Integer minP = Integer::Power2(pbits-1);
01061 Integer maxP = Integer::Power2(pbits) - 1;
01062
01063 do
01064 {
01065 q.Randomize(rng, minQ, maxQ, Integer::PRIME);
01066 } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
01067
01068
01069 if (delta==1)
01070 {
01071 do
01072 {
01073 Integer h(rng, 2, p-2, Integer::ANY);
01074 g = a_exp_b_mod_c(h, (p-1)/q, p);
01075 } while (g <= 1);
01076 assert(a_exp_b_mod_c(g, q, p)==1);
01077 }
01078 else
01079 {
01080 assert(delta==-1);
01081 do
01082 {
01083 Integer h(rng, 3, p-1, Integer::ANY);
01084 if (Jacobi(h*h-4, p)==1)
01085 continue;
01086 g = Lucas((p+1)/q, h, p);
01087 } while (g <= 2);
01088 assert(Lucas(q, g, p) == 2);
01089 }
01090 }
01091 }
01092
01093 NAMESPACE_END
01094
01095 #endif