nbtheory.cpp

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

Generated on Fri Dec 16 03:04:16 2005 for Crypto++ by  doxygen 1.4.5