Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members

nbtheory.cpp

00001 // nbtheory.cpp - written and placed in the public domain by Wei Dai 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; // last prime 32719 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 // calculate a = largest power of 2 that divides (n-1) 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()) // avoid infinite loop if n is a square 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()) // avoid infinite loop if n is a square 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 // calculate a = largest power of 2 that divides n1 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 // delta == 1 or -1 means double sieve with p = 2*q + delta 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 // if the first multiple of p is p, skip it 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 // the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv) 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 // the following two functions are based on code and comments provided by Preda Mihailescu 00477 static bool ProvePrime(const Integer &p, const Integer &q) 00478 { 00479 assert(p < q*q*q); 00480 assert(p % q == 1); 00481 00482 // this is the Quisquater test. Numbers p having passed the Lucas - Lehmer test 00483 // for q and verifying p < q^3 can only be built up of two factors, both = 1 mod q, 00484 // or be prime. The next two lines build the discriminant of a quadratic equation 00485 // which holds iff p is built up of two factors (excercise ... ) 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 // Randomize() will generate a prime provable by trial division 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 // this initializes the sieve to search in the arithmetic 00521 // progression p = p_0 + \lambda * q2 = p_0 + 2 * \lambda * q, 00522 // with q the recursively generated prime above. We will be able 00523 // to use Lucas tets for proving primality. A trick of Quisquater 00524 // allows taking q > cubic_root(p) rather then square_root: this 00525 // decreases the recursion. 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 // not reached 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 // isn't operator overloading great? 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); // not reached 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 Integer GCDI(const Integer &x, const Integer &y) 00687 { 00688 Integer a=x, b=y; 00689 unsigned k=0; 00690 00691 assert(!!a && !!b); 00692 00693 while (a[0]==0 && b[0]==0) 00694 { 00695 a >>= 1; 00696 b >>= 1; 00697 k++; 00698 } 00699 00700 while (a[0]==0) 00701 a >>= 1; 00702 00703 while (b[0]==0) 00704 b >>= 1; 00705 00706 while (1) 00707 { 00708 switch (a.Compare(b)) 00709 { 00710 case -1: 00711 b -= a; 00712 while (b[0]==0) 00713 b >>= 1; 00714 break; 00715 00716 case 0: 00717 return (a <<= k); 00718 00719 case 1: 00720 a -= b; 00721 while (a[0]==0) 00722 a >>= 1; 00723 break; 00724 00725 default: 00726 assert(false); 00727 } 00728 } 00729 } 00730 00731 Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b) 00732 { 00733 assert(b.Positive()); 00734 00735 if (a.Negative()) 00736 return EuclideanMultiplicativeInverse(a%b, b); 00737 00738 if (b[0]==0) 00739 { 00740 if (!b || a[0]==0) 00741 return Integer::Zero(); // no inverse 00742 if (a==1) 00743 return 1; 00744 Integer u = EuclideanMultiplicativeInverse(b, a); 00745 if (!u) 00746 return Integer::Zero(); // no inverse 00747 else 00748 return (b*(a-u)+1)/a; 00749 } 00750 00751 Integer u=1, d=a, v1=b, v3=b, t1, t3, b2=(b+1)>>1; 00752 00753 if (a[0]) 00754 { 00755 t1 = Integer::Zero(); 00756 t3 = -b; 00757 } 00758 else 00759 { 00760 t1 = b2; 00761 t3 = a>>1; 00762 } 00763 00764 while (!!t3) 00765 { 00766 while (t3[0]==0) 00767 { 00768 t3 >>= 1; 00769 if (t1[0]==0) 00770 t1 >>= 1; 00771 else 00772 { 00773 t1 >>= 1; 00774 t1 += b2; 00775 } 00776 } 00777 if (t3.Positive()) 00778 { 00779 u = t1; 00780 d = t3; 00781 } 00782 else 00783 { 00784 v1 = b-t1; 00785 v3 = -t3; 00786 } 00787 t1 = u-v1; 00788 t3 = d-v3; 00789 if (t1.Negative()) 00790 t1 += b; 00791 } 00792 if (d==1) 00793 return u; 00794 else 00795 return Integer::Zero(); // no inverse 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 // v = (v*v1 - p) % m; 00842 v = m.Subtract(m.Multiply(v,v1), p); 00843 // v1 = (v1*v1 - 2) % m; 00844 v1 = m.Subtract(m.Square(v1), two); 00845 } 00846 else 00847 { 00848 // v1 = (v*v1 - p) % m; 00849 v1 = m.Subtract(m.Multiply(v,v1), p); 00850 // v = (v*v - 2) % m; 00851 v = m.Subtract(m.Square(v), two); 00852 } 00853 } 00854 return m.ConvertOut(v); 00855 } 00856 00857 // This is Peter Montgomery's unpublished Lucas sequence evalutation algorithm. 00858 // The total number of multiplies and squares used is less than the binary 00859 // algorithm (see above). Unfortunately I can't get it to run as fast as 00860 // the binary algorithm because of the extra overhead. 00861 /* 00862 Integer Lucas(const Integer &n, const Integer &P, const Integer &modulus) 00863 { 00864 if (!n) 00865 return 2; 00866 00867 #define f(A, B, C) m.Subtract(m.Multiply(A, B), C) 00868 #define X2(A) m.Subtract(m.Square(A), two) 00869 #define X3(A) m.Multiply(A, m.Subtract(m.Square(A), three)) 00870 00871 MontgomeryRepresentation m(modulus); 00872 Integer two=m.ConvertIn(2), three=m.ConvertIn(3); 00873 Integer A=m.ConvertIn(P), B, C, p, d=n, e, r, t, T, U; 00874 00875 while (d!=1) 00876 { 00877 p = d; 00878 unsigned int b = WORD_BITS * p.WordCount(); 00879 Integer alpha = (Integer(5)<<(2*b-2)).SquareRoot() - Integer::Power2(b-1); 00880 r = (p*alpha)>>b; 00881 e = d-r; 00882 B = A; 00883 C = two; 00884 d = r; 00885 00886 while (d!=e) 00887 { 00888 if (d<e) 00889 { 00890 swap(d, e); 00891 swap(A, B); 00892 } 00893 00894 unsigned int dm2 = d[0], em2 = e[0]; 00895 unsigned int dm3 = d%3, em3 = e%3; 00896 00897 // if ((dm6+em6)%3 == 0 && d <= e + (e>>2)) 00898 if ((dm3+em3==0 || dm3+em3==3) && (t = e, t >>= 2, t += e, d <= t)) 00899 { 00900 // #1 00901 // t = (d+d-e)/3; 00902 // t = d; t += d; t -= e; t /= 3; 00903 // e = (e+e-d)/3; 00904 // e += e; e -= d; e /= 3; 00905 // d = t; 00906 00907 // t = (d+e)/3 00908 t = d; t += e; t /= 3; 00909 e -= t; 00910 d -= t; 00911 00912 T = f(A, B, C); 00913 U = f(T, A, B); 00914 B = f(T, B, A); 00915 A = U; 00916 continue; 00917 } 00918 00919 // if (dm6 == em6 && d <= e + (e>>2)) 00920 if (dm3 == em3 && dm2 == em2 && (t = e, t >>= 2, t += e, d <= t)) 00921 { 00922 // #2 00923 // d = (d-e)>>1; 00924 d -= e; d >>= 1; 00925 B = f(A, B, C); 00926 A = X2(A); 00927 continue; 00928 } 00929 00930 // if (d <= (e<<2)) 00931 if (d <= (t = e, t <<= 2)) 00932 { 00933 // #3 00934 d -= e; 00935 C = f(A, B, C); 00936 swap(B, C); 00937 continue; 00938 } 00939 00940 if (dm2 == em2) 00941 { 00942 // #4 00943 // d = (d-e)>>1; 00944 d -= e; d >>= 1; 00945 B = f(A, B, C); 00946 A = X2(A); 00947 continue; 00948 } 00949 00950 if (dm2 == 0) 00951 { 00952 // #5 00953 d >>= 1; 00954 C = f(A, C, B); 00955 A = X2(A); 00956 continue; 00957 } 00958 00959 if (dm3 == 0) 00960 { 00961 // #6 00962 // d = d/3 - e; 00963 d /= 3; d -= e; 00964 T = X2(A); 00965 C = f(T, f(A, B, C), C); 00966 swap(B, C); 00967 A = f(T, A, A); 00968 continue; 00969 } 00970 00971 if (dm3+em3==0 || dm3+em3==3) 00972 { 00973 // #7 00974 // d = (d-e-e)/3; 00975 d -= e; d -= e; d /= 3; 00976 T = f(A, B, C); 00977 B = f(T, A, B); 00978 A = X3(A); 00979 continue; 00980 } 00981 00982 if (dm3 == em3) 00983 { 00984 // #8 00985 // d = (d-e)/3; 00986 d -= e; d /= 3; 00987 T = f(A, B, C); 00988 C = f(A, C, B); 00989 B = T; 00990 A = X3(A); 00991 continue; 00992 } 00993 00994 assert(em2 == 0); 00995 // #9 00996 e >>= 1; 00997 C = f(C, B, A); 00998 B = X2(B); 00999 } 01000 01001 A = f(A, B, C); 01002 } 01003 01004 #undef f 01005 #undef X2 01006 #undef X3 01007 01008 return m.ConvertOut(A); 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 // extrapolated from the table in Odlyzko's "The Future of Integer Factorization" 01028 // updated to reflect the factoring of RSA-130 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 // assuming discrete log takes about the same time as factoring 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 // no prime exists for delta = -1, qbits = 4, and pbits = 5 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 // find g such that g is a quadratic residue mod p, then g has order q 01075 // g=4 always works, but this way we get the smallest quadratic residue (other than 1) 01076 for (g=2; Jacobi(g, p) != 1; ++g) {} 01077 // contributed by Walt Tuvell: g should be the following according to the Law of Quadratic Reciprocity 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 // find g such that g*g-4 is a quadratic non-residue, 01084 // and such that g has order q 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 // find a random g of order q 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

Generated on Fri Aug 13 09:56:54 2004 for Crypto++ by doxygen 1.3.7