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