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 unsigned int maxPrimeTableSize = 3511;
00017 const word lastSmallPrime = 32719;
00018 unsigned int primeTableSize=552;
00019
00020 word primeTable[maxPrimeTableSize] =
00021 {2, 3, 5, 7, 11, 13, 17, 19,
00022 23, 29, 31, 37, 41, 43, 47, 53,
00023 59, 61, 67, 71, 73, 79, 83, 89,
00024 97, 101, 103, 107, 109, 113, 127, 131,
00025 137, 139, 149, 151, 157, 163, 167, 173,
00026 179, 181, 191, 193, 197, 199, 211, 223,
00027 227, 229, 233, 239, 241, 251, 257, 263,
00028 269, 271, 277, 281, 283, 293, 307, 311,
00029 313, 317, 331, 337, 347, 349, 353, 359,
00030 367, 373, 379, 383, 389, 397, 401, 409,
00031 419, 421, 431, 433, 439, 443, 449, 457,
00032 461, 463, 467, 479, 487, 491, 499, 503,
00033 509, 521, 523, 541, 547, 557, 563, 569,
00034 571, 577, 587, 593, 599, 601, 607, 613,
00035 617, 619, 631, 641, 643, 647, 653, 659,
00036 661, 673, 677, 683, 691, 701, 709, 719,
00037 727, 733, 739, 743, 751, 757, 761, 769,
00038 773, 787, 797, 809, 811, 821, 823, 827,
00039 829, 839, 853, 857, 859, 863, 877, 881,
00040 883, 887, 907, 911, 919, 929, 937, 941,
00041 947, 953, 967, 971, 977, 983, 991, 997,
00042 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049,
00043 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097,
00044 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163,
00045 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
00046 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283,
00047 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
00048 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423,
00049 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459,
00050 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
00051 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571,
00052 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619,
00053 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693,
00054 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747,
00055 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
00056 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877,
00057 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949,
00058 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
00059 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069,
00060 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
00061 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203,
00062 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267,
00063 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311,
00064 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377,
00065 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
00066 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503,
00067 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579,
00068 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657,
00069 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,
00070 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
00071 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801,
00072 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861,
00073 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939,
00074 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011,
00075 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
00076 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167,
00077 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221,
00078 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301,
00079 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347,
00080 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
00081 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491,
00082 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541,
00083 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607,
00084 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671,
00085 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
00086 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797,
00087 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863,
00088 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923,
00089 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003};
00090
00091 void BuildPrimeTable()
00092 {
00093 unsigned int p=primeTable[primeTableSize-1];
00094 for (unsigned int i=primeTableSize; i<maxPrimeTableSize; i++)
00095 {
00096 int j;
00097 do
00098 {
00099 p+=2;
00100 for (j=1; j<54; j++)
00101 if (p%primeTable[j] == 0)
00102 break;
00103 } while (j!=54);
00104 primeTable[i] = p;
00105 }
00106 primeTableSize = maxPrimeTableSize;
00107 assert(primeTable[primeTableSize-1] == lastSmallPrime);
00108 }
00109
00110 bool IsSmallPrime(const Integer &p)
00111 {
00112 BuildPrimeTable();
00113
00114 if (p.IsPositive() && p <= primeTable[primeTableSize-1])
00115 return std::binary_search(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
00116 else
00117 return false;
00118 }
00119
00120 bool TrialDivision(const Integer &p, unsigned bound)
00121 {
00122 assert(primeTable[primeTableSize-1] >= bound);
00123
00124 unsigned int i;
00125 for (i = 0; primeTable[i]<bound; i++)
00126 if ((p % primeTable[i]) == 0)
00127 return true;
00128
00129 if (bound == primeTable[i])
00130 return (p % bound == 0);
00131 else
00132 return false;
00133 }
00134
00135 bool SmallDivisorsTest(const Integer &p)
00136 {
00137 BuildPrimeTable();
00138 return !TrialDivision(p, primeTable[primeTableSize-1]);
00139 }
00140
00141 bool IsFermatProbablePrime(const Integer &n, const Integer &b)
00142 {
00143 if (n <= 3)
00144 return n==2 || n==3;
00145
00146 assert(n>3 && b>1 && b<n-1);
00147 return a_exp_b_mod_c(b, n-1, n)==1;
00148 }
00149
00150 bool IsStrongProbablePrime(const Integer &n, const Integer &b)
00151 {
00152 if (n <= 3)
00153 return n==2 || n==3;
00154
00155 assert(n>3 && b>1 && b<n-1);
00156
00157 if ((n.IsEven() && n!=2) || GCD(b, n) != 1)
00158 return false;
00159
00160 Integer nminus1 = (n-1);
00161 unsigned int a;
00162
00163
00164 for (a=0; ; a++)
00165 if (nminus1.GetBit(a))
00166 break;
00167 Integer m = nminus1>>a;
00168
00169 Integer z = a_exp_b_mod_c(b, m, n);
00170 if (z==1 || z==nminus1)
00171 return true;
00172 for (unsigned j=1; j<a; j++)
00173 {
00174 z = z.Squared()%n;
00175 if (z==nminus1)
00176 return true;
00177 if (z==1)
00178 return false;
00179 }
00180 return false;
00181 }
00182
00183 bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds)
00184 {
00185 if (n <= 3)
00186 return n==2 || n==3;
00187
00188 assert(n>3);
00189
00190 Integer b;
00191 for (unsigned int i=0; i<rounds; i++)
00192 {
00193 b.Randomize(rng, 2, n-2);
00194 if (!IsStrongProbablePrime(n, b))
00195 return false;
00196 }
00197 return true;
00198 }
00199
00200 bool IsLucasProbablePrime(const Integer &n)
00201 {
00202 if (n <= 1)
00203 return false;
00204
00205 if (n.IsEven())
00206 return n==2;
00207
00208 assert(n>2);
00209
00210 Integer b=3;
00211 unsigned int i=0;
00212 int j;
00213
00214 while ((j=Jacobi(b.Squared()-4, n)) == 1)
00215 {
00216 if (++i==64 && n.IsSquare())
00217 return false;
00218 ++b; ++b;
00219 }
00220
00221 if (j==0)
00222 return false;
00223 else
00224 return Lucas(n+1, b, n)==2;
00225 }
00226
00227 bool IsStrongLucasProbablePrime(const Integer &n)
00228 {
00229 if (n <= 1)
00230 return false;
00231
00232 if (n.IsEven())
00233 return n==2;
00234
00235 assert(n>2);
00236
00237 Integer b=3;
00238 unsigned int i=0;
00239 int j;
00240
00241 while ((j=Jacobi(b.Squared()-4, n)) == 1)
00242 {
00243 if (++i==64 && n.IsSquare())
00244 return false;
00245 ++b; ++b;
00246 }
00247
00248 if (j==0)
00249 return false;
00250
00251 Integer n1 = n+1;
00252 unsigned int a;
00253
00254
00255 for (a=0; ; a++)
00256 if (n1.GetBit(a))
00257 break;
00258 Integer m = n1>>a;
00259
00260 Integer z = Lucas(m, b, n);
00261 if (z==2 || z==n-2)
00262 return true;
00263 for (i=1; i<a; i++)
00264 {
00265 z = (z.Squared()-2)%n;
00266 if (z==n-2)
00267 return true;
00268 if (z==2)
00269 return false;
00270 }
00271 return false;
00272 }
00273
00274 bool IsPrime(const Integer &p)
00275 {
00276 static const Integer lastSmallPrimeSquared = Integer(lastSmallPrime).Squared();
00277
00278 if (p <= lastSmallPrime)
00279 return IsSmallPrime(p);
00280 else if (p <= lastSmallPrimeSquared)
00281 return SmallDivisorsTest(p);
00282 else
00283 return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p);
00284 }
00285
00286 bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level)
00287 {
00288 bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1);
00289 if (level >= 1)
00290 pass = pass && RabinMillerTest(rng, p, 10);
00291 return pass;
00292 }
00293
00294 unsigned int PrimeSearchInterval(const Integer &max)
00295 {
00296 return max.BitCount();
00297 }
00298
00299 static inline bool FastProbablePrimeTest(const Integer &n)
00300 {
00301 return IsStrongProbablePrime(n,2);
00302 }
00303
00304 AlgorithmParameters<AlgorithmParameters<AlgorithmParameters<NullNameValuePairs, Integer::RandomNumberType>, Integer>, Integer>
00305 MakeParametersForTwoPrimesOfEqualSize(unsigned int productBitLength)
00306 {
00307 if (productBitLength < 16)
00308 throw InvalidArgument("invalid bit length");
00309
00310 Integer minP, maxP;
00311
00312 if (productBitLength%2==0)
00313 {
00314 minP = Integer(182) << (productBitLength/2-8);
00315 maxP = Integer::Power2(productBitLength/2)-1;
00316 }
00317 else
00318 {
00319 minP = Integer::Power2((productBitLength-1)/2);
00320 maxP = Integer(181) << ((productBitLength+1)/2-8);
00321 }
00322
00323 return MakeParameters("RandomNumberType", Integer::PRIME)("Min", minP)("Max", maxP);
00324 }
00325
00326 class PrimeSieve
00327 {
00328 public:
00329
00330 PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta=0);
00331 bool NextCandidate(Integer &c);
00332
00333 void DoSieve();
00334 static void SieveSingle(std::vector<bool> &sieve, word p, const Integer &first, const Integer &step, word stepInv);
00335
00336 Integer m_first, m_last, m_step;
00337 signed int m_delta;
00338 word m_next;
00339 std::vector<bool> m_sieve;
00340 };
00341
00342 PrimeSieve::PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta)
00343 : m_first(first), m_last(last), m_step(step), m_delta(delta), m_next(0)
00344 {
00345 DoSieve();
00346 }
00347
00348 bool PrimeSieve::NextCandidate(Integer &c)
00349 {
00350 m_next = std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin();
00351 if (m_next == m_sieve.size())
00352 {
00353 m_first += m_sieve.size()*m_step;
00354 if (m_first > m_last)
00355 return false;
00356 else
00357 {
00358 m_next = 0;
00359 DoSieve();
00360 return NextCandidate(c);
00361 }
00362 }
00363 else
00364 {
00365 c = m_first + m_next*m_step;
00366 ++m_next;
00367 return true;
00368 }
00369 }
00370
00371 void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word p, const Integer &first, const Integer &step, word stepInv)
00372 {
00373 if (stepInv)
00374 {
00375 unsigned int sieveSize = sieve.size();
00376 word j = word((dword(p-(first%p))*stepInv) % p);
00377
00378 if (first.WordCount() <= 1 && first + step*j == p)
00379 j += p;
00380 for (; j < sieveSize; j += p)
00381 sieve[j] = true;
00382 }
00383 }
00384
00385 void PrimeSieve::DoSieve()
00386 {
00387 BuildPrimeTable();
00388
00389 const unsigned int maxSieveSize = 32768;
00390 unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong();
00391
00392 m_sieve.clear();
00393 m_sieve.resize(sieveSize, false);
00394
00395 if (m_delta == 0)
00396 {
00397 for (unsigned int i = 0; i < primeTableSize; ++i)
00398 SieveSingle(m_sieve, primeTable[i], m_first, m_step, m_step.InverseMod(primeTable[i]));
00399 }
00400 else
00401 {
00402 assert(m_step%2==0);
00403 Integer qFirst = (m_first-m_delta) >> 1;
00404 Integer halfStep = m_step >> 1;
00405 for (unsigned int i = 0; i < primeTableSize; ++i)
00406 {
00407 word p = primeTable[i];
00408 word stepInv = m_step.InverseMod(p);
00409 SieveSingle(m_sieve, p, m_first, m_step, stepInv);
00410
00411 word halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p;
00412 SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
00413 }
00414 }
00415 }
00416
00417 bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector)
00418 {
00419 assert(!equiv.IsNegative() && equiv < mod);
00420
00421 Integer gcd = GCD(equiv, mod);
00422 if (gcd != Integer::One())
00423 {
00424
00425 if (p <= gcd && gcd <= max && IsPrime(gcd))
00426 {
00427 p = gcd;
00428 return true;
00429 }
00430 else
00431 return false;
00432 }
00433
00434 BuildPrimeTable();
00435
00436 if (p <= primeTable[primeTableSize-1])
00437 {
00438 word *pItr;
00439
00440 --p;
00441 if (p.IsPositive())
00442 pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
00443 else
00444 pItr = primeTable;
00445
00446 while (pItr < primeTable+primeTableSize && *pItr%mod != equiv)
00447 ++pItr;
00448
00449 if (pItr < primeTable+primeTableSize)
00450 {
00451 p = *pItr;
00452 return p <= max;
00453 }
00454
00455 p = primeTable[primeTableSize-1]+1;
00456 }
00457
00458 assert(p > primeTable[primeTableSize-1]);
00459
00460 if (mod.IsOdd())
00461 return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector);
00462
00463 p += (equiv-p)%mod;
00464
00465 if (p>max)
00466 return false;
00467
00468 PrimeSieve sieve(p, max, mod);
00469
00470 while (sieve.NextCandidate(p))
00471 {
00472 if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p))
00473 return true;
00474 }
00475
00476 return false;
00477 }
00478
00479
00480 static bool ProvePrime(const Integer &p, const Integer &q)
00481 {
00482 assert(p < q*q*q);
00483 assert(p % q == 1);
00484
00485
00486
00487
00488
00489
00490 Integer r = (p-1)/q;
00491 if (((r%q).Squared()-4*(r/q)).IsSquare())
00492 return false;
00493
00494 assert(primeTableSize >= 50);
00495 for (int i=0; i<50; i++)
00496 {
00497 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
00498 if (b != 1)
00499 return a_exp_b_mod_c(b, q, p) == 1;
00500 }
00501 return false;
00502 }
00503
00504 Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits)
00505 {
00506 Integer p;
00507 Integer minP = Integer::Power2(pbits-1);
00508 Integer maxP = Integer::Power2(pbits) - 1;
00509
00510 if (maxP <= Integer(lastSmallPrime).Squared())
00511 {
00512
00513 p.Randomize(rng, minP, maxP, Integer::PRIME);
00514 return p;
00515 }
00516
00517 unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36);
00518 Integer q = MihailescuProvablePrime(rng, qbits);
00519 Integer q2 = q<<1;
00520
00521 while (true)
00522 {
00523
00524
00525
00526
00527
00528
00529
00530 p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2);
00531 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
00532
00533 while (sieve.NextCandidate(p))
00534 {
00535 if (FastProbablePrimeTest(p) && ProvePrime(p, q))
00536 return p;
00537 }
00538 }
00539
00540
00541 return p;
00542 }
00543
00544 Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
00545 {
00546 const unsigned smallPrimeBound = 29, c_opt=10;
00547 Integer p;
00548
00549 BuildPrimeTable();
00550 if (bits < smallPrimeBound)
00551 {
00552 do
00553 p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2);
00554 while (TrialDivision(p, 1 << ((bits+1)/2)));
00555 }
00556 else
00557 {
00558 const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
00559 double relativeSize;
00560 do
00561 relativeSize = pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1);
00562 while (bits * relativeSize >= bits - margin);
00563
00564 Integer a,b;
00565 Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize));
00566 Integer I = Integer::Power2(bits-2)/q;
00567 Integer I2 = I << 1;
00568 unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt);
00569 bool success = false;
00570 while (!success)
00571 {
00572 p.Randomize(rng, I, I2, Integer::ANY);
00573 p *= q; p <<= 1; ++p;
00574 if (!TrialDivision(p, trialDivisorBound))
00575 {
00576 a.Randomize(rng, 2, p-1, Integer::ANY);
00577 b = a_exp_b_mod_c(a, (p-1)/q, p);
00578 success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1);
00579 }
00580 }
00581 }
00582 return p;
00583 }
00584
00585 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u)
00586 {
00587
00588 return p * (u * (xq-xp) % q) + xp;
00589 }
00590
00591 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q)
00592 {
00593 return CRT(xp, p, xq, q, EuclideanMultiplicativeInverse(p, q));
00594 }
00595
00596 Integer ModularSquareRoot(const Integer &a, const Integer &p)
00597 {
00598 if (p%4 == 3)
00599 return a_exp_b_mod_c(a, (p+1)/4, p);
00600
00601 Integer q=p-1;
00602 unsigned int r=0;
00603 while (q.IsEven())
00604 {
00605 r++;
00606 q >>= 1;
00607 }
00608
00609 Integer n=2;
00610 while (Jacobi(n, p) != -1)
00611 ++n;
00612
00613 Integer y = a_exp_b_mod_c(n, q, p);
00614 Integer x = a_exp_b_mod_c(a, (q-1)/2, p);
00615 Integer b = (x.Squared()%p)*a%p;
00616 x = a*x%p;
00617 Integer tempb, t;
00618
00619 while (b != 1)
00620 {
00621 unsigned m=0;
00622 tempb = b;
00623 do
00624 {
00625 m++;
00626 b = b.Squared()%p;
00627 if (m==r)
00628 return Integer::Zero();
00629 }
00630 while (b != 1);
00631
00632 t = y;
00633 for (unsigned i=0; i<r-m-1; i++)
00634 t = t.Squared()%p;
00635 y = t.Squared()%p;
00636 r = m;
00637 x = x*t%p;
00638 b = tempb*y%p;
00639 }
00640
00641 assert(x.Squared()%p == a);
00642 return x;
00643 }
00644
00645 bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p)
00646 {
00647 Integer D = (b.Squared() - 4*a*c) % p;
00648 switch (Jacobi(D, p))
00649 {
00650 default:
00651 assert(false);
00652 return false;
00653 case -1:
00654 return false;
00655 case 0:
00656 r1 = r2 = (-b*(a+a).InverseMod(p)) % p;
00657 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
00658 return true;
00659 case 1:
00660 Integer s = ModularSquareRoot(D, p);
00661 Integer t = (a+a).InverseMod(p);
00662 r1 = (s-b)*t % p;
00663 r2 = (-s-b)*t % p;
00664 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
00665 assert(((r2.Squared()*a + r2*b + c) % p).IsZero());
00666 return true;
00667 }
00668 }
00669
00670 Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq,
00671 const Integer &p, const Integer &q, const Integer &u)
00672 {
00673 Integer p2 = ModularExponentiation((a % p), dp, p);
00674 Integer q2 = ModularExponentiation((a % q), dq, q);
00675 return CRT(p2, p, q2, q, u);
00676 }
00677
00678 Integer ModularRoot(const Integer &a, const Integer &e,
00679 const Integer &p, const Integer &q)
00680 {
00681 Integer dp = EuclideanMultiplicativeInverse(e, p-1);
00682 Integer dq = EuclideanMultiplicativeInverse(e, q-1);
00683 Integer u = EuclideanMultiplicativeInverse(p, q);
00684 assert(!!dp && !!dq && !!u);
00685 return ModularRoot(a, dp, dq, p, q, u);
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
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802 int Jacobi(const Integer &aIn, const Integer &bIn)
00803 {
00804 assert(bIn.IsOdd());
00805
00806 Integer b = bIn, a = aIn%bIn;
00807 int result = 1;
00808
00809 while (!!a)
00810 {
00811 unsigned i=0;
00812 while (a.GetBit(i)==0)
00813 i++;
00814 a>>=i;
00815
00816 if (i%2==1 && (b%8==3 || b%8==5))
00817 result = -result;
00818
00819 if (a%4==3 && b%4==3)
00820 result = -result;
00821
00822 std::swap(a, b);
00823 a %= b;
00824 }
00825
00826 return (b==1) ? result : 0;
00827 }
00828
00829 Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n)
00830 {
00831 unsigned i = e.BitCount();
00832 if (i==0)
00833 return Integer::Two();
00834
00835 MontgomeryRepresentation m(n);
00836 Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two());
00837 Integer v=p, v1=m.Subtract(m.Square(p), two);
00838
00839 i--;
00840 while (i--)
00841 {
00842 if (e.GetBit(i))
00843 {
00844
00845 v = m.Subtract(m.Multiply(v,v1), p);
00846
00847 v1 = m.Subtract(m.Square(v1), two);
00848 }
00849 else
00850 {
00851
00852 v1 = m.Subtract(m.Multiply(v,v1), p);
00853
00854 v = m.Subtract(m.Square(v), two);
00855 }
00856 }
00857 return m.ConvertOut(v);
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
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
01016 {
01017 Integer d = (m*m-4);
01018 Integer p2 = p-Jacobi(d,p);
01019 Integer q2 = q-Jacobi(d,q);
01020 return CRT(Lucas(EuclideanMultiplicativeInverse(e,p2), m, p), p, Lucas(EuclideanMultiplicativeInverse(e,q2), m, q), q, u);
01021 }
01022
01023 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q)
01024 {
01025 return InverseLucas(e, m, p, q, EuclideanMultiplicativeInverse(p, q));
01026 }
01027
01028 unsigned int FactoringWorkFactor(unsigned int n)
01029 {
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 unsigned int DiscreteLogWorkFactor(unsigned int n)
01037 {
01038
01039 if (n<5) return 0;
01040 else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
01041 }
01042
01043
01044
01045 void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits)
01046 {
01047
01048 assert(qbits > 4);
01049 assert(pbits > qbits);
01050
01051 if (qbits+1 == pbits)
01052 {
01053 Integer minP = Integer::Power2(pbits-1);
01054 Integer maxP = Integer::Power2(pbits) - 1;
01055 bool success = false;
01056
01057 while (!success)
01058 {
01059 p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
01060 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
01061
01062 while (sieve.NextCandidate(p))
01063 {
01064 assert(IsSmallPrime(p) || SmallDivisorsTest(p));
01065 q = (p-delta) >> 1;
01066 assert(IsSmallPrime(q) || SmallDivisorsTest(q));
01067 if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
01068 {
01069 success = true;
01070 break;
01071 }
01072 }
01073 }
01074
01075 if (delta == 1)
01076 {
01077
01078
01079 for (g=2; Jacobi(g, p) != 1; ++g) {}
01080
01081 assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
01082 }
01083 else
01084 {
01085 assert(delta == -1);
01086
01087
01088 for (g=3; ; ++g)
01089 if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
01090 break;
01091 }
01092 }
01093 else
01094 {
01095 Integer minQ = Integer::Power2(qbits-1);
01096 Integer maxQ = Integer::Power2(qbits) - 1;
01097 Integer minP = Integer::Power2(pbits-1);
01098 Integer maxP = Integer::Power2(pbits) - 1;
01099
01100 do
01101 {
01102 q.Randomize(rng, minQ, maxQ, Integer::PRIME);
01103 } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
01104
01105
01106 if (delta==1)
01107 {
01108 do
01109 {
01110 Integer h(rng, 2, p-2, Integer::ANY);
01111 g = a_exp_b_mod_c(h, (p-1)/q, p);
01112 } while (g <= 1);
01113 assert(a_exp_b_mod_c(g, q, p)==1);
01114 }
01115 else
01116 {
01117 assert(delta==-1);
01118 do
01119 {
01120 Integer h(rng, 3, p-1, Integer::ANY);
01121 if (Jacobi(h*h-4, p)==1)
01122 continue;
01123 g = Lucas((p+1)/q, h, p);
01124 } while (g <= 2);
01125 assert(Lucas(q, g, p) == 2);
01126 }
01127 }
01128 }
01129
01130 NAMESPACE_END
01131
01132 #endif