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 #ifdef _OPENMP
00015
00016 #include <omp.h>
00017 #endif
00018
00019 NAMESPACE_BEGIN(CryptoPP)
00020
00021 const word s_lastSmallPrime = 32719;
00022
00023 struct NewPrimeTable
00024 {
00025 std::vector<word16> * operator()() const
00026 {
00027 const unsigned int maxPrimeTableSize = 3511;
00028
00029 std::auto_ptr<std::vector<word16> > pPrimeTable(new std::vector<word16>);
00030 std::vector<word16> &primeTable = *pPrimeTable;
00031 primeTable.reserve(maxPrimeTableSize);
00032
00033 primeTable.push_back(2);
00034 unsigned int testEntriesEnd = 1;
00035
00036 for (unsigned int p=3; p<=s_lastSmallPrime; p+=2)
00037 {
00038 unsigned int j;
00039 for (j=1; j<testEntriesEnd; j++)
00040 if (p%primeTable[j] == 0)
00041 break;
00042 if (j == testEntriesEnd)
00043 {
00044 primeTable.push_back(p);
00045 testEntriesEnd = UnsignedMin(54U, primeTable.size());
00046 }
00047 }
00048
00049 return pPrimeTable.release();
00050 }
00051 };
00052
00053 const word16 * GetPrimeTable(unsigned int &size)
00054 {
00055 const std::vector<word16> &primeTable = Singleton<std::vector<word16>, NewPrimeTable>().Ref();
00056 size = (unsigned int)primeTable.size();
00057 return &primeTable[0];
00058 }
00059
00060 bool IsSmallPrime(const Integer &p)
00061 {
00062 unsigned int primeTableSize;
00063 const word16 * primeTable = GetPrimeTable(primeTableSize);
00064
00065 if (p.IsPositive() && p <= primeTable[primeTableSize-1])
00066 return std::binary_search(primeTable, primeTable+primeTableSize, (word16)p.ConvertToLong());
00067 else
00068 return false;
00069 }
00070
00071 bool TrialDivision(const Integer &p, unsigned bound)
00072 {
00073 unsigned int primeTableSize;
00074 const word16 * primeTable = GetPrimeTable(primeTableSize);
00075
00076 assert(primeTable[primeTableSize-1] >= bound);
00077
00078 unsigned int i;
00079 for (i = 0; primeTable[i]<bound; i++)
00080 if ((p % primeTable[i]) == 0)
00081 return true;
00082
00083 if (bound == primeTable[i])
00084 return (p % bound == 0);
00085 else
00086 return false;
00087 }
00088
00089 bool SmallDivisorsTest(const Integer &p)
00090 {
00091 unsigned int primeTableSize;
00092 const word16 * primeTable = GetPrimeTable(primeTableSize);
00093 return !TrialDivision(p, primeTable[primeTableSize-1]);
00094 }
00095
00096 bool IsFermatProbablePrime(const Integer &n, const Integer &b)
00097 {
00098 if (n <= 3)
00099 return n==2 || n==3;
00100
00101 assert(n>3 && b>1 && b<n-1);
00102 return a_exp_b_mod_c(b, n-1, n)==1;
00103 }
00104
00105 bool IsStrongProbablePrime(const Integer &n, const Integer &b)
00106 {
00107 if (n <= 3)
00108 return n==2 || n==3;
00109
00110 assert(n>3 && b>1 && b<n-1);
00111
00112 if ((n.IsEven() && n!=2) || GCD(b, n) != 1)
00113 return false;
00114
00115 Integer nminus1 = (n-1);
00116 unsigned int a;
00117
00118
00119 for (a=0; ; a++)
00120 if (nminus1.GetBit(a))
00121 break;
00122 Integer m = nminus1>>a;
00123
00124 Integer z = a_exp_b_mod_c(b, m, n);
00125 if (z==1 || z==nminus1)
00126 return true;
00127 for (unsigned j=1; j<a; j++)
00128 {
00129 z = z.Squared()%n;
00130 if (z==nminus1)
00131 return true;
00132 if (z==1)
00133 return false;
00134 }
00135 return false;
00136 }
00137
00138 bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds)
00139 {
00140 if (n <= 3)
00141 return n==2 || n==3;
00142
00143 assert(n>3);
00144
00145 Integer b;
00146 for (unsigned int i=0; i<rounds; i++)
00147 {
00148 b.Randomize(rng, 2, n-2);
00149 if (!IsStrongProbablePrime(n, b))
00150 return false;
00151 }
00152 return true;
00153 }
00154
00155 bool IsLucasProbablePrime(const Integer &n)
00156 {
00157 if (n <= 1)
00158 return false;
00159
00160 if (n.IsEven())
00161 return n==2;
00162
00163 assert(n>2);
00164
00165 Integer b=3;
00166 unsigned int i=0;
00167 int j;
00168
00169 while ((j=Jacobi(b.Squared()-4, n)) == 1)
00170 {
00171 if (++i==64 && n.IsSquare())
00172 return false;
00173 ++b; ++b;
00174 }
00175
00176 if (j==0)
00177 return false;
00178 else
00179 return Lucas(n+1, b, n)==2;
00180 }
00181
00182 bool IsStrongLucasProbablePrime(const Integer &n)
00183 {
00184 if (n <= 1)
00185 return false;
00186
00187 if (n.IsEven())
00188 return n==2;
00189
00190 assert(n>2);
00191
00192 Integer b=3;
00193 unsigned int i=0;
00194 int j;
00195
00196 while ((j=Jacobi(b.Squared()-4, n)) == 1)
00197 {
00198 if (++i==64 && n.IsSquare())
00199 return false;
00200 ++b; ++b;
00201 }
00202
00203 if (j==0)
00204 return false;
00205
00206 Integer n1 = n+1;
00207 unsigned int a;
00208
00209
00210 for (a=0; ; a++)
00211 if (n1.GetBit(a))
00212 break;
00213 Integer m = n1>>a;
00214
00215 Integer z = Lucas(m, b, n);
00216 if (z==2 || z==n-2)
00217 return true;
00218 for (i=1; i<a; i++)
00219 {
00220 z = (z.Squared()-2)%n;
00221 if (z==n-2)
00222 return true;
00223 if (z==2)
00224 return false;
00225 }
00226 return false;
00227 }
00228
00229 struct NewLastSmallPrimeSquared
00230 {
00231 Integer * operator()() const
00232 {
00233 return new Integer(Integer(s_lastSmallPrime).Squared());
00234 }
00235 };
00236
00237 bool IsPrime(const Integer &p)
00238 {
00239 if (p <= s_lastSmallPrime)
00240 return IsSmallPrime(p);
00241 else if (p <= Singleton<Integer, NewLastSmallPrimeSquared>().Ref())
00242 return SmallDivisorsTest(p);
00243 else
00244 return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p);
00245 }
00246
00247 bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level)
00248 {
00249 bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1);
00250 if (level >= 1)
00251 pass = pass && RabinMillerTest(rng, p, 10);
00252 return pass;
00253 }
00254
00255 unsigned int PrimeSearchInterval(const Integer &max)
00256 {
00257 return max.BitCount();
00258 }
00259
00260 static inline bool FastProbablePrimeTest(const Integer &n)
00261 {
00262 return IsStrongProbablePrime(n,2);
00263 }
00264
00265 AlgorithmParameters MakeParametersForTwoPrimesOfEqualSize(unsigned int productBitLength)
00266 {
00267 if (productBitLength < 16)
00268 throw InvalidArgument("invalid bit length");
00269
00270 Integer minP, maxP;
00271
00272 if (productBitLength%2==0)
00273 {
00274 minP = Integer(182) << (productBitLength/2-8);
00275 maxP = Integer::Power2(productBitLength/2)-1;
00276 }
00277 else
00278 {
00279 minP = Integer::Power2((productBitLength-1)/2);
00280 maxP = Integer(181) << ((productBitLength+1)/2-8);
00281 }
00282
00283 return MakeParameters("RandomNumberType", Integer::PRIME)("Min", minP)("Max", maxP);
00284 }
00285
00286 class PrimeSieve
00287 {
00288 public:
00289
00290 PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta=0);
00291 bool NextCandidate(Integer &c);
00292
00293 void DoSieve();
00294 static void SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv);
00295
00296 Integer m_first, m_last, m_step;
00297 signed int m_delta;
00298 word m_next;
00299 std::vector<bool> m_sieve;
00300 };
00301
00302 PrimeSieve::PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta)
00303 : m_first(first), m_last(last), m_step(step), m_delta(delta), m_next(0)
00304 {
00305 DoSieve();
00306 }
00307
00308 bool PrimeSieve::NextCandidate(Integer &c)
00309 {
00310 bool safe = SafeConvert(std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin(), m_next);
00311 assert(safe);
00312 if (m_next == m_sieve.size())
00313 {
00314 m_first += long(m_sieve.size())*m_step;
00315 if (m_first > m_last)
00316 return false;
00317 else
00318 {
00319 m_next = 0;
00320 DoSieve();
00321 return NextCandidate(c);
00322 }
00323 }
00324 else
00325 {
00326 c = m_first + long(m_next)*m_step;
00327 ++m_next;
00328 return true;
00329 }
00330 }
00331
00332 void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv)
00333 {
00334 if (stepInv)
00335 {
00336 size_t sieveSize = sieve.size();
00337 size_t j = (word32(p-(first%p))*stepInv) % p;
00338
00339 if (first.WordCount() <= 1 && first + step*long(j) == p)
00340 j += p;
00341 for (; j < sieveSize; j += p)
00342 sieve[j] = true;
00343 }
00344 }
00345
00346 void PrimeSieve::DoSieve()
00347 {
00348 unsigned int primeTableSize;
00349 const word16 * primeTable = GetPrimeTable(primeTableSize);
00350
00351 const unsigned int maxSieveSize = 32768;
00352 unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong();
00353
00354 m_sieve.clear();
00355 m_sieve.resize(sieveSize, false);
00356
00357 if (m_delta == 0)
00358 {
00359 for (unsigned int i = 0; i < primeTableSize; ++i)
00360 SieveSingle(m_sieve, primeTable[i], m_first, m_step, (word16)m_step.InverseMod(primeTable[i]));
00361 }
00362 else
00363 {
00364 assert(m_step%2==0);
00365 Integer qFirst = (m_first-m_delta) >> 1;
00366 Integer halfStep = m_step >> 1;
00367 for (unsigned int i = 0; i < primeTableSize; ++i)
00368 {
00369 word16 p = primeTable[i];
00370 word16 stepInv = (word16)m_step.InverseMod(p);
00371 SieveSingle(m_sieve, p, m_first, m_step, stepInv);
00372
00373 word16 halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p;
00374 SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
00375 }
00376 }
00377 }
00378
00379 bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector)
00380 {
00381 assert(!equiv.IsNegative() && equiv < mod);
00382
00383 Integer gcd = GCD(equiv, mod);
00384 if (gcd != Integer::One())
00385 {
00386
00387 if (p <= gcd && gcd <= max && IsPrime(gcd) && (!pSelector || pSelector->IsAcceptable(gcd)))
00388 {
00389 p = gcd;
00390 return true;
00391 }
00392 else
00393 return false;
00394 }
00395
00396 unsigned int primeTableSize;
00397 const word16 * primeTable = GetPrimeTable(primeTableSize);
00398
00399 if (p <= primeTable[primeTableSize-1])
00400 {
00401 const word16 *pItr;
00402
00403 --p;
00404 if (p.IsPositive())
00405 pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
00406 else
00407 pItr = primeTable;
00408
00409 while (pItr < primeTable+primeTableSize && !(*pItr%mod == equiv && (!pSelector || pSelector->IsAcceptable(*pItr))))
00410 ++pItr;
00411
00412 if (pItr < primeTable+primeTableSize)
00413 {
00414 p = *pItr;
00415 return p <= max;
00416 }
00417
00418 p = primeTable[primeTableSize-1]+1;
00419 }
00420
00421 assert(p > primeTable[primeTableSize-1]);
00422
00423 if (mod.IsOdd())
00424 return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector);
00425
00426 p += (equiv-p)%mod;
00427
00428 if (p>max)
00429 return false;
00430
00431 PrimeSieve sieve(p, max, mod);
00432
00433 while (sieve.NextCandidate(p))
00434 {
00435 if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p))
00436 return true;
00437 }
00438
00439 return false;
00440 }
00441
00442
00443 static bool ProvePrime(const Integer &p, const Integer &q)
00444 {
00445 assert(p < q*q*q);
00446 assert(p % q == 1);
00447
00448
00449
00450
00451
00452
00453 Integer r = (p-1)/q;
00454 if (((r%q).Squared()-4*(r/q)).IsSquare())
00455 return false;
00456
00457 unsigned int primeTableSize;
00458 const word16 * primeTable = GetPrimeTable(primeTableSize);
00459
00460 assert(primeTableSize >= 50);
00461 for (int i=0; i<50; i++)
00462 {
00463 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
00464 if (b != 1)
00465 return a_exp_b_mod_c(b, q, p) == 1;
00466 }
00467 return false;
00468 }
00469
00470 Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits)
00471 {
00472 Integer p;
00473 Integer minP = Integer::Power2(pbits-1);
00474 Integer maxP = Integer::Power2(pbits) - 1;
00475
00476 if (maxP <= Integer(s_lastSmallPrime).Squared())
00477 {
00478
00479 p.Randomize(rng, minP, maxP, Integer::PRIME);
00480 return p;
00481 }
00482
00483 unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36);
00484 Integer q = MihailescuProvablePrime(rng, qbits);
00485 Integer q2 = q<<1;
00486
00487 while (true)
00488 {
00489
00490
00491
00492
00493
00494
00495
00496 p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2);
00497 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
00498
00499 while (sieve.NextCandidate(p))
00500 {
00501 if (FastProbablePrimeTest(p) && ProvePrime(p, q))
00502 return p;
00503 }
00504 }
00505
00506
00507 return p;
00508 }
00509
00510 Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
00511 {
00512 const unsigned smallPrimeBound = 29, c_opt=10;
00513 Integer p;
00514
00515 unsigned int primeTableSize;
00516 const word16 * primeTable = GetPrimeTable(primeTableSize);
00517
00518 if (bits < smallPrimeBound)
00519 {
00520 do
00521 p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2);
00522 while (TrialDivision(p, 1 << ((bits+1)/2)));
00523 }
00524 else
00525 {
00526 const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
00527 double relativeSize;
00528 do
00529 relativeSize = pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1);
00530 while (bits * relativeSize >= bits - margin);
00531
00532 Integer a,b;
00533 Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize));
00534 Integer I = Integer::Power2(bits-2)/q;
00535 Integer I2 = I << 1;
00536 unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt);
00537 bool success = false;
00538 while (!success)
00539 {
00540 p.Randomize(rng, I, I2, Integer::ANY);
00541 p *= q; p <<= 1; ++p;
00542 if (!TrialDivision(p, trialDivisorBound))
00543 {
00544 a.Randomize(rng, 2, p-1, Integer::ANY);
00545 b = a_exp_b_mod_c(a, (p-1)/q, p);
00546 success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1);
00547 }
00548 }
00549 }
00550 return p;
00551 }
00552
00553 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u)
00554 {
00555
00556 return p * (u * (xq-xp) % q) + xp;
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570 }
00571
00572 Integer ModularSquareRoot(const Integer &a, const Integer &p)
00573 {
00574 if (p%4 == 3)
00575 return a_exp_b_mod_c(a, (p+1)/4, p);
00576
00577 Integer q=p-1;
00578 unsigned int r=0;
00579 while (q.IsEven())
00580 {
00581 r++;
00582 q >>= 1;
00583 }
00584
00585 Integer n=2;
00586 while (Jacobi(n, p) != -1)
00587 ++n;
00588
00589 Integer y = a_exp_b_mod_c(n, q, p);
00590 Integer x = a_exp_b_mod_c(a, (q-1)/2, p);
00591 Integer b = (x.Squared()%p)*a%p;
00592 x = a*x%p;
00593 Integer tempb, t;
00594
00595 while (b != 1)
00596 {
00597 unsigned m=0;
00598 tempb = b;
00599 do
00600 {
00601 m++;
00602 b = b.Squared()%p;
00603 if (m==r)
00604 return Integer::Zero();
00605 }
00606 while (b != 1);
00607
00608 t = y;
00609 for (unsigned i=0; i<r-m-1; i++)
00610 t = t.Squared()%p;
00611 y = t.Squared()%p;
00612 r = m;
00613 x = x*t%p;
00614 b = tempb*y%p;
00615 }
00616
00617 assert(x.Squared()%p == a);
00618 return x;
00619 }
00620
00621 bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p)
00622 {
00623 Integer D = (b.Squared() - 4*a*c) % p;
00624 switch (Jacobi(D, p))
00625 {
00626 default:
00627 assert(false);
00628 return false;
00629 case -1:
00630 return false;
00631 case 0:
00632 r1 = r2 = (-b*(a+a).InverseMod(p)) % p;
00633 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
00634 return true;
00635 case 1:
00636 Integer s = ModularSquareRoot(D, p);
00637 Integer t = (a+a).InverseMod(p);
00638 r1 = (s-b)*t % p;
00639 r2 = (-s-b)*t % p;
00640 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
00641 assert(((r2.Squared()*a + r2*b + c) % p).IsZero());
00642 return true;
00643 }
00644 }
00645
00646 Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq,
00647 const Integer &p, const Integer &q, const Integer &u)
00648 {
00649 Integer p2, q2;
00650 #pragma omp parallel
00651 #pragma omp sections
00652 {
00653 #pragma omp section
00654 p2 = ModularExponentiation((a % p), dp, p);
00655 #pragma omp section
00656 q2 = ModularExponentiation((a % q), dq, q);
00657 }
00658 return CRT(p2, p, q2, q, u);
00659 }
00660
00661 Integer ModularRoot(const Integer &a, const Integer &e,
00662 const Integer &p, const Integer &q)
00663 {
00664 Integer dp = EuclideanMultiplicativeInverse(e, p-1);
00665 Integer dq = EuclideanMultiplicativeInverse(e, q-1);
00666 Integer u = EuclideanMultiplicativeInverse(p, q);
00667 assert(!!dp && !!dq && !!u);
00668 return ModularRoot(a, dp, dq, p, q, u);
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
00780
00781
00782
00783
00784
00785 int Jacobi(const Integer &aIn, const Integer &bIn)
00786 {
00787 assert(bIn.IsOdd());
00788
00789 Integer b = bIn, a = aIn%bIn;
00790 int result = 1;
00791
00792 while (!!a)
00793 {
00794 unsigned i=0;
00795 while (a.GetBit(i)==0)
00796 i++;
00797 a>>=i;
00798
00799 if (i%2==1 && (b%8==3 || b%8==5))
00800 result = -result;
00801
00802 if (a%4==3 && b%4==3)
00803 result = -result;
00804
00805 std::swap(a, b);
00806 a %= b;
00807 }
00808
00809 return (b==1) ? result : 0;
00810 }
00811
00812 Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n)
00813 {
00814 unsigned i = e.BitCount();
00815 if (i==0)
00816 return Integer::Two();
00817
00818 MontgomeryRepresentation m(n);
00819 Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two());
00820 Integer v=p, v1=m.Subtract(m.Square(p), two);
00821
00822 i--;
00823 while (i--)
00824 {
00825 if (e.GetBit(i))
00826 {
00827
00828 v = m.Subtract(m.Multiply(v,v1), p);
00829
00830 v1 = m.Subtract(m.Square(v1), two);
00831 }
00832 else
00833 {
00834
00835 v1 = m.Subtract(m.Multiply(v,v1), p);
00836
00837 v = m.Subtract(m.Square(v), two);
00838 }
00839 }
00840 return m.ConvertOut(v);
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
00993
00994
00995
00996
00997
00998 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
00999 {
01000 Integer d = (m*m-4);
01001 Integer p2, q2;
01002 #pragma omp parallel
01003 #pragma omp sections
01004 {
01005 #pragma omp section
01006 {
01007 p2 = p-Jacobi(d,p);
01008 p2 = Lucas(EuclideanMultiplicativeInverse(e,p2), m, p);
01009 }
01010 #pragma omp section
01011 {
01012 q2 = q-Jacobi(d,q);
01013 q2 = Lucas(EuclideanMultiplicativeInverse(e,q2), m, q);
01014 }
01015 }
01016 return CRT(p2, p, q2, q, u);
01017 }
01018
01019 unsigned int FactoringWorkFactor(unsigned int n)
01020 {
01021
01022
01023 if (n<5) return 0;
01024 else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
01025 }
01026
01027 unsigned int DiscreteLogWorkFactor(unsigned int n)
01028 {
01029
01030 if (n<5) return 0;
01031 else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
01032 }
01033
01034
01035
01036 void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits)
01037 {
01038
01039 assert(qbits > 4);
01040 assert(pbits > qbits);
01041
01042 if (qbits+1 == pbits)
01043 {
01044 Integer minP = Integer::Power2(pbits-1);
01045 Integer maxP = Integer::Power2(pbits) - 1;
01046 bool success = false;
01047
01048 while (!success)
01049 {
01050 p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
01051 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
01052
01053 while (sieve.NextCandidate(p))
01054 {
01055 assert(IsSmallPrime(p) || SmallDivisorsTest(p));
01056 q = (p-delta) >> 1;
01057 assert(IsSmallPrime(q) || SmallDivisorsTest(q));
01058 if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
01059 {
01060 success = true;
01061 break;
01062 }
01063 }
01064 }
01065
01066 if (delta == 1)
01067 {
01068
01069
01070 for (g=2; Jacobi(g, p) != 1; ++g) {}
01071
01072 assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
01073 }
01074 else
01075 {
01076 assert(delta == -1);
01077
01078
01079 for (g=3; ; ++g)
01080 if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
01081 break;
01082 }
01083 }
01084 else
01085 {
01086 Integer minQ = Integer::Power2(qbits-1);
01087 Integer maxQ = Integer::Power2(qbits) - 1;
01088 Integer minP = Integer::Power2(pbits-1);
01089 Integer maxP = Integer::Power2(pbits) - 1;
01090
01091 do
01092 {
01093 q.Randomize(rng, minQ, maxQ, Integer::PRIME);
01094 } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
01095
01096
01097 if (delta==1)
01098 {
01099 do
01100 {
01101 Integer h(rng, 2, p-2, Integer::ANY);
01102 g = a_exp_b_mod_c(h, (p-1)/q, p);
01103 } while (g <= 1);
01104 assert(a_exp_b_mod_c(g, q, p)==1);
01105 }
01106 else
01107 {
01108 assert(delta==-1);
01109 do
01110 {
01111 Integer h(rng, 3, p-1, Integer::ANY);
01112 if (Jacobi(h*h-4, p)==1)
01113 continue;
01114 g = Lucas((p+1)/q, h, p);
01115 } while (g <= 2);
01116 assert(Lucas(q, g, p) == 2);
01117 }
01118 }
01119 }
01120
01121 NAMESPACE_END
01122
01123 #endif