Crypto++
|
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