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