- Main Page
- Namespaces
- Classes
- Files

- File List
- 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 word s_lastSmallPrime = 32719;
    00017
    00018 struct NewPrimeTable
    00019 {
    00020         std::vector<word16> * operator()() const
    00021         {
    00022                 const unsigned int maxPrimeTableSize = 3511;
    00023
    00024                 std::auto_ptr<std::vector<word16> > pPrimeTable(new std::vector<word16>);
    00025                 std::vector<word16> &primeTable = *pPrimeTable;
    00026                 primeTable.reserve(maxPrimeTableSize);
    00027
    00028                 primeTable.push_back(2);
    00029                 unsigned int testEntriesEnd = 1;
    00030
    00031                 for (unsigned int p=3; p<=s_lastSmallPrime; p+=2)
    00032                 {
    00033                         unsigned int j;
    00034                         for (j=1; j<testEntriesEnd; j++)
    00035                                 if (p%primeTable[j] == 0)
    00036                                         break;
    00037                         if (j == testEntriesEnd)
    00038                         {
    00039                                 primeTable.push_back(p);
    00040                                 testEntriesEnd = UnsignedMin(54U, primeTable.size());
    00041                         }
    00042                 }
    00043
    00044                 return pPrimeTable.release();
    00045         }
    00046 };
    00047
    00048 const word16 * GetPrimeTable(unsigned int &size)
    00049 {
    00050         const std::vector<word16> &primeTable = Singleton<std::vector<word16>, NewPrimeTable>().Ref();
    00051         size = (unsigned int)primeTable.size();
    00052         return &primeTable[0];
    00053 }
    00054
    00055 bool IsSmallPrime(const Integer &p)
    00056 {
    00057         unsigned int primeTableSize;
    00058         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00059
    00060         if (p.IsPositive() && p <= primeTable[primeTableSize-1])
    00061                 return std::binary_search(primeTable, primeTable+primeTableSize, (word16)p.ConvertToLong());
    00062         else
    00063                 return false;
    00064 }
    00065
    00066 bool TrialDivision(const Integer &p, unsigned bound)
    00067 {
    00068         unsigned int primeTableSize;
    00069         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00070
    00071         assert(primeTable[primeTableSize-1] >= bound);
    00072
    00073         unsigned int i;
    00074         for (i = 0; primeTable[i]<bound; i++)
    00075                 if ((p % primeTable[i]) == 0)
    00076                         return true;
    00077
    00078         if (bound == primeTable[i])
    00079                 return (p % bound == 0);
    00080         else
    00081                 return false;
    00082 }
    00083
    00084 bool SmallDivisorsTest(const Integer &p)
    00085 {
    00086         unsigned int primeTableSize;
    00087         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00088         return !TrialDivision(p, primeTable[primeTableSize-1]);
    00089 }
    00090
    00091 bool IsFermatProbablePrime(const Integer &n, const Integer &b)
    00092 {
    00093         if (n <= 3)
    00094                 return n==2 || n==3;
    00095
    00096         assert(n>3 && b>1 && b<n-1);
    00097         return a_exp_b_mod_c(b, n-1, n)==1;
    00098 }
    00099
    00100 bool IsStrongProbablePrime(const Integer &n, const Integer &b)
    00101 {
    00102         if (n <= 3)
    00103                 return n==2 || n==3;
    00104
    00105         assert(n>3 && b>1 && b<n-1);
    00106
    00107         if ((n.IsEven() && n!=2) || GCD(b, n) != 1)
    00108                 return false;
    00109
    00110         Integer nminus1 = (n-1);
    00111         unsigned int a;
    00112
    00113         // calculate a = largest power of 2 that divides (n-1)
    00114         for (a=0; ; a++)
    00115                 if (nminus1.GetBit(a))
    00116                         break;
    00117         Integer m = nminus1>>a;
    00118
    00119         Integer z = a_exp_b_mod_c(b, m, n);
    00120         if (z==1 || z==nminus1)
    00121                 return true;
    00122         for (unsigned j=1; j<a; j++)
    00123         {
    00124                 z = z.Squared()%n;
    00125                 if (z==nminus1)
    00126                         return true;
    00127                 if (z==1)
    00128                         return false;
    00129         }
    00130         return false;
    00131 }
    00132
    00133 bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds)
    00134 {
    00135         if (n <= 3)
    00136                 return n==2 || n==3;
    00137
    00138         assert(n>3);
    00139
    00140         Integer b;
    00141         for (unsigned int i=0; i<rounds; i++)
    00142         {
    00143                 b.Randomize(rng, 2, n-2);
    00144                 if (!IsStrongProbablePrime(n, b))
    00145                         return false;
    00146         }
    00147         return true;
    00148 }
    00149
    00150 bool IsLucasProbablePrime(const Integer &n)
    00151 {
    00152         if (n <= 1)
    00153                 return false;
    00154
    00155         if (n.IsEven())
    00156                 return n==2;
    00157
    00158         assert(n>2);
    00159
    00160         Integer b=3;
    00161         unsigned int i=0;
    00162         int j;
    00163
    00164         while ((j=Jacobi(b.Squared()-4, n)) == 1)
    00165         {
    00166                 if (++i==64 && n.IsSquare())    // avoid infinite loop if n is a square
    00167                         return false;
    00168                 ++b; ++b;
    00169         }
    00170
    00171         if (j==0)
    00172                 return false;
    00173         else
    00174                 return Lucas(n+1, b, n)==2;
    00175 }
    00176
    00177 bool IsStrongLucasProbablePrime(const Integer &n)
    00178 {
    00179         if (n <= 1)
    00180                 return false;
    00181
    00182         if (n.IsEven())
    00183                 return n==2;
    00184
    00185         assert(n>2);
    00186
    00187         Integer b=3;
    00188         unsigned int i=0;
    00189         int j;
    00190
    00191         while ((j=Jacobi(b.Squared()-4, n)) == 1)
    00192         {
    00193                 if (++i==64 && n.IsSquare())    // avoid infinite loop if n is a square
    00194                         return false;
    00195                 ++b; ++b;
    00196         }
    00197
    00198         if (j==0)
    00199                 return false;
    00200
    00201         Integer n1 = n+1;
    00202         unsigned int a;
    00203
    00204         // calculate a = largest power of 2 that divides n1
    00205         for (a=0; ; a++)
    00206                 if (n1.GetBit(a))
    00207                         break;
    00208         Integer m = n1>>a;
    00209
    00210         Integer z = Lucas(m, b, n);
    00211         if (z==2 || z==n-2)
    00212                 return true;
    00213         for (i=1; i<a; i++)
    00214         {
    00215                 z = (z.Squared()-2)%n;
    00216                 if (z==n-2)
    00217                         return true;
    00218                 if (z==2)
    00219                         return false;
    00220         }
    00221         return false;
    00222 }
    00223
    00224 struct NewLastSmallPrimeSquared
    00225 {
    00226         Integer * operator()() const
    00227         {
    00228                 return new Integer(Integer(s_lastSmallPrime).Squared());
    00229         }
    00230 };
    00231
    00232 bool IsPrime(const Integer &p)
    00233 {
    00234         if (p <= s_lastSmallPrime)
    00235                 return IsSmallPrime(p);
    00236         else if (p <= Singleton<Integer, NewLastSmallPrimeSquared>().Ref())
    00237                 return SmallDivisorsTest(p);
    00238         else
    00239                 return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p);
    00240 }
    00241
    00242 bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level)
    00243 {
    00244         bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1);
    00245         if (level >= 1)
    00246                 pass = pass && RabinMillerTest(rng, p, 10);
    00247         return pass;
    00248 }
    00249
    00250 unsigned int PrimeSearchInterval(const Integer &max)
    00251 {
    00252         return max.BitCount();
    00253 }
    00254
    00255 static inline bool FastProbablePrimeTest(const Integer &n)
    00256 {
    00257         return IsStrongProbablePrime(n,2);
    00258 }
    00259
    00260 AlgorithmParameters<AlgorithmParameters<AlgorithmParameters<NullNameValuePairs, Integer::RandomNumberType>, Integer>, Integer>
    00261         MakeParametersForTwoPrimesOfEqualSize(unsigned int productBitLength)
    00262 {
    00263         if (productBitLength < 16)
    00264                 throw InvalidArgument("invalid bit length");
    00265
    00266         Integer minP, maxP;
    00267
    00268         if (productBitLength%2==0)
    00269         {
    00270                 minP = Integer(182) << (productBitLength/2-8);
    00271                 maxP = Integer::Power2(productBitLength/2)-1;
    00272         }
    00273         else
    00274         {
    00275                 minP = Integer::Power2((productBitLength-1)/2);
    00276                 maxP = Integer(181) << ((productBitLength+1)/2-8);
    00277         }
    00278
    00279         return MakeParameters("RandomNumberType", Integer::PRIME)("Min", minP)("Max", maxP);
    00280 }
    00281
    00282 class PrimeSieve
    00283 {
    00284 public:
    00285         // delta == 1 or -1 means double sieve with p = 2*q + delta
    00286         PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta=0);
    00287         bool NextCandidate(Integer &c);
    00288
    00289         void DoSieve();
    00290         static void SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv);
    00291
    00292         Integer m_first, m_last, m_step;
    00293         signed int m_delta;
    00294         word m_next;
    00295         std::vector<bool> m_sieve;
    00296 };
    00297
    00298 PrimeSieve::PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta)
    00299         : m_first(first), m_last(last), m_step(step), m_delta(delta), m_next(0)
    00300 {
    00301         DoSieve();
    00302 }
    00303
    00304 bool PrimeSieve::NextCandidate(Integer &c)
    00305 {
    00306         bool safe = SafeConvert(std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin(), m_next);
    00307         assert(safe);
    00308         if (m_next == m_sieve.size())
    00309         {
    00310                 m_first += long(m_sieve.size())*m_step;
    00311                 if (m_first > m_last)
    00312                         return false;
    00313                 else
    00314                 {
    00315                         m_next = 0;
    00316                         DoSieve();
    00317                         return NextCandidate(c);
    00318                 }
    00319         }
    00320         else
    00321         {
    00322                 c = m_first + long(m_next)*m_step;
    00323                 ++m_next;
    00324                 return true;
    00325         }
    00326 }
    00327
    00328 void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv)
    00329 {
    00330         if (stepInv)
    00331         {
    00332                 size_t sieveSize = sieve.size();
    00333                 size_t j = (word32(p-(first%p))*stepInv) % p;
    00334                 // if the first multiple of p is p, skip it
    00335                 if (first.WordCount() <= 1 && first + step*long(j) == p)
    00336                         j += p;
    00337                 for (; j < sieveSize; j += p)
    00338                         sieve[j] = true;
    00339         }
    00340 }
    00341
    00342 void PrimeSieve::DoSieve()
    00343 {
    00344         unsigned int primeTableSize;
    00345         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00346
    00347         const unsigned int maxSieveSize = 32768;
    00348         unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong();
    00349
    00350         m_sieve.clear();
    00351         m_sieve.resize(sieveSize, false);
    00352
    00353         if (m_delta == 0)
    00354         {
    00355                 for (unsigned int i = 0; i < primeTableSize; ++i)
    00356                         SieveSingle(m_sieve, primeTable[i], m_first, m_step, (word16)m_step.InverseMod(primeTable[i]));
    00357         }
    00358         else
    00359         {
    00360                 assert(m_step%2==0);
    00361                 Integer qFirst = (m_first-m_delta) >> 1;
    00362                 Integer halfStep = m_step >> 1;
    00363                 for (unsigned int i = 0; i < primeTableSize; ++i)
    00364                 {
    00365                         word16 p = primeTable[i];
    00366                         word16 stepInv = (word16)m_step.InverseMod(p);
    00367                         SieveSingle(m_sieve, p, m_first, m_step, stepInv);
    00368
    00369                         word16 halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p;
    00370                         SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
    00371                 }
    00372         }
    00373 }
    00374
    00375 bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector)
    00376 {
    00377         assert(!equiv.IsNegative() && equiv < mod);
    00378
    00379         Integer gcd = GCD(equiv, mod);
    00380         if (gcd != Integer::One())
    00381         {
    00382                 // the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv)
    00383                 if (p <= gcd && gcd <= max && IsPrime(gcd) && (!pSelector || pSelector->IsAcceptable(gcd)))
    00384                 {
    00385                         p = gcd;
    00386                         return true;
    00387                 }
    00388                 else
    00389                         return false;
    00390         }
    00391
    00392         unsigned int primeTableSize;
    00393         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00394
    00395         if (p <= primeTable[primeTableSize-1])
    00396         {
    00397                 const word16 *pItr;
    00398
    00399                 --p;
    00400                 if (p.IsPositive())
    00401                         pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
    00402                 else
    00403                         pItr = primeTable;
    00404
    00405                 while (pItr < primeTable+primeTableSize && !(*pItr%mod == equiv && (!pSelector || pSelector->IsAcceptable(*pItr))))
    00406                         ++pItr;
    00407
    00408                 if (pItr < primeTable+primeTableSize)
    00409                 {
    00410                         p = *pItr;
    00411                         return p <= max;
    00412                 }
    00413
    00414                 p = primeTable[primeTableSize-1]+1;
    00415         }
    00416
    00417         assert(p > primeTable[primeTableSize-1]);
    00418
    00419         if (mod.IsOdd())
    00420                 return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector);
    00421
    00422         p += (equiv-p)%mod;
    00423
    00424         if (p>max)
    00425                 return false;
    00426
    00427         PrimeSieve sieve(p, max, mod);
    00428
    00429         while (sieve.NextCandidate(p))
    00430         {
    00431                 if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p))
    00432                         return true;
    00433         }
    00434
    00435         return false;
    00436 }
    00437
    00438 // the following two functions are based on code and comments provided by Preda Mihailescu
    00439 static bool ProvePrime(const Integer &p, const Integer &q)
    00440 {
    00441         assert(p < q*q*q);
    00442         assert(p % q == 1);
    00443
    00444 // this is the Quisquater test. Numbers p having passed the Lucas - Lehmer test
    00445 // for q and verifying p < q^3 can only be built up of two factors, both = 1 mod q,
    00446 // or be prime. The next two lines build the discriminant of a quadratic equation
    00447 // which holds iff p is built up of two factors (excercise ... )
    00448
    00449         Integer r = (p-1)/q;
    00450         if (((r%q).Squared()-4*(r/q)).IsSquare())
    00451                 return false;
    00452
    00453         unsigned int primeTableSize;
    00454         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00455
    00456         assert(primeTableSize >= 50);
    00457         for (int i=0; i<50; i++)
    00458         {
    00459                 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
    00460                 if (b != 1)
    00461                         return a_exp_b_mod_c(b, q, p) == 1;
    00462         }
    00463         return false;
    00464 }
    00465
    00466 Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits)
    00467 {
    00468         Integer p;
    00469         Integer minP = Integer::Power2(pbits-1);
    00470         Integer maxP = Integer::Power2(pbits) - 1;
    00471
    00472         if (maxP <= Integer(s_lastSmallPrime).Squared())
    00473         {
    00474                 // Randomize() will generate a prime provable by trial division
    00475                 p.Randomize(rng, minP, maxP, Integer::PRIME);
    00476                 return p;
    00477         }
    00478
    00479         unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36);
    00480         Integer q = MihailescuProvablePrime(rng, qbits);
    00481         Integer q2 = q<<1;
    00482
    00483         while (true)
    00484         {
    00485                 // this initializes the sieve to search in the arithmetic
    00486                 // progression p = p_0 + \lambda * q2 = p_0 + 2 * \lambda * q,
    00487                 // with q the recursively generated prime above. We will be able
    00488                 // to use Lucas tets for proving primality. A trick of Quisquater
    00489                 // allows taking q > cubic_root(p) rather then square_root: this
    00490                 // decreases the recursion.
    00491
    00492                 p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2);
    00493                 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
    00494
    00495                 while (sieve.NextCandidate(p))
    00496                 {
    00497                         if (FastProbablePrimeTest(p) && ProvePrime(p, q))
    00498                                 return p;
    00499                 }
    00500         }
    00501
    00502         // not reached
    00503         return p;
    00504 }
    00505
    00506 Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
    00507 {
    00508         const unsigned smallPrimeBound = 29, c_opt=10;
    00509         Integer p;
    00510
    00511         unsigned int primeTableSize;
    00512         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00513
    00514         if (bits < smallPrimeBound)
    00515         {
    00516                 do
    00517                         p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2);
    00518                 while (TrialDivision(p, 1 << ((bits+1)/2)));
    00519         }
    00520         else
    00521         {
    00522                 const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
    00523                 double relativeSize;
    00524                 do
    00525                         relativeSize = pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1);
    00526                 while (bits * relativeSize >= bits - margin);
    00527
    00528                 Integer a,b;
    00529                 Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize));
    00530                 Integer I = Integer::Power2(bits-2)/q;
    00531                 Integer I2 = I << 1;
    00532                 unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt);
    00533                 bool success = false;
    00534                 while (!success)
    00535                 {
    00536                         p.Randomize(rng, I, I2, Integer::ANY);
    00537                         p *= q; p <<= 1; ++p;
    00538                         if (!TrialDivision(p, trialDivisorBound))
    00539                         {
    00540                                 a.Randomize(rng, 2, p-1, Integer::ANY);
    00541                                 b = a_exp_b_mod_c(a, (p-1)/q, p);
    00542                                 success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1);
    00543                         }
    00544                 }
    00545         }
    00546         return p;
    00547 }
    00548
    00549 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u)
    00550 {
    00551         // isn't operator overloading great?
    00552         return p * (u * (xq-xp) % q) + xp;
    00553 /*
    00554         Integer t1 = xq-xp;
    00555         cout << hex << t1 << endl;
    00556         Integer t2 = u * t1;
    00557         cout << hex << t2 << endl;
    00558         Integer t3 = t2 % q;
    00559         cout << hex << t3 << endl;
    00560         Integer t4 = p * t3;
    00561         cout << hex << t4 << endl;
    00562         Integer t5 = t4 + xp;
    00563         cout << hex << t5 << endl;
    00564         return t5;
    00565 */
    00566 }
    00567
    00568 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q)
    00569 {
    00570         return CRT(xp, p, xq, q, EuclideanMultiplicativeInverse(p, q));
    00571 }
    00572
    00573 Integer ModularSquareRoot(const Integer &a, const Integer &p)
    00574 {
    00575         if (p%4 == 3)
    00576                 return a_exp_b_mod_c(a, (p+1)/4, p);
    00577
    00578         Integer q=p-1;
    00579         unsigned int r=0;
    00580         while (q.IsEven())
    00581         {
    00582                 r++;
    00583                 q >>= 1;
    00584         }
    00585
    00586         Integer n=2;
    00587         while (Jacobi(n, p) != -1)
    00588                 ++n;
    00589
    00590         Integer y = a_exp_b_mod_c(n, q, p);
    00591         Integer x = a_exp_b_mod_c(a, (q-1)/2, p);
    00592         Integer b = (x.Squared()%p)*a%p;
    00593         x = a*x%p;
    00594         Integer tempb, t;
    00595
    00596         while (b != 1)
    00597         {
    00598                 unsigned m=0;
    00599                 tempb = b;
    00600                 do
    00601                 {
    00602                         m++;
    00603                         b = b.Squared()%p;
    00604                         if (m==r)
    00605                                 return Integer::Zero();
    00606                 }
    00607                 while (b != 1);
    00608
    00609                 t = y;
    00610                 for (unsigned i=0; i<r-m-1; i++)
    00611                         t = t.Squared()%p;
    00612                 y = t.Squared()%p;
    00613                 r = m;
    00614                 x = x*t%p;
    00615                 b = tempb*y%p;
    00616         }
    00617
    00618         assert(x.Squared()%p == a);
    00619         return x;
    00620 }
    00621
    00622 bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p)
    00623 {
    00624         Integer D = (b.Squared() - 4*a*c) % p;
    00625         switch (Jacobi(D, p))
    00626         {
    00627         default:
    00628                 assert(false);  // not reached
    00629                 return false;
    00630         case -1:
    00631                 return false;
    00632         case 0:
    00633                 r1 = r2 = (-b*(a+a).InverseMod(p)) % p;
    00634                 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
    00635                 return true;
    00636         case 1:
    00637                 Integer s = ModularSquareRoot(D, p);
    00638                 Integer t = (a+a).InverseMod(p);
    00639                 r1 = (s-b)*t % p;
    00640                 r2 = (-s-b)*t % p;
    00641                 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
    00642                 assert(((r2.Squared()*a + r2*b + c) % p).IsZero());
    00643                 return true;
    00644         }
    00645 }
    00646
    00647 Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq,
    00648                                         const Integer &p, const Integer &q, const Integer &u)
    00649 {
    00650         Integer p2 = ModularExponentiation((a % p), dp, p);
    00651         Integer q2 = ModularExponentiation((a % q), dq, q);
    00652         return CRT(p2, p, q2, q, u);
    00653 }
    00654
    00655 Integer ModularRoot(const Integer &a, const Integer &e,
    00656                                         const Integer &p, const Integer &q)
    00657 {
    00658         Integer dp = EuclideanMultiplicativeInverse(e, p-1);
    00659         Integer dq = EuclideanMultiplicativeInverse(e, q-1);
    00660         Integer u = EuclideanMultiplicativeInverse(p, q);
    00661         assert(!!dp && !!dq && !!u);
    00662         return ModularRoot(a, dp, dq, p, q, u);
    00663 }
    00664
    00665 /*
    00666 Integer GCDI(const Integer &x, const Integer &y)
    00667 {
    00668         Integer a=x, b=y;
    00669         unsigned k=0;
    00670
    00671         assert(!!a && !!b);
    00672
    00673         while (a[0]==0 && b[0]==0)
    00674         {
    00675                 a >>= 1;
    00676                 b >>= 1;
    00677                 k++;
    00678         }
    00679
    00680         while (a[0]==0)
    00681                 a >>= 1;
    00682
    00683         while (b[0]==0)
    00684                 b >>= 1;
    00685
    00686         while (1)
    00687         {
    00688                 switch (a.Compare(b))
    00689                 {
    00690                         case -1:
    00691                                 b -= a;
    00692                                 while (b[0]==0)
    00693                                         b >>= 1;
    00694                                 break;
    00695
    00696                         case 0:
    00697                                 return (a <<= k);
    00698
    00699                         case 1:
    00700                                 a -= b;
    00701                                 while (a[0]==0)
    00702                                         a >>= 1;
    00703                                 break;
    00704
    00705                         default:
    00706                                 assert(false);
    00707                 }
    00708         }
    00709 }
    00710
    00711 Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b)
    00712 {
    00713         assert(b.Positive());
    00714
    00715         if (a.Negative())
    00716                 return EuclideanMultiplicativeInverse(a%b, b);
    00717
    00718         if (b[0]==0)
    00719         {
    00720                 if (!b || a[0]==0)
    00721                         return Integer::Zero();       // no inverse
    00722                 if (a==1)
    00723                         return 1;
    00724                 Integer u = EuclideanMultiplicativeInverse(b, a);
    00725                 if (!u)
    00726                         return Integer::Zero();       // no inverse
    00727                 else
    00728                         return (b*(a-u)+1)/a;
    00729         }
    00730
    00731         Integer u=1, d=a, v1=b, v3=b, t1, t3, b2=(b+1)>>1;
    00732
    00733         if (a[0])
    00734         {
    00735                 t1 = Integer::Zero();
    00736                 t3 = -b;
    00737         }
    00738         else
    00739         {
    00740                 t1 = b2;
    00741                 t3 = a>>1;
    00742         }
    00743
    00744         while (!!t3)
    00745         {
    00746                 while (t3[0]==0)
    00747                 {
    00748                         t3 >>= 1;
    00749                         if (t1[0]==0)
    00750                                 t1 >>= 1;
    00751                         else
    00752                         {
    00753                                 t1 >>= 1;
    00754                                 t1 += b2;
    00755                         }
    00756                 }
    00757                 if (t3.Positive())
    00758                 {
    00759                         u = t1;
    00760                         d = t3;
    00761                 }
    00762                 else
    00763                 {
    00764                         v1 = b-t1;
    00765                         v3 = -t3;
    00766                 }
    00767                 t1 = u-v1;
    00768                 t3 = d-v3;
    00769                 if (t1.Negative())
    00770                         t1 += b;
    00771         }
    00772         if (d==1)
    00773                 return u;
    00774         else
    00775                 return Integer::Zero();   // no inverse
    00776 }
    00777 */
    00778
    00779 int Jacobi(const Integer &aIn, const Integer &bIn)
    00780 {
    00781         assert(bIn.IsOdd());
    00782
    00783         Integer b = bIn, a = aIn%bIn;
    00784         int result = 1;
    00785
    00786         while (!!a)
    00787         {
    00788                 unsigned i=0;
    00789                 while (a.GetBit(i)==0)
    00790                         i++;
    00791                 a>>=i;
    00792
    00793                 if (i%2==1 && (b%8==3 || b%8==5))
    00794                         result = -result;
    00795
    00796                 if (a%4==3 && b%4==3)
    00797                         result = -result;
    00798
    00799                 std::swap(a, b);
    00800                 a %= b;
    00801         }
    00802
    00803         return (b==1) ? result : 0;
    00804 }
    00805
    00806 Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n)
    00807 {
    00808         unsigned i = e.BitCount();
    00809         if (i==0)
    00810                 return Integer::Two();
    00811
    00812         MontgomeryRepresentation m(n);
    00813         Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two());
    00814         Integer v=p, v1=m.Subtract(m.Square(p), two);
    00815
    00816         i--;
    00817         while (i--)
    00818         {
    00819                 if (e.GetBit(i))
    00820                 {
    00821                         // v = (v*v1 - p) % m;
    00822                         v = m.Subtract(m.Multiply(v,v1), p);
    00823                         // v1 = (v1*v1 - 2) % m;
    00824                         v1 = m.Subtract(m.Square(v1), two);
    00825                 }
    00826                 else
    00827                 {
    00828                         // v1 = (v*v1 - p) % m;
    00829                         v1 = m.Subtract(m.Multiply(v,v1), p);
    00830                         // v = (v*v - 2) % m;
    00831                         v = m.Subtract(m.Square(v), two);
    00832                 }
    00833         }
    00834         return m.ConvertOut(v);
    00835 }
    00836
    00837 // This is Peter Montgomery's unpublished Lucas sequence evalutation algorithm.
    00838 // The total number of multiplies and squares used is less than the binary
    00839 // algorithm (see above).  Unfortunately I can't get it to run as fast as
    00840 // the binary algorithm because of the extra overhead.
    00841 /*
    00842 Integer Lucas(const Integer &n, const Integer &P, const Integer &modulus)
    00843 {
    00844         if (!n)
    00845                 return 2;
    00846
    00847 #define f(A, B, C)      m.Subtract(m.Multiply(A, B), C)
    00848 #define X2(A) m.Subtract(m.Square(A), two)
    00849 #define X3(A) m.Multiply(A, m.Subtract(m.Square(A), three))
    00850
    00851         MontgomeryRepresentation m(modulus);
    00852         Integer two=m.ConvertIn(2), three=m.ConvertIn(3);
    00853         Integer A=m.ConvertIn(P), B, C, p, d=n, e, r, t, T, U;
    00854
    00855         while (d!=1)
    00856         {
    00857                 p = d;
    00858                 unsigned int b = WORD_BITS * p.WordCount();
    00859                 Integer alpha = (Integer(5)<<(2*b-2)).SquareRoot() - Integer::Power2(b-1);
    00860                 r = (p*alpha)>>b;
    00861                 e = d-r;
    00862                 B = A;
    00863                 C = two;
    00864                 d = r;
    00865
    00866                 while (d!=e)
    00867                 {
    00868                         if (d<e)
    00869                         {
    00870                                 swap(d, e);
    00871                                 swap(A, B);
    00872                         }
    00873
    00874                         unsigned int dm2 = d[0], em2 = e[0];
    00875                         unsigned int dm3 = d%3, em3 = e%3;
    00876
    00877 //                      if ((dm6+em6)%3 == 0 && d <= e + (e>>2))
    00878                         if ((dm3+em3==0 || dm3+em3==3) && (t = e, t >>= 2, t += e, d <= t))
    00879                         {
    00880                                 // #1
    00881 //                              t = (d+d-e)/3;
    00882 //                              t = d; t += d; t -= e; t /= 3;
    00883 //                              e = (e+e-d)/3;
    00884 //                              e += e; e -= d; e /= 3;
    00885 //                              d = t;
    00886
    00887 //                              t = (d+e)/3
    00888                                 t = d; t += e; t /= 3;
    00889                                 e -= t;
    00890                                 d -= t;
    00891
    00892                                 T = f(A, B, C);
    00893                                 U = f(T, A, B);
    00894                                 B = f(T, B, A);
    00895                                 A = U;
    00896                                 continue;
    00897                         }
    00898
    00899 //                      if (dm6 == em6 && d <= e + (e>>2))
    00900                         if (dm3 == em3 && dm2 == em2 && (t = e, t >>= 2, t += e, d <= t))
    00901                         {
    00902                                 // #2
    00903 //                              d = (d-e)>>1;
    00904                                 d -= e; d >>= 1;
    00905                                 B = f(A, B, C);
    00906                                 A = X2(A);
    00907                                 continue;
    00908                         }
    00909
    00910 //                      if (d <= (e<<2))
    00911                         if (d <= (t = e, t <<= 2))
    00912                         {
    00913                                 // #3
    00914                                 d -= e;
    00915                                 C = f(A, B, C);
    00916                                 swap(B, C);
    00917                                 continue;
    00918                         }
    00919
    00920                         if (dm2 == em2)
    00921                         {
    00922                                 // #4
    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 (dm2 == 0)
    00931                         {
    00932                                 // #5
    00933                                 d >>= 1;
    00934                                 C = f(A, C, B);
    00935                                 A = X2(A);
    00936                                 continue;
    00937                         }
    00938
    00939                         if (dm3 == 0)
    00940                         {
    00941                                 // #6
    00942 //                              d = d/3 - e;
    00943                                 d /= 3; d -= e;
    00944                                 T = X2(A);
    00945                                 C = f(T, f(A, B, C), C);
    00946                                 swap(B, C);
    00947                                 A = f(T, A, A);
    00948                                 continue;
    00949                         }
    00950
    00951                         if (dm3+em3==0 || dm3+em3==3)
    00952                         {
    00953                                 // #7
    00954 //                              d = (d-e-e)/3;
    00955                                 d -= e; d -= e; d /= 3;
    00956                                 T = f(A, B, C);
    00957                                 B = f(T, A, B);
    00958                                 A = X3(A);
    00959                                 continue;
    00960                         }
    00961
    00962                         if (dm3 == em3)
    00963                         {
    00964                                 // #8
    00965 //                              d = (d-e)/3;
    00966                                 d -= e; d /= 3;
    00967                                 T = f(A, B, C);
    00968                                 C = f(A, C, B);
    00969                                 B = T;
    00970                                 A = X3(A);
    00971                                 continue;
    00972                         }
    00973
    00974                         assert(em2 == 0);
    00975                         // #9
    00976                         e >>= 1;
    00977                         C = f(C, B, A);
    00978                         B = X2(B);
    00979                 }
    00980
    00981                 A = f(A, B, C);
    00982         }
    00983
    00984 #undef f
    00985 #undef X2
    00986 #undef X3
    00987
    00988         return m.ConvertOut(A);
    00989 }
    00990 */
    00991
    00992 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
    00993 {
    00994         Integer d = (m*m-4);
    00995         Integer p2 = p-Jacobi(d,p);
    00996         Integer q2 = q-Jacobi(d,q);
    00997         return CRT(Lucas(EuclideanMultiplicativeInverse(e,p2), m, p), p, Lucas(EuclideanMultiplicativeInverse(e,q2), m, q), q, u);
    00998 }
    00999
    01000 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q)
    01001 {
    01002         return InverseLucas(e, m, p, q, EuclideanMultiplicativeInverse(p, q));
    01003 }
    01004
    01005 unsigned int FactoringWorkFactor(unsigned int n)
    01006 {
    01007         // extrapolated from the table in Odlyzko's "The Future of Integer Factorization"
    01008         // updated to reflect the factoring of RSA-130
    01009         if (n<5) return 0;
    01010         else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
    01011 }
    01012
    01013 unsigned int DiscreteLogWorkFactor(unsigned int n)
    01014 {
    01015         // assuming discrete log takes about the same time as factoring
    01016         if (n<5) return 0;
    01017         else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
    01018 }
    01019
    01020 // ********************************************************
    01021
    01022 void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits)
    01023 {
    01024         // no prime exists for delta = -1, qbits = 4, and pbits = 5
    01025         assert(qbits > 4);
    01026         assert(pbits > qbits);
    01027
    01028         if (qbits+1 == pbits)
    01029         {
    01030                 Integer minP = Integer::Power2(pbits-1);
    01031                 Integer maxP = Integer::Power2(pbits) - 1;
    01032                 bool success = false;
    01033
    01034                 while (!success)
    01035                 {
    01036                         p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
    01037                         PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
    01038
    01039                         while (sieve.NextCandidate(p))
    01040                         {
    01041                                 assert(IsSmallPrime(p) || SmallDivisorsTest(p));
    01042                                 q = (p-delta) >> 1;
    01043                                 assert(IsSmallPrime(q) || SmallDivisorsTest(q));
    01044                                 if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
    01045                                 {
    01046                                         success = true;
    01047                                         break;
    01048                                 }
    01049                         }
    01050                 }
    01051
    01052                 if (delta == 1)
    01053                 {
    01054                         // find g such that g is a quadratic residue mod p, then g has order q
    01055                         // g=4 always works, but this way we get the smallest quadratic residue (other than 1)
    01056                         for (g=2; Jacobi(g, p) != 1; ++g) {}
    01057                         // contributed by Walt Tuvell: g should be the following according to the Law of Quadratic Reciprocity
    01058                         assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
    01059                 }
    01060                 else
    01061                 {
    01062                         assert(delta == -1);
    01063                         // find g such that g*g-4 is a quadratic non-residue,
    01064                         // and such that g has order q
    01065                         for (g=3; ; ++g)
    01066                                 if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
    01067                                         break;
    01068                 }
    01069         }
    01070         else
    01071         {
    01072                 Integer minQ = Integer::Power2(qbits-1);
    01073                 Integer maxQ = Integer::Power2(qbits) - 1;
    01074                 Integer minP = Integer::Power2(pbits-1);
    01075                 Integer maxP = Integer::Power2(pbits) - 1;
    01076
    01077                 do
    01078                 {
    01079                         q.Randomize(rng, minQ, maxQ, Integer::PRIME);
    01080                 } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
    01081
    01082                 // find a random g of order q
    01083                 if (delta==1)
    01084                 {
    01085                         do
    01086                         {
    01087                                 Integer h(rng, 2, p-2, Integer::ANY);
    01088                                 g = a_exp_b_mod_c(h, (p-1)/q, p);
    01089                         } while (g <= 1);
    01090                         assert(a_exp_b_mod_c(g, q, p)==1);
    01091                 }
    01092                 else
    01093                 {
    01094                         assert(delta==-1);
    01095                         do
    01096                         {
    01097                                 Integer h(rng, 3, p-1, Integer::ANY);
    01098                                 if (Jacobi(h*h-4, p)==1)
    01099                                         continue;
    01100                                 g = Lucas((p+1)/q, h, p);
    01101                         } while (g <= 2);
    01102                         assert(Lucas(q, g, p) == 2);
    01103                 }
    01104         }
    01105 }
    01106
    01107 NAMESPACE_END
    01108
    01109 #endif