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 Sun Mar 14 20:44:27 2004 for Crypto++ by doxygen 1.3.6-20040222