Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members

nbtheory.cpp

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

Generated on Fri Aug 27 13:36:34 2004 for Crypto++ by doxygen 1.3.8