integer.cpp

来自「各种加密算法的集合」· C++ 代码 · 共 2,363 行 · 第 1/4 页

			SetWords(T0, 0, N); 
			carry = 0; 
		} 
 
		RecursiveMultiply(T2, R0, A1, B1, N2); 
 
		// now T[01] holds (A1-A0)*(B0-B1), T[23] holds A1*B1 
 
		CopyWords(R0, L+N2, N2); 
		word c2 = Subtract(R0, R0, L, N2); 
		c2 += Subtract(R0, R0, T0, N2); 
		word t = (Compare(R0, T2, N2) == -1); 
 
		carry += t; 
		carry += Increment(R0, N2, c2+t); 
		carry += Add(R0, R0, T1, N2); 
		carry += Add(R0, R0, T3, N2); 
 
		CopyWords(R1, T3, N2); 
		assert (carry >= 0 && carry <= 2); 
		Increment(R1, N2, carry); 
	} 
} 
 
// R[NA+NB] - result = A*B 
// T[NA+NB] - temporary work space 
// A[NA] ---- multiplier 
// B[NB] ---- multiplicant 
 
void AsymmetricMultiply(word *R, word *T, const word *A, unsigned int NA, const word *B, unsigned int NB) 
{ 
	if (NA == NB) 
	{ 
		if (A == B) 
			RecursiveSquare(R, T, A, NA); 
		else 
			RecursiveMultiply(R, T, A, B, NA); 
 
		return; 
	} 
 
	if (NA > NB) 
	{ 
		std::swap(A, B); 
		std::swap(NA, NB); 
	} 
 
	assert(NB % NA == 0); 
	assert((NB/NA)%2 == 0);		// NB is an even multiple of NA 
 
	if (NA==2 && !A[1]) 
	{ 
		switch (A[0]) 
		{ 
		case 0: 
			SetWords(R, 0, NB+2); 
			return; 
		case 1: 
			CopyWords(R, B, NB); 
			R[NB] = R[NB+1] = 0; 
			return; 
		default: 
			R[NB] = LinearMultiply(R, B, A[0], NB); 
			R[NB+1] = 0; 
			return; 
		} 
	} 
 
	RecursiveMultiply(R, T, A, B, NA); 
	CopyWords(T+2*NA, R+NA, NA); 
 
	unsigned i; 
 
	for (i=2*NA; i<NB; i+=2*NA) 
		RecursiveMultiply(T+NA+i, T, A, B+i, NA); 
	for (i=NA; i<NB; i+=2*NA) 
		RecursiveMultiply(R+i, T, A, B+i, NA); 
 
	if (Add(R+NA, R+NA, T+2*NA, NB-NA)) 
		Increment(R+NB, NA); 
} 
 
// R[N] ----- result = A inverse mod 2**(WORD_BITS*N) 
// T[3*N/2] - temporary work space 
// A[N] ----- an odd number as input 
 
void RecursiveInverseModPower2(word *R, word *T, const word *A, unsigned int N) 
{ 
	if (N==2) 
		AtomicInverseModPower2(R, A[0], A[1]); 
	else 
	{ 
		const unsigned int N2 = N/2; 
		RecursiveInverseModPower2(R0, T0, A0, N2); 
		T0[0] = 1; 
		SetWords(T0+1, 0, N2-1); 
		RecursiveMultiplyTop(R1, T1, T0, R0, A0, N2); 
		RecursiveMultiplyBottom(T0, T1, R0, A1, N2); 
		Add(T0, R1, T0, N2); 
		TwosComplement(T0, N2); 
		RecursiveMultiplyBottom(R1, T1, R0, T0, N2); 
	} 
} 
 
// R[N] --- result = X/(2**(WORD_BITS*N)) mod M 
// T[3*N] - temporary work space 
// X[2*N] - number to be reduced 
// M[N] --- modulus 
// U[N] --- multiplicative inverse of M mod 2**(WORD_BITS*N) 
 
void MontgomeryReduce(word *R, word *T, const word *X, const word *M, const word *U, unsigned int N) 
{ 
	RecursiveMultiplyBottom(R, T, X, U, N); 
	RecursiveMultiplyTop(T, T+N, X, R, M, N); 
	if (Subtract(R, X+N, T, N)) 
	{ 
		word carry = Add(R, R, M, N); 
		assert(carry); 
	} 
} 
 
// R[N] --- result = X/(2**(WORD_BITS*N/2)) mod M 
// T[2*N] - temporary work space 
// X[2*N] - number to be reduced 
// M[N] --- modulus 
// U[N/2] - multiplicative inverse of M mod 2**(WORD_BITS*N/2) 
// V[N] --- 2**(WORD_BITS*3*N/2) mod M 
 
void HalfMontgomeryReduce(word *R, word *T, const word *X, const word *M, const word *U, const word *V, unsigned int N) 
{ 
	assert(N%2==0 && N>=4); 
 
#define M0		M 
#define M1		(M+N2) 
#define V0		V 
#define V1		(V+N2) 
 
#define X0		X 
#define X1		(X+N2) 
#define X2		(X+N) 
#define X3		(X+N+N2) 
 
	const unsigned int N2 = N/2; 
	RecursiveMultiply(T0, T2, V0, X3, N2); 
	int c2 = Add(T0, T0, X0, N); 
	RecursiveMultiplyBottom(T3, T2, T0, U, N2); 
	RecursiveMultiplyTop(T2, R, T0, T3, M0, N2); 
	c2 -= Subtract(T2, T1, T2, N2); 
	RecursiveMultiply(T0, R, T3, M1, N2); 
	c2 -= Subtract(T0, T2, T0, N2); 
	int c3 = -(int)Subtract(T1, X2, T1, N2); 
	RecursiveMultiply(R0, T2, V1, X3, N2); 
	c3 += Add(R, R, T, N); 
 
	if (c2>0) 
		c3 += Increment(R1, N2); 
	else if (c2<0) 
		c3 -= Decrement(R1, N2, -c2); 
 
	assert(c3>=-1 && c3<=1); 
	if (c3>0) 
		Subtract(R, R, M, N); 
	else if (c3<0) 
		Add(R, R, M, N); 
 
#undef M0 
#undef M1 
#undef V0 
#undef V1 
 
#undef X0 
#undef X1 
#undef X2 
#undef X3 
} 
 
#undef A0 
#undef A1 
#undef B0 
#undef B1 
 
#undef T0 
#undef T1 
#undef T2 
#undef T3 
 
#undef R0 
#undef R1 
#undef R2 
#undef R3 
 
// do a 3 word by 2 word divide, returns quotient and leaves remainder in A 
static word SubatomicDivide(word *A, word B0, word B1) 
{ 
	// assert {A[2],A[1]} < {B1,B0}, so quotient can fit in a word 
	assert(A[2] < B1 || (A[2]==B1 && A[1] < B0)); 
 
	dword p, u; 
	word Q; 
 
	// estimate the quotient: do a 2 word by 1 word divide 
	if (B1+1 == 0) 
		Q = A[2]; 
	else 
		Q = word(MAKE_DWORD(A[1], A[2]) / (B1+1)); 
 
	// now subtract Q*B from A 
	p = (dword) B0*Q; 
	u = (dword) A[0] - LOW_WORD(p); 
	A[0] = LOW_WORD(u); 
	u = (dword) A[1] - HIGH_WORD(p) - (word)(0-HIGH_WORD(u)) - (dword)B1*Q; 
	A[1] = LOW_WORD(u); 
	A[2] += HIGH_WORD(u); 
 
	// Q <= actual quotient, so fix it 
	while (A[2] || A[1] > B1 || (A[1]==B1 && A[0]>=B0)) 
	{ 
		u = (dword) A[0] - B0; 
		A[0] = LOW_WORD(u); 
		u = (dword) A[1] - B1 - (word)(0-HIGH_WORD(u)); 
		A[1] = LOW_WORD(u); 
		A[2] += HIGH_WORD(u); 
		Q++; 
		assert(Q);  // shouldn't overflow 
	} 
 
	return Q; 
} 
 
// do a 4 word by 2 word divide, returns 2 word quotient in Q0 and Q1 
static inline void AtomicDivide(word &Q0, word &Q1, const word *A, word B0, word B1) 
{ 
	if (!B0 && !B1) // if divisor is 0, we assume divisor==2**(2*WORD_BITS) 
	{ 
		Q0 = A[2]; 
		Q1 = A[3]; 
	} 
	else 
	{ 
		word T[4]; 
		T[0] = A[0]; T[1] = A[1]; T[2] = A[2]; T[3] = A[3]; 
		Q1 = SubatomicDivide(T+1, B0, B1); 
		Q0 = SubatomicDivide(T, B0, B1); 
 
#ifndef NDEBUG 
		// multiply quotient and divisor and add remainder, make sure it equals dividend 
		assert(!T[2] && !T[3] && (T[1] < B1 || (T[1]==B1 && T[0]<B0))); 
		word P[4]; 
		AtomicMultiply(P, Q0, Q1, B0, B1); 
		Add(P, P, T, 4); 
		assert(memcmp(P, A, 4*WORD_SIZE)==0); 
#endif 
	} 
} 
 
// for use by Divide(), corrects the underestimated quotient {Q1,Q0} 
static void CorrectQuotientEstimate(word *R, word *T, word &Q0, word &Q1, const word *B, unsigned int N) 
{ 
	assert(N && N%2==0); 
 
	if (Q1) 
	{ 
		T[N] = T[N+1] = 0; 
		unsigned i; 
		for (i=0; i<N; i+=4) 
			AtomicMultiply(T+i, Q0, Q1, B[i], B[i+1]); 
		for (i=2; i<N; i+=4) 
			if (AtomicMultiplyAdd(T+i, Q0, Q1, B[i], B[i+1])) 
				T[i+5] += (++T[i+4]==0); 
	} 
	else 
	{ 
		T[N] = LinearMultiply(T, B, Q0, N); 
		T[N+1] = 0; 
	} 
 
	word borrow = Subtract(R, R, T, N+2); 
	assert(!borrow && !R[N+1]); 
 
	while (R[N] || Compare(R, B, N) >= 0) 
	{ 
		R[N] -= Subtract(R, R, B, N); 
		Q1 += (++Q0==0); 
		assert(Q0 || Q1); // no overflow 
	} 
} 
 
// R[NB] -------- remainder = A%B 
// Q[NA-NB+2] --- quotient  = A/B 
// T[NA+2*NB+4] - temp work space 
// A[NA] -------- dividend 
// B[NB] -------- divisor 
 
void Divide(word *R, word *Q, word *T, const word *A, unsigned int NA, const word *B, unsigned int NB) 
{ 
	assert(NA && NB && NA%2==0 && NB%2==0); 
	assert(B[NB-1] || B[NB-2]); 
	assert(NB <= NA); 
 
	// set up temporary work space 
	word *const TA=T; 
	word *const TB=T+NA+2; 
	word *const TP=T+NA+2+NB; 
 
	// copy B into TB and normalize it so that TB has highest bit set to 1 
	unsigned shiftWords = (B[NB-1]==0); 
	TB[0] = TB[NB-1] = 0; 
	CopyWords(TB+shiftWords, B, NB-shiftWords); 
	unsigned shiftBits = WORD_BITS - BitPrecision(TB[NB-1]); 
	assert(shiftBits < WORD_BITS); 
	ShiftWordsLeftByBits(TB, NB, shiftBits); 
 
	// copy A into TA and normalize it 
	TA[0] = TA[NA] = TA[NA+1] = 0; 
	CopyWords(TA+shiftWords, A, NA); 
	ShiftWordsLeftByBits(TA, NA+2, shiftBits); 
 
	if (TA[NA+1]==0 && TA[NA] <= 1) 
	{ 
		Q[NA-NB+1] = Q[NA-NB] = 0; 
		while (TA[NA] || Compare(TA+NA-NB, TB, NB) >= 0) 
		{ 
			TA[NA] -= Subtract(TA+NA-NB, TA+NA-NB, TB, NB); 
			++Q[NA-NB]; 
		} 
	} 
	else 
	{ 
		NA+=2; 
		assert(Compare(TA+NA-NB, TB, NB) < 0); 
	} 
 
	word B0 = TB[NB-2] + 1; 
	word B1 = TB[NB-1] + (B0==0); 
 
	// start reducing TA mod TB, 2 words at a time 
	for (unsigned i=NA-2; i>=NB; i-=2) 
	{ 
		AtomicDivide(Q[i-NB], Q[i-NB+1], TA+i-2, B0, B1); 
		CorrectQuotientEstimate(TA+i-NB, TP, Q[i-NB], Q[i-NB+1], TB, NB); 
	} 
 
	// copy TA into R, and denormalize it 
	CopyWords(R, TA+shiftWords, NB); 
	ShiftWordsRightByBits(R, NB, shiftBits); 
} 
 
static inline unsigned int EvenWordCount(const word *X, unsigned int N) 
{ 
	while (N && X[N-2]==0 && X[N-1]==0) 
		N-=2; 
	return N; 
} 
 
// return k 
// R[N] --- result = A^(-1) * 2^k mod M 
// T[4*N] - temporary work space 
// A[NA] -- number to take inverse of 
// M[N] --- modulus 
 
unsigned int AlmostInverse(word *R, word *T, const word *A, unsigned int NA, const word *M, unsigned int N) 
{ 
	assert(NA<=N && N && N%2==0); 
 
	word *b = T; 
	word *c = T+N; 
	word *f = T+2*N; 
	word *g = T+3*N; 
	unsigned int bcLen=2, fgLen=EvenWordCount(M, N); 
	unsigned int k=0, s=0; 
 
	SetWords(T, 0, 3*N); 
	b[0]=1; 
	CopyWords(f, A, NA); 
	CopyWords(g, M, N); 
 
	while (1) 
	{ 
		word t=f[0]; 
		while (!t) 
		{ 
			if (EvenWordCount(f, fgLen)==0) 
			{ 
				SetWords(R, 0, N); 
				return 0; 
			} 
 
			ShiftWordsRightByWords(f, fgLen, 1); 
			if (c[bcLen-1]) bcLen+=2; 
			assert(bcLen <= N); 
			ShiftWordsLeftByWords(c, bcLen, 1); 
			k+=WORD_BITS; 
			t=f[0]; 
		} 
 
		unsigned int i=0; 
		while (t%2 == 0) 
		{ 
			t>>=1; 
			i++; 
		} 
		k+=i; 
 
		if (t==1 && f[1]==0 && EvenWordCount(f, fgLen)==2) 
		{ 
			if (s%2==0) 
				CopyWords(R, b, N); 
			else 
				Subtract(R, M, b, N); 
			return k; 
		} 
 
		ShiftWordsRightByBits(f, fgLen, i); 
		t=ShiftWordsLeftByBits(c, bcLen, i); 
		if (t) 
		{ 
			c[bcLen] = t; 
			bcLen+=2; 
			assert(bcLen <= N); 
		} 
 
		if (f[fgLen-2]==0 && g[fgLen-2]==0 && f[fgLen-1]==0 && g[fgLen-1]==0) 
			fgLen-=2; 
 
		if (Compare(f, g, fgLen)==-1) 
		{ 
			std::swap(f, g); 
			std::swap(b, c); 
			s++; 
		} 
 
		Subtract(f, f, g, fgLen); 
 
		if (Add(b, b, c, bcLen)) 
		{ 
			b[bcLen] = 1; 
			bcLen+=2; 
			assert(bcLen <= N); 
		} 
	} 
} 
 
// R[N] - result = A/(2^k) mod M 
// A[N] - input 
// M[N] - modulus 
 
void DivideByPower2Mod(word *R, const word *A, unsigned int k, const word *M, unsigned int N) 
{ 
	CopyWords(R, A, N); 
 
	while (k--) 
	{ 
		if (R[0]%2==0) 
			ShiftWordsRightByBits(R, N, 1); 
		else 
		{ 
			word carry = Add(R, R, M, N); 
			ShiftWordsRightByBits(R, N, 1); 
			R[N-1] += carry<<(WORD_BITS-1); 
		} 
	} 
} 
 
// R[N] - result = A*(2^k) mod M 
// A[N] - input 
// M[N] - modulus 
 
void MultiplyByPower2Mod(word *R, const word *A, unsigned int k, const word *M, unsigned int N) 
{ 
	CopyWords(R, A, N); 
 
	while (k--) 
		if (ShiftWordsLeftByBits(R, N, 1) || Compare(R, M, N)>=0) 
			Subtract(R, R, M, N); 
} 
 
// ****************************************************************** 
 
static const unsigned int RoundupSizeTable[] = {2, 2, 2, 4, 4, 8, 8, 8, 8}; 
 
static inline unsigned int RoundupSize(unsigned int n) 
{ 
	if (n<=8) 
		return RoundupSizeTable[n]; 
	else if (n<=16) 
		return 16; 
	else if (n<=32) 
		return 32; 
	else if (n<=64) 
		return 64; 
	else return 1U << BitPrecision(n-1); 
} 
 
Integer::Integer() 
	: reg(2), sign(POSITIVE) 
{ 
	reg[0] = reg[1] = 0; 
} 
 
Integer::Integer(const Integer& t) 
	: reg(RoundupSize(t.WordCount())), sign(t.sign) 
{ 
	CopyWords(reg, t.reg, reg.size); 
} 
 
Integer::Integer(signed long value) 
	: reg(2) 
{ 
	if (value >= 0) 
		sign = POSITIVE; 
	else 
	{ 
		sign = NEGATIVE; 
		value = -value; 
	} 
	reg[0] = word(value); 
	reg[1] = sizeof(value)>WORD_SIZE ? word(value>>WORD_BITS) : 0; 
} 
 
signed long Integer::ConvertToLong() const 
{ 
	unsigned long value = reg[0]; 
	value += sizeof(value)>WORD_SIZE ? ((unsigned long)reg[1]<<WORD_BITS) : 0; 
	return sign==POSITIVE ? value : -long(value); 
} 
 
Integer::Integer(const byte *encodedInteger, unsigned int byteCount, Signedness s) 
{ 
	Decode(encodedInteger, byteCount, s); 
} 
 
Integer::Integer(const byte *BEREncodedInteger) 
{ 
	BERDecode(BEREncodedInteger); 
} 
 
Integer::Integer(BufferedTransformation &bt) 
{ 
	BERDecode(bt); 
} 
 
Integer::Integer(RandomNumberGenerator &rng, unsigned int bitcount) 
{ 
	Randomize(rng, bitcount); 
} 
 
Integer::Integer(RandomNumberGenerator &rng, const Integer &min, const Integer &max, RandomNumberType rnType, const Integer &equiv, const Integer &mod) 
{ 
	if (!Randomize(rng, min, max, rnType, equiv, mod)) 
		throw Integer::RandomNumberNotFound(); 
} 
 
Integer Integer::Power2(unsigned int e) 
{ 
	Integer r((word)0, bitsToWords(e+1));  
	r.SetBit(e);  
	return r; 
} 
 
const Integer &Integer::Zero() 
{ 
	static const Integer zero; 
	return zero; 
} 
 
const Integer &Integer::One() 
{ 
	static const Integer one(1,2); 
	return one; 
} 
 
bool Integer::operator!() const 
{ 
	return IsNegative() ? false : (reg[0]==0 && WordCount()==0); 
} 
 
Integer& Integer::operator=(const Integer& t) 
{ 
	if (this != &t) 
	{ 
		reg.New(RoundupSize(t.WordCount())); 
		CopyWords(reg, t.reg, reg.size); 
		sign = t.sign; 
	} 
	return *this; 
} 
 
bool Integer::GetBit(unsigned int n) const 
{ 
	if (n/WORD_BITS >= reg.size) 
		return 0;