- Main Page
- Namespaces
- Classes
- Files

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