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

nbtheory.cpp

00001 // nbtheory.cpp - written and placed in the public domain by Wei Dai
00002 
00003 #include "pch.h"
00004 
00005 #ifndef CRYPTOPP_IMPORTS
00006 
00007 #include "nbtheory.h"
00008 #include "modarith.h"
00009 #include "algparam.h"
00010 
00011 #include <math.h>
00012 #include <vector>
00013 
00014 NAMESPACE_BEGIN(CryptoPP)
00015 
00016 const unsigned int maxPrimeTableSize = 3511;    // last prime 32719
00017 const word lastSmallPrime = 32719;
00018 unsigned int primeTableSize=552;
00019 
00020 word primeTable[maxPrimeTableSize] =
00021         {2, 3, 5, 7, 11, 13, 17, 19,
00022         23, 29, 31, 37, 41, 43, 47, 53,
00023         59, 61, 67, 71, 73, 79, 83, 89,
00024         97, 101, 103, 107, 109, 113, 127, 131,
00025         137, 139, 149, 151, 157, 163, 167, 173,
00026         179, 181, 191, 193, 197, 199, 211, 223,
00027         227, 229, 233, 239, 241, 251, 257, 263,
00028         269, 271, 277, 281, 283, 293, 307, 311,
00029         313, 317, 331, 337, 347, 349, 353, 359,
00030         367, 373, 379, 383, 389, 397, 401, 409,
00031         419, 421, 431, 433, 439, 443, 449, 457,
00032         461, 463, 467, 479, 487, 491, 499, 503,
00033         509, 521, 523, 541, 547, 557, 563, 569,
00034         571, 577, 587, 593, 599, 601, 607, 613,
00035         617, 619, 631, 641, 643, 647, 653, 659,
00036         661, 673, 677, 683, 691, 701, 709, 719,
00037         727, 733, 739, 743, 751, 757, 761, 769,
00038         773, 787, 797, 809, 811, 821, 823, 827,
00039         829, 839, 853, 857, 859, 863, 877, 881,
00040         883, 887, 907, 911, 919, 929, 937, 941,
00041         947, 953, 967, 971, 977, 983, 991, 997,
00042         1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049,
00043         1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097,
00044         1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163,
00045         1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
00046         1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283,
00047         1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
00048         1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423,
00049         1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459,
00050         1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
00051         1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571,
00052         1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619,
00053         1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693,
00054         1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747,
00055         1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
00056         1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877,
00057         1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949,
00058         1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
00059         2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069,
00060         2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
00061         2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203,
00062         2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267,
00063         2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311,
00064         2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377,
00065         2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
00066         2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503,
00067         2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579,
00068         2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657,
00069         2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,
00070         2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
00071         2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801,
00072         2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861,
00073         2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939,
00074         2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011,
00075         3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
00076         3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167,
00077         3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221,
00078         3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301,
00079         3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347,
00080         3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
00081         3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491,
00082         3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541,
00083         3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607,
00084         3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671,
00085         3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
00086         3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797,
00087         3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863,
00088         3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923,
00089         3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003};
00090 
00091 void BuildPrimeTable()
00092 {
00093         unsigned int p=primeTable[primeTableSize-1];
00094         for (unsigned int i=primeTableSize; i<maxPrimeTableSize; i++)
00095         {
00096                 int j;
00097                 do
00098                 {
00099                         p+=2;
00100                         for (j=1; j<54; j++)
00101                                 if (p%primeTable[j] == 0)
00102                                         break;
00103                 } while (j!=54);
00104                 primeTable[i] = p;
00105         }
00106         primeTableSize = maxPrimeTableSize;
00107         assert(primeTable[primeTableSize-1] == lastSmallPrime);
00108 }
00109 
00110 bool IsSmallPrime(const Integer &p)
00111 {
00112         BuildPrimeTable();
00113 
00114         if (p.IsPositive() && p <= primeTable[primeTableSize-1])
00115                 return std::binary_search(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
00116         else
00117                 return false;
00118 }
00119 
00120 bool TrialDivision(const Integer &p, unsigned bound)
00121 {
00122         assert(primeTable[primeTableSize-1] >= bound);
00123 
00124         unsigned int i;
00125         for (i = 0; primeTable[i]<bound; i++)
00126                 if ((p % primeTable[i]) == 0)
00127                         return true;
00128 
00129         if (bound == primeTable[i])
00130                 return (p % bound == 0);
00131         else
00132                 return false;
00133 }
00134 
00135 bool SmallDivisorsTest(const Integer &p)
00136 {
00137         BuildPrimeTable();
00138         return !TrialDivision(p, primeTable[primeTableSize-1]);
00139 }
00140 
00141 bool IsFermatProbablePrime(const Integer &n, const Integer &b)
00142 {
00143         if (n <= 3)
00144                 return n==2 || n==3;
00145 
00146         assert(n>3 && b>1 && b<n-1);
00147         return a_exp_b_mod_c(b, n-1, n)==1;
00148 }
00149 
00150 bool IsStrongProbablePrime(const Integer &n, const Integer &b)
00151 {
00152         if (n <= 3)
00153                 return n==2 || n==3;
00154 
00155         assert(n>3 && b>1 && b<n-1);
00156 
00157         if ((n.IsEven() && n!=2) || GCD(b, n) != 1)
00158                 return false;
00159 
00160         Integer nminus1 = (n-1);
00161         unsigned int a;
00162 
00163         // calculate a = largest power of 2 that divides (n-1)
00164         for (a=0; ; a++)
00165                 if (nminus1.GetBit(a))
00166                         break;
00167         Integer m = nminus1>>a;
00168 
00169         Integer z = a_exp_b_mod_c(b, m, n);
00170         if (z==1 || z==nminus1)
00171                 return true;
00172         for (unsigned j=1; j<a; j++)
00173         {
00174                 z = z.Squared()%n;
00175                 if (z==nminus1)
00176                         return true;
00177                 if (z==1)
00178                         return false;
00179         }
00180         return false;
00181 }
00182 
00183 bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds)
00184 {
00185         if (n <= 3)
00186                 return n==2 || n==3;
00187 
00188         assert(n>3);
00189 
00190         Integer b;
00191         for (unsigned int i=0; i<rounds; i++)
00192         {
00193                 b.Randomize(rng, 2, n-2);
00194                 if (!IsStrongProbablePrime(n, b))
00195                         return false;
00196         }
00197         return true;
00198 }
00199 
00200 bool IsLucasProbablePrime(const Integer &n)
00201 {
00202         if (n <= 1)
00203                 return false;
00204 
00205         if (n.IsEven())
00206                 return n==2;
00207 
00208         assert(n>2);
00209 
00210         Integer b=3;
00211         unsigned int i=0;
00212         int j;
00213 
00214         while ((j=Jacobi(b.Squared()-4, n)) == 1)
00215         {
00216                 if (++i==64 && n.IsSquare())    // avoid infinite loop if n is a square
00217                         return false;
00218                 ++b; ++b;
00219         }
00220 
00221         if (j==0) 
00222                 return false;
00223         else
00224                 return Lucas(n+1, b, n)==2;
00225 }
00226 
00227 bool IsStrongLucasProbablePrime(const Integer &n)
00228 {
00229         if (n <= 1)
00230                 return false;
00231 
00232         if (n.IsEven())
00233                 return n==2;
00234 
00235         assert(n>2);
00236 
00237         Integer b=3;
00238         unsigned int i=0;
00239         int j;
00240 
00241         while ((j=Jacobi(b.Squared()-4, n)) == 1)
00242         {
00243                 if (++i==64 && n.IsSquare())    // avoid infinite loop if n is a square
00244                         return false;
00245                 ++b; ++b;
00246         }
00247 
00248         if (j==0) 
00249                 return false;
00250 
00251         Integer n1 = n+1;
00252         unsigned int a;
00253 
00254         // calculate a = largest power of 2 that divides n1
00255         for (a=0; ; a++)
00256                 if (n1.GetBit(a))
00257                         break;
00258         Integer m = n1>>a;
00259 
00260         Integer z = Lucas(m, b, n);
00261         if (z==2 || z==n-2)
00262                 return true;
00263         for (i=1; i<a; i++)
00264         {
00265                 z = (z.Squared()-2)%n;
00266                 if (z==n-2)
00267                         return true;
00268                 if (z==2)
00269                         return false;
00270         }
00271         return false;
00272 }
00273 
00274 bool IsPrime(const Integer &p)
00275 {
00276         static const Integer lastSmallPrimeSquared = Integer(lastSmallPrime).Squared();
00277 
00278         if (p <= lastSmallPrime)
00279                 return IsSmallPrime(p);
00280         else if (p <= lastSmallPrimeSquared)
00281                 return SmallDivisorsTest(p);
00282         else
00283                 return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p);
00284 }
00285 
00286 bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level)
00287 {
00288         bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1);
00289         if (level >= 1)
00290                 pass = pass && RabinMillerTest(rng, p, 10);
00291         return pass;
00292 }
00293 
00294 unsigned int PrimeSearchInterval(const Integer &max)
00295 {
00296         return max.BitCount();
00297 }
00298 
00299 static inline bool FastProbablePrimeTest(const Integer &n)
00300 {
00301         return IsStrongProbablePrime(n,2);
00302 }
00303 
00304 AlgorithmParameters<AlgorithmParameters<AlgorithmParameters<NullNameValuePairs, Integer::RandomNumberType>, Integer>, Integer>
00305         MakeParametersForTwoPrimesOfEqualSize(unsigned int productBitLength)
00306 {
00307         if (productBitLength < 16)
00308                 throw InvalidArgument("invalid bit length");
00309 
00310         Integer minP, maxP;
00311 
00312         if (productBitLength%2==0)
00313         {
00314                 minP = Integer(182) << (productBitLength/2-8);
00315                 maxP = Integer::Power2(productBitLength/2)-1;
00316         }
00317         else
00318         {
00319                 minP = Integer::Power2((productBitLength-1)/2);
00320                 maxP = Integer(181) << ((productBitLength+1)/2-8);
00321         }
00322 
00323         return MakeParameters("RandomNumberType", Integer::PRIME)("Min", minP)("Max", maxP);
00324 }
00325 
00326 class PrimeSieve
00327 {
00328 public:
00329         // delta == 1 or -1 means double sieve with p = 2*q + delta
00330         PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta=0);
00331         bool NextCandidate(Integer &c);
00332 
00333         void DoSieve();
00334         static void SieveSingle(std::vector<bool> &sieve, word p, const Integer &first, const Integer &step, word stepInv);
00335 
00336         Integer m_first, m_last, m_step;
00337         signed int m_delta;
00338         word m_next;
00339         std::vector<bool> m_sieve;
00340 };
00341 
00342 PrimeSieve::PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta)
00343         : m_first(first), m_last(last), m_step(step), m_delta(delta), m_next(0)
00344 {
00345         DoSieve();
00346 }
00347 
00348 bool PrimeSieve::NextCandidate(Integer &c)
00349 {
00350         m_next = std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin();
00351         if (m_next == m_sieve.size())
00352         {
00353                 m_first += m_sieve.size()*m_step;
00354                 if (m_first > m_last)
00355                         return false;
00356                 else
00357                 {
00358                         m_next = 0;
00359                         DoSieve();
00360                         return NextCandidate(c);
00361                 }
00362         }
00363         else
00364         {
00365                 c = m_first + m_next*m_step;
00366                 ++m_next;
00367                 return true;
00368         }
00369 }
00370 
00371 void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word p, const Integer &first, const Integer &step, word stepInv)
00372 {
00373         if (stepInv)
00374         {
00375                 unsigned int sieveSize = sieve.size();
00376                 word j = word((dword(p-(first%p))*stepInv) % p);
00377                 // if the first multiple of p is p, skip it
00378                 if (first.WordCount() <= 1 && first + step*j == p)
00379                         j += p;
00380                 for (; j < sieveSize; j += p)
00381                         sieve[j] = true;
00382         }
00383 }
00384 
00385 void PrimeSieve::DoSieve()
00386 {
00387         BuildPrimeTable();
00388 
00389         const unsigned int maxSieveSize = 32768;
00390         unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong();
00391 
00392         m_sieve.clear();
00393         m_sieve.resize(sieveSize, false);
00394 
00395         if (m_delta == 0)
00396         {
00397                 for (unsigned int i = 0; i < primeTableSize; ++i)
00398                         SieveSingle(m_sieve, primeTable[i], m_first, m_step, m_step.InverseMod(primeTable[i]));
00399         }
00400         else
00401         {
00402                 assert(m_step%2==0);
00403                 Integer qFirst = (m_first-m_delta) >> 1;
00404                 Integer halfStep = m_step >> 1;
00405                 for (unsigned int i = 0; i < primeTableSize; ++i)
00406                 {
00407                         word p = primeTable[i];
00408                         word stepInv = m_step.InverseMod(p);
00409                         SieveSingle(m_sieve, p, m_first, m_step, stepInv);
00410 
00411                         word halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p;
00412                         SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
00413                 }
00414         }
00415 }
00416 
00417 bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector)
00418 {
00419         assert(!equiv.IsNegative() && equiv < mod);
00420 
00421         Integer gcd = GCD(equiv, mod);
00422         if (gcd != Integer::One())
00423         {
00424                 // the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv)
00425                 if (p <= gcd && gcd <= max && IsPrime(gcd))
00426                 {
00427                         p = gcd;
00428                         return true;
00429                 }
00430                 else
00431                         return false;
00432         }
00433 
00434         BuildPrimeTable();
00435 
00436         if (p <= primeTable[primeTableSize-1])
00437         {
00438                 word *pItr;
00439 
00440                 --p;
00441                 if (p.IsPositive())
00442                         pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
00443                 else
00444                         pItr = primeTable;
00445 
00446                 while (pItr < primeTable+primeTableSize && *pItr%mod != equiv)
00447                         ++pItr;
00448 
00449                 if (pItr < primeTable+primeTableSize)
00450                 {
00451                         p = *pItr;
00452                         return p <= max;
00453                 }
00454 
00455                 p = primeTable[primeTableSize-1]+1;
00456         }
00457 
00458         assert(p > primeTable[primeTableSize-1]);
00459 
00460         if (mod.IsOdd())
00461                 return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector);
00462 
00463         p += (equiv-p)%mod;
00464 
00465         if (p>max)
00466                 return false;
00467 
00468         PrimeSieve sieve(p, max, mod);
00469 
00470         while (sieve.NextCandidate(p))
00471         {
00472                 if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p))
00473                         return true;
00474         }
00475 
00476         return false;
00477 }
00478 
00479 // the following two functions are based on code and comments provided by Preda Mihailescu
00480 static bool ProvePrime(const Integer &p, const Integer &q)
00481 {
00482         assert(p < q*q*q);
00483         assert(p % q == 1);
00484 
00485 // this is the Quisquater test. Numbers p having passed the Lucas - Lehmer test
00486 // for q and verifying p < q^3 can only be built up of two factors, both = 1 mod q,
00487 // or be prime. The next two lines build the discriminant of a quadratic equation
00488 // which holds iff p is built up of two factors (excercise ... )
00489 
00490         Integer r = (p-1)/q;
00491         if (((r%q).Squared()-4*(r/q)).IsSquare())
00492                 return false;
00493 
00494         assert(primeTableSize >= 50);
00495         for (int i=0; i<50; i++) 
00496         {
00497                 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
00498                 if (b != 1) 
00499                         return a_exp_b_mod_c(b, q, p) == 1;
00500         }
00501         return false;
00502 }
00503 
00504 Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits)
00505 {
00506         Integer p;
00507         Integer minP = Integer::Power2(pbits-1);
00508         Integer maxP = Integer::Power2(pbits) - 1;
00509 
00510         if (maxP <= Integer(lastSmallPrime).Squared())
00511         {
00512                 // Randomize() will generate a prime provable by trial division
00513                 p.Randomize(rng, minP, maxP, Integer::PRIME);
00514                 return p;
00515         }
00516 
00517         unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36);
00518         Integer q = MihailescuProvablePrime(rng, qbits);
00519         Integer q2 = q<<1;
00520 
00521         while (true)
00522         {
00523                 // this initializes the sieve to search in the arithmetic
00524                 // progression p = p_0 + \lambda * q2 = p_0 + 2 * \lambda * q,
00525                 // with q the recursively generated prime above. We will be able
00526                 // to use Lucas tets for proving primality. A trick of Quisquater
00527                 // allows taking q > cubic_root(p) rather then square_root: this
00528                 // decreases the recursion.
00529 
00530                 p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2);
00531                 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
00532 
00533                 while (sieve.NextCandidate(p))
00534                 {
00535                         if (FastProbablePrimeTest(p) && ProvePrime(p, q))
00536                                 return p;
00537                 }
00538         }
00539 
00540         // not reached
00541         return p;
00542 }
00543 
00544 Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
00545 {
00546         const unsigned smallPrimeBound = 29, c_opt=10;
00547         Integer p;
00548 
00549         BuildPrimeTable();
00550         if (bits < smallPrimeBound)
00551         {
00552                 do
00553                         p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2);
00554                 while (TrialDivision(p, 1 << ((bits+1)/2)));
00555         }
00556         else
00557         {
00558                 const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
00559                 double relativeSize;
00560                 do
00561                         relativeSize = pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1);
00562                 while (bits * relativeSize >= bits - margin);
00563 
00564                 Integer a,b;
00565                 Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize));
00566                 Integer I = Integer::Power2(bits-2)/q;
00567                 Integer I2 = I << 1;
00568                 unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt);
00569                 bool success = false;
00570                 while (!success)
00571                 {
00572                         p.Randomize(rng, I, I2, Integer::ANY);
00573                         p *= q; p <<= 1; ++p;
00574                         if (!TrialDivision(p, trialDivisorBound))
00575                         {
00576                                 a.Randomize(rng, 2, p-1, Integer::ANY);
00577                                 b = a_exp_b_mod_c(a, (p-1)/q, p);
00578                                 success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1);
00579                         }
00580                 }
00581         }
00582         return p;
00583 }
00584 
00585 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u)
00586 {
00587         // isn't operator overloading great?
00588         return p * (u * (xq-xp) % q) + xp;
00589 }
00590 
00591 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q)
00592 {
00593         return CRT(xp, p, xq, q, EuclideanMultiplicativeInverse(p, q));
00594 }
00595 
00596 Integer ModularSquareRoot(const Integer &a, const Integer &p)
00597 {
00598         if (p%4 == 3)
00599                 return a_exp_b_mod_c(a, (p+1)/4, p);
00600 
00601         Integer q=p-1;
00602         unsigned int r=0;
00603         while (q.IsEven())
00604         {
00605                 r++;
00606                 q >>= 1;
00607         }
00608 
00609         Integer n=2;
00610         while (Jacobi(n, p) != -1)
00611                 ++n;
00612 
00613         Integer y = a_exp_b_mod_c(n, q, p);
00614         Integer x = a_exp_b_mod_c(a, (q-1)/2, p);
00615         Integer b = (x.Squared()%p)*a%p;
00616         x = a*x%p;
00617         Integer tempb, t;
00618 
00619         while (b != 1)
00620         {
00621                 unsigned m=0;
00622                 tempb = b;
00623                 do
00624                 {
00625                         m++;
00626                         b = b.Squared()%p;
00627                         if (m==r)
00628                                 return Integer::Zero();
00629                 }
00630                 while (b != 1);
00631 
00632                 t = y;
00633                 for (unsigned i=0; i<r-m-1; i++)
00634                         t = t.Squared()%p;
00635                 y = t.Squared()%p;
00636                 r = m;
00637                 x = x*t%p;
00638                 b = tempb*y%p;
00639         }
00640 
00641         assert(x.Squared()%p == a);
00642         return x;
00643 }
00644 
00645 bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p)
00646 {
00647         Integer D = (b.Squared() - 4*a*c) % p;
00648         switch (Jacobi(D, p))
00649         {
00650         default:
00651                 assert(false);  // not reached
00652                 return false;
00653         case -1:
00654                 return false;
00655         case 0:
00656                 r1 = r2 = (-b*(a+a).InverseMod(p)) % p;
00657                 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
00658                 return true;
00659         case 1:
00660                 Integer s = ModularSquareRoot(D, p);
00661                 Integer t = (a+a).InverseMod(p);
00662                 r1 = (s-b)*t % p;
00663                 r2 = (-s-b)*t % p;
00664                 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
00665                 assert(((r2.Squared()*a + r2*b + c) % p).IsZero());
00666                 return true;
00667         }
00668 }
00669 
00670 Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq,
00671                                         const Integer &p, const Integer &q, const Integer &u)
00672 {
00673         Integer p2 = ModularExponentiation((a % p), dp, p);
00674         Integer q2 = ModularExponentiation((a % q), dq, q);
00675         return CRT(p2, p, q2, q, u);
00676 }
00677 
00678 Integer ModularRoot(const Integer &a, const Integer &e,
00679                                         const Integer &p, const Integer &q)
00680 {
00681         Integer dp = EuclideanMultiplicativeInverse(e, p-1);
00682         Integer dq = EuclideanMultiplicativeInverse(e, q-1);
00683         Integer u = EuclideanMultiplicativeInverse(p, q);
00684         assert(!!dp && !!dq && !!u);
00685         return ModularRoot(a, dp, dq, p, q, u);
00686 }
00687 
00688 /*
00689 Integer GCDI(const Integer &x, const Integer &y)
00690 {
00691         Integer a=x, b=y;
00692         unsigned k=0;
00693 
00694         assert(!!a && !!b);
00695 
00696         while (a[0]==0 && b[0]==0)
00697         {
00698                 a >>= 1;
00699                 b >>= 1;
00700                 k++;
00701         }
00702 
00703         while (a[0]==0)
00704                 a >>= 1;
00705 
00706         while (b[0]==0)
00707                 b >>= 1;
00708 
00709         while (1)
00710         {
00711                 switch (a.Compare(b))
00712                 {
00713                         case -1:
00714                                 b -= a;
00715                                 while (b[0]==0)
00716                                         b >>= 1;
00717                                 break;
00718 
00719                         case 0:
00720                                 return (a <<= k);
00721 
00722                         case 1:
00723                                 a -= b;
00724                                 while (a[0]==0)
00725                                         a >>= 1;
00726                                 break;
00727 
00728                         default:
00729                                 assert(false);
00730                 }
00731         }
00732 }
00733 
00734 Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b)
00735 {
00736         assert(b.Positive());
00737 
00738         if (a.Negative())
00739                 return EuclideanMultiplicativeInverse(a%b, b);
00740 
00741         if (b[0]==0)
00742         {
00743                 if (!b || a[0]==0)
00744                         return Integer::Zero();       // no inverse
00745                 if (a==1)
00746                         return 1;
00747                 Integer u = EuclideanMultiplicativeInverse(b, a);
00748                 if (!u)
00749                         return Integer::Zero();       // no inverse
00750                 else
00751                         return (b*(a-u)+1)/a;
00752         }
00753 
00754         Integer u=1, d=a, v1=b, v3=b, t1, t3, b2=(b+1)>>1;
00755 
00756         if (a[0])
00757         {
00758                 t1 = Integer::Zero();
00759                 t3 = -b;
00760         }
00761         else
00762         {
00763                 t1 = b2;
00764                 t3 = a>>1;
00765         }
00766 
00767         while (!!t3)
00768         {
00769                 while (t3[0]==0)
00770                 {
00771                         t3 >>= 1;
00772                         if (t1[0]==0)
00773                                 t1 >>= 1;
00774                         else
00775                         {
00776                                 t1 >>= 1;
00777                                 t1 += b2;
00778                         }
00779                 }
00780                 if (t3.Positive())
00781                 {
00782                         u = t1;
00783                         d = t3;
00784                 }
00785                 else
00786                 {
00787                         v1 = b-t1;
00788                         v3 = -t3;
00789                 }
00790                 t1 = u-v1;
00791                 t3 = d-v3;
00792                 if (t1.Negative())
00793                         t1 += b;
00794         }
00795         if (d==1)
00796                 return u;
00797         else
00798                 return Integer::Zero();   // no inverse
00799 }
00800 */
00801 
00802 int Jacobi(const Integer &aIn, const Integer &bIn)
00803 {
00804         assert(bIn.IsOdd());
00805 
00806         Integer b = bIn, a = aIn%bIn;
00807         int result = 1;
00808 
00809         while (!!a)
00810         {
00811                 unsigned i=0;
00812                 while (a.GetBit(i)==0)
00813                         i++;
00814                 a>>=i;
00815 
00816                 if (i%2==1 && (b%8==3 || b%8==5))
00817                         result = -result;
00818 
00819                 if (a%4==3 && b%4==3)
00820                         result = -result;
00821 
00822                 std::swap(a, b);
00823                 a %= b;
00824         }
00825 
00826         return (b==1) ? result : 0;
00827 }
00828 
00829 Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n)
00830 {
00831         unsigned i = e.BitCount();
00832         if (i==0)
00833                 return Integer::Two();
00834 
00835         MontgomeryRepresentation m(n);
00836         Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two());
00837         Integer v=p, v1=m.Subtract(m.Square(p), two);
00838 
00839         i--;
00840         while (i--)
00841         {
00842                 if (e.GetBit(i))
00843                 {
00844                         // v = (v*v1 - p) % m;
00845                         v = m.Subtract(m.Multiply(v,v1), p);
00846                         // v1 = (v1*v1 - 2) % m;
00847                         v1 = m.Subtract(m.Square(v1), two);
00848                 }
00849                 else
00850                 {
00851                         // v1 = (v*v1 - p) % m;
00852                         v1 = m.Subtract(m.Multiply(v,v1), p);
00853                         // v = (v*v - 2) % m;
00854                         v = m.Subtract(m.Square(v), two);
00855                 }
00856         }
00857         return m.ConvertOut(v);
00858 }
00859 
00860 // This is Peter Montgomery's unpublished Lucas sequence evalutation algorithm.
00861 // The total number of multiplies and squares used is less than the binary
00862 // algorithm (see above).  Unfortunately I can't get it to run as fast as
00863 // the binary algorithm because of the extra overhead.
00864 /*
00865 Integer Lucas(const Integer &n, const Integer &P, const Integer &modulus)
00866 {
00867         if (!n)
00868                 return 2;
00869 
00870 #define f(A, B, C)      m.Subtract(m.Multiply(A, B), C)
00871 #define X2(A) m.Subtract(m.Square(A), two)
00872 #define X3(A) m.Multiply(A, m.Subtract(m.Square(A), three))
00873 
00874         MontgomeryRepresentation m(modulus);
00875         Integer two=m.ConvertIn(2), three=m.ConvertIn(3);
00876         Integer A=m.ConvertIn(P), B, C, p, d=n, e, r, t, T, U;
00877 
00878         while (d!=1)
00879         {
00880                 p = d;
00881                 unsigned int b = WORD_BITS * p.WordCount();
00882                 Integer alpha = (Integer(5)<<(2*b-2)).SquareRoot() - Integer::Power2(b-1);
00883                 r = (p*alpha)>>b;
00884                 e = d-r;
00885                 B = A;
00886                 C = two;
00887                 d = r;
00888 
00889                 while (d!=e)
00890                 {
00891                         if (d<e)
00892                         {
00893                                 swap(d, e);
00894                                 swap(A, B);
00895                         }
00896 
00897                         unsigned int dm2 = d[0], em2 = e[0];
00898                         unsigned int dm3 = d%3, em3 = e%3;
00899 
00900 //                      if ((dm6+em6)%3 == 0 && d <= e + (e>>2))
00901                         if ((dm3+em3==0 || dm3+em3==3) && (t = e, t >>= 2, t += e, d <= t))
00902                         {
00903                                 // #1
00904 //                              t = (d+d-e)/3;
00905 //                              t = d; t += d; t -= e; t /= 3;
00906 //                              e = (e+e-d)/3;
00907 //                              e += e; e -= d; e /= 3;
00908 //                              d = t;
00909 
00910 //                              t = (d+e)/3
00911                                 t = d; t += e; t /= 3;
00912                                 e -= t;
00913                                 d -= t;
00914 
00915                                 T = f(A, B, C);
00916                                 U = f(T, A, B);
00917                                 B = f(T, B, A);
00918                                 A = U;
00919                                 continue;
00920                         }
00921 
00922 //                      if (dm6 == em6 && d <= e + (e>>2))
00923                         if (dm3 == em3 && dm2 == em2 && (t = e, t >>= 2, t += e, d <= t))
00924                         {
00925                                 // #2
00926 //                              d = (d-e)>>1;
00927                                 d -= e; d >>= 1;
00928                                 B = f(A, B, C);
00929                                 A = X2(A);
00930                                 continue;
00931                         }
00932 
00933 //                      if (d <= (e<<2))
00934                         if (d <= (t = e, t <<= 2))
00935                         {
00936                                 // #3
00937                                 d -= e;
00938                                 C = f(A, B, C);
00939                                 swap(B, C);
00940                                 continue;
00941                         }
00942 
00943                         if (dm2 == em2)
00944                         {
00945                                 // #4
00946 //                              d = (d-e)>>1;
00947                                 d -= e; d >>= 1;
00948                                 B = f(A, B, C);
00949                                 A = X2(A);
00950                                 continue;
00951                         }
00952 
00953                         if (dm2 == 0)
00954                         {
00955                                 // #5
00956                                 d >>= 1;
00957                                 C = f(A, C, B);
00958                                 A = X2(A);
00959                                 continue;
00960                         }
00961 
00962                         if (dm3 == 0)
00963                         {
00964                                 // #6
00965 //                              d = d/3 - e;
00966                                 d /= 3; d -= e;
00967                                 T = X2(A);
00968                                 C = f(T, f(A, B, C), C);
00969                                 swap(B, C);
00970                                 A = f(T, A, A);
00971                                 continue;
00972                         }
00973 
00974                         if (dm3+em3==0 || dm3+em3==3)
00975                         {
00976                                 // #7
00977 //                              d = (d-e-e)/3;
00978                                 d -= e; d -= e; d /= 3;
00979                                 T = f(A, B, C);
00980                                 B = f(T, A, B);
00981                                 A = X3(A);
00982                                 continue;
00983                         }
00984 
00985                         if (dm3 == em3)
00986                         {
00987                                 // #8
00988 //                              d = (d-e)/3;
00989                                 d -= e; d /= 3;
00990                                 T = f(A, B, C);
00991                                 C = f(A, C, B);
00992                                 B = T;
00993                                 A = X3(A);
00994                                 continue;
00995                         }
00996 
00997                         assert(em2 == 0);
00998                         // #9
00999                         e >>= 1;
01000                         C = f(C, B, A);
01001                         B = X2(B);
01002                 }
01003 
01004                 A = f(A, B, C);
01005         }
01006 
01007 #undef f
01008 #undef X2
01009 #undef X3
01010 
01011         return m.ConvertOut(A);
01012 }
01013 */
01014 
01015 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
01016 {
01017         Integer d = (m*m-4);
01018         Integer p2 = p-Jacobi(d,p);
01019         Integer q2 = q-Jacobi(d,q);
01020         return CRT(Lucas(EuclideanMultiplicativeInverse(e,p2), m, p), p, Lucas(EuclideanMultiplicativeInverse(e,q2), m, q), q, u);
01021 }
01022 
01023 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q)
01024 {
01025         return InverseLucas(e, m, p, q, EuclideanMultiplicativeInverse(p, q));
01026 }
01027 
01028 unsigned int FactoringWorkFactor(unsigned int n)
01029 {
01030         // extrapolated from the table in Odlyzko's "The Future of Integer Factorization"
01031         // updated to reflect the factoring of RSA-130
01032         if (n<5) return 0;
01033         else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
01034 }
01035 
01036 unsigned int DiscreteLogWorkFactor(unsigned int n)
01037 {
01038         // assuming discrete log takes about the same time as factoring
01039         if (n<5) return 0;
01040         else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
01041 }
01042 
01043 // ********************************************************
01044 
01045 void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits)
01046 {
01047         // no prime exists for delta = -1, qbits = 4, and pbits = 5
01048         assert(qbits > 4);
01049         assert(pbits > qbits);
01050 
01051         if (qbits+1 == pbits)
01052         {
01053                 Integer minP = Integer::Power2(pbits-1);
01054                 Integer maxP = Integer::Power2(pbits) - 1;
01055                 bool success = false;
01056 
01057                 while (!success)
01058                 {
01059                         p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
01060                         PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
01061 
01062                         while (sieve.NextCandidate(p))
01063                         {
01064                                 assert(IsSmallPrime(p) || SmallDivisorsTest(p));
01065                                 q = (p-delta) >> 1;
01066                                 assert(IsSmallPrime(q) || SmallDivisorsTest(q));
01067                                 if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
01068                                 {
01069                                         success = true;
01070                                         break;
01071                                 }
01072                         }
01073                 }
01074 
01075                 if (delta == 1)
01076                 {
01077                         // find g such that g is a quadratic residue mod p, then g has order q
01078                         // g=4 always works, but this way we get the smallest quadratic residue (other than 1)
01079                         for (g=2; Jacobi(g, p) != 1; ++g) {}
01080                         // contributed by Walt Tuvell: g should be the following according to the Law of Quadratic Reciprocity
01081                         assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
01082                 }
01083                 else
01084                 {
01085                         assert(delta == -1);
01086                         // find g such that g*g-4 is a quadratic non-residue, 
01087                         // and such that g has order q
01088                         for (g=3; ; ++g)
01089                                 if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
01090                                         break;
01091                 }
01092         }
01093         else
01094         {
01095                 Integer minQ = Integer::Power2(qbits-1);
01096                 Integer maxQ = Integer::Power2(qbits) - 1;
01097                 Integer minP = Integer::Power2(pbits-1);
01098                 Integer maxP = Integer::Power2(pbits) - 1;
01099 
01100                 do
01101                 {
01102                         q.Randomize(rng, minQ, maxQ, Integer::PRIME);
01103                 } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
01104 
01105                 // find a random g of order q
01106                 if (delta==1)
01107                 {
01108                         do
01109                         {
01110                                 Integer h(rng, 2, p-2, Integer::ANY);
01111                                 g = a_exp_b_mod_c(h, (p-1)/q, p);
01112                         } while (g <= 1);
01113                         assert(a_exp_b_mod_c(g, q, p)==1);
01114                 }
01115                 else
01116                 {
01117                         assert(delta==-1);
01118                         do
01119                         {
01120                                 Integer h(rng, 3, p-1, Integer::ANY);
01121                                 if (Jacobi(h*h-4, p)==1)
01122                                         continue;
01123                                 g = Lucas((p+1)/q, h, p);
01124                         } while (g <= 2);
01125                         assert(Lucas(q, g, p) == 2);
01126                 }
01127         }
01128 }
01129 
01130 NAMESPACE_END
01131 
01132 #endif

Generated on Tue Jul 8 23:34:21 2003 for Crypto++ by doxygen 1.3.2