From 539364846a89987ac2679988653f50332cb91d26 Mon Sep 17 00:00:00 2001 From: "madmaxoft@gmail.com" Date: Thu, 30 Aug 2012 21:06:13 +0000 Subject: Implemented 1.3.2 protocol encryption using CryptoPP, up to Client Status packet (http://wiki.vg/Protocol_FAQ step 14) git-svn-id: http://mc-server.googlecode.com/svn/trunk@808 0a769ca7-a7f5-676a-18bf-c427514a06d6 --- CryptoPP/nbtheory.cpp | 1123 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1123 insertions(+) create mode 100644 CryptoPP/nbtheory.cpp (limited to 'CryptoPP/nbtheory.cpp') diff --git a/CryptoPP/nbtheory.cpp b/CryptoPP/nbtheory.cpp new file mode 100644 index 000000000..3fdea4e69 --- /dev/null +++ b/CryptoPP/nbtheory.cpp @@ -0,0 +1,1123 @@ +// 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 +#include + +#ifdef _OPENMP +// needed in MSVC 2005 to generate correct manifest +#include +#endif + +NAMESPACE_BEGIN(CryptoPP) + +const word s_lastSmallPrime = 32719; + +struct NewPrimeTable +{ + std::vector * operator()() const + { + const unsigned int maxPrimeTableSize = 3511; + + std::auto_ptr > pPrimeTable(new std::vector); + std::vector &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 &primeTable = Singleton, 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]3 && b>1 && b3 && b>1 && b>a; + + Integer z = a_exp_b_mod_c(b, m, n); + if (z==1 || z==nminus1) + return true; + for (unsigned j=1; j3); + + Integer b; + for (unsigned int i=0; i2); + + 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().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 &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 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 &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>= 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>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 -- cgit v1.2.3