1// nbtheory.cpp - written and placed in the public domain by Wei Dai 2 3#include "pch.h" 4 5#ifndef CRYPTOPP_IMPORTS 6 7#include "nbtheory.h" 8#include "modarith.h" 9#include "algparam.h" 10 11#include <math.h> 12#include <vector> 13 14#ifdef _OPENMP 15// needed in MSVC 2005 to generate correct manifest 16#include <omp.h> 17#endif 18 19NAMESPACE_BEGIN(CryptoPP) 20 21const word s_lastSmallPrime = 32719; 22 23struct NewPrimeTable 24{ 25 std::vector<word16> * operator()() const 26 { 27 const unsigned int maxPrimeTableSize = 3511; 28 29 std::auto_ptr<std::vector<word16> > pPrimeTable(new std::vector<word16>); 30 std::vector<word16> &primeTable = *pPrimeTable; 31 primeTable.reserve(maxPrimeTableSize); 32 33 primeTable.push_back(2); 34 unsigned int testEntriesEnd = 1; 35 36 for (unsigned int p=3; p<=s_lastSmallPrime; p+=2) 37 { 38 unsigned int j; 39 for (j=1; j<testEntriesEnd; j++) 40 if (p%primeTable[j] == 0) 41 break; 42 if (j == testEntriesEnd) 43 { 44 primeTable.push_back(p); 45 testEntriesEnd = UnsignedMin(54U, primeTable.size()); 46 } 47 } 48 49 return pPrimeTable.release(); 50 } 51}; 52 53const word16 * GetPrimeTable(unsigned int &size) 54{ 55 const std::vector<word16> &primeTable = Singleton<std::vector<word16>, NewPrimeTable>().Ref(); 56 size = (unsigned int)primeTable.size(); 57 return &primeTable[0]; 58} 59 60bool IsSmallPrime(const Integer &p) 61{ 62 unsigned int primeTableSize; 63 const word16 * primeTable = GetPrimeTable(primeTableSize); 64 65 if (p.IsPositive() && p <= primeTable[primeTableSize-1]) 66 return std::binary_search(primeTable, primeTable+primeTableSize, (word16)p.ConvertToLong()); 67 else 68 return false; 69} 70 71bool TrialDivision(const Integer &p, unsigned bound) 72{ 73 unsigned int primeTableSize; 74 const word16 * primeTable = GetPrimeTable(primeTableSize); 75 76 assert(primeTable[primeTableSize-1] >= bound); 77 78 unsigned int i; 79 for (i = 0; primeTable[i]<bound; i++) 80 if ((p % primeTable[i]) == 0) 81 return true; 82 83 if (bound == primeTable[i]) 84 return (p % bound == 0); 85 else 86 return false; 87} 88 89bool SmallDivisorsTest(const Integer &p) 90{ 91 unsigned int primeTableSize; 92 const word16 * primeTable = GetPrimeTable(primeTableSize); 93 return !TrialDivision(p, primeTable[primeTableSize-1]); 94} 95 96bool IsFermatProbablePrime(const Integer &n, const Integer &b) 97{ 98 if (n <= 3) 99 return n==2 || n==3; 100 101 assert(n>3 && b>1 && b<n-1); 102 return a_exp_b_mod_c(b, n-1, n)==1; 103} 104 105bool IsStrongProbablePrime(const Integer &n, const Integer &b) 106{ 107 if (n <= 3) 108 return n==2 || n==3; 109 110 assert(n>3 && b>1 && b<n-1); 111 112 if ((n.IsEven() && n!=2) || GCD(b, n) != 1) 113 return false; 114 115 Integer nminus1 = (n-1); 116 unsigned int a; 117 118 // calculate a = largest power of 2 that divides (n-1) 119 for (a=0; ; a++) 120 if (nminus1.GetBit(a)) 121 break; 122 Integer m = nminus1>>a; 123 124 Integer z = a_exp_b_mod_c(b, m, n); 125 if (z==1 || z==nminus1) 126 return true; 127 for (unsigned j=1; j<a; j++) 128 { 129 z = z.Squared()%n; 130 if (z==nminus1) 131 return true; 132 if (z==1) 133 return false; 134 } 135 return false; 136} 137 138bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds) 139{ 140 if (n <= 3) 141 return n==2 || n==3; 142 143 assert(n>3); 144 145 Integer b; 146 for (unsigned int i=0; i<rounds; i++) 147 { 148 b.Randomize(rng, 2, n-2); 149 if (!IsStrongProbablePrime(n, b)) 150 return false; 151 } 152 return true; 153} 154 155bool IsLucasProbablePrime(const Integer &n) 156{ 157 if (n <= 1) 158 return false; 159 160 if (n.IsEven()) 161 return n==2; 162 163 assert(n>2); 164 165 Integer b=3; 166 unsigned int i=0; 167 int j; 168 169 while ((j=Jacobi(b.Squared()-4, n)) == 1) 170 { 171 if (++i==64 && n.IsSquare()) // avoid infinite loop if n is a square 172 return false; 173 ++b; ++b; 174 } 175 176 if (j==0) 177 return false; 178 else 179 return Lucas(n+1, b, n)==2; 180} 181 182bool IsStrongLucasProbablePrime(const Integer &n) 183{ 184 if (n <= 1) 185 return false; 186 187 if (n.IsEven()) 188 return n==2; 189 190 assert(n>2); 191 192 Integer b=3; 193 unsigned int i=0; 194 int j; 195 196 while ((j=Jacobi(b.Squared()-4, n)) == 1) 197 { 198 if (++i==64 && n.IsSquare()) // avoid infinite loop if n is a square 199 return false; 200 ++b; ++b; 201 } 202 203 if (j==0) 204 return false; 205 206 Integer n1 = n+1; 207 unsigned int a; 208 209 // calculate a = largest power of 2 that divides n1 210 for (a=0; ; a++) 211 if (n1.GetBit(a)) 212 break; 213 Integer m = n1>>a; 214 215 Integer z = Lucas(m, b, n); 216 if (z==2 || z==n-2) 217 return true; 218 for (i=1; i<a; i++) 219 { 220 z = (z.Squared()-2)%n; 221 if (z==n-2) 222 return true; 223 if (z==2) 224 return false; 225 } 226 return false; 227} 228 229struct NewLastSmallPrimeSquared 230{ 231 Integer * operator()() const 232 { 233 return new Integer(Integer(s_lastSmallPrime).Squared()); 234 } 235}; 236 237bool IsPrime(const Integer &p) 238{ 239 if (p <= s_lastSmallPrime) 240 return IsSmallPrime(p); 241 else if (p <= Singleton<Integer, NewLastSmallPrimeSquared>().Ref()) 242 return SmallDivisorsTest(p); 243 else 244 return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p); 245} 246 247bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level) 248{ 249 bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1); 250 if (level >= 1) 251 pass = pass && RabinMillerTest(rng, p, 10); 252 return pass; 253} 254 255unsigned int PrimeSearchInterval(const Integer &max) 256{ 257 return max.BitCount(); 258} 259 260static inline bool FastProbablePrimeTest(const Integer &n) 261{ 262 return IsStrongProbablePrime(n,2); 263} 264 265AlgorithmParameters MakeParametersForTwoPrimesOfEqualSize(unsigned int productBitLength) 266{ 267 if (productBitLength < 16) 268 throw InvalidArgument("invalid bit length"); 269 270 Integer minP, maxP; 271 272 if (productBitLength%2==0) 273 { 274 minP = Integer(182) << (productBitLength/2-8); 275 maxP = Integer::Power2(productBitLength/2)-1; 276 } 277 else 278 { 279 minP = Integer::Power2((productBitLength-1)/2); 280 maxP = Integer(181) << ((productBitLength+1)/2-8); 281 } 282 283 return MakeParameters("RandomNumberType", Integer::PRIME)("Min", minP)("Max", maxP); 284} 285 286class PrimeSieve 287{ 288public: 289 // delta == 1 or -1 means double sieve with p = 2*q + delta 290 PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta=0); 291 bool NextCandidate(Integer &c); 292 293 void DoSieve(); 294 static void SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv); 295 296 Integer m_first, m_last, m_step; 297 signed int m_delta; 298 word m_next; 299 std::vector<bool> m_sieve; 300}; 301 302PrimeSieve::PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta) 303 : m_first(first), m_last(last), m_step(step), m_delta(delta), m_next(0) 304{ 305 DoSieve(); 306} 307 308bool PrimeSieve::NextCandidate(Integer &c) 309{ 310 bool safe = SafeConvert(std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin(), m_next); 311 assert(safe); 312 if (m_next == m_sieve.size()) 313 { 314 m_first += long(m_sieve.size())*m_step; 315 if (m_first > m_last) 316 return false; 317 else 318 { 319 m_next = 0; 320 DoSieve(); 321 return NextCandidate(c); 322 } 323 } 324 else 325 { 326 c = m_first + long(m_next)*m_step; 327 ++m_next; 328 return true; 329 } 330} 331 332void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv) 333{ 334 if (stepInv) 335 { 336 size_t sieveSize = sieve.size(); 337 size_t j = (word32(p-(first%p))*stepInv) % p; 338 // if the first multiple of p is p, skip it 339 if (first.WordCount() <= 1 && first + step*long(j) == p) 340 j += p; 341 for (; j < sieveSize; j += p) 342 sieve[j] = true; 343 } 344} 345 346void PrimeSieve::DoSieve() 347{ 348 unsigned int primeTableSize; 349 const word16 * primeTable = GetPrimeTable(primeTableSize); 350 351 const unsigned int maxSieveSize = 32768; 352 unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong(); 353 354 m_sieve.clear(); 355 m_sieve.resize(sieveSize, false); 356 357 if (m_delta == 0) 358 { 359 for (unsigned int i = 0; i < primeTableSize; ++i) 360 SieveSingle(m_sieve, primeTable[i], m_first, m_step, (word16)m_step.InverseMod(primeTable[i])); 361 } 362 else 363 { 364 assert(m_step%2==0); 365 Integer qFirst = (m_first-m_delta) >> 1; 366 Integer halfStep = m_step >> 1; 367 for (unsigned int i = 0; i < primeTableSize; ++i) 368 { 369 word16 p = primeTable[i]; 370 word16 stepInv = (word16)m_step.InverseMod(p); 371 SieveSingle(m_sieve, p, m_first, m_step, stepInv); 372 373 word16 halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p; 374 SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv); 375 } 376 } 377} 378 379bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector) 380{ 381 assert(!equiv.IsNegative() && equiv < mod); 382 383 Integer gcd = GCD(equiv, mod); 384 if (gcd != Integer::One()) 385 { 386 // the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv) 387 if (p <= gcd && gcd <= max && IsPrime(gcd) && (!pSelector || pSelector->IsAcceptable(gcd))) 388 { 389 p = gcd; 390 return true; 391 } 392 else 393 return false; 394 } 395 396 unsigned int primeTableSize; 397 const word16 * primeTable = GetPrimeTable(primeTableSize); 398 399 if (p <= primeTable[primeTableSize-1]) 400 { 401 const word16 *pItr; 402 403 --p; 404 if (p.IsPositive()) 405 pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong()); 406 else 407 pItr = primeTable; 408 409 while (pItr < primeTable+primeTableSize && !(*pItr%mod == equiv && (!pSelector || pSelector->IsAcceptable(*pItr)))) 410 ++pItr; 411 412 if (pItr < primeTable+primeTableSize) 413 { 414 p = *pItr; 415 return p <= max; 416 } 417 418 p = primeTable[primeTableSize-1]+1; 419 } 420 421 assert(p > primeTable[primeTableSize-1]); 422 423 if (mod.IsOdd()) 424 return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector); 425 426 p += (equiv-p)%mod; 427 428 if (p>max) 429 return false; 430 431 PrimeSieve sieve(p, max, mod); 432 433 while (sieve.NextCandidate(p)) 434 { 435 if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p)) 436 return true; 437 } 438 439 return false; 440} 441 442// the following two functions are based on code and comments provided by Preda Mihailescu 443static bool ProvePrime(const Integer &p, const Integer &q) 444{ 445 assert(p < q*q*q); 446 assert(p % q == 1); 447 448// this is the Quisquater test. Numbers p having passed the Lucas - Lehmer test 449// for q and verifying p < q^3 can only be built up of two factors, both = 1 mod q, 450// or be prime. The next two lines build the discriminant of a quadratic equation 451// which holds iff p is built up of two factors (excercise ... ) 452 453 Integer r = (p-1)/q; 454 if (((r%q).Squared()-4*(r/q)).IsSquare()) 455 return false; 456 457 unsigned int primeTableSize; 458 const word16 * primeTable = GetPrimeTable(primeTableSize); 459 460 assert(primeTableSize >= 50); 461 for (int i=0; i<50; i++) 462 { 463 Integer b = a_exp_b_mod_c(primeTable[i], r, p); 464 if (b != 1) 465 return a_exp_b_mod_c(b, q, p) == 1; 466 } 467 return false; 468} 469 470Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits) 471{ 472 Integer p; 473 Integer minP = Integer::Power2(pbits-1); 474 Integer maxP = Integer::Power2(pbits) - 1; 475 476 if (maxP <= Integer(s_lastSmallPrime).Squared()) 477 { 478 // Randomize() will generate a prime provable by trial division 479 p.Randomize(rng, minP, maxP, Integer::PRIME); 480 return p; 481 } 482 483 unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36); 484 Integer q = MihailescuProvablePrime(rng, qbits); 485 Integer q2 = q<<1; 486 487 while (true) 488 { 489 // this initializes the sieve to search in the arithmetic 490 // progression p = p_0 + \lambda * q2 = p_0 + 2 * \lambda * q, 491 // with q the recursively generated prime above. We will be able 492 // to use Lucas tets for proving primality. A trick of Quisquater 493 // allows taking q > cubic_root(p) rather then square_root: this 494 // decreases the recursion. 495 496 p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2); 497 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2); 498 499 while (sieve.NextCandidate(p)) 500 { 501 if (FastProbablePrimeTest(p) && ProvePrime(p, q)) 502 return p; 503 } 504 } 505 506 // not reached 507 return p; 508} 509 510Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits) 511{ 512 const unsigned smallPrimeBound = 29, c_opt=10; 513 Integer p; 514 515 unsigned int primeTableSize; 516 const word16 * primeTable = GetPrimeTable(primeTableSize); 517 518 if (bits < smallPrimeBound) 519 { 520 do 521 p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2); 522 while (TrialDivision(p, 1 << ((bits+1)/2))); 523 } 524 else 525 { 526 const unsigned margin = bits > 50 ? 20 : (bits-10)/2; 527 double relativeSize; 528 do 529 relativeSize = pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1); 530 while (bits * relativeSize >= bits - margin); 531 532 Integer a,b; 533 Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize)); 534 Integer I = Integer::Power2(bits-2)/q; 535 Integer I2 = I << 1; 536 unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt); 537 bool success = false; 538 while (!success) 539 { 540 p.Randomize(rng, I, I2, Integer::ANY); 541 p *= q; p <<= 1; ++p; 542 if (!TrialDivision(p, trialDivisorBound)) 543 { 544 a.Randomize(rng, 2, p-1, Integer::ANY); 545 b = a_exp_b_mod_c(a, (p-1)/q, p); 546 success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1); 547 } 548 } 549 } 550 return p; 551} 552 553Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u) 554{ 555 // isn't operator overloading great? 556 return p * (u * (xq-xp) % q) + xp; 557/* 558 Integer t1 = xq-xp; 559 cout << hex << t1 << endl; 560 Integer t2 = u * t1; 561 cout << hex << t2 << endl; 562 Integer t3 = t2 % q; 563 cout << hex << t3 << endl; 564 Integer t4 = p * t3; 565 cout << hex << t4 << endl; 566 Integer t5 = t4 + xp; 567 cout << hex << t5 << endl; 568 return t5; 569*/ 570} 571 572Integer ModularSquareRoot(const Integer &a, const Integer &p) 573{ 574 if (p%4 == 3) 575 return a_exp_b_mod_c(a, (p+1)/4, p); 576 577 Integer q=p-1; 578 unsigned int r=0; 579 while (q.IsEven()) 580 { 581 r++; 582 q >>= 1; 583 } 584 585 Integer n=2; 586 while (Jacobi(n, p) != -1) 587 ++n; 588 589 Integer y = a_exp_b_mod_c(n, q, p); 590 Integer x = a_exp_b_mod_c(a, (q-1)/2, p); 591 Integer b = (x.Squared()%p)*a%p; 592 x = a*x%p; 593 Integer tempb, t; 594 595 while (b != 1) 596 { 597 unsigned m=0; 598 tempb = b; 599 do 600 { 601 m++; 602 b = b.Squared()%p; 603 if (m==r) 604 return Integer::Zero(); 605 } 606 while (b != 1); 607 608 t = y; 609 for (unsigned i=0; i<r-m-1; i++) 610 t = t.Squared()%p; 611 y = t.Squared()%p; 612 r = m; 613 x = x*t%p; 614 b = tempb*y%p; 615 } 616 617 assert(x.Squared()%p == a); 618 return x; 619} 620 621bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p) 622{ 623 Integer D = (b.Squared() - 4*a*c) % p; 624 switch (Jacobi(D, p)) 625 { 626 default: 627 assert(false); // not reached 628 return false; 629 case -1: 630 return false; 631 case 0: 632 r1 = r2 = (-b*(a+a).InverseMod(p)) % p; 633 assert(((r1.Squared()*a + r1*b + c) % p).IsZero()); 634 return true; 635 case 1: 636 Integer s = ModularSquareRoot(D, p); 637 Integer t = (a+a).InverseMod(p); 638 r1 = (s-b)*t % p; 639 r2 = (-s-b)*t % p; 640 assert(((r1.Squared()*a + r1*b + c) % p).IsZero()); 641 assert(((r2.Squared()*a + r2*b + c) % p).IsZero()); 642 return true; 643 } 644} 645 646Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq, 647 const Integer &p, const Integer &q, const Integer &u) 648{ 649 Integer p2, q2; 650 #pragma omp parallel 651 #pragma omp sections 652 { 653 #pragma omp section 654 p2 = ModularExponentiation((a % p), dp, p); 655 #pragma omp section 656 q2 = ModularExponentiation((a % q), dq, q); 657 } 658 return CRT(p2, p, q2, q, u); 659} 660 661Integer ModularRoot(const Integer &a, const Integer &e, 662 const Integer &p, const Integer &q) 663{ 664 Integer dp = EuclideanMultiplicativeInverse(e, p-1); 665 Integer dq = EuclideanMultiplicativeInverse(e, q-1); 666 Integer u = EuclideanMultiplicativeInverse(p, q); 667 assert(!!dp && !!dq && !!u); 668 return ModularRoot(a, dp, dq, p, q, u); 669} 670 671/* 672Integer GCDI(const Integer &x, const Integer &y) 673{ 674 Integer a=x, b=y; 675 unsigned k=0; 676 677 assert(!!a && !!b); 678 679 while (a[0]==0 && b[0]==0) 680 { 681 a >>= 1; 682 b >>= 1; 683 k++; 684 } 685 686 while (a[0]==0) 687 a >>= 1; 688 689 while (b[0]==0) 690 b >>= 1; 691 692 while (1) 693 { 694 switch (a.Compare(b)) 695 { 696 case -1: 697 b -= a; 698 while (b[0]==0) 699 b >>= 1; 700 break; 701 702 case 0: 703 return (a <<= k); 704 705 case 1: 706 a -= b; 707 while (a[0]==0) 708 a >>= 1; 709 break; 710 711 default: 712 assert(false); 713 } 714 } 715} 716 717Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b) 718{ 719 assert(b.Positive()); 720 721 if (a.Negative()) 722 return EuclideanMultiplicativeInverse(a%b, b); 723 724 if (b[0]==0) 725 { 726 if (!b || a[0]==0) 727 return Integer::Zero(); // no inverse 728 if (a==1) 729 return 1; 730 Integer u = EuclideanMultiplicativeInverse(b, a); 731 if (!u) 732 return Integer::Zero(); // no inverse 733 else 734 return (b*(a-u)+1)/a; 735 } 736 737 Integer u=1, d=a, v1=b, v3=b, t1, t3, b2=(b+1)>>1; 738 739 if (a[0]) 740 { 741 t1 = Integer::Zero(); 742 t3 = -b; 743 } 744 else 745 { 746 t1 = b2; 747 t3 = a>>1; 748 } 749 750 while (!!t3) 751 { 752 while (t3[0]==0) 753 { 754 t3 >>= 1; 755 if (t1[0]==0) 756 t1 >>= 1; 757 else 758 { 759 t1 >>= 1; 760 t1 += b2; 761 } 762 } 763 if (t3.Positive()) 764 { 765 u = t1; 766 d = t3; 767 } 768 else 769 { 770 v1 = b-t1; 771 v3 = -t3; 772 } 773 t1 = u-v1; 774 t3 = d-v3; 775 if (t1.Negative()) 776 t1 += b; 777 } 778 if (d==1) 779 return u; 780 else 781 return Integer::Zero(); // no inverse 782} 783*/ 784 785int Jacobi(const Integer &aIn, const Integer &bIn) 786{ 787 assert(bIn.IsOdd()); 788 789 Integer b = bIn, a = aIn%bIn; 790 int result = 1; 791 792 while (!!a) 793 { 794 unsigned i=0; 795 while (a.GetBit(i)==0) 796 i++; 797 a>>=i; 798 799 if (i%2==1 && (b%8==3 || b%8==5)) 800 result = -result; 801 802 if (a%4==3 && b%4==3) 803 result = -result; 804 805 std::swap(a, b); 806 a %= b; 807 } 808 809 return (b==1) ? result : 0; 810} 811 812Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n) 813{ 814 unsigned i = e.BitCount(); 815 if (i==0) 816 return Integer::Two(); 817 818 MontgomeryRepresentation m(n); 819 Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two()); 820 Integer v=p, v1=m.Subtract(m.Square(p), two); 821 822 i--; 823 while (i--) 824 { 825 if (e.GetBit(i)) 826 { 827 // v = (v*v1 - p) % m; 828 v = m.Subtract(m.Multiply(v,v1), p); 829 // v1 = (v1*v1 - 2) % m; 830 v1 = m.Subtract(m.Square(v1), two); 831 } 832 else 833 { 834 // v1 = (v*v1 - p) % m; 835 v1 = m.Subtract(m.Multiply(v,v1), p); 836 // v = (v*v - 2) % m; 837 v = m.Subtract(m.Square(v), two); 838 } 839 } 840 return m.ConvertOut(v); 841} 842 843// This is Peter Montgomery's unpublished Lucas sequence evalutation algorithm. 844// The total number of multiplies and squares used is less than the binary 845// algorithm (see above). Unfortunately I can't get it to run as fast as 846// the binary algorithm because of the extra overhead. 847/* 848Integer Lucas(const Integer &n, const Integer &P, const Integer &modulus) 849{ 850 if (!n) 851 return 2; 852 853#define f(A, B, C) m.Subtract(m.Multiply(A, B), C) 854#define X2(A) m.Subtract(m.Square(A), two) 855#define X3(A) m.Multiply(A, m.Subtract(m.Square(A), three)) 856 857 MontgomeryRepresentation m(modulus); 858 Integer two=m.ConvertIn(2), three=m.ConvertIn(3); 859 Integer A=m.ConvertIn(P), B, C, p, d=n, e, r, t, T, U; 860 861 while (d!=1) 862 { 863 p = d; 864 unsigned int b = WORD_BITS * p.WordCount(); 865 Integer alpha = (Integer(5)<<(2*b-2)).SquareRoot() - Integer::Power2(b-1); 866 r = (p*alpha)>>b; 867 e = d-r; 868 B = A; 869 C = two; 870 d = r; 871 872 while (d!=e) 873 { 874 if (d<e) 875 { 876 swap(d, e); 877 swap(A, B); 878 } 879 880 unsigned int dm2 = d[0], em2 = e[0]; 881 unsigned int dm3 = d%3, em3 = e%3; 882 883// if ((dm6+em6)%3 == 0 && d <= e + (e>>2)) 884 if ((dm3+em3==0 || dm3+em3==3) && (t = e, t >>= 2, t += e, d <= t)) 885 { 886 // #1 887// t = (d+d-e)/3; 888// t = d; t += d; t -= e; t /= 3; 889// e = (e+e-d)/3; 890// e += e; e -= d; e /= 3; 891// d = t; 892 893// t = (d+e)/3 894 t = d; t += e; t /= 3; 895 e -= t; 896 d -= t; 897 898 T = f(A, B, C); 899 U = f(T, A, B); 900 B = f(T, B, A); 901 A = U; 902 continue; 903 } 904 905// if (dm6 == em6 && d <= e + (e>>2)) 906 if (dm3 == em3 && dm2 == em2 && (t = e, t >>= 2, t += e, d <= t)) 907 { 908 // #2 909// d = (d-e)>>1; 910 d -= e; d >>= 1; 911 B = f(A, B, C); 912 A = X2(A); 913 continue; 914 } 915 916// if (d <= (e<<2)) 917 if (d <= (t = e, t <<= 2)) 918 { 919 // #3 920 d -= e; 921 C = f(A, B, C); 922 swap(B, C); 923 continue; 924 } 925 926 if (dm2 == em2) 927 { 928 // #4 929// d = (d-e)>>1; 930 d -= e; d >>= 1; 931 B = f(A, B, C); 932 A = X2(A); 933 continue; 934 } 935 936 if (dm2 == 0) 937 { 938 // #5 939 d >>= 1; 940 C = f(A, C, B); 941 A = X2(A); 942 continue; 943 } 944 945 if (dm3 == 0) 946 { 947 // #6 948// d = d/3 - e; 949 d /= 3; d -= e; 950 T = X2(A); 951 C = f(T, f(A, B, C), C); 952 swap(B, C); 953 A = f(T, A, A); 954 continue; 955 } 956 957 if (dm3+em3==0 || dm3+em3==3) 958 { 959 // #7 960// d = (d-e-e)/3; 961 d -= e; d -= e; d /= 3; 962 T = f(A, B, C); 963 B = f(T, A, B); 964 A = X3(A); 965 continue; 966 } 967 968 if (dm3 == em3) 969 { 970 // #8 971// d = (d-e)/3; 972 d -= e; d /= 3; 973 T = f(A, B, C); 974 C = f(A, C, B); 975 B = T; 976 A = X3(A); 977 continue; 978 } 979 980 assert(em2 == 0); 981 // #9 982 e >>= 1; 983 C = f(C, B, A); 984 B = X2(B); 985 } 986 987 A = f(A, B, C); 988 } 989 990#undef f 991#undef X2 992#undef X3 993 994 return m.ConvertOut(A); 995} 996*/ 997 998Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u) 999{ 1000 Integer d = (m*m-4); 1001 Integer p2, q2; 1002 #pragma omp parallel 1003 #pragma omp sections 1004 { 1005 #pragma omp section 1006 { 1007 p2 = p-Jacobi(d,p); 1008 p2 = Lucas(EuclideanMultiplicativeInverse(e,p2), m, p); 1009 } 1010 #pragma omp section 1011 { 1012 q2 = q-Jacobi(d,q); 1013 q2 = Lucas(EuclideanMultiplicativeInverse(e,q2), m, q); 1014 } 1015 } 1016 return CRT(p2, p, q2, q, u); 1017} 1018 1019unsigned int FactoringWorkFactor(unsigned int n) 1020{ 1021 // extrapolated from the table in Odlyzko's "The Future of Integer Factorization" 1022 // updated to reflect the factoring of RSA-130 1023 if (n<5) return 0; 1024 else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5); 1025} 1026 1027unsigned int DiscreteLogWorkFactor(unsigned int n) 1028{ 1029 // assuming discrete log takes about the same time as factoring 1030 if (n<5) return 0; 1031 else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5); 1032} 1033 1034// ******************************************************** 1035 1036void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits) 1037{ 1038 // no prime exists for delta = -1, qbits = 4, and pbits = 5 1039 assert(qbits > 4); 1040 assert(pbits > qbits); 1041 1042 if (qbits+1 == pbits) 1043 { 1044 Integer minP = Integer::Power2(pbits-1); 1045 Integer maxP = Integer::Power2(pbits) - 1; 1046 bool success = false; 1047 1048 while (!success) 1049 { 1050 p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12); 1051 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta); 1052 1053 while (sieve.NextCandidate(p)) 1054 { 1055 assert(IsSmallPrime(p) || SmallDivisorsTest(p)); 1056 q = (p-delta) >> 1; 1057 assert(IsSmallPrime(q) || SmallDivisorsTest(q)); 1058 if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p)) 1059 { 1060 success = true; 1061 break; 1062 } 1063 } 1064 } 1065 1066 if (delta == 1) 1067 { 1068 // find g such that g is a quadratic residue mod p, then g has order q 1069 // g=4 always works, but this way we get the smallest quadratic residue (other than 1) 1070 for (g=2; Jacobi(g, p) != 1; ++g) {} 1071 // contributed by Walt Tuvell: g should be the following according to the Law of Quadratic Reciprocity 1072 assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4); 1073 } 1074 else 1075 { 1076 assert(delta == -1); 1077 // find g such that g*g-4 is a quadratic non-residue, 1078 // and such that g has order q 1079 for (g=3; ; ++g) 1080 if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2) 1081 break; 1082 } 1083 } 1084 else 1085 { 1086 Integer minQ = Integer::Power2(qbits-1); 1087 Integer maxQ = Integer::Power2(qbits) - 1; 1088 Integer minP = Integer::Power2(pbits-1); 1089 Integer maxP = Integer::Power2(pbits) - 1; 1090 1091 do 1092 { 1093 q.Randomize(rng, minQ, maxQ, Integer::PRIME); 1094 } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q)); 1095 1096 // find a random g of order q 1097 if (delta==1) 1098 { 1099 do 1100 { 1101 Integer h(rng, 2, p-2, Integer::ANY); 1102 g = a_exp_b_mod_c(h, (p-1)/q, p); 1103 } while (g <= 1); 1104 assert(a_exp_b_mod_c(g, q, p)==1); 1105 } 1106 else 1107 { 1108 assert(delta==-1); 1109 do 1110 { 1111 Integer h(rng, 3, p-1, Integer::ANY); 1112 if (Jacobi(h*h-4, p)==1) 1113 continue; 1114 g = Lucas((p+1)/q, h, p); 1115 } while (g <= 2); 1116 assert(Lucas(q, g, p) == 2); 1117 } 1118 } 1119} 1120 1121NAMESPACE_END 1122 1123#endif 1124