diff options
Diffstat (limited to 'embeddedcryptopp/nbtheory.cpp')
| -rw-r--r-- | embeddedcryptopp/nbtheory.cpp | 1123 |
1 files changed, 0 insertions, 1123 deletions
diff --git a/embeddedcryptopp/nbtheory.cpp b/embeddedcryptopp/nbtheory.cpp deleted file mode 100644 index 3fdea4e..0000000 --- a/embeddedcryptopp/nbtheory.cpp +++ /dev/null @@ -1,1123 +0,0 @@ -// nbtheory.cpp - written and placed in the public domain by Wei Dai - -#include "pch.h" - -#ifndef CRYPTOPP_IMPORTS - -#include "nbtheory.h" -#include "modarith.h" -#include "algparam.h" - -#include <math.h> -#include <vector> - -#ifdef _OPENMP -// needed in MSVC 2005 to generate correct manifest -#include <omp.h> -#endif - -NAMESPACE_BEGIN(CryptoPP) - -const word s_lastSmallPrime = 32719; - -struct NewPrimeTable -{ - std::vector<word16> * operator()() const - { - const unsigned int maxPrimeTableSize = 3511; - - std::auto_ptr<std::vector<word16> > pPrimeTable(new std::vector<word16>); - std::vector<word16> &primeTable = *pPrimeTable; - primeTable.reserve(maxPrimeTableSize); - - primeTable.push_back(2); - unsigned int testEntriesEnd = 1; - - for (unsigned int p=3; p<=s_lastSmallPrime; p+=2) - { - unsigned int j; - for (j=1; j<testEntriesEnd; j++) - if (p%primeTable[j] == 0) - break; - if (j == testEntriesEnd) - { - primeTable.push_back(p); - testEntriesEnd = UnsignedMin(54U, primeTable.size()); - } - } - - return pPrimeTable.release(); - } -}; - -const word16 * GetPrimeTable(unsigned int &size) -{ - const std::vector<word16> &primeTable = Singleton<std::vector<word16>, NewPrimeTable>().Ref(); - size = (unsigned int)primeTable.size(); - return &primeTable[0]; -} - -bool IsSmallPrime(const Integer &p) -{ - unsigned int primeTableSize; - const word16 * primeTable = GetPrimeTable(primeTableSize); - - if (p.IsPositive() && p <= primeTable[primeTableSize-1]) - return std::binary_search(primeTable, primeTable+primeTableSize, (word16)p.ConvertToLong()); - else - return false; -} - -bool TrialDivision(const Integer &p, unsigned bound) -{ - unsigned int primeTableSize; - const word16 * primeTable = GetPrimeTable(primeTableSize); - - assert(primeTable[primeTableSize-1] >= bound); - - unsigned int i; - for (i = 0; primeTable[i]<bound; i++) - if ((p % primeTable[i]) == 0) - return true; - - if (bound == primeTable[i]) - return (p % bound == 0); - else - return false; -} - -bool SmallDivisorsTest(const Integer &p) -{ - unsigned int primeTableSize; - const word16 * primeTable = GetPrimeTable(primeTableSize); - return !TrialDivision(p, primeTable[primeTableSize-1]); -} - -bool IsFermatProbablePrime(const Integer &n, const Integer &b) -{ - if (n <= 3) - return n==2 || n==3; - - assert(n>3 && b>1 && b<n-1); - return a_exp_b_mod_c(b, n-1, n)==1; -} - -bool IsStrongProbablePrime(const Integer &n, const Integer &b) -{ - if (n <= 3) - return n==2 || n==3; - - assert(n>3 && b>1 && b<n-1); - - if ((n.IsEven() && n!=2) || GCD(b, n) != 1) - return false; - - Integer nminus1 = (n-1); - unsigned int a; - - // calculate a = largest power of 2 that divides (n-1) - for (a=0; ; a++) - if (nminus1.GetBit(a)) - break; - Integer m = nminus1>>a; - - Integer z = a_exp_b_mod_c(b, m, n); - if (z==1 || z==nminus1) - return true; - for (unsigned j=1; j<a; j++) - { - z = z.Squared()%n; - if (z==nminus1) - return true; - if (z==1) - return false; - } - return false; -} - -bool RabinMillerTest(RandomNumberGenerator &rng, const Integer &n, unsigned int rounds) -{ - if (n <= 3) - return n==2 || n==3; - - assert(n>3); - - Integer b; - for (unsigned int i=0; i<rounds; i++) - { - b.Randomize(rng, 2, n-2); - if (!IsStrongProbablePrime(n, b)) - return false; - } - return true; -} - -bool IsLucasProbablePrime(const Integer &n) -{ - if (n <= 1) - return false; - - if (n.IsEven()) - return n==2; - - assert(n>2); - - Integer b=3; - unsigned int i=0; - int j; - - while ((j=Jacobi(b.Squared()-4, n)) == 1) - { - if (++i==64 && n.IsSquare()) // avoid infinite loop if n is a square - return false; - ++b; ++b; - } - - if (j==0) - return false; - else - return Lucas(n+1, b, n)==2; -} - -bool IsStrongLucasProbablePrime(const Integer &n) -{ - if (n <= 1) - return false; - - if (n.IsEven()) - return n==2; - - assert(n>2); - - Integer b=3; - unsigned int i=0; - int j; - - while ((j=Jacobi(b.Squared()-4, n)) == 1) - { - if (++i==64 && n.IsSquare()) // avoid infinite loop if n is a square - return false; - ++b; ++b; - } - - if (j==0) - return false; - - Integer n1 = n+1; - unsigned int a; - - // calculate a = largest power of 2 that divides n1 - for (a=0; ; a++) - if (n1.GetBit(a)) - break; - Integer m = n1>>a; - - Integer z = Lucas(m, b, n); - if (z==2 || z==n-2) - return true; - for (i=1; i<a; i++) - { - z = (z.Squared()-2)%n; - if (z==n-2) - return true; - if (z==2) - return false; - } - return false; -} - -struct NewLastSmallPrimeSquared -{ - Integer * operator()() const - { - return new Integer(Integer(s_lastSmallPrime).Squared()); - } -}; - -bool IsPrime(const Integer &p) -{ - if (p <= s_lastSmallPrime) - return IsSmallPrime(p); - else if (p <= Singleton<Integer, NewLastSmallPrimeSquared>().Ref()) - return SmallDivisorsTest(p); - else - return SmallDivisorsTest(p) && IsStrongProbablePrime(p, 3) && IsStrongLucasProbablePrime(p); -} - -bool VerifyPrime(RandomNumberGenerator &rng, const Integer &p, unsigned int level) -{ - bool pass = IsPrime(p) && RabinMillerTest(rng, p, 1); - if (level >= 1) - pass = pass && RabinMillerTest(rng, p, 10); - return pass; -} - -unsigned int PrimeSearchInterval(const Integer &max) -{ - return max.BitCount(); -} - -static inline bool FastProbablePrimeTest(const Integer &n) -{ - return IsStrongProbablePrime(n,2); -} - -AlgorithmParameters MakeParametersForTwoPrimesOfEqualSize(unsigned int productBitLength) -{ - if (productBitLength < 16) - throw InvalidArgument("invalid bit length"); - - Integer minP, maxP; - - if (productBitLength%2==0) - { - minP = Integer(182) << (productBitLength/2-8); - maxP = Integer::Power2(productBitLength/2)-1; - } - else - { - minP = Integer::Power2((productBitLength-1)/2); - maxP = Integer(181) << ((productBitLength+1)/2-8); - } - - return MakeParameters("RandomNumberType", Integer::PRIME)("Min", minP)("Max", maxP); -} - -class PrimeSieve -{ -public: - // delta == 1 or -1 means double sieve with p = 2*q + delta - PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta=0); - bool NextCandidate(Integer &c); - - void DoSieve(); - static void SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv); - - Integer m_first, m_last, m_step; - signed int m_delta; - word m_next; - std::vector<bool> m_sieve; -}; - -PrimeSieve::PrimeSieve(const Integer &first, const Integer &last, const Integer &step, signed int delta) - : m_first(first), m_last(last), m_step(step), m_delta(delta), m_next(0) -{ - DoSieve(); -} - -bool PrimeSieve::NextCandidate(Integer &c) -{ - bool safe = SafeConvert(std::find(m_sieve.begin()+m_next, m_sieve.end(), false) - m_sieve.begin(), m_next); - assert(safe); - if (m_next == m_sieve.size()) - { - m_first += long(m_sieve.size())*m_step; - if (m_first > m_last) - return false; - else - { - m_next = 0; - DoSieve(); - return NextCandidate(c); - } - } - else - { - c = m_first + long(m_next)*m_step; - ++m_next; - return true; - } -} - -void PrimeSieve::SieveSingle(std::vector<bool> &sieve, word16 p, const Integer &first, const Integer &step, word16 stepInv) -{ - if (stepInv) - { - size_t sieveSize = sieve.size(); - size_t j = (word32(p-(first%p))*stepInv) % p; - // if the first multiple of p is p, skip it - if (first.WordCount() <= 1 && first + step*long(j) == p) - j += p; - for (; j < sieveSize; j += p) - sieve[j] = true; - } -} - -void PrimeSieve::DoSieve() -{ - unsigned int primeTableSize; - const word16 * primeTable = GetPrimeTable(primeTableSize); - - const unsigned int maxSieveSize = 32768; - unsigned int sieveSize = STDMIN(Integer(maxSieveSize), (m_last-m_first)/m_step+1).ConvertToLong(); - - m_sieve.clear(); - m_sieve.resize(sieveSize, false); - - if (m_delta == 0) - { - for (unsigned int i = 0; i < primeTableSize; ++i) - SieveSingle(m_sieve, primeTable[i], m_first, m_step, (word16)m_step.InverseMod(primeTable[i])); - } - else - { - assert(m_step%2==0); - Integer qFirst = (m_first-m_delta) >> 1; - Integer halfStep = m_step >> 1; - for (unsigned int i = 0; i < primeTableSize; ++i) - { - word16 p = primeTable[i]; - word16 stepInv = (word16)m_step.InverseMod(p); - SieveSingle(m_sieve, p, m_first, m_step, stepInv); - - word16 halfStepInv = 2*stepInv < p ? 2*stepInv : 2*stepInv-p; - SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv); - } - } -} - -bool FirstPrime(Integer &p, const Integer &max, const Integer &equiv, const Integer &mod, const PrimeSelector *pSelector) -{ - assert(!equiv.IsNegative() && equiv < mod); - - Integer gcd = GCD(equiv, mod); - if (gcd != Integer::One()) - { - // the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv) - if (p <= gcd && gcd <= max && IsPrime(gcd) && (!pSelector || pSelector->IsAcceptable(gcd))) - { - p = gcd; - return true; - } - else - return false; - } - - unsigned int primeTableSize; - const word16 * primeTable = GetPrimeTable(primeTableSize); - - if (p <= primeTable[primeTableSize-1]) - { - const word16 *pItr; - - --p; - if (p.IsPositive()) - pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.ConvertToLong()); - else - pItr = primeTable; - - while (pItr < primeTable+primeTableSize && !(*pItr%mod == equiv && (!pSelector || pSelector->IsAcceptable(*pItr)))) - ++pItr; - - if (pItr < primeTable+primeTableSize) - { - p = *pItr; - return p <= max; - } - - p = primeTable[primeTableSize-1]+1; - } - - assert(p > primeTable[primeTableSize-1]); - - if (mod.IsOdd()) - return FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod<<1, pSelector); - - p += (equiv-p)%mod; - - if (p>max) - return false; - - PrimeSieve sieve(p, max, mod); - - while (sieve.NextCandidate(p)) - { - if ((!pSelector || pSelector->IsAcceptable(p)) && FastProbablePrimeTest(p) && IsPrime(p)) - return true; - } - - return false; -} - -// the following two functions are based on code and comments provided by Preda Mihailescu -static bool ProvePrime(const Integer &p, const Integer &q) -{ - assert(p < q*q*q); - assert(p % q == 1); - -// this is the Quisquater test. Numbers p having passed the Lucas - Lehmer test -// for q and verifying p < q^3 can only be built up of two factors, both = 1 mod q, -// or be prime. The next two lines build the discriminant of a quadratic equation -// which holds iff p is built up of two factors (excercise ... ) - - Integer r = (p-1)/q; - if (((r%q).Squared()-4*(r/q)).IsSquare()) - return false; - - unsigned int primeTableSize; - const word16 * primeTable = GetPrimeTable(primeTableSize); - - assert(primeTableSize >= 50); - for (int i=0; i<50; i++) - { - Integer b = a_exp_b_mod_c(primeTable[i], r, p); - if (b != 1) - return a_exp_b_mod_c(b, q, p) == 1; - } - return false; -} - -Integer MihailescuProvablePrime(RandomNumberGenerator &rng, unsigned int pbits) -{ - Integer p; - Integer minP = Integer::Power2(pbits-1); - Integer maxP = Integer::Power2(pbits) - 1; - - if (maxP <= Integer(s_lastSmallPrime).Squared()) - { - // Randomize() will generate a prime provable by trial division - p.Randomize(rng, minP, maxP, Integer::PRIME); - return p; - } - - unsigned int qbits = (pbits+2)/3 + 1 + rng.GenerateWord32(0, pbits/36); - Integer q = MihailescuProvablePrime(rng, qbits); - Integer q2 = q<<1; - - while (true) - { - // this initializes the sieve to search in the arithmetic - // progression p = p_0 + \lambda * q2 = p_0 + 2 * \lambda * q, - // with q the recursively generated prime above. We will be able - // to use Lucas tets for proving primality. A trick of Quisquater - // allows taking q > cubic_root(p) rather then square_root: this - // decreases the recursion. - - p.Randomize(rng, minP, maxP, Integer::ANY, 1, q2); - PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2); - - while (sieve.NextCandidate(p)) - { - if (FastProbablePrimeTest(p) && ProvePrime(p, q)) - return p; - } - } - - // not reached - return p; -} - -Integer MaurerProvablePrime(RandomNumberGenerator &rng, unsigned int bits) -{ - const unsigned smallPrimeBound = 29, c_opt=10; - Integer p; - - unsigned int primeTableSize; - const word16 * primeTable = GetPrimeTable(primeTableSize); - - if (bits < smallPrimeBound) - { - do - p.Randomize(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2); - while (TrialDivision(p, 1 << ((bits+1)/2))); - } - else - { - const unsigned margin = bits > 50 ? 20 : (bits-10)/2; - double relativeSize; - do - relativeSize = pow(2.0, double(rng.GenerateWord32())/0xffffffff - 1); - while (bits * relativeSize >= bits - margin); - - Integer a,b; - Integer q = MaurerProvablePrime(rng, unsigned(bits*relativeSize)); - Integer I = Integer::Power2(bits-2)/q; - Integer I2 = I << 1; - unsigned int trialDivisorBound = (unsigned int)STDMIN((unsigned long)primeTable[primeTableSize-1], (unsigned long)bits*bits/c_opt); - bool success = false; - while (!success) - { - p.Randomize(rng, I, I2, Integer::ANY); - p *= q; p <<= 1; ++p; - if (!TrialDivision(p, trialDivisorBound)) - { - a.Randomize(rng, 2, p-1, Integer::ANY); - b = a_exp_b_mod_c(a, (p-1)/q, p); - success = (GCD(b-1, p) == 1) && (a_exp_b_mod_c(b, q, p) == 1); - } - } - } - return p; -} - -Integer CRT(const Integer &xp, const Integer &p, const Integer &xq, const Integer &q, const Integer &u) -{ - // isn't operator overloading great? - return p * (u * (xq-xp) % q) + xp; -/* - Integer t1 = xq-xp; - cout << hex << t1 << endl; - Integer t2 = u * t1; - cout << hex << t2 << endl; - Integer t3 = t2 % q; - cout << hex << t3 << endl; - Integer t4 = p * t3; - cout << hex << t4 << endl; - Integer t5 = t4 + xp; - cout << hex << t5 << endl; - return t5; -*/ -} - -Integer ModularSquareRoot(const Integer &a, const Integer &p) -{ - if (p%4 == 3) - return a_exp_b_mod_c(a, (p+1)/4, p); - - Integer q=p-1; - unsigned int r=0; - while (q.IsEven()) - { - r++; - q >>= 1; - } - - Integer n=2; - while (Jacobi(n, p) != -1) - ++n; - - Integer y = a_exp_b_mod_c(n, q, p); - Integer x = a_exp_b_mod_c(a, (q-1)/2, p); - Integer b = (x.Squared()%p)*a%p; - x = a*x%p; - Integer tempb, t; - - while (b != 1) - { - unsigned m=0; - tempb = b; - do - { - m++; - b = b.Squared()%p; - if (m==r) - return Integer::Zero(); - } - while (b != 1); - - t = y; - for (unsigned i=0; i<r-m-1; i++) - t = t.Squared()%p; - y = t.Squared()%p; - r = m; - x = x*t%p; - b = tempb*y%p; - } - - assert(x.Squared()%p == a); - return x; -} - -bool SolveModularQuadraticEquation(Integer &r1, Integer &r2, const Integer &a, const Integer &b, const Integer &c, const Integer &p) -{ - Integer D = (b.Squared() - 4*a*c) % p; - switch (Jacobi(D, p)) - { - default: - assert(false); // not reached - return false; - case -1: - return false; - case 0: - r1 = r2 = (-b*(a+a).InverseMod(p)) % p; - assert(((r1.Squared()*a + r1*b + c) % p).IsZero()); - return true; - case 1: - Integer s = ModularSquareRoot(D, p); - Integer t = (a+a).InverseMod(p); - r1 = (s-b)*t % p; - r2 = (-s-b)*t % p; - assert(((r1.Squared()*a + r1*b + c) % p).IsZero()); - assert(((r2.Squared()*a + r2*b + c) % p).IsZero()); - return true; - } -} - -Integer ModularRoot(const Integer &a, const Integer &dp, const Integer &dq, - const Integer &p, const Integer &q, const Integer &u) -{ - Integer p2, q2; - #pragma omp parallel - #pragma omp sections - { - #pragma omp section - p2 = ModularExponentiation((a % p), dp, p); - #pragma omp section - q2 = ModularExponentiation((a % q), dq, q); - } - return CRT(p2, p, q2, q, u); -} - -Integer ModularRoot(const Integer &a, const Integer &e, - const Integer &p, const Integer &q) -{ - Integer dp = EuclideanMultiplicativeInverse(e, p-1); - Integer dq = EuclideanMultiplicativeInverse(e, q-1); - Integer u = EuclideanMultiplicativeInverse(p, q); - assert(!!dp && !!dq && !!u); - return ModularRoot(a, dp, dq, p, q, u); -} - -/* -Integer GCDI(const Integer &x, const Integer &y) -{ - Integer a=x, b=y; - unsigned k=0; - - assert(!!a && !!b); - - while (a[0]==0 && b[0]==0) - { - a >>= 1; - b >>= 1; - k++; - } - - while (a[0]==0) - a >>= 1; - - while (b[0]==0) - b >>= 1; - - while (1) - { - switch (a.Compare(b)) - { - case -1: - b -= a; - while (b[0]==0) - b >>= 1; - break; - - case 0: - return (a <<= k); - - case 1: - a -= b; - while (a[0]==0) - a >>= 1; - break; - - default: - assert(false); - } - } -} - -Integer EuclideanMultiplicativeInverse(const Integer &a, const Integer &b) -{ - assert(b.Positive()); - - if (a.Negative()) - return EuclideanMultiplicativeInverse(a%b, b); - - if (b[0]==0) - { - if (!b || a[0]==0) - return Integer::Zero(); // no inverse - if (a==1) - return 1; - Integer u = EuclideanMultiplicativeInverse(b, a); - if (!u) - return Integer::Zero(); // no inverse - else - return (b*(a-u)+1)/a; - } - - Integer u=1, d=a, v1=b, v3=b, t1, t3, b2=(b+1)>>1; - - if (a[0]) - { - t1 = Integer::Zero(); - t3 = -b; - } - else - { - t1 = b2; - t3 = a>>1; - } - - while (!!t3) - { - while (t3[0]==0) - { - t3 >>= 1; - if (t1[0]==0) - t1 >>= 1; - else - { - t1 >>= 1; - t1 += b2; - } - } - if (t3.Positive()) - { - u = t1; - d = t3; - } - else - { - v1 = b-t1; - v3 = -t3; - } - t1 = u-v1; - t3 = d-v3; - if (t1.Negative()) - t1 += b; - } - if (d==1) - return u; - else - return Integer::Zero(); // no inverse -} -*/ - -int Jacobi(const Integer &aIn, const Integer &bIn) -{ - assert(bIn.IsOdd()); - - Integer b = bIn, a = aIn%bIn; - int result = 1; - - while (!!a) - { - unsigned i=0; - while (a.GetBit(i)==0) - i++; - a>>=i; - - if (i%2==1 && (b%8==3 || b%8==5)) - result = -result; - - if (a%4==3 && b%4==3) - result = -result; - - std::swap(a, b); - a %= b; - } - - return (b==1) ? result : 0; -} - -Integer Lucas(const Integer &e, const Integer &pIn, const Integer &n) -{ - unsigned i = e.BitCount(); - if (i==0) - return Integer::Two(); - - MontgomeryRepresentation m(n); - Integer p=m.ConvertIn(pIn%n), two=m.ConvertIn(Integer::Two()); - Integer v=p, v1=m.Subtract(m.Square(p), two); - - i--; - while (i--) - { - if (e.GetBit(i)) - { - // v = (v*v1 - p) % m; - v = m.Subtract(m.Multiply(v,v1), p); - // v1 = (v1*v1 - 2) % m; - v1 = m.Subtract(m.Square(v1), two); - } - else - { - // v1 = (v*v1 - p) % m; - v1 = m.Subtract(m.Multiply(v,v1), p); - // v = (v*v - 2) % m; - v = m.Subtract(m.Square(v), two); - } - } - return m.ConvertOut(v); -} - -// This is Peter Montgomery's unpublished Lucas sequence evalutation algorithm. -// The total number of multiplies and squares used is less than the binary -// algorithm (see above). Unfortunately I can't get it to run as fast as -// the binary algorithm because of the extra overhead. -/* -Integer Lucas(const Integer &n, const Integer &P, const Integer &modulus) -{ - if (!n) - return 2; - -#define f(A, B, C) m.Subtract(m.Multiply(A, B), C) -#define X2(A) m.Subtract(m.Square(A), two) -#define X3(A) m.Multiply(A, m.Subtract(m.Square(A), three)) - - MontgomeryRepresentation m(modulus); - Integer two=m.ConvertIn(2), three=m.ConvertIn(3); - Integer A=m.ConvertIn(P), B, C, p, d=n, e, r, t, T, U; - - while (d!=1) - { - p = d; - unsigned int b = WORD_BITS * p.WordCount(); - Integer alpha = (Integer(5)<<(2*b-2)).SquareRoot() - Integer::Power2(b-1); - r = (p*alpha)>>b; - e = d-r; - B = A; - C = two; - d = r; - - while (d!=e) - { - if (d<e) - { - swap(d, e); - swap(A, B); - } - - unsigned int dm2 = d[0], em2 = e[0]; - unsigned int dm3 = d%3, em3 = e%3; - -// if ((dm6+em6)%3 == 0 && d <= e + (e>>2)) - if ((dm3+em3==0 || dm3+em3==3) && (t = e, t >>= 2, t += e, d <= t)) - { - // #1 -// t = (d+d-e)/3; -// t = d; t += d; t -= e; t /= 3; -// e = (e+e-d)/3; -// e += e; e -= d; e /= 3; -// d = t; - -// t = (d+e)/3 - t = d; t += e; t /= 3; - e -= t; - d -= t; - - T = f(A, B, C); - U = f(T, A, B); - B = f(T, B, A); - A = U; - continue; - } - -// if (dm6 == em6 && d <= e + (e>>2)) - if (dm3 == em3 && dm2 == em2 && (t = e, t >>= 2, t += e, d <= t)) - { - // #2 -// d = (d-e)>>1; - d -= e; d >>= 1; - B = f(A, B, C); - A = X2(A); - continue; - } - -// if (d <= (e<<2)) - if (d <= (t = e, t <<= 2)) - { - // #3 - d -= e; - C = f(A, B, C); - swap(B, C); - continue; - } - - if (dm2 == em2) - { - // #4 -// d = (d-e)>>1; - d -= e; d >>= 1; - B = f(A, B, C); - A = X2(A); - continue; - } - - if (dm2 == 0) - { - // #5 - d >>= 1; - C = f(A, C, B); - A = X2(A); - continue; - } - - if (dm3 == 0) - { - // #6 -// d = d/3 - e; - d /= 3; d -= e; - T = X2(A); - C = f(T, f(A, B, C), C); - swap(B, C); - A = f(T, A, A); - continue; - } - - if (dm3+em3==0 || dm3+em3==3) - { - // #7 -// d = (d-e-e)/3; - d -= e; d -= e; d /= 3; - T = f(A, B, C); - B = f(T, A, B); - A = X3(A); - continue; - } - - if (dm3 == em3) - { - // #8 -// d = (d-e)/3; - d -= e; d /= 3; - T = f(A, B, C); - C = f(A, C, B); - B = T; - A = X3(A); - continue; - } - - assert(em2 == 0); - // #9 - e >>= 1; - C = f(C, B, A); - B = X2(B); - } - - A = f(A, B, C); - } - -#undef f -#undef X2 -#undef X3 - - return m.ConvertOut(A); -} -*/ - -Integer InverseLucas(const Integer &e, const Integer &m, const Integer &p, const Integer &q, const Integer &u) -{ - Integer d = (m*m-4); - Integer p2, q2; - #pragma omp parallel - #pragma omp sections - { - #pragma omp section - { - p2 = p-Jacobi(d,p); - p2 = Lucas(EuclideanMultiplicativeInverse(e,p2), m, p); - } - #pragma omp section - { - q2 = q-Jacobi(d,q); - q2 = Lucas(EuclideanMultiplicativeInverse(e,q2), m, q); - } - } - return CRT(p2, p, q2, q, u); -} - -unsigned int FactoringWorkFactor(unsigned int n) -{ - // extrapolated from the table in Odlyzko's "The Future of Integer Factorization" - // updated to reflect the factoring of RSA-130 - if (n<5) return 0; - else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5); -} - -unsigned int DiscreteLogWorkFactor(unsigned int n) -{ - // assuming discrete log takes about the same time as factoring - if (n<5) return 0; - else return (unsigned int)(2.4 * pow((double)n, 1.0/3.0) * pow(log(double(n)), 2.0/3.0) - 5); -} - -// ******************************************************** - -void PrimeAndGenerator::Generate(signed int delta, RandomNumberGenerator &rng, unsigned int pbits, unsigned int qbits) -{ - // no prime exists for delta = -1, qbits = 4, and pbits = 5 - assert(qbits > 4); - assert(pbits > qbits); - - if (qbits+1 == pbits) - { - Integer minP = Integer::Power2(pbits-1); - Integer maxP = Integer::Power2(pbits) - 1; - bool success = false; - - while (!success) - { - p.Randomize(rng, minP, maxP, Integer::ANY, 6+5*delta, 12); - PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*12, maxP), 12, delta); - - while (sieve.NextCandidate(p)) - { - assert(IsSmallPrime(p) || SmallDivisorsTest(p)); - q = (p-delta) >> 1; - assert(IsSmallPrime(q) || SmallDivisorsTest(q)); - if (FastProbablePrimeTest(q) && FastProbablePrimeTest(p) && IsPrime(q) && IsPrime(p)) - { - success = true; - break; - } - } - } - - if (delta == 1) - { - // find g such that g is a quadratic residue mod p, then g has order q - // g=4 always works, but this way we get the smallest quadratic residue (other than 1) - for (g=2; Jacobi(g, p) != 1; ++g) {} - // contributed by Walt Tuvell: g should be the following according to the Law of Quadratic Reciprocity - assert((p%8==1 || p%8==7) ? g==2 : (p%12==1 || p%12==11) ? g==3 : g==4); - } - else - { - assert(delta == -1); - // find g such that g*g-4 is a quadratic non-residue, - // and such that g has order q - for (g=3; ; ++g) - if (Jacobi(g*g-4, p)==-1 && Lucas(q, g, p)==2) - break; - } - } - else - { - Integer minQ = Integer::Power2(qbits-1); - Integer maxQ = Integer::Power2(qbits) - 1; - Integer minP = Integer::Power2(pbits-1); - Integer maxP = Integer::Power2(pbits) - 1; - - do - { - q.Randomize(rng, minQ, maxQ, Integer::PRIME); - } while (!p.Randomize(rng, minP, maxP, Integer::PRIME, delta%q, q)); - - // find a random g of order q - if (delta==1) - { - do - { - Integer h(rng, 2, p-2, Integer::ANY); - g = a_exp_b_mod_c(h, (p-1)/q, p); - } while (g <= 1); - assert(a_exp_b_mod_c(g, q, p)==1); - } - else - { - assert(delta==-1); - do - { - Integer h(rng, 3, p-1, Integer::ANY); - if (Jacobi(h*h-4, p)==1) - continue; - g = Lucas((p+1)/q, h, p); - } while (g <= 2); - assert(Lucas(q, g, p) == 2); - } - } -} - -NAMESPACE_END - -#endif |