00001
00002
00003
#include "pch.h"
00004
#include "nbtheory.h"
00005
#include "modarith.h"
00006
#include "algparam.h"
00007
00008
#include <math.h>
00009
#include <vector>
00010
00011 NAMESPACE_BEGIN(CryptoPP)
00012
00013 const
unsigned int maxPrimeTableSize = 3511;
00014 const word lastSmallPrime = 32719;
00015
unsigned int primeTableSize=552;
00016
00017 word primeTable[maxPrimeTableSize] =
00018 {2, 3, 5, 7, 11, 13, 17, 19,
00019 23, 29, 31, 37, 41, 43, 47, 53,
00020 59, 61, 67, 71, 73, 79, 83, 89,
00021 97, 101, 103, 107, 109, 113, 127, 131,
00022 137, 139, 149, 151, 157, 163, 167, 173,
00023 179, 181, 191, 193, 197, 199, 211, 223,
00024 227, 229, 233, 239, 241, 251, 257, 263,
00025 269, 271, 277, 281, 283, 293, 307, 311,
00026 313, 317, 331, 337, 347, 349, 353, 359,
00027 367, 373, 379, 383, 389, 397, 401, 409,
00028 419, 421, 431, 433, 439, 443, 449, 457,
00029 461, 463, 467, 479, 487, 491, 499, 503,
00030 509, 521, 523, 541, 547, 557, 563, 569,
00031 571, 577, 587, 593, 599, 601, 607, 613,
00032 617, 619, 631, 641, 643, 647, 653, 659,
00033 661, 673, 677, 683, 691, 701, 709, 719,
00034 727, 733, 739, 743, 751, 757, 761, 769,
00035 773, 787, 797, 809, 811, 821, 823, 827,
00036 829, 839, 853, 857, 859, 863, 877, 881,
00037 883, 887, 907, 911, 919, 929, 937, 941,
00038 947, 953, 967, 971, 977, 983, 991, 997,
00039 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049,
00040 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097,
00041 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163,
00042 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
00043 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283,
00044 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
00045 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423,
00046 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459,
00047 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
00048 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571,
00049 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619,
00050 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693,
00051 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747,
00052 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
00053 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877,
00054 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949,
00055 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
00056 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069,
00057 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
00058 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203,
00059 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267,
00060 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311,
00061 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377,
00062 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
00063 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503,
00064 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579,
00065 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657,
00066 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,
00067 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
00068 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801,
00069 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861,
00070 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939,
00071 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011,
00072 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
00073 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167,
00074 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221,
00075 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301,
00076 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347,
00077 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
00078 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491,
00079 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541,
00080 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607,
00081 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671,
00082 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
00083 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797,
00084 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863,
00085 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923,
00086 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003};
00087
00088
void BuildPrimeTable()
00089 {
00090
unsigned int p=primeTable[primeTableSize-1];
00091
for (
unsigned int i=primeTableSize; i<maxPrimeTableSize; i++)
00092 {
00093
int j;
00094
do
00095 {
00096 p+=2;
00097
for (j=1; j<54; j++)
00098
if (p%primeTable[j] == 0)
00099
break;
00100 }
while (j!=54);
00101 primeTable[i] = p;
00102 }
00103 primeTableSize = maxPrimeTableSize;
00104 assert(primeTable[primeTableSize-1] == lastSmallPrime);
00105 }
00106
00107
bool IsSmallPrime(
const Integer &p)
00108 {
00109 BuildPrimeTable();
00110
00111
if (p.
IsPositive() && p <= primeTable[primeTableSize-1])
00112
return std::binary_search(primeTable, primeTable+primeTableSize, (word)p.
ConvertToLong());
00113
else
00114
return false;
00115 }
00116
00117
bool TrialDivision(
const Integer &p,
unsigned bound)
00118 {
00119 assert(primeTable[primeTableSize-1] >= bound);
00120
00121
unsigned int i;
00122
for (i = 0; primeTable[i]<bound; i++)
00123
if ((p % primeTable[i]) == 0)
00124
return true;
00125
00126
if (bound == primeTable[i])
00127
return (p % bound == 0);
00128
else
00129
return false;
00130 }
00131
00132
bool SmallDivisorsTest(
const Integer &p)
00133 {
00134 BuildPrimeTable();
00135
return !TrialDivision(p, primeTable[primeTableSize-1]);
00136 }
00137
00138
bool IsFermatProbablePrime(
const Integer &n,
const Integer &b)
00139 {
00140
if (n <= 3)
00141
return n==2 || n==3;
00142
00143 assert(n>3 && b>1 && b<n-1);
00144
return a_exp_b_mod_c(b, n-1, n)==1;
00145 }
00146
00147
bool IsStrongProbablePrime(
const Integer &n,
const Integer &b)
00148 {
00149
if (n <= 3)
00150
return n==2 || n==3;
00151
00152 assert(n>3 && b>1 && b<n-1);
00153
00154
if ((n.
IsEven() && n!=2) || GCD(b, n) != 1)
00155
return false;
00156
00157
Integer nminus1 = (n-1);
00158
unsigned int a;
00159
00160
00161
for (a=0; ; a++)
00162
if (nminus1.
GetBit(a))
00163
break;
00164
Integer m = nminus1>>a;
00165
00166
Integer z = a_exp_b_mod_c(b, m, n);
00167
if (z==1 || z==nminus1)
00168
return true;
00169
for (
unsigned j=1; j<a; j++)
00170 {
00171 z = z.Squared()%n;
00172
if (z==nminus1)
00173
return true;
00174
if (z==1)
00175
return false;
00176 }
00177
return false;
00178 }
00179
00180
bool RabinMillerTest(
RandomNumberGenerator &rng,
const Integer &n,
unsigned int rounds)
00181 {
00182
if (n <= 3)
00183
return n==2 || n==3;
00184
00185 assert(n>3);
00186
00187
Integer b;
00188
for (
unsigned int i=0; i<rounds; i++)
00189 {
00190 b.Randomize(rng, 2, n-2);
00191
if (!IsStrongProbablePrime(n, b))
00192
return false;
00193 }
00194
return true;
00195 }
00196
00197
bool IsLucasProbablePrime(
const Integer &n)
00198 {
00199
if (n <= 1)
00200
return false;
00201
00202
if (n.
IsEven())
00203
return n==2;
00204
00205 assert(n>2);
00206
00207
Integer b=3;
00208
unsigned int i=0;
00209
int j;
00210
00211
while ((j=Jacobi(b.Squared()-4, n)) == 1)
00212 {
00213
if (++i==64 && n.
IsSquare())
00214
return false;
00215 ++b; ++b;
00216 }
00217
00218
if (j==0)
00219
return false;
00220
else
00221
return Lucas(n+1, b, n)==2;
00222 }
00223
00224
bool IsStrongLucasProbablePrime(
const Integer &n)
00225 {
00226
if (n <= 1)
00227
return false;
00228
00229
if (n.
IsEven())
00230
return n==2;
00231
00232 assert(n>2);
00233
00234
Integer b=3;
00235
unsigned int i=0;
00236
int j;
00237
00238
while ((j=Jacobi(b.Squared()-4, n)) == 1)
00239 {
00240
if (++i==64 && n.
IsSquare())
00241
return false;
00242 ++b; ++b;
00243 }
00244
00245
if (j==0)
00246
return false;
00247
00248
Integer n1 = n+1;
00249
unsigned int a;
00250
00251
00252
for (a=0; ; a++)
00253
if (n1.
GetBit(a))
00254
break;
00255
Integer m = n1>>a;
00256
00257
Integer z = Lucas(m, b, n);
00258
if (z==2 || z==n-2)
00259
return true;
00260
for (i=1; i<a; i++)
00261 {
00262 z = (z.Squared()-2)%n;
00263
if (z==n-2)
00264
return true;
00265
if (z==2)
00266
return false;
00267 }
00268
return false;
00269 }
00270
00271
bool IsPrime(
const Integer &p)
00272 {
00273
static const Integer lastSmallPrimeSquared =
Integer(lastSmallPrime).Squared();
00274
00275
if (p <= lastSmallPrime)
00276
return IsSmallPrime(p);
00277
else if (p <= lastSmallPrimeSquared)
00278
return SmallDivisorsTest(p);
00279
else
00280
return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p);
00281 }
00282
00283
bool VerifyPrime(
RandomNumberGenerator &rng,
const Integer &p,
unsigned int level)
00284 {
00285
bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1);
00286
if (level >= 1)
00287 pass = pass && RabinMillerTest(rng, p, 10);
00288
return pass;
00289 }
00290
00291
unsigned int PrimeSearchInterval(
const Integer &max)
00292 {
00293
return max.
BitCount();
00294 }
00295
00296
static inline bool FastProbablePrimeTest(
const Integer &n)
00297 {
00298
return IsStrongProbablePrime(n,2);
00299 }
00300
00301 AlgorithmParameters<AlgorithmParameters<AlgorithmParameters<NullNameValuePairs, Integer::RandomNumberType>,
Integer>, Integer>
00302 MakeParametersForTwoPrimesOfEqualSize(
unsigned int productBitLength)
00303 {
00304
if (productBitLength < 16)
00305
throw InvalidArgument(
"invalid bit length");
00306
00307 Integer minP, maxP;
00308
00309
if (productBitLength%2==0)
00310 {
00311 minP = Integer(182) << (productBitLength/2-8);
00312 maxP =
Integer::Power2(productBitLength/2)-1;
00313 }
00314
else
00315 {
00316 minP =
Integer::Power2((productBitLength-1)/2);
00317 maxP = Integer(181) << ((productBitLength+1)/2-8);
00318 }
00319
00320
return MakeParameters(
"RandomNumberType", Integer::PRIME)(
"Min", minP)(
"Max", maxP);
00321 }
00322
00323
class PrimeSieve
00324 {
00325
public:
00326
00327 PrimeSieve(
const Integer &first,
const Integer &last,
const Integer &step,
signed int delta=0);
00328
bool NextCandidate(Integer &c);
00329
00330
void DoSieve();
00331
static void SieveSingle(std::vector<bool> &sieve, word p,
const Integer &first,
const Integer &step, word stepInv);
00332
00333 Integer m_first, m_last, m_step;
00334
signed int m_delta;
00335 word m_next;
00336 std::vector<bool> m_sieve;
00337 };
00338
00339 PrimeSieve::PrimeSieve(
const Integer &first,
const Integer &last,
const Integer &step,
signed int delta)
00340 : m_first(first), m_last(last), m_step(step), m_delta(delta), m_next(0)
00341 {
00342 DoSieve();
00343 }
00344
00345
bool PrimeSieve::NextCandidate(Integer &c)
00346 {
00347 m_next = std::find(m_sieve.begin()+m_next, m_sieve.end(),
false) - m_sieve.begin();
00348
if (m_next == m_sieve.size())
00349 {
00350 m_first += m_sieve.size()*m_step;
00351
if (m_first > m_last)
00352
return false;
00353
else
00354 {
00355 m_next = 0;
00356 DoSieve();
00357
return NextCandidate(c);
00358 }
00359 }
00360
else
00361 {
00362 c = m_first + m_next*m_step;
00363 ++m_next;
00364
return true;
00365 }
00366 }
00367
00368
void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word p,
const Integer &first,
const Integer &step, word stepInv)
00369 {
00370
if (stepInv)
00371 {
00372
unsigned int sieveSize = sieve.size();
00373 word j = word((dword(p-(first%p))*stepInv) % p);
00374
00375
if (first.
WordCount() <= 1 && first + step*j == p)
00376 j += p;
00377
for (; j < sieveSize; j += p)
00378 sieve[j] =
true;
00379 }
00380 }
00381
00382
void PrimeSieve::DoSieve()
00383 {
00384 BuildPrimeTable();
00385
00386
const unsigned int maxSieveSize = 32768;
00387
unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong();
00388
00389 m_sieve.clear();
00390 m_sieve.resize(sieveSize,
false);
00391
00392
if (m_delta == 0)
00393 {
00394
for (
unsigned int i = 0; i < primeTableSize; ++i)
00395 SieveSingle(m_sieve, primeTable[i], m_first, m_step, m_step.InverseMod(primeTable[i]));
00396 }
00397
else
00398 {
00399 assert(m_step%2==0);
00400 Integer qFirst = (m_first-m_delta) >> 1;
00401 Integer halfStep = m_step >> 1;
00402
for (
unsigned int i = 0; i < primeTableSize; ++i)
00403 {
00404 word p = primeTable[i];
00405 word stepInv = m_step.
InverseMod(p);
00406 SieveSingle(m_sieve, p, m_first, m_step, stepInv);
00407
00408 word halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p;
00409 SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
00410 }
00411 }
00412 }
00413
00414
bool FirstPrime(Integer &p,
const Integer &max,
const Integer &equiv,
const Integer &mod,
const PrimeSelector *pSelector)
00415 {
00416 assert(!equiv.
IsNegative() && equiv < mod);
00417
00418 Integer gcd = GCD(equiv, mod);
00419
if (gcd !=
Integer::One())
00420 {
00421
00422
if (p <= gcd && gcd <= max && IsPrime(gcd))
00423 {
00424 p = gcd;
00425
return true;
00426 }
00427
else
00428
return false;
00429 }
00430
00431 BuildPrimeTable();
00432
00433
if (p <= primeTable[primeTableSize-1])
00434 {
00435 word *pItr;
00436
00437 --p;
00438
if (p.
IsPositive())
00439 pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.
ConvertToLong());
00440
else
00441 pItr = primeTable;
00442
00443
while (pItr < primeTable+primeTableSize && *pItr%mod != equiv)
00444 ++pItr;
00445
00446
if (pItr < primeTable+primeTableSize)
00447 {
00448 p = *pItr;
00449
return p <= max;
00450 }
00451
00452 p = primeTable[primeTableSize-1]+1;
00453 }
00454
00455 assert(p > primeTable[primeTableSize-1]);
00456
00457
if (mod.
IsOdd())
00458
return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector);
00459
00460 p += (equiv-p)%mod;
00461
00462
if (p>max)
00463
return false;
00464
00465 PrimeSieve sieve(p, max, mod);
00466
00467
while (sieve.NextCandidate(p))
00468 {
00469
if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p))
00470
return true;
00471 }
00472
00473
return false;
00474 }
00475
00476
00477
static bool ProvePrime(
const Integer &p,
const Integer &q)
00478 {
00479 assert(p < q*q*q);
00480 assert(p % q == 1);
00481
00482
00483
00484
00485
00486
00487 Integer r = (p-1)/q;
00488
if (((r%q).Squared()-4*(r/q)).IsSquare())
00489
return false;
00490
00491 assert(primeTableSize >= 50);
00492
for (
int i=0; i<50; i++)
00493 {
00494 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
00495
if (b != 1)
00496
return a_exp_b_mod_c(b, q, p) == 1;
00497 }
00498
return false;
00499 }
00500
00501 Integer MihailescuProvablePrime(
RandomNumberGenerator &rng,
unsigned int pbits)
00502 {
00503 Integer p;
00504 Integer minP =
Integer::Power2(pbits-1);
00505 Integer maxP =
Integer::Power2(pbits) - 1;
00506
00507
if (maxP <= Integer(lastSmallPrime).
Squared())
00508 {
00509
00510 p.
Randomize(rng, minP, maxP, Integer::PRIME);
00511
return p;
00512 }
00513
00514
unsigned int qbits = (pbits+2)/3 + 1 + rng.
GenerateWord32(0, pbits/36);
00515 Integer q = MihailescuProvablePrime(rng, qbits);
00516 Integer q2 = q<<1;
00517
00518
while (
true)
00519 {
00520
00521
00522
00523
00524
00525
00526
00527 p.
Randomize(rng, minP, maxP, Integer::ANY, 1, q2);
00528 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
00529
00530
while (sieve.NextCandidate(p))
00531 {
00532
if (FastProbablePrimeTest(p) && ProvePrime(p, q))
00533
return p;
00534 }
00535 }
00536
00537
00538
return p;
00539 }
00540
00541 Integer MaurerProvablePrime(
RandomNumberGenerator &rng,
unsigned int bits)
00542 {
00543
const unsigned smallPrimeBound = 29, c_opt=10;
00544 Integer p;
00545
00546 BuildPrimeTable();
00547
if (bits < smallPrimeBound)
00548 {
00549
do
00550 p.
Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2);
00551
while (TrialDivision(p, 1 << ((bits+1)/2)));
00552 }
00553
else
00554 {
00555
const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
00556
double relativeSize;
00557
do
00558 relativeSize = pow(2.0,
double(rng.
GenerateWord32())/0xffffffff - 1);
00559
while (bits * relativeSize >= bits - margin);
00560
00561 Integer a,b;
00562 Integer q = MaurerProvablePrime(rng,
unsigned(bits*relativeSize));
00563 Integer I =
Integer::Power2(bits-2)/q;
00564 Integer I2 = I << 1;
00565
unsigned int trialDivisorBound = (
unsigned int)STDMIN((
unsigned long)primeTable[primeTableSize-1], (
unsigned long)bits*bits/c_opt);
00566
bool success =
false;
00567
while (!success)
00568 {
00569 p.
Randomize(rng, I, I2, Integer::ANY);
00570 p *= q; p <<= 1; ++p;
00571
if (!TrialDivision(p, trialDivisorBound))
00572 {
00573 a.Randomize(rng, 2, p-1, Integer::ANY);
00574 b = a_exp_b_mod_c(a, (p-1)/q, p);
00575 success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1);
00576 }
00577 }
00578 }
00579
return p;
00580 }
00581
00582 Integer CRT(
const Integer &xp,
const Integer &p,
const Integer &xq,
const Integer &q,
const Integer &u)
00583 {
00584
00585
return p * (u * (xq-xp) % q) + xp;
00586 }
00587
00588 Integer CRT(
const Integer &xp,
const Integer &p,
const Integer &xq,
const Integer &q)
00589 {
00590
return CRT(xp, p, xq, q, EuclideanMultiplicativeInverse(p, q));
00591 }
00592
00593 Integer ModularSquareRoot(
const Integer &a,
const Integer &p)
00594 {
00595
if (p%4 == 3)
00596
return a_exp_b_mod_c(a, (p+1)/4, p);
00597
00598 Integer q=p-1;
00599
unsigned int r=0;
00600
while (q.
IsEven())
00601 {
00602 r++;
00603 q >>= 1;
00604 }
00605
00606 Integer n=2;
00607
while (Jacobi(n, p) != -1)
00608 ++n;
00609
00610 Integer y = a_exp_b_mod_c(n, q, p);
00611 Integer x = a_exp_b_mod_c(a, (q-1)/2, p);
00612 Integer b = (x.Squared()%p)*a%p;
00613 x = a*x%p;
00614 Integer tempb, t;
00615
00616
while (b != 1)
00617 {
00618
unsigned m=0;
00619 tempb = b;
00620
do
00621 {
00622 m++;
00623 b = b.Squared()%p;
00624
if (m==r)
00625
return Integer::Zero();
00626 }
00627
while (b != 1);
00628
00629 t = y;
00630
for (
unsigned i=0; i<r-m-1; i++)
00631 t = t.Squared()%p;
00632 y = t.Squared()%p;
00633 r = m;
00634 x = x*t%p;
00635 b = tempb*y%p;
00636 }
00637
00638 assert(x.Squared()%p == a);
00639
return x;
00640 }
00641
00642
bool SolveModularQuadraticEquation(Integer &r1, Integer &r2,
const Integer &a,
const Integer &b,
const Integer &c,
const Integer &p)
00643 {
00644 Integer D = (b.Squared() - 4*a*c) % p;
00645
switch (Jacobi(D, p))
00646 {
00647
default:
00648 assert(
false);
00649
return false;
00650
case -1:
00651
return false;
00652
case 0:
00653 r1 = r2 = (-b*(a+a).InverseMod(p)) % p;
00654 assert(((r1.
Squared()*a + r1*b + c) % p).IsZero());
00655
return true;
00656
case 1:
00657 Integer s = ModularSquareRoot(D, p);
00658 Integer t = (a+a).InverseMod(p);
00659 r1 = (s-b)*t % p;
00660 r2 = (-s-b)*t % p;
00661 assert(((r1.
Squared()*a + r1*b + c) % p).IsZero());
00662 assert(((r2.Squared()*a + r2*b + c) % p).IsZero());
00663
return true;
00664 }
00665 }
00666
00667 Integer ModularRoot(
const Integer &a,
const Integer &dp,
const Integer &dq,
00668
const Integer &p,
const Integer &q,
const Integer &u)
00669 {
00670 Integer p2 = ModularExponentiation((a % p), dp, p);
00671 Integer q2 = ModularExponentiation((a % q), dq, q);
00672
return CRT(p2, p, q2, q, u);
00673 }
00674
00675 Integer ModularRoot(
const Integer &a,
const Integer &e,
00676
const Integer &p,
const Integer &q)
00677 {
00678 Integer dp = EuclideanMultiplicativeInverse(e, p-1);
00679 Integer dq = EuclideanMultiplicativeInverse(e, q-1);
00680 Integer u = EuclideanMultiplicativeInverse(p, q);
00681 assert(!!dp && !!dq && !!u);
00682
return ModularRoot(a, dp, dq, p, q, u);
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
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
int Jacobi(
const Integer &aIn,
const Integer &bIn)
00800 {
00801 assert(bIn.
IsOdd());
00802
00803 Integer b = bIn, a = aIn%bIn;
00804
int result = 1;
00805
00806
while (!!a)
00807 {
00808
unsigned i=0;
00809
while (a.GetBit(i)==0)
00810 i++;
00811 a>>=i;
00812
00813
if (i%2==1 && (b%8==3 || b%8==5))
00814 result = -result;
00815
00816
if (a%4==3 && b%4==3)
00817 result = -result;
00818
00819 std::swap(a, b);
00820 a %= b;
00821 }
00822
00823
return (b==1) ? result : 0;
00824 }
00825
00826 Integer Lucas(
const Integer &e,
const Integer &pIn,
const Integer &n)
00827 {
00828
unsigned i = e.BitCount();
00829
if (i==0)
00830
return Integer::Two();
00831
00832
MontgomeryRepresentation m(n);
00833 Integer p=m.
ConvertIn(pIn%n), two=m.
ConvertIn(Integer::Two());
00834 Integer v=p, v1=m.
Subtract(m.
Square(p), two);
00835
00836 i--;
00837
while (i--)
00838 {
00839
if (e.GetBit(i))
00840 {
00841
00842 v = m.
Subtract(m.
Multiply(v,v1), p);
00843
00844 v1 = m.
Subtract(m.
Square(v1), two);
00845 }
00846
else
00847 {
00848
00849 v1 = m.
Subtract(m.
Multiply(v,v1), p);
00850
00851 v = m.
Subtract(m.
Square(v), two);
00852 }
00853 }
00854
return m.
ConvertOut(v);
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
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012 Integer InverseLucas(
const Integer &e,
const Integer &m,
const Integer &p,
const Integer &q,
const Integer &u)
01013 {
01014 Integer d = (m*m-4);
01015 Integer p2 = p-Jacobi(d,p);
01016 Integer q2 = q-Jacobi(d,q);
01017
return CRT(Lucas(EuclideanMultiplicativeInverse(e,p2), m, p), p, Lucas(EuclideanMultiplicativeInverse(e,q2), m, q), q, u);
01018 }
01019
01020 Integer InverseLucas(
const Integer &e,
const Integer &m,
const Integer &p,
const Integer &q)
01021 {
01022
return InverseLucas(e, m, p, q, EuclideanMultiplicativeInverse(p, q));
01023 }
01024
01025
unsigned int FactoringWorkFactor(
unsigned int n)
01026 {
01027
01028
01029
if (n<5)
return 0;
01030
else return (
unsigned int)(2.4 * pow((
double)n, 1.0/3.0) * pow(log(
double(n)), 2.0/3.0) - 5);
01031 }
01032
01033
unsigned int DiscreteLogWorkFactor(
unsigned int n)
01034 {
01035
01036
if (n<5)
return 0;
01037
else return (
unsigned int)(2.4 * pow((
double)n, 1.0/3.0) * pow(log(
double(n)), 2.0/3.0) - 5);
01038 }
01039
01040
01041
01042
void PrimeAndGenerator::Generate(
signed int delta,
RandomNumberGenerator &rng,
unsigned int pbits,
unsigned int qbits)
01043 {
01044
01045 assert(qbits > 4);
01046 assert(pbits > qbits);
01047
01048
if (qbits+1 == pbits)
01049 {
01050 Integer minP =
Integer::Power2(pbits-1);
01051 Integer maxP =
Integer::Power2(pbits) - 1;
01052
bool success =
false;
01053
01054
while (!success)
01055 {
01056 p.
Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
01057 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
01058
01059
while (sieve.NextCandidate(p))
01060 {
01061 assert(IsSmallPrime(p) || SmallDivisorsTest(p));
01062 q = (p-delta) >> 1;
01063 assert(IsSmallPrime(q) || SmallDivisorsTest(q));
01064
if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
01065 {
01066 success =
true;
01067
break;
01068 }
01069 }
01070 }
01071
01072
if (delta == 1)
01073 {
01074
01075
01076
for (g=2; Jacobi(g, p) != 1; ++g) {}
01077
01078 assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
01079 }
01080
else
01081 {
01082 assert(delta == -1);
01083
01084
01085
for (g=3; ; ++g)
01086
if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
01087
break;
01088 }
01089 }
01090
else
01091 {
01092 Integer minQ =
Integer::Power2(qbits-1);
01093 Integer maxQ =
Integer::Power2(qbits) - 1;
01094 Integer minP =
Integer::Power2(pbits-1);
01095 Integer maxP =
Integer::Power2(pbits) - 1;
01096
01097
do
01098 {
01099 q.
Randomize(rng, minQ, maxQ, Integer::PRIME);
01100 }
while (!p.
Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
01101
01102
01103
if (delta==1)
01104 {
01105
do
01106 {
01107 Integer h(rng, 2, p-2, Integer::ANY);
01108 g = a_exp_b_mod_c(h, (p-1)/q, p);
01109 }
while (g <= 1);
01110 assert(a_exp_b_mod_c(g, q, p)==1);
01111 }
01112
else
01113 {
01114 assert(delta==-1);
01115
do
01116 {
01117 Integer h(rng, 3, p-1, Integer::ANY);
01118
if (Jacobi(h*h-4, p)==1)
01119
continue;
01120 g = Lucas((p+1)/q, h, p);
01121 }
while (g <= 2);
01122 assert(Lucas(q, g, p) == 2);
01123 }
01124 }
01125 }
01126
01127 NAMESPACE_END