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