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 = UnsignedMin(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 = (unsigned int)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 bool safe = SafeConvert(std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin(), m_next);
00307 assert(safe);
00308 if (m_next == m_sieve.size())
00309 {
00310 m_first += long(m_sieve.size())*m_step;
00311 if (m_first > m_last)
00312 return false;
00313 else
00314 {
00315 m_next = 0;
00316 DoSieve();
00317 return NextCandidate(c);
00318 }
00319 }
00320 else
00321 {
00322 c = m_first + long(m_next)*m_step;
00323 ++m_next;
00324 return true;
00325 }
00326 }
00327
00328 void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv)
00329 {
00330 if (stepInv)
00331 {
00332 size_t sieveSize = sieve.size();
00333 size_t j = (word32(p-(first%p))*stepInv) % p;
00334
00335 if (first.WordCount() <= 1 && first + step*long(j) == p)
00336 j += p;
00337 for (; j < sieveSize; j += p)
00338 sieve[j] = true;
00339 }
00340 }
00341
00342 void PrimeSieve::DoSieve()
00343 {
00344 unsigned int primeTableSize;
00345 const word16 * primeTable = GetPrimeTable(primeTableSize);
00346
00347 const unsigned int maxSieveSize = 32768;
00348 unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong();
00349
00350 m_sieve.clear();
00351 m_sieve.resize(sieveSize, false);
00352
00353 if (m_delta == 0)
00354 {
00355 for (unsigned int i = 0; i < primeTableSize; ++i)
00356 SieveSingle(m_sieve, primeTable[i], m_first, m_step, (word16)m_step.InverseMod(primeTable[i]));
00357 }
00358 else
00359 {
00360 assert(m_step%2==0);
00361 Integer qFirst = (m_first-m_delta) >> 1;
00362 Integer halfStep = m_step >> 1;
00363 for (unsigned int i = 0; i < primeTableSize; ++i)
00364 {
00365 word16 p = primeTable[i];
00366 word16 stepInv = (word16)m_step.InverseMod(p);
00367 SieveSingle(m_sieve, p, m_first, m_step, stepInv);
00368
00369 word16 halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p;
00370 SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
00371 }
00372 }
00373 }
00374
00375 bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector)
00376 {
00377 assert(!equiv.IsNegative() && equiv < mod);
00378
00379 Integer gcd = GCD(equiv, mod);
00380 if (gcd != Integer::One())
00381 {
00382
00383 if (p <= gcd && gcd <= max && IsPrime(gcd) && (!pSelector || pSelector->IsAcceptable(gcd)))
00384 {
00385 p = gcd;
00386 return true;
00387 }
00388 else
00389 return false;
00390 }
00391
00392 unsigned int primeTableSize;
00393 const word16 * primeTable = GetPrimeTable(primeTableSize);
00394
00395 if (p <= primeTable[primeTableSize-1])
00396 {
00397 const word16 *pItr;
00398
00399 --p;
00400 if (p.IsPositive())
00401 pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
00402 else
00403 pItr = primeTable;
00404
00405 while (pItr < primeTable+primeTableSize && !(*pItr%mod == equiv && (!pSelector || pSelector->IsAcceptable(*pItr))))
00406 ++pItr;
00407
00408 if (pItr < primeTable+primeTableSize)
00409 {
00410 p = *pItr;
00411 return p <= max;
00412 }
00413
00414 p = primeTable[primeTableSize-1]+1;
00415 }
00416
00417 assert(p > primeTable[primeTableSize-1]);
00418
00419 if (mod.IsOdd())
00420 return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector);
00421
00422 p += (equiv-p)%mod;
00423
00424 if (p>max)
00425 return false;
00426
00427 PrimeSieve sieve(p, max, mod);
00428
00429 while (sieve.NextCandidate(p))
00430 {
00431 if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p))
00432 return true;
00433 }
00434
00435 return false;
00436 }
00437
00438
00439 static bool ProvePrime(const Integer &p, const Integer &q)
00440 {
00441 assert(p < q*q*q);
00442 assert(p % q == 1);
00443
00444
00445
00446
00447
00448
00449 Integer r = (p-1)/q;
00450 if (((r%q).Squared()-4*(r/q)).IsSquare())
00451 return false;
00452
00453 unsigned int primeTableSize;
00454 const word16 * primeTable = GetPrimeTable(primeTableSize);
00455
00456 assert(primeTableSize >= 50);
00457 for (int i=0; i<50; i++)
00458 {
00459 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
00460 if (b != 1)
00461 return a_exp_b_mod_c(b, q, p) == 1;
00462 }
00463 return false;
00464 }
00465
00466 Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits)
00467 {
00468 Integer p;
00469 Integer minP = Integer::Power2(pbits-1);
00470 Integer maxP = Integer::Power2(pbits) - 1;
00471
00472 if (maxP <= Integer(s_lastSmallPrime).Squared())
00473 {
00474
00475 p.Randomize(rng, minP, maxP, Integer::PRIME);
00476 return p;
00477 }
00478
00479 unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36);
00480 Integer q = MihailescuProvablePrime(rng, qbits);
00481 Integer q2 = q<<1;
00482
00483 while (true)
00484 {
00485
00486
00487
00488
00489
00490
00491
00492 p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2);
00493 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
00494
00495 while (sieve.NextCandidate(p))
00496 {
00497 if (FastProbablePrimeTest(p) && ProvePrime(p, q))
00498 return p;
00499 }
00500 }
00501
00502
00503 return p;
00504 }
00505
00506 Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
00507 {
00508 const unsigned smallPrimeBound = 29, c_opt=10;
00509 Integer p;
00510
00511 unsigned int primeTableSize;
00512 const word16 * primeTable = GetPrimeTable(primeTableSize);
00513
00514 if (bits < smallPrimeBound)
00515 {
00516 do
00517 p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2);
00518 while (TrialDivision(p, 1 << ((bits+1)/2)));
00519 }
00520 else
00521 {
00522 const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
00523 double relativeSize;
00524 do
00525 relativeSize = pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1);
00526 while (bits * relativeSize >= bits - margin);
00527
00528 Integer a,b;
00529 Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize));
00530 Integer I = Integer::Power2(bits-2)/q;
00531 Integer I2 = I << 1;
00532 unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt);
00533 bool success = false;
00534 while (!success)
00535 {
00536 p.Randomize(rng, I, I2, Integer::ANY);
00537 p *= q; p <<= 1; ++p;
00538 if (!TrialDivision(p, trialDivisorBound))
00539 {
00540 a.Randomize(rng, 2, p-1, Integer::ANY);
00541 b = a_exp_b_mod_c(a, (p-1)/q, p);
00542 success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1);
00543 }
00544 }
00545 }
00546 return p;
00547 }
00548
00549 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u)
00550 {
00551
00552 return p * (u * (xq-xp) % q) + xp;
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566 }
00567
00568 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q)
00569 {
00570 return CRT(xp, p, xq, q, EuclideanMultiplicativeInverse(p, q));
00571 }
00572
00573 Integer ModularSquareRoot(const Integer &a, const Integer &p)
00574 {
00575 if (p%4 == 3)
00576 return a_exp_b_mod_c(a, (p+1)/4, p);
00577
00578 Integer q=p-1;
00579 unsigned int r=0;
00580 while (q.IsEven())
00581 {
00582 r++;
00583 q >>= 1;
00584 }
00585
00586 Integer n=2;
00587 while (Jacobi(n, p) != -1)
00588 ++n;
00589
00590 Integer y = a_exp_b_mod_c(n, q, p);
00591 Integer x = a_exp_b_mod_c(a, (q-1)/2, p);
00592 Integer b = (x.Squared()%p)*a%p;
00593 x = a*x%p;
00594 Integer tempb, t;
00595
00596 while (b != 1)
00597 {
00598 unsigned m=0;
00599 tempb = b;
00600 do
00601 {
00602 m++;
00603 b = b.Squared()%p;
00604 if (m==r)
00605 return Integer::Zero();
00606 }
00607 while (b != 1);
00608
00609 t = y;
00610 for (unsigned i=0; i<r-m-1; i++)
00611 t = t.Squared()%p;
00612 y = t.Squared()%p;
00613 r = m;
00614 x = x*t%p;
00615 b = tempb*y%p;
00616 }
00617
00618 assert(x.Squared()%p == a);
00619 return x;
00620 }
00621
00622 bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p)
00623 {
00624 Integer D = (b.Squared() - 4*a*c) % p;
00625 switch (Jacobi(D, p))
00626 {
00627 default:
00628 assert(false);
00629 return false;
00630 case -1:
00631 return false;
00632 case 0:
00633 r1 = r2 = (-b*(a+a).InverseMod(p)) % p;
00634 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
00635 return true;
00636 case 1:
00637 Integer s = ModularSquareRoot(D, p);
00638 Integer t = (a+a).InverseMod(p);
00639 r1 = (s-b)*t % p;
00640 r2 = (-s-b)*t % p;
00641 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
00642 assert(((r2.Squared()*a + r2*b + c) % p).IsZero());
00643 return true;
00644 }
00645 }
00646
00647 Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq,
00648 const Integer &p, const Integer &q, const Integer &u)
00649 {
00650 Integer p2 = ModularExponentiation((a % p), dp, p);
00651 Integer q2 = ModularExponentiation((a % q), dq, q);
00652 return CRT(p2, p, q2, q, u);
00653 }
00654
00655 Integer ModularRoot(const Integer &a, const Integer &e,
00656 const Integer &p, const Integer &q)
00657 {
00658 Integer dp = EuclideanMultiplicativeInverse(e, p-1);
00659 Integer dq = EuclideanMultiplicativeInverse(e, q-1);
00660 Integer u = EuclideanMultiplicativeInverse(p, q);
00661 assert(!!dp && !!dq && !!u);
00662 return ModularRoot(a, dp, dq, p, q, u);
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
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779 int Jacobi(const Integer &aIn, const Integer &bIn)
00780 {
00781 assert(bIn.IsOdd());
00782
00783 Integer b = bIn, a = aIn%bIn;
00784 int result = 1;
00785
00786 while (!!a)
00787 {
00788 unsigned i=0;
00789 while (a.GetBit(i)==0)
00790 i++;
00791 a>>=i;
00792
00793 if (i%2==1 && (b%8==3 || b%8==5))
00794 result = -result;
00795
00796 if (a%4==3 && b%4==3)
00797 result = -result;
00798
00799 std::swap(a, b);
00800 a %= b;
00801 }
00802
00803 return (b==1) ? result : 0;
00804 }
00805
00806 Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n)
00807 {
00808 unsigned i = e.BitCount();
00809 if (i==0)
00810 return Integer::Two();
00811
00812 MontgomeryRepresentation m(n);
00813 Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two());
00814 Integer v=p, v1=m.Subtract(m.Square(p), two);
00815
00816 i--;
00817 while (i--)
00818 {
00819 if (e.GetBit(i))
00820 {
00821
00822 v = m.Subtract(m.Multiply(v,v1), p);
00823
00824 v1 = m.Subtract(m.Square(v1), two);
00825 }
00826 else
00827 {
00828
00829 v1 = m.Subtract(m.Multiply(v,v1), p);
00830
00831 v = m.Subtract(m.Square(v), two);
00832 }
00833 }
00834 return m.ConvertOut(v);
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
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
00993 {
00994 Integer d = (m*m-4);
00995 Integer p2 = p-Jacobi(d,p);
00996 Integer q2 = q-Jacobi(d,q);
00997 return CRT(Lucas(EuclideanMultiplicativeInverse(e,p2), m, p), p, Lucas(EuclideanMultiplicativeInverse(e,q2), m, q), q, u);
00998 }
00999
01000 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q)
01001 {
01002 return InverseLucas(e, m, p, q, EuclideanMultiplicativeInverse(p, q));
01003 }
01004
01005 unsigned int FactoringWorkFactor(unsigned int n)
01006 {
01007
01008
01009 if (n<5) return 0;
01010 else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
01011 }
01012
01013 unsigned int DiscreteLogWorkFactor(unsigned int n)
01014 {
01015
01016 if (n<5) return 0;
01017 else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
01018 }
01019
01020
01021
01022 void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits)
01023 {
01024
01025 assert(qbits > 4);
01026 assert(pbits > qbits);
01027
01028 if (qbits+1 == pbits)
01029 {
01030 Integer minP = Integer::Power2(pbits-1);
01031 Integer maxP = Integer::Power2(pbits) - 1;
01032 bool success = false;
01033
01034 while (!success)
01035 {
01036 p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
01037 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
01038
01039 while (sieve.NextCandidate(p))
01040 {
01041 assert(IsSmallPrime(p) || SmallDivisorsTest(p));
01042 q = (p-delta) >> 1;
01043 assert(IsSmallPrime(q) || SmallDivisorsTest(q));
01044 if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
01045 {
01046 success = true;
01047 break;
01048 }
01049 }
01050 }
01051
01052 if (delta == 1)
01053 {
01054
01055
01056 for (g=2; Jacobi(g, p) != 1; ++g) {}
01057
01058 assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
01059 }
01060 else
01061 {
01062 assert(delta == -1);
01063
01064
01065 for (g=3; ; ++g)
01066 if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
01067 break;
01068 }
01069 }
01070 else
01071 {
01072 Integer minQ = Integer::Power2(qbits-1);
01073 Integer maxQ = Integer::Power2(qbits) - 1;
01074 Integer minP = Integer::Power2(pbits-1);
01075 Integer maxP = Integer::Power2(pbits) - 1;
01076
01077 do
01078 {
01079 q.Randomize(rng, minQ, maxQ, Integer::PRIME);
01080 } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
01081
01082
01083 if (delta==1)
01084 {
01085 do
01086 {
01087 Integer h(rng, 2, p-2, Integer::ANY);
01088 g = a_exp_b_mod_c(h, (p-1)/q, p);
01089 } while (g <= 1);
01090 assert(a_exp_b_mod_c(g, q, p)==1);
01091 }
01092 else
01093 {
01094 assert(delta==-1);
01095 do
01096 {
01097 Integer h(rng, 3, p-1, Integer::ANY);
01098 if (Jacobi(h*h-4, p)==1)
01099 continue;
01100 g = Lucas((p+1)/q, h, p);
01101 } while (g <= 2);
01102 assert(Lucas(q, g, p) == 2);
01103 }
01104 }
01105 }
01106
01107 NAMESPACE_END
01108
01109 #endif