From 5494fc310acf0aabb9d828451331e44483eb21c7 Mon Sep 17 00:00:00 2001 From: Malfurious Date: Mon, 21 Oct 2024 11:09:00 -0400 Subject: Remove Crypto++ library The tracked version of Crypto++ is going on 10 years old and doesn't always compile properly on modern tooling. This removes the entire subdirectory as well as references to files in the build script. Due to the number of files touched by this commit, I opt to add its replacement in the next commit. Signed-off-by: Malfurious --- cryptopp562/nbtheory.cpp | 1123 ---------------------------------------------- 1 file changed, 1123 deletions(-) delete mode 100644 cryptopp562/nbtheory.cpp (limited to 'cryptopp562/nbtheory.cpp') diff --git a/cryptopp562/nbtheory.cpp b/cryptopp562/nbtheory.cpp deleted file mode 100644 index 3fdea4e..0000000 --- a/cryptopp562/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 -#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