00001
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
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())
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())
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
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
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
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
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
00438
static bool ProvePrime(
const Integer &p,
const Integer &q)
00439 {
00440 assert(p < q*q*q);
00441 assert(p % q == 1);
00442
00443
00444
00445
00446
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
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
00485
00486
00487
00488
00489
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
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
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);
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
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
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
00808 v = m.Subtract(m.Multiply(v,v1), p);
00809
00810 v1 = m.Subtract(m.Square(v1), two);
00811 }
00812
else
00813 {
00814
00815 v1 = m.Subtract(m.Multiply(v,v1), p);
00816
00817 v = m.Subtract(m.Square(v), two);
00818 }
00819 }
00820
return m.ConvertOut(v);
00821 }
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
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
00994
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
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
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
01041
01042
for (g=2; Jacobi(g, p) != 1; ++g) {}
01043
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
01050
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
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