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
    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 = STDMIN((size_t)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 = 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         m_next = std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin();
    00307         if (m_next == m_sieve.size())
    00308         {
    00309                 m_first += m_sieve.size()*m_step;
    00310                 if (m_first > m_last)
    00311                         return false;
    00312                 else
    00313                 {
    00314                         m_next = 0;
    00315                         DoSieve();
    00316                         return NextCandidate(c);
    00317                 }
    00318         }
    00319         else
    00320         {
    00321                 c = m_first + m_next*m_step;
    00322                 ++m_next;
    00323                 return true;
    00324         }
    00325 }
    00326
    00327 void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv)
    00328 {
    00329         if (stepInv)
    00330         {
    00331                 unsigned int sieveSize = sieve.size();
    00332                 word j = word((word32(p-(first%p))*stepInv) % p);
    00333                 // if the first multiple of p is p, skip it
    00334                 if (first.WordCount() <= 1 && first + step*j == p)
    00335                         j += p;
    00336                 for (; j < sieveSize; j += p)
    00337                         sieve[j] = true;
    00338         }
    00339 }
    00340
    00341 void PrimeSieve::DoSieve()
    00342 {
    00343         unsigned int primeTableSize;
    00344         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00345
    00346         const unsigned int maxSieveSize = 32768;
    00347         unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong();
    00348
    00349         m_sieve.clear();
    00350         m_sieve.resize(sieveSize, false);
    00351
    00352         if (m_delta == 0)
    00353         {
    00354                 for (unsigned int i = 0; i < primeTableSize; ++i)
    00355                         SieveSingle(m_sieve, primeTable[i], m_first, m_step, m_step.InverseMod(primeTable[i]));
    00356         }
    00357         else
    00358         {
    00359                 assert(m_step%2==0);
    00360                 Integer qFirst = (m_first-m_delta) >> 1;
    00361                 Integer halfStep = m_step >> 1;
    00362                 for (unsigned int i = 0; i < primeTableSize; ++i)
    00363                 {
    00364                         word16 p = primeTable[i];
    00365                         word16 stepInv = m_step.InverseMod(p);
    00366                         SieveSingle(m_sieve, p, m_first, m_step, stepInv);
    00367
    00368                         word16 halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p;
    00369                         SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
    00370                 }
    00371         }
    00372 }
    00373
    00374 bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector)
    00375 {
    00376         assert(!equiv.IsNegative() && equiv < mod);
    00377
    00378         Integer gcd = GCD(equiv, mod);
    00379         if (gcd != Integer::One())
    00380         {
    00381                 // the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv)
    00382                 if (p <= gcd && gcd <= max && IsPrime(gcd) && (!pSelector || pSelector->IsAcceptable(gcd)))
    00383                 {
    00384                         p = gcd;
    00385                         return true;
    00386                 }
    00387                 else
    00388                         return false;
    00389         }
    00390
    00391         unsigned int primeTableSize;
    00392         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00393
    00394         if (p <= primeTable[primeTableSize-1])
    00395         {
    00396                 const word16 *pItr;
    00397
    00398                 --p;
    00399                 if (p.IsPositive())
    00400                         pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
    00401                 else
    00402                         pItr = primeTable;
    00403
    00404                 while (pItr < primeTable+primeTableSize && !(*pItr%mod == equiv && (!pSelector || pSelector->IsAcceptable(*pItr))))
    00405                         ++pItr;
    00406
    00407                 if (pItr < primeTable+primeTableSize)
    00408                 {
    00409                         p = *pItr;
    00410                         return p <= max;
    00411                 }
    00412
    00413                 p = primeTable[primeTableSize-1]+1;
    00414         }
    00415
    00416         assert(p > primeTable[primeTableSize-1]);
    00417
    00418         if (mod.IsOdd())
    00419                 return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector);
    00420
    00421         p += (equiv-p)%mod;
    00422
    00423         if (p>max)
    00424                 return false;
    00425
    00426         PrimeSieve sieve(p, max, mod);
    00427
    00428         while (sieve.NextCandidate(p))
    00429         {
    00430                 if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p))
    00431                         return true;
    00432         }
    00433
    00434         return false;
    00435 }
    00436
    00437 // the following two functions are based on code and comments provided by Preda Mihailescu
    00438 static bool ProvePrime(const Integer &p, const Integer &q)
    00439 {
    00440         assert(p < q*q*q);
    00441         assert(p % q == 1);
    00442
    00443 // this is the Quisquater test. Numbers p having passed the Lucas - Lehmer test
    00444 // for q and verifying p < q^3 can only be built up of two factors, both = 1 mod q,
    00445 // or be prime. The next two lines build the discriminant of a quadratic equation
    00446 // which holds iff p is built up of two factors (excercise ... )
    00447
    00448         Integer r = (p-1)/q;
    00449         if (((r%q).Squared()-4*(r/q)).IsSquare())
    00450                 return false;
    00451
    00452         unsigned int primeTableSize;
    00453         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00454
    00455         assert(primeTableSize >= 50);
    00456         for (int i=0; i<50; i++)
    00457         {
    00458                 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
    00459                 if (b != 1)
    00460                         return a_exp_b_mod_c(b, q, p) == 1;
    00461         }
    00462         return false;
    00463 }
    00464
    00465 Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits)
    00466 {
    00467         Integer p;
    00468         Integer minP = Integer::Power2(pbits-1);
    00469         Integer maxP = Integer::Power2(pbits) - 1;
    00470
    00471         if (maxP <= Integer(s_lastSmallPrime).Squared())
    00472         {
    00473                 // Randomize() will generate a prime provable by trial division
    00474                 p.Randomize(rng, minP, maxP, Integer::PRIME);
    00475                 return p;
    00476         }
    00477
    00478         unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36);
    00479         Integer q = MihailescuProvablePrime(rng, qbits);
    00480         Integer q2 = q<<1;
    00481
    00482         while (true)
    00483         {
    00484                 // this initializes the sieve to search in the arithmetic
    00485                 // progression p = p_0 + \lambda * q2 = p_0 + 2 * \lambda * q,
    00486                 // with q the recursively generated prime above. We will be able
    00487                 // to use Lucas tets for proving primality. A trick of Quisquater
    00488                 // allows taking q > cubic_root(p) rather then square_root: this
    00489                 // decreases the recursion.
    00490
    00491                 p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2);
    00492                 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
    00493
    00494                 while (sieve.NextCandidate(p))
    00495                 {
    00496                         if (FastProbablePrimeTest(p) && ProvePrime(p, q))
    00497                                 return p;
    00498                 }
    00499         }
    00500
    00501         // not reached
    00502         return p;
    00503 }
    00504
    00505 Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
    00506 {
    00507         const unsigned smallPrimeBound = 29, c_opt=10;
    00508         Integer p;
    00509
    00510         unsigned int primeTableSize;
    00511         const word16 * primeTable = GetPrimeTable(primeTableSize);
    00512
    00513         if (bits < smallPrimeBound)
    00514         {
    00515                 do
    00516                         p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2);
    00517                 while (TrialDivision(p, 1 << ((bits+1)/2)));
    00518         }
    00519         else
    00520         {
    00521                 const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
    00522                 double relativeSize;
    00523                 do
    00524                         relativeSize = pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1);
    00525                 while (bits * relativeSize >= bits - margin);
    00526
    00527                 Integer a,b;
    00528                 Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize));
    00529                 Integer I = Integer::Power2(bits-2)/q;
    00530                 Integer I2 = I << 1;
    00531                 unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt);
    00532                 bool success = false;
    00533                 while (!success)
    00534                 {
    00535                         p.Randomize(rng, I, I2, Integer::ANY);
    00536                         p *= q; p <<= 1; ++p;
    00537                         if (!TrialDivision(p, trialDivisorBound))
    00538                         {
    00539                                 a.Randomize(rng, 2, p-1, Integer::ANY);
    00540                                 b = a_exp_b_mod_c(a, (p-1)/q, p);
    00541                                 success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1);
    00542                         }
    00543                 }
    00544         }
    00545         return p;
    00546 }
    00547
    00548 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u)
    00549 {
    00550         // isn't operator overloading great?
    00551         return p * (u * (xq-xp) % q) + xp;
    00552 }
    00553
    00554 Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q)
    00555 {
    00556         return CRT(xp, p, xq, q, EuclideanMultiplicativeInverse(p, q));
    00557 }
    00558
    00559 Integer ModularSquareRoot(const Integer &a, const Integer &p)
    00560 {
    00561         if (p%4 == 3)
    00562                 return a_exp_b_mod_c(a, (p+1)/4, p);
    00563
    00564         Integer q=p-1;
    00565         unsigned int r=0;
    00566         while (q.IsEven())
    00567         {
    00568                 r++;
    00569                 q >>= 1;
    00570         }
    00571
    00572         Integer n=2;
    00573         while (Jacobi(n, p) != -1)
    00574                 ++n;
    00575
    00576         Integer y = a_exp_b_mod_c(n, q, p);
    00577         Integer x = a_exp_b_mod_c(a, (q-1)/2, p);
    00578         Integer b = (x.Squared()%p)*a%p;
    00579         x = a*x%p;
    00580         Integer tempb, t;
    00581
    00582         while (b != 1)
    00583         {
    00584                 unsigned m=0;
    00585                 tempb = b;
    00586                 do
    00587                 {
    00588                         m++;
    00589                         b = b.Squared()%p;
    00590                         if (m==r)
    00591                                 return Integer::Zero();
    00592                 }
    00593                 while (b != 1);
    00594
    00595                 t = y;
    00596                 for (unsigned i=0; i<r-m-1; i++)
    00597                         t = t.Squared()%p;
    00598                 y = t.Squared()%p;
    00599                 r = m;
    00600                 x = x*t%p;
    00601                 b = tempb*y%p;
    00602         }
    00603
    00604         assert(x.Squared()%p == a);
    00605         return x;
    00606 }
    00607
    00608 bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p)
    00609 {
    00610         Integer D = (b.Squared() - 4*a*c) % p;
    00611         switch (Jacobi(D, p))
    00612         {
    00613         default:
    00614                 assert(false);  // not reached
    00615                 return false;
    00616         case -1:
    00617                 return false;
    00618         case 0:
    00619                 r1 = r2 = (-b*(a+a).InverseMod(p)) % p;
    00620                 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
    00621                 return true;
    00622         case 1:
    00623                 Integer s = ModularSquareRoot(D, p);
    00624                 Integer t = (a+a).InverseMod(p);
    00625                 r1 = (s-b)*t % p;
    00626                 r2 = (-s-b)*t % p;
    00627                 assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
    00628                 assert(((r2.Squared()*a + r2*b + c) % p).IsZero());
    00629                 return true;
    00630         }
    00631 }
    00632
    00633 Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq,
    00634                                         const Integer &p, const Integer &q, const Integer &u)
    00635 {
    00636         Integer p2 = ModularExponentiation((a % p), dp, p);
    00637         Integer q2 = ModularExponentiation((a % q), dq, q);
    00638         return CRT(p2, p, q2, q, u);
    00639 }
    00640
    00641 Integer ModularRoot(const Integer &a, const Integer &e,
    00642                                         const Integer &p, const Integer &q)
    00643 {
    00644         Integer dp = EuclideanMultiplicativeInverse(e, p-1);
    00645         Integer dq = EuclideanMultiplicativeInverse(e, q-1);
    00646         Integer u = EuclideanMultiplicativeInverse(p, q);
    00647         assert(!!dp && !!dq && !!u);
    00648         return ModularRoot(a, dp, dq, p, q, u);
    00649 }
    00650
    00651 /*
    00652 Integer GCDI(const Integer &x, const Integer &y)
    00653 {
    00654         Integer a=x, b=y;
    00655         unsigned k=0;
    00656
    00657         assert(!!a && !!b);
    00658
    00659         while (a[0]==0 && b[0]==0)
    00660         {
    00661                 a >>= 1;
    00662                 b >>= 1;
    00663                 k++;
    00664         }
    00665
    00666         while (a[0]==0)
    00667                 a >>= 1;
    00668
    00669         while (b[0]==0)
    00670                 b >>= 1;
    00671
    00672         while (1)
    00673         {
    00674                 switch (a.Compare(b))
    00675                 {
    00676                         case -1:
    00677                                 b -= a;
    00678                                 while (b[0]==0)
    00679                                         b >>= 1;
    00680                                 break;
    00681
    00682                         case 0:
    00683                                 return (a <<= k);
    00684
    00685                         case 1:
    00686                                 a -= b;
    00687                                 while (a[0]==0)
    00688                                         a >>= 1;
    00689                                 break;
    00690
    00691                         default:
    00692                                 assert(false);
    00693                 }
    00694         }
    00695 }
    00696
    00697 Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b)
    00698 {
    00699         assert(b.Positive());
    00700
    00701         if (a.Negative())
    00702                 return EuclideanMultiplicativeInverse(a%b, b);
    00703
    00704         if (b[0]==0)
    00705         {
    00706                 if (!b || a[0]==0)
    00707                         return Integer::Zero();       // no inverse
    00708                 if (a==1)
    00709                         return 1;
    00710                 Integer u = EuclideanMultiplicativeInverse(b, a);
    00711                 if (!u)
    00712                         return Integer::Zero();       // no inverse
    00713                 else
    00714                         return (b*(a-u)+1)/a;
    00715         }
    00716
    00717         Integer u=1, d=a, v1=b, v3=b, t1, t3, b2=(b+1)>>1;
    00718
    00719         if (a[0])
    00720         {
    00721                 t1 = Integer::Zero();
    00722                 t3 = -b;
    00723         }
    00724         else
    00725         {
    00726                 t1 = b2;
    00727                 t3 = a>>1;
    00728         }
    00729
    00730         while (!!t3)
    00731         {
    00732                 while (t3[0]==0)
    00733                 {
    00734                         t3 >>= 1;
    00735                         if (t1[0]==0)
    00736                                 t1 >>= 1;
    00737                         else
    00738                         {
    00739                                 t1 >>= 1;
    00740                                 t1 += b2;
    00741                         }
    00742                 }
    00743                 if (t3.Positive())
    00744                 {
    00745                         u = t1;
    00746                         d = t3;
    00747                 }
    00748                 else
    00749                 {
    00750                         v1 = b-t1;
    00751                         v3 = -t3;
    00752                 }
    00753                 t1 = u-v1;
    00754                 t3 = d-v3;
    00755                 if (t1.Negative())
    00756                         t1 += b;
    00757         }
    00758         if (d==1)
    00759                 return u;
    00760         else
    00761                 return Integer::Zero();   // no inverse
    00762 }
    00763 */
    00764
    00765 int Jacobi(const Integer &aIn, const Integer &bIn)
    00766 {
    00767         assert(bIn.IsOdd());
    00768
    00769         Integer b = bIn, a = aIn%bIn;
    00770         int result = 1;
    00771
    00772         while (!!a)
    00773         {
    00774                 unsigned i=0;
    00775                 while (a.GetBit(i)==0)
    00776                         i++;
    00777                 a>>=i;
    00778
    00779                 if (i%2==1 && (b%8==3 || b%8==5))
    00780                         result = -result;
    00781
    00782                 if (a%4==3 && b%4==3)
    00783                         result = -result;
    00784
    00785                 std::swap(a, b);
    00786                 a %= b;
    00787         }
    00788
    00789         return (b==1) ? result : 0;
    00790 }
    00791
    00792 Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n)
    00793 {
    00794         unsigned i = e.BitCount();
    00795         if (i==0)
    00796                 return Integer::Two();
    00797
    00798         MontgomeryRepresentation m(n);
    00799         Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two());
    00800         Integer v=p, v1=m.Subtract(m.Square(p), two);
    00801
    00802         i--;
    00803         while (i--)
    00804         {
    00805                 if (e.GetBit(i))
    00806                 {
    00807                         // v = (v*v1 - p) % m;
    00808                         v = m.Subtract(m.Multiply(v,v1), p);
    00809                         // v1 = (v1*v1 - 2) % m;
    00810                         v1 = m.Subtract(m.Square(v1), two);
    00811                 }
    00812                 else
    00813                 {
    00814                         // v1 = (v*v1 - p) % m;
    00815                         v1 = m.Subtract(m.Multiply(v,v1), p);
    00816                         // v = (v*v - 2) % m;
    00817                         v = m.Subtract(m.Square(v), two);
    00818                 }
    00819         }
    00820         return m.ConvertOut(v);
    00821 }
    00822
    00823 // This is Peter Montgomery's unpublished Lucas sequence evalutation algorithm.
    00824 // The total number of multiplies and squares used is less than the binary
    00825 // algorithm (see above).  Unfortunately I can't get it to run as fast as
    00826 // the binary algorithm because of the extra overhead.
    00827 /*
    00828 Integer Lucas(const Integer &n, const Integer &P, const Integer &modulus)
    00829 {
    00830         if (!n)
    00831                 return 2;
    00832
    00833 #define f(A, B, C)      m.Subtract(m.Multiply(A, B), C)
    00834 #define X2(A) m.Subtract(m.Square(A), two)
    00835 #define X3(A) m.Multiply(A, m.Subtract(m.Square(A), three))
    00836
    00837         MontgomeryRepresentation m(modulus);
    00838         Integer two=m.ConvertIn(2), three=m.ConvertIn(3);
    00839         Integer A=m.ConvertIn(P), B, C, p, d=n, e, r, t, T, U;
    00840
    00841         while (d!=1)
    00842         {
    00843                 p = d;
    00844                 unsigned int b = WORD_BITS * p.WordCount();
    00845                 Integer alpha = (Integer(5)<<(2*b-2)).SquareRoot() - Integer::Power2(b-1);
    00846                 r = (p*alpha)>>b;
    00847                 e = d-r;
    00848                 B = A;
    00849                 C = two;
    00850                 d = r;
    00851
    00852                 while (d!=e)
    00853                 {
    00854                         if (d<e)
    00855                         {
    00856                                 swap(d, e);
    00857                                 swap(A, B);
    00858                         }
    00859
    00860                         unsigned int dm2 = d[0], em2 = e[0];
    00861                         unsigned int dm3 = d%3, em3 = e%3;
    00862
    00863 //                      if ((dm6+em6)%3 == 0 && d <= e + (e>>2))
    00864                         if ((dm3+em3==0 || dm3+em3==3) && (t = e, t >>= 2, t += e, d <= t))
    00865                         {
    00866                                 // #1
    00867 //                              t = (d+d-e)/3;
    00868 //                              t = d; t += d; t -= e; t /= 3;
    00869 //                              e = (e+e-d)/3;
    00870 //                              e += e; e -= d; e /= 3;
    00871 //                              d = t;
    00872
    00873 //                              t = (d+e)/3
    00874                                 t = d; t += e; t /= 3;
    00875                                 e -= t;
    00876                                 d -= t;
    00877
    00878                                 T = f(A, B, C);
    00879                                 U = f(T, A, B);
    00880                                 B = f(T, B, A);
    00881                                 A = U;
    00882                                 continue;
    00883                         }
    00884
    00885 //                      if (dm6 == em6 && d <= e + (e>>2))
    00886                         if (dm3 == em3 && dm2 == em2 && (t = e, t >>= 2, t += e, d <= t))
    00887                         {
    00888                                 // #2
    00889 //                              d = (d-e)>>1;
    00890                                 d -= e; d >>= 1;
    00891                                 B = f(A, B, C);
    00892                                 A = X2(A);
    00893                                 continue;
    00894                         }
    00895
    00896 //                      if (d <= (e<<2))
    00897                         if (d <= (t = e, t <<= 2))
    00898                         {
    00899                                 // #3
    00900                                 d -= e;
    00901                                 C = f(A, B, C);
    00902                                 swap(B, C);
    00903                                 continue;
    00904                         }
    00905
    00906                         if (dm2 == em2)
    00907                         {
    00908                                 // #4
    00909 //                              d = (d-e)>>1;
    00910                                 d -= e; d >>= 1;
    00911                                 B = f(A, B, C);
    00912                                 A = X2(A);
    00913                                 continue;
    00914                         }
    00915
    00916                         if (dm2 == 0)
    00917                         {
    00918                                 // #5
    00919                                 d >>= 1;
    00920                                 C = f(A, C, B);
    00921                                 A = X2(A);
    00922                                 continue;
    00923                         }
    00924
    00925                         if (dm3 == 0)
    00926                         {
    00927                                 // #6
    00928 //                              d = d/3 - e;
    00929                                 d /= 3; d -= e;
    00930                                 T = X2(A);
    00931                                 C = f(T, f(A, B, C), C);
    00932                                 swap(B, C);
    00933                                 A = f(T, A, A);
    00934                                 continue;
    00935                         }
    00936
    00937                         if (dm3+em3==0 || dm3+em3==3)
    00938                         {
    00939                                 // #7
    00940 //                              d = (d-e-e)/3;
    00941                                 d -= e; d -= e; d /= 3;
    00942                                 T = f(A, B, C);
    00943                                 B = f(T, A, B);
    00944                                 A = X3(A);
    00945                                 continue;
    00946                         }
    00947
    00948                         if (dm3 == em3)
    00949                         {
    00950                                 // #8
    00951 //                              d = (d-e)/3;
    00952                                 d -= e; d /= 3;
    00953                                 T = f(A, B, C);
    00954                                 C = f(A, C, B);
    00955                                 B = T;
    00956                                 A = X3(A);
    00957                                 continue;
    00958                         }
    00959
    00960                         assert(em2 == 0);
    00961                         // #9
    00962                         e >>= 1;
    00963                         C = f(C, B, A);
    00964                         B = X2(B);
    00965                 }
    00966
    00967                 A = f(A, B, C);
    00968         }
    00969
    00970 #undef f
    00971 #undef X2
    00972 #undef X3
    00973
    00974         return m.ConvertOut(A);
    00975 }
    00976 */
    00977
    00978 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
    00979 {
    00980         Integer d = (m*m-4);
    00981         Integer p2 = p-Jacobi(d,p);
    00982         Integer q2 = q-Jacobi(d,q);
    00983         return CRT(Lucas(EuclideanMultiplicativeInverse(e,p2), m, p), p, Lucas(EuclideanMultiplicativeInverse(e,q2), m, q), q, u);
    00984 }
    00985
    00986 Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q)
    00987 {
    00988         return InverseLucas(e, m, p, q, EuclideanMultiplicativeInverse(p, q));
    00989 }
    00990
    00991 unsigned int FactoringWorkFactor(unsigned int n)
    00992 {
    00993         // extrapolated from the table in Odlyzko's "The Future of Integer Factorization"
    00994         // updated to reflect the factoring of RSA-130
    00995         if (n<5) return 0;
    00996         else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
    00997 }
    00998
    00999 unsigned int DiscreteLogWorkFactor(unsigned int n)
    01000 {
    01001         // assuming discrete log takes about the same time as factoring
    01002         if (n<5) return 0;
    01003         else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
    01004 }
    01005
    01006 // ********************************************************
    01007
    01008 void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits)
    01009 {
    01010         // no prime exists for delta = -1, qbits = 4, and pbits = 5
    01011         assert(qbits > 4);
    01012         assert(pbits > qbits);
    01013
    01014         if (qbits+1 == pbits)
    01015         {
    01016                 Integer minP = Integer::Power2(pbits-1);
    01017                 Integer maxP = Integer::Power2(pbits) - 1;
    01018                 bool success = false;
    01019
    01020                 while (!success)
    01021                 {
    01022                         p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
    01023                         PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
    01024
    01025                         while (sieve.NextCandidate(p))
    01026                         {
    01027                                 assert(IsSmallPrime(p) || SmallDivisorsTest(p));
    01028                                 q = (p-delta) >> 1;
    01029                                 assert(IsSmallPrime(q) || SmallDivisorsTest(q));
    01030                                 if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
    01031                                 {
    01032                                         success = true;
    01033                                         break;
    01034                                 }
    01035                         }
    01036                 }
    01037
    01038                 if (delta == 1)
    01039                 {
    01040                         // find g such that g is a quadratic residue mod p, then g has order q
    01041                         // g=4 always works, but this way we get the smallest quadratic residue (other than 1)
    01042                         for (g=2; Jacobi(g, p) != 1; ++g) {}
    01043                         // contributed by Walt Tuvell: g should be the following according to the Law of Quadratic Reciprocity
    01044                         assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
    01045                 }
    01046                 else
    01047                 {
    01048                         assert(delta == -1);
    01049                         // find g such that g*g-4 is a quadratic non-residue,
    01050                         // and such that g has order q
    01051                         for (g=3; ; ++g)
    01052                                 if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
    01053                                         break;
    01054                 }
    01055         }
    01056         else
    01057         {
    01058                 Integer minQ = Integer::Power2(qbits-1);
    01059                 Integer maxQ = Integer::Power2(qbits) - 1;
    01060                 Integer minP = Integer::Power2(pbits-1);
    01061                 Integer maxP = Integer::Power2(pbits) - 1;
    01062
    01063                 do
    01064                 {
    01065                         q.Randomize(rng, minQ, maxQ, Integer::PRIME);
    01066                 } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
    01067
    01068                 // find a random g of order q
    01069                 if (delta==1)
    01070                 {
    01071                         do
    01072                         {
    01073                                 Integer h(rng, 2, p-2, Integer::ANY);
    01074                                 g = a_exp_b_mod_c(h, (p-1)/q, p);
    01075                         } while (g <= 1);
    01076                         assert(a_exp_b_mod_c(g, q, p)==1);
    01077                 }
    01078                 else
    01079                 {
    01080                         assert(delta==-1);
    01081                         do
    01082                         {
    01083                                 Integer h(rng, 3, p-1, Integer::ANY);
    01084                                 if (Jacobi(h*h-4, p)==1)
    01085                                         continue;
    01086                                 g = Lucas((p+1)/q, h, p);
    01087                         } while (g <= 2);
    01088                         assert(Lucas(q, g, p) == 2);
    01089                 }
    01090         }
    01091 }
    01092
    01093 NAMESPACE_END
    01094
    01095 #endif