source: git/src-cryptopp/nbtheory.cpp

Last change on this file was e230cb0, checked in by David Stainton <dstainton415@…>, at 2016-10-12T13:27:29Z

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