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