monotone

monotone Mtn Source Tree

Root/cryptopp/nbtheory.cpp

1// nbtheory.cpp - written and placed in the public domain by Wei Dai
2
3#include "pch.h"
4#include "nbtheory.h"
5#include "modarith.h"
6#include "algparam.h"
7
8#include <math.h>
9#include <vector>
10
11NAMESPACE_BEGIN(CryptoPP)
12
13const unsigned int maxPrimeTableSize = 3511;// last prime 32719
14const word lastSmallPrime = 32719;
15unsigned int primeTableSize=552;
16
17word primeTable[maxPrimeTableSize] =
18{2, 3, 5, 7, 11, 13, 17, 19,
1923, 29, 31, 37, 41, 43, 47, 53,
2059, 61, 67, 71, 73, 79, 83, 89,
2197, 101, 103, 107, 109, 113, 127, 131,
22137, 139, 149, 151, 157, 163, 167, 173,
23179, 181, 191, 193, 197, 199, 211, 223,
24227, 229, 233, 239, 241, 251, 257, 263,
25269, 271, 277, 281, 283, 293, 307, 311,
26313, 317, 331, 337, 347, 349, 353, 359,
27367, 373, 379, 383, 389, 397, 401, 409,
28419, 421, 431, 433, 439, 443, 449, 457,
29461, 463, 467, 479, 487, 491, 499, 503,
30509, 521, 523, 541, 547, 557, 563, 569,
31571, 577, 587, 593, 599, 601, 607, 613,
32617, 619, 631, 641, 643, 647, 653, 659,
33661, 673, 677, 683, 691, 701, 709, 719,
34727, 733, 739, 743, 751, 757, 761, 769,
35773, 787, 797, 809, 811, 821, 823, 827,
36829, 839, 853, 857, 859, 863, 877, 881,
37883, 887, 907, 911, 919, 929, 937, 941,
38947, 953, 967, 971, 977, 983, 991, 997,
391009, 1013, 1019, 1021, 1031, 1033, 1039, 1049,
401051, 1061, 1063, 1069, 1087, 1091, 1093, 1097,
411103, 1109, 1117, 1123, 1129, 1151, 1153, 1163,
421171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
431229, 1231, 1237, 1249, 1259, 1277, 1279, 1283,
441289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
451327, 1361, 1367, 1373, 1381, 1399, 1409, 1423,
461427, 1429, 1433, 1439, 1447, 1451, 1453, 1459,
471471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
481523, 1531, 1543, 1549, 1553, 1559, 1567, 1571,
491579, 1583, 1597, 1601, 1607, 1609, 1613, 1619,
501621, 1627, 1637, 1657, 1663, 1667, 1669, 1693,
511697, 1699, 1709, 1721, 1723, 1733, 1741, 1747,
521753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
531823, 1831, 1847, 1861, 1867, 1871, 1873, 1877,
541879, 1889, 1901, 1907, 1913, 1931, 1933, 1949,
551951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
562011, 2017, 2027, 2029, 2039, 2053, 2063, 2069,
572081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
582131, 2137, 2141, 2143, 2153, 2161, 2179, 2203,
592207, 2213, 2221, 2237, 2239, 2243, 2251, 2267,
602269, 2273, 2281, 2287, 2293, 2297, 2309, 2311,
612333, 2339, 2341, 2347, 2351, 2357, 2371, 2377,
622381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
632437, 2441, 2447, 2459, 2467, 2473, 2477, 2503,
642521, 2531, 2539, 2543, 2549, 2551, 2557, 2579,
652591, 2593, 2609, 2617, 2621, 2633, 2647, 2657,
662659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,
672699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
682749, 2753, 2767, 2777, 2789, 2791, 2797, 2801,
692803, 2819, 2833, 2837, 2843, 2851, 2857, 2861,
702879, 2887, 2897, 2903, 2909, 2917, 2927, 2939,
712953, 2957, 2963, 2969, 2971, 2999, 3001, 3011,
723019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
733083, 3089, 3109, 3119, 3121, 3137, 3163, 3167,
743169, 3181, 3187, 3191, 3203, 3209, 3217, 3221,
753229, 3251, 3253, 3257, 3259, 3271, 3299, 3301,
763307, 3313, 3319, 3323, 3329, 3331, 3343, 3347,
773359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
783433, 3449, 3457, 3461, 3463, 3467, 3469, 3491,
793499, 3511, 3517, 3527, 3529, 3533, 3539, 3541,
803547, 3557, 3559, 3571, 3581, 3583, 3593, 3607,
813613, 3617, 3623, 3631, 3637, 3643, 3659, 3671,
823673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
833733, 3739, 3761, 3767, 3769, 3779, 3793, 3797,
843803, 3821, 3823, 3833, 3847, 3851, 3853, 3863,
853877, 3881, 3889, 3907, 3911, 3917, 3919, 3923,
863929, 3931, 3943, 3947, 3967, 3989, 4001, 4003};
87
88void BuildPrimeTable()
89{
90unsigned int p=primeTable[primeTableSize-1];
91for (unsigned int i=primeTableSize; i<maxPrimeTableSize; i++)
92{
93int j;
94do
95{
96p+=2;
97for (j=1; j<54; j++)
98if (p%primeTable[j] == 0)
99break;
100} while (j!=54);
101primeTable[i] = p;
102}
103primeTableSize = maxPrimeTableSize;
104assert(primeTable[primeTableSize-1] == lastSmallPrime);
105}
106
107bool IsSmallPrime(const Integer &p)
108{
109BuildPrimeTable();
110
111if (p.IsPositive() && p <= primeTable[primeTableSize-1])
112return std::binary_search(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
113else
114return false;
115}
116
117bool TrialDivision(const Integer &p, unsigned bound)
118{
119assert(primeTable[primeTableSize-1] >= bound);
120
121unsigned int i;
122for (i = 0; primeTable[i]<bound; i++)
123if ((p % primeTable[i]) == 0)
124return true;
125
126if (bound == primeTable[i])
127return (p % bound == 0);
128else
129return false;
130}
131
132bool SmallDivisorsTest(const Integer &p)
133{
134BuildPrimeTable();
135return !TrialDivision(p, primeTable[primeTableSize-1]);
136}
137
138bool IsFermatProbablePrime(const Integer &n, const Integer &b)
139{
140if (n <= 3)
141return n==2 || n==3;
142
143assert(n>3 && b>1 && b<n-1);
144return a_exp_b_mod_c(b, n-1, n)==1;
145}
146
147bool IsStrongProbablePrime(const Integer &n, const Integer &b)
148{
149if (n <= 3)
150return n==2 || n==3;
151
152assert(n>3 && b>1 && b<n-1);
153
154if ((n.IsEven() && n!=2) || GCD(b, n) != 1)
155return false;
156
157Integer nminus1 = (n-1);
158unsigned int a;
159
160// calculate a = largest power of 2 that divides (n-1)
161for (a=0; ; a++)
162if (nminus1.GetBit(a))
163break;
164Integer m = nminus1>>a;
165
166Integer z = a_exp_b_mod_c(b, m, n);
167if (z==1 || z==nminus1)
168return true;
169for (unsigned j=1; j<a; j++)
170{
171z = z.Squared()%n;
172if (z==nminus1)
173return true;
174if (z==1)
175return false;
176}
177return false;
178}
179
180bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds)
181{
182if (n <= 3)
183return n==2 || n==3;
184
185assert(n>3);
186
187Integer b;
188for (unsigned int i=0; i<rounds; i++)
189{
190b.Randomize(rng, 2, n-2);
191if (!IsStrongProbablePrime(n, b))
192return false;
193}
194return true;
195}
196
197bool IsLucasProbablePrime(const Integer &n)
198{
199if (n <= 1)
200return false;
201
202if (n.IsEven())
203return n==2;
204
205assert(n>2);
206
207Integer b=3;
208unsigned int i=0;
209int j;
210
211while ((j=Jacobi(b.Squared()-4, n)) == 1)
212{
213if (++i==64 && n.IsSquare())// avoid infinite loop if n is a square
214return false;
215++b; ++b;
216}
217
218if (j==0)
219return false;
220else
221return Lucas(n+1, b, n)==2;
222}
223
224bool IsStrongLucasProbablePrime(const Integer &n)
225{
226if (n <= 1)
227return false;
228
229if (n.IsEven())
230return n==2;
231
232assert(n>2);
233
234Integer b=3;
235unsigned int i=0;
236int j;
237
238while ((j=Jacobi(b.Squared()-4, n)) == 1)
239{
240if (++i==64 && n.IsSquare())// avoid infinite loop if n is a square
241return false;
242++b; ++b;
243}
244
245if (j==0)
246return false;
247
248Integer n1 = n+1;
249unsigned int a;
250
251// calculate a = largest power of 2 that divides n1
252for (a=0; ; a++)
253if (n1.GetBit(a))
254break;
255Integer m = n1>>a;
256
257Integer z = Lucas(m, b, n);
258if (z==2 || z==n-2)
259return true;
260for (i=1; i<a; i++)
261{
262z = (z.Squared()-2)%n;
263if (z==n-2)
264return true;
265if (z==2)
266return false;
267}
268return false;
269}
270
271bool IsPrime(const Integer &p)
272{
273static const Integer lastSmallPrimeSquared = Integer(lastSmallPrime).Squared();
274
275if (p <= lastSmallPrime)
276return IsSmallPrime(p);
277else if (p <= lastSmallPrimeSquared)
278return SmallDivisorsTest(p);
279else
280return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p);
281}
282
283bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level)
284{
285bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1);
286if (level >= 1)
287pass = pass && RabinMillerTest(rng, p, 10);
288return pass;
289}
290
291unsigned int PrimeSearchInterval(const Integer &max)
292{
293return max.BitCount();
294}
295
296static inline bool FastProbablePrimeTest(const Integer &n)
297{
298return IsStrongProbablePrime(n,2);
299}
300
301AlgorithmParameters<AlgorithmParameters<AlgorithmParameters<NullNameValuePairs, Integer::RandomNumberType>, Integer>, Integer>
302MakeParametersForTwoPrimesOfEqualSize(unsigned int productBitLength)
303{
304if (productBitLength < 16)
305throw InvalidArgument("invalid bit length");
306
307Integer minP, maxP;
308
309if (productBitLength%2==0)
310{
311minP = Integer(182) << (productBitLength/2-8);
312maxP = Integer::Power2(productBitLength/2)-1;
313}
314else
315{
316minP = Integer::Power2((productBitLength-1)/2);
317maxP = Integer(181) << ((productBitLength+1)/2-8);
318}
319
320return MakeParameters("RandomNumberType", Integer::PRIME)("Min", minP)("Max", maxP);
321}
322
323class PrimeSieve
324{
325public:
326// delta == 1 or -1 means double sieve with p = 2*q + delta
327PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta=0);
328bool NextCandidate(Integer &c);
329
330void DoSieve();
331static void SieveSingle(std::vector<bool> &sieve, word p, const Integer &first, const Integer &step, word stepInv);
332
333Integer m_first, m_last, m_step;
334signed int m_delta;
335word m_next;
336std::vector<bool> m_sieve;
337};
338
339PrimeSieve::PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta)
340: m_first(first), m_last(last), m_step(step), m_delta(delta), m_next(0)
341{
342DoSieve();
343}
344
345bool PrimeSieve::NextCandidate(Integer &c)
346{
347m_next = std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin();
348if (m_next == m_sieve.size())
349{
350m_first += m_sieve.size()*m_step;
351if (m_first > m_last)
352return false;
353else
354{
355m_next = 0;
356DoSieve();
357return NextCandidate(c);
358}
359}
360else
361{
362c = m_first + m_next*m_step;
363++m_next;
364return true;
365}
366}
367
368void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word p, const Integer &first, const Integer &step, word stepInv)
369{
370if (stepInv)
371{
372unsigned int sieveSize = sieve.size();
373word j = word((dword(p-(first%p))*stepInv) % p);
374// if the first multiple of p is p, skip it
375if (first.WordCount() <= 1 && first + step*j == p)
376j += p;
377for (; j < sieveSize; j += p)
378sieve[j] = true;
379}
380}
381
382void PrimeSieve::DoSieve()
383{
384BuildPrimeTable();
385
386const unsigned int maxSieveSize = 32768;
387unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong();
388
389m_sieve.clear();
390m_sieve.resize(sieveSize, false);
391
392if (m_delta == 0)
393{
394for (unsigned int i = 0; i < primeTableSize; ++i)
395SieveSingle(m_sieve, primeTable[i], m_first, m_step, m_step.InverseMod(primeTable[i]));
396}
397else
398{
399assert(m_step%2==0);
400Integer qFirst = (m_first-m_delta) >> 1;
401Integer halfStep = m_step >> 1;
402for (unsigned int i = 0; i < primeTableSize; ++i)
403{
404word p = primeTable[i];
405word stepInv = m_step.InverseMod(p);
406SieveSingle(m_sieve, p, m_first, m_step, stepInv);
407
408word halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p;
409SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
410}
411}
412}
413
414bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector)
415{
416assert(!equiv.IsNegative() && equiv < mod);
417
418Integer gcd = GCD(equiv, mod);
419if (gcd != Integer::One())
420{
421// the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv)
422if (p <= gcd && gcd <= max && IsPrime(gcd))
423{
424p = gcd;
425return true;
426}
427else
428return false;
429}
430
431BuildPrimeTable();
432
433if (p <= primeTable[primeTableSize-1])
434{
435word *pItr;
436
437--p;
438if (p.IsPositive())
439pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong());
440else
441pItr = primeTable;
442
443while (pItr < primeTable+primeTableSize && *pItr%mod != equiv)
444++pItr;
445
446if (pItr < primeTable+primeTableSize)
447{
448p = *pItr;
449return p <= max;
450}
451
452p = primeTable[primeTableSize-1]+1;
453}
454
455assert(p > primeTable[primeTableSize-1]);
456
457if (mod.IsOdd())
458return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector);
459
460p += (equiv-p)%mod;
461
462if (p>max)
463return false;
464
465PrimeSieve sieve(p, max, mod);
466
467while (sieve.NextCandidate(p))
468{
469if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p))
470return true;
471}
472
473return false;
474}
475
476// the following two functions are based on code and comments provided by Preda Mihailescu
477static bool ProvePrime(const Integer &p, const Integer &q)
478{
479assert(p < q*q*q);
480assert(p % q == 1);
481
482// this is the Quisquater test. Numbers p having passed the Lucas - Lehmer test
483// for q and verifying p < q^3 can only be built up of two factors, both = 1 mod q,
484// or be prime. The next two lines build the discriminant of a quadratic equation
485// which holds iff p is built up of two factors (excercise ... )
486
487Integer r = (p-1)/q;
488if (((r%q).Squared()-4*(r/q)).IsSquare())
489return false;
490
491assert(primeTableSize >= 50);
492for (int i=0; i<50; i++)
493{
494Integer b = a_exp_b_mod_c(primeTable[i], r, p);
495if (b != 1)
496return a_exp_b_mod_c(b, q, p) == 1;
497}
498return false;
499}
500
501Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits)
502{
503Integer p;
504Integer minP = Integer::Power2(pbits-1);
505Integer maxP = Integer::Power2(pbits) - 1;
506
507if (maxP <= Integer(lastSmallPrime).Squared())
508{
509// Randomize() will generate a prime provable by trial division
510p.Randomize(rng, minP, maxP, Integer::PRIME);
511return p;
512}
513
514unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36);
515Integer q = MihailescuProvablePrime(rng, qbits);
516Integer q2 = q<<1;
517
518while (true)
519{
520// this initializes the sieve to search in the arithmetic
521// progression p = p_0 + \lambda * q2 = p_0 + 2 * \lambda * q,
522// with q the recursively generated prime above. We will be able
523// to use Lucas tets for proving primality. A trick of Quisquater
524// allows taking q > cubic_root(p) rather then square_root: this
525// decreases the recursion.
526
527p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2);
528PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
529
530while (sieve.NextCandidate(p))
531{
532if (FastProbablePrimeTest(p) && ProvePrime(p, q))
533return p;
534}
535}
536
537// not reached
538return p;
539}
540
541Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits)
542{
543const unsigned smallPrimeBound = 29, c_opt=10;
544Integer p;
545
546BuildPrimeTable();
547if (bits < smallPrimeBound)
548{
549do
550p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2);
551while (TrialDivision(p, 1 << ((bits+1)/2)));
552}
553else
554{
555const unsigned margin = bits > 50 ? 20 : (bits-10)/2;
556double relativeSize;
557do
558relativeSize = pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1);
559while (bits * relativeSize >= bits - margin);
560
561Integer a,b;
562Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize));
563Integer I = Integer::Power2(bits-2)/q;
564Integer I2 = I << 1;
565unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt);
566bool success = false;
567while (!success)
568{
569p.Randomize(rng, I, I2, Integer::ANY);
570p *= q; p <<= 1; ++p;
571if (!TrialDivision(p, trialDivisorBound))
572{
573a.Randomize(rng, 2, p-1, Integer::ANY);
574b = a_exp_b_mod_c(a, (p-1)/q, p);
575success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1);
576}
577}
578}
579return p;
580}
581
582Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u)
583{
584// isn't operator overloading great?
585return p * (u * (xq-xp) % q) + xp;
586}
587
588Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q)
589{
590return CRT(xp, p, xq, q, EuclideanMultiplicativeInverse(p, q));
591}
592
593Integer ModularSquareRoot(const Integer &a, const Integer &p)
594{
595if (p%4 == 3)
596return a_exp_b_mod_c(a, (p+1)/4, p);
597
598Integer q=p-1;
599unsigned int r=0;
600while (q.IsEven())
601{
602r++;
603q >>= 1;
604}
605
606Integer n=2;
607while (Jacobi(n, p) != -1)
608++n;
609
610Integer y = a_exp_b_mod_c(n, q, p);
611Integer x = a_exp_b_mod_c(a, (q-1)/2, p);
612Integer b = (x.Squared()%p)*a%p;
613x = a*x%p;
614Integer tempb, t;
615
616while (b != 1)
617{
618unsigned m=0;
619tempb = b;
620do
621{
622m++;
623b = b.Squared()%p;
624if (m==r)
625return Integer::Zero();
626}
627while (b != 1);
628
629t = y;
630for (unsigned i=0; i<r-m-1; i++)
631t = t.Squared()%p;
632y = t.Squared()%p;
633r = m;
634x = x*t%p;
635b = tempb*y%p;
636}
637
638assert(x.Squared()%p == a);
639return x;
640}
641
642bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p)
643{
644Integer D = (b.Squared() - 4*a*c) % p;
645switch (Jacobi(D, p))
646{
647default:
648assert(false);// not reached
649return false;
650case -1:
651return false;
652case 0:
653r1 = r2 = (-b*(a+a).InverseMod(p)) % p;
654assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
655return true;
656case 1:
657Integer s = ModularSquareRoot(D, p);
658Integer t = (a+a).InverseMod(p);
659r1 = (s-b)*t % p;
660r2 = (-s-b)*t % p;
661assert(((r1.Squared()*a + r1*b + c) % p).IsZero());
662assert(((r2.Squared()*a + r2*b + c) % p).IsZero());
663return true;
664}
665}
666
667Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq,
668const Integer &p, const Integer &q, const Integer &u)
669{
670Integer p2 = ModularExponentiation((a % p), dp, p);
671Integer q2 = ModularExponentiation((a % q), dq, q);
672return CRT(p2, p, q2, q, u);
673}
674
675Integer ModularRoot(const Integer &a, const Integer &e,
676const Integer &p, const Integer &q)
677{
678Integer dp = EuclideanMultiplicativeInverse(e, p-1);
679Integer dq = EuclideanMultiplicativeInverse(e, q-1);
680Integer u = EuclideanMultiplicativeInverse(p, q);
681assert(!!dp && !!dq && !!u);
682return ModularRoot(a, dp, dq, p, q, u);
683}
684
685/*
686Integer GCDI(const Integer &x, const Integer &y)
687{
688Integer a=x, b=y;
689unsigned k=0;
690
691assert(!!a && !!b);
692
693while (a[0]==0 && b[0]==0)
694{
695a >>= 1;
696b >>= 1;
697k++;
698}
699
700while (a[0]==0)
701a >>= 1;
702
703while (b[0]==0)
704b >>= 1;
705
706while (1)
707{
708switch (a.Compare(b))
709{
710case -1:
711b -= a;
712while (b[0]==0)
713b >>= 1;
714break;
715
716case 0:
717return (a <<= k);
718
719case 1:
720a -= b;
721while (a[0]==0)
722a >>= 1;
723break;
724
725default:
726assert(false);
727}
728}
729}
730
731Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b)
732{
733assert(b.Positive());
734
735if (a.Negative())
736return EuclideanMultiplicativeInverse(a%b, b);
737
738if (b[0]==0)
739{
740if (!b || a[0]==0)
741return Integer::Zero(); // no inverse
742if (a==1)
743return 1;
744Integer u = EuclideanMultiplicativeInverse(b, a);
745if (!u)
746return Integer::Zero(); // no inverse
747else
748return (b*(a-u)+1)/a;
749}
750
751Integer u=1, d=a, v1=b, v3=b, t1, t3, b2=(b+1)>>1;
752
753if (a[0])
754{
755t1 = Integer::Zero();
756t3 = -b;
757}
758else
759{
760t1 = b2;
761t3 = a>>1;
762}
763
764while (!!t3)
765{
766while (t3[0]==0)
767{
768t3 >>= 1;
769if (t1[0]==0)
770t1 >>= 1;
771else
772{
773t1 >>= 1;
774t1 += b2;
775}
776}
777if (t3.Positive())
778{
779u = t1;
780d = t3;
781}
782else
783{
784v1 = b-t1;
785v3 = -t3;
786}
787t1 = u-v1;
788t3 = d-v3;
789if (t1.Negative())
790t1 += b;
791}
792if (d==1)
793return u;
794else
795return Integer::Zero(); // no inverse
796}
797*/
798
799int Jacobi(const Integer &aIn, const Integer &bIn)
800{
801assert(bIn.IsOdd());
802
803Integer b = bIn, a = aIn%bIn;
804int result = 1;
805
806while (!!a)
807{
808unsigned i=0;
809while (a.GetBit(i)==0)
810i++;
811a>>=i;
812
813if (i%2==1 && (b%8==3 || b%8==5))
814result = -result;
815
816if (a%4==3 && b%4==3)
817result = -result;
818
819std::swap(a, b);
820a %= b;
821}
822
823return (b==1) ? result : 0;
824}
825
826Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n)
827{
828unsigned i = e.BitCount();
829if (i==0)
830return Integer::Two();
831
832MontgomeryRepresentation m(n);
833Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two());
834Integer v=p, v1=m.Subtract(m.Square(p), two);
835
836i--;
837while (i--)
838{
839if (e.GetBit(i))
840{
841// v = (v*v1 - p) % m;
842v = m.Subtract(m.Multiply(v,v1), p);
843// v1 = (v1*v1 - 2) % m;
844v1 = m.Subtract(m.Square(v1), two);
845}
846else
847{
848// v1 = (v*v1 - p) % m;
849v1 = m.Subtract(m.Multiply(v,v1), p);
850// v = (v*v - 2) % m;
851v = m.Subtract(m.Square(v), two);
852}
853}
854return m.ConvertOut(v);
855}
856
857// This is Peter Montgomery's unpublished Lucas sequence evalutation algorithm.
858// The total number of multiplies and squares used is less than the binary
859// algorithm (see above). Unfortunately I can't get it to run as fast as
860// the binary algorithm because of the extra overhead.
861/*
862Integer Lucas(const Integer &n, const Integer &P, const Integer &modulus)
863{
864if (!n)
865return 2;
866
867#define f(A, B, C)m.Subtract(m.Multiply(A, B), C)
868#define X2(A) m.Subtract(m.Square(A), two)
869#define X3(A) m.Multiply(A, m.Subtract(m.Square(A), three))
870
871MontgomeryRepresentation m(modulus);
872Integer two=m.ConvertIn(2), three=m.ConvertIn(3);
873Integer A=m.ConvertIn(P), B, C, p, d=n, e, r, t, T, U;
874
875while (d!=1)
876{
877p = d;
878unsigned int b = WORD_BITS * p.WordCount();
879Integer alpha = (Integer(5)<<(2*b-2)).SquareRoot() - Integer::Power2(b-1);
880r = (p*alpha)>>b;
881e = d-r;
882B = A;
883C = two;
884d = r;
885
886while (d!=e)
887{
888if (d<e)
889{
890swap(d, e);
891swap(A, B);
892}
893
894unsigned int dm2 = d[0], em2 = e[0];
895unsigned int dm3 = d%3, em3 = e%3;
896
897//if ((dm6+em6)%3 == 0 && d <= e + (e>>2))
898if ((dm3+em3==0 || dm3+em3==3) && (t = e, t >>= 2, t += e, d <= t))
899{
900// #1
901//t = (d+d-e)/3;
902//t = d; t += d; t -= e; t /= 3;
903//e = (e+e-d)/3;
904//e += e; e -= d; e /= 3;
905//d = t;
906
907//t = (d+e)/3
908t = d; t += e; t /= 3;
909e -= t;
910d -= t;
911
912T = f(A, B, C);
913U = f(T, A, B);
914B = f(T, B, A);
915A = U;
916continue;
917}
918
919//if (dm6 == em6 && d <= e + (e>>2))
920if (dm3 == em3 && dm2 == em2 && (t = e, t >>= 2, t += e, d <= t))
921{
922// #2
923//d = (d-e)>>1;
924d -= e; d >>= 1;
925B = f(A, B, C);
926A = X2(A);
927continue;
928}
929
930//if (d <= (e<<2))
931if (d <= (t = e, t <<= 2))
932{
933// #3
934d -= e;
935C = f(A, B, C);
936swap(B, C);
937continue;
938}
939
940if (dm2 == em2)
941{
942// #4
943//d = (d-e)>>1;
944d -= e; d >>= 1;
945B = f(A, B, C);
946A = X2(A);
947continue;
948}
949
950if (dm2 == 0)
951{
952// #5
953d >>= 1;
954C = f(A, C, B);
955A = X2(A);
956continue;
957}
958
959if (dm3 == 0)
960{
961// #6
962//d = d/3 - e;
963d /= 3; d -= e;
964T = X2(A);
965C = f(T, f(A, B, C), C);
966swap(B, C);
967A = f(T, A, A);
968continue;
969}
970
971if (dm3+em3==0 || dm3+em3==3)
972{
973// #7
974//d = (d-e-e)/3;
975d -= e; d -= e; d /= 3;
976T = f(A, B, C);
977B = f(T, A, B);
978A = X3(A);
979continue;
980}
981
982if (dm3 == em3)
983{
984// #8
985//d = (d-e)/3;
986d -= e; d /= 3;
987T = f(A, B, C);
988C = f(A, C, B);
989B = T;
990A = X3(A);
991continue;
992}
993
994assert(em2 == 0);
995// #9
996e >>= 1;
997C = f(C, B, A);
998B = X2(B);
999}
1000
1001A = f(A, B, C);
1002}
1003
1004#undef f
1005#undef X2
1006#undef X3
1007
1008return m.ConvertOut(A);
1009}
1010*/
1011
1012Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u)
1013{
1014Integer d = (m*m-4);
1015Integer p2 = p-Jacobi(d,p);
1016Integer q2 = q-Jacobi(d,q);
1017return CRT(Lucas(EuclideanMultiplicativeInverse(e,p2), m, p), p, Lucas(EuclideanMultiplicativeInverse(e,q2), m, q), q, u);
1018}
1019
1020Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q)
1021{
1022return InverseLucas(e, m, p, q, EuclideanMultiplicativeInverse(p, q));
1023}
1024
1025unsigned int FactoringWorkFactor(unsigned int n)
1026{
1027// extrapolated from the table in Odlyzko's "The Future of Integer Factorization"
1028// updated to reflect the factoring of RSA-130
1029if (n<5) return 0;
1030else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
1031}
1032
1033unsigned int DiscreteLogWorkFactor(unsigned int n)
1034{
1035// assuming discrete log takes about the same time as factoring
1036if (n<5) return 0;
1037else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5);
1038}
1039
1040// ********************************************************
1041
1042void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits)
1043{
1044// no prime exists for delta = -1, qbits = 4, and pbits = 5
1045assert(qbits > 4);
1046assert(pbits > qbits);
1047
1048if (qbits+1 == pbits)
1049{
1050Integer minP = Integer::Power2(pbits-1);
1051Integer maxP = Integer::Power2(pbits) - 1;
1052bool success = false;
1053
1054while (!success)
1055{
1056p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12);
1057PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta);
1058
1059while (sieve.NextCandidate(p))
1060{
1061assert(IsSmallPrime(p) || SmallDivisorsTest(p));
1062q = (p-delta) >> 1;
1063assert(IsSmallPrime(q) || SmallDivisorsTest(q));
1064if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p))
1065{
1066success = true;
1067break;
1068}
1069}
1070}
1071
1072if (delta == 1)
1073{
1074// find g such that g is a quadratic residue mod p, then g has order q
1075// g=4 always works, but this way we get the smallest quadratic residue (other than 1)
1076for (g=2; Jacobi(g, p) != 1; ++g) {}
1077// contributed by Walt Tuvell: g should be the following according to the Law of Quadratic Reciprocity
1078assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4);
1079}
1080else
1081{
1082assert(delta == -1);
1083// find g such that g*g-4 is a quadratic non-residue,
1084// and such that g has order q
1085for (g=3; ; ++g)
1086if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2)
1087break;
1088}
1089}
1090else
1091{
1092Integer minQ = Integer::Power2(qbits-1);
1093Integer maxQ = Integer::Power2(qbits) - 1;
1094Integer minP = Integer::Power2(pbits-1);
1095Integer maxP = Integer::Power2(pbits) - 1;
1096
1097do
1098{
1099q.Randomize(rng, minQ, maxQ, Integer::PRIME);
1100} while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q));
1101
1102// find a random g of order q
1103if (delta==1)
1104{
1105do
1106{
1107Integer h(rng, 2, p-2, Integer::ANY);
1108g = a_exp_b_mod_c(h, (p-1)/q, p);
1109} while (g <= 1);
1110assert(a_exp_b_mod_c(g, q, p)==1);
1111}
1112else
1113{
1114assert(delta==-1);
1115do
1116{
1117Integer h(rng, 3, p-1, Integer::ANY);
1118if (Jacobi(h*h-4, p)==1)
1119continue;
1120g = Lucas((p+1)/q, h, p);
1121} while (g <= 2);
1122assert(Lucas(q, g, p) == 2);
1123}
1124}
1125}
1126
1127NAMESPACE_END

Archive Download this file

Branches

Tags

Quick Links:     www.monotone.ca    -     Downloads    -     Documentation    -     Wiki    -     Code Forge    -     Build Status