using System; using System.Collections; using System.Globalization; using org.bouncycastle.security; namespace org.bouncycastle.math { public class BigInteger { private int sign; // -1 means -ve; +1 means +ve; 0 means 0; private int[] magnitude; // array of ints with [0] being the most significant private int nBits = -1; // cache bitCount() value private int nBitLength = -1; // cache bitLength() value private static readonly long IMASK = 0xffffffffL; private long mQuote = -1L; // -m^(-1) mod b, b = 2^32 (see Montgomery mult.) private BigInteger() { } private BigInteger(int nWords) { sign = 1; magnitude = new int[nWords]; } private BigInteger(int signum, int[] mag) { sign = signum; if (mag.Length > 0) { int i = 0; while (i < mag.Length && mag[i] == 0) { i++; } if (i == 0) { magnitude = mag; } else { // strip leading 0 bytes int[] newMag = new int[mag.Length - i]; Array.Copy(mag, i, newMag, 0, newMag.Length); magnitude = newMag; if (newMag.Length == 0) sign = 0; } } else { magnitude = mag; sign = 0; } } public BigInteger(String sval) //throws FormatException : this(sval, 10) { } public BigInteger(String sval, int rdx) //throws FormatException { if (sval.Length == 0) { throw new FormatException("Zero length BigInteger"); } NumberStyles style; switch (rdx) { case 10: style = NumberStyles.Integer; break; case 16: style = NumberStyles.AllowHexSpecifier; break; default: throw new FormatException("Only base 10 or 16 alllowed"); } int index = 0; sign = 1; if (sval[0] == '-') { if (sval.Length == 1) { throw new FormatException("Zero length BigInteger"); } sign = -1; index = 1; } // strip leading zeros from the string value while (index < sval.Length && Int32.Parse(sval[index].ToString(), style) == 0) { index++; } if (index >= sval.Length) { // zero value - we're done sign = 0; magnitude = new int[0]; return; } ////// // could we work out the max number of ints required to store // sval.length digits in the given base, then allocate that // storage in one hit?, then generate the magnitude in one hit too? ////// BigInteger b = BigInteger.ZERO; BigInteger r = valueOf(rdx); while (index < sval.Length) { // (optimise this by taking chunks of digits instead?) b = b.multiply(r).add(valueOf(Int32.Parse(sval[index].ToString(), style))); index++; } magnitude = b.magnitude; return; } public BigInteger(byte[] bval) //throws FormatException { if (bval.Length == 0) { throw new FormatException("Zero length BigInteger"); } sign = 1; if (bval[0] < 0) { // FIXME: int iBval; sign = -1; // strip leading sign bytes for (iBval = 0; iBval < bval.Length && ((sbyte)bval[iBval] == -1); iBval++) ; magnitude = new int[(bval.Length - iBval) / 2 + 1]; // copy bytes to magnitude // invert bytes then add one to find magnitude of value } else { // strip leading zero bytes and return magnitude bytes magnitude = makeMagnitude(bval); } } private int[] makeMagnitude(byte[] bval) { int i; int[] mag; int firstSignificant; // strip leading zeros for (firstSignificant = 0; firstSignificant < bval.Length && bval[firstSignificant] == 0; firstSignificant++) ; if (firstSignificant >= bval.Length) { return new int[0]; } int nInts = (bval.Length - firstSignificant + 3) / 4; int bCount = (bval.Length - firstSignificant) % 4; if (bCount == 0) bCount = 4; mag = new int[nInts]; int v = 0; int magnitudeIndex = 0; for (i = firstSignificant; i < bval.Length; i++) { v <<= 8; v |= bval[i] & 0xff; bCount--; if (bCount <= 0) { mag[magnitudeIndex] = v; magnitudeIndex++; bCount = 4; v = 0; } } if (magnitudeIndex < mag.Length) { mag[magnitudeIndex] = v; } return mag; } public BigInteger(int sign, byte[] mag) //throws FormatException { if (sign < -1 || sign > 1) { throw new FormatException("Invalid sign value"); } if (sign == 0) { this.sign = 0; this.magnitude = new int[0]; return; } // copy bytes this.magnitude = makeMagnitude(mag); this.sign = sign; } public BigInteger(int numBits, Random rnd) //throws ArgumentException { if (numBits < 0) { throw new ArgumentException("numBits must be non-negative"); } int nBytes = (numBits + 7) / 8; byte[] b = new byte[nBytes]; if (nBytes > 0) { nextRndBytes(rnd, b); // strip off any excess bits in the MSB b[0] &= rndMask[8 * nBytes - numBits]; } this.magnitude = makeMagnitude(b); this.sign = 1; this.nBits = -1; this.nBitLength = -1; } private static readonly int BITS_PER_BYTE = 8; private static readonly int BYTES_PER_INT = 4; /** * strictly speaking this is a little dodgey from a compliance * point of view as it forces people to be using SecureRandom as * well, that being said - this implementation is for a crypto * library and you do have the source! */ private void nextRndBytes(Random rnd, byte[] bytes) { int numRequested = bytes.Length; int numGot = 0, r = 0; if (typeof(SecureRandom).IsInstanceOfType(rnd)) { ((SecureRandom)rnd).nextBytes(bytes); } else { for (; ; ) { for (int i = 0; i < BYTES_PER_INT; i++) { if (numGot == numRequested) { return; } r = (i == 0 ? rnd.Next() : r >> BITS_PER_BYTE); bytes[numGot++] = (byte)r; } } } } private static readonly byte[] rndMask = { (byte)255, 127, 63, 31, 15, 7, 3, 1 }; public BigInteger(int bitLength, int certainty, Random rnd) //throws ArithmeticException { int nBytes = (bitLength + 7) / 8; byte[] b = new byte[nBytes]; do { if (nBytes > 0) { nextRndBytes(rnd, b); // strip off any excess bits in the MSB b[0] &= rndMask[8 * nBytes - bitLength]; } this.magnitude = makeMagnitude(b); this.sign = 1; this.nBits = -1; this.nBitLength = -1; this.mQuote = -1L; if (certainty > 0 && bitLength > 2) { this.magnitude[this.magnitude.Length - 1] |= 1; } } while (this.bitLength() != bitLength || !this.isProbablePrime(certainty)); } public BigInteger abs() { return (sign >= 0) ? this : this.negate(); } /** * return a = a + b - b preserved. */ private int[] add(int[] a, int[] b) { int tI = a.Length - 1; int vI = b.Length - 1; long m = 0; while (vI >= 0) { m += (((long)a[tI]) & IMASK) + (((long)b[vI--]) & IMASK); a[tI--] = (int)m; m = (long)((ulong)m >> 32); } while (tI >= 0 && m != 0) { m += (((long)a[tI]) & IMASK); a[tI--] = (int)m; m = (long)((ulong)m >> 32); } return a; } public BigInteger add(BigInteger val) //throws ArithmeticException { if (val.sign == 0 || val.magnitude.Length == 0) return this; if (this.sign == 0 || this.magnitude.Length == 0) return val; if (val.sign < 0) { if (this.sign > 0) return this.subtract(val.negate()); } else { if (this.sign < 0) return val.subtract(this.negate()); } // both BigIntegers are either +ve or -ve; set the sign later int[] mag, op; if (this.magnitude.Length < val.magnitude.Length) { mag = new int[val.magnitude.Length + 1]; Array.Copy(val.magnitude, 0, mag, 1, val.magnitude.Length); op = this.magnitude; } else { mag = new int[this.magnitude.Length + 1]; Array.Copy(this.magnitude, 0, mag, 1, this.magnitude.Length); op = val.magnitude; } return new BigInteger(this.sign, add(mag, op)); } public int bitCount() { if (nBits == -1) { nBits = 0; for (int i = 0; i < magnitude.Length; i++) { nBits += bitCounts[magnitude[i] & 0xff]; nBits += bitCounts[(magnitude[i] >> 8) & 0xff]; nBits += bitCounts[(magnitude[i] >> 16) & 0xff]; nBits += bitCounts[(magnitude[i] >> 24) & 0xff]; } } return nBits; } private readonly static byte[] bitCounts = {0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8}; private int bitLength(int indx, int[] mag) { int bitLength; if (mag.Length == 0) { return 0; } else { while (indx != mag.Length && mag[indx] == 0) { indx++; } if (indx == mag.Length) { return 0; } // bit length for everything after the first int bitLength = 32 * ((mag.Length - indx) - 1); // and determine bitlength of first int bitLength += bitLen(mag[indx]); if (sign < 0) { // Check if magnitude is a power of two bool pow2 = ((bitCounts[mag[indx] & 0xff]) + (bitCounts[(mag[indx] >> 8) & 0xff]) + (bitCounts[(mag[indx] >> 16) & 0xff]) + (bitCounts[(mag[indx] >> 24) & 0xff])) == 1; for (int i = indx + 1; i < mag.Length && pow2; i++) { pow2 = (mag[i] == 0); } bitLength -= (pow2 ? 1 : 0); } } return bitLength; } public int bitLength() { if (nBitLength == -1) { if (sign == 0) { nBitLength = 0; } else { nBitLength = bitLength(0, magnitude); } } return nBitLength; } // // bitLen(val) is the number of bits in val. // static int bitLen(int w) { // Binary search - decision tree (5 tests, rarely 6) return (w < 1 << 15 ? (w < 1 << 7 ? (w < 1 << 3 ? (w < 1 << 1 ? (w < 1 << 0 ? (w < 0 ? 32 : 0) : 1) : (w < 1 << 2 ? 2 : 3)) : (w < 1 << 5 ? (w < 1 << 4 ? 4 : 5) : (w < 1 << 6 ? 6 : 7))) : (w < 1 << 11 ? (w < 1 << 9 ? (w < 1 << 8 ? 8 : 9) : (w < 1 << 10 ? 10 : 11)) : (w < 1 << 13 ? (w < 1 << 12 ? 12 : 13) : (w < 1 << 14 ? 14 : 15)))) : (w < 1 << 23 ? (w < 1 << 19 ? (w < 1 << 17 ? (w < 1 << 16 ? 16 : 17) : (w < 1 << 18 ? 18 : 19)) : (w < 1 << 21 ? (w < 1 << 20 ? 20 : 21) : (w < 1 << 22 ? 22 : 23))) : (w < 1 << 27 ? (w < 1 << 25 ? (w < 1 << 24 ? 24 : 25) : (w < 1 << 26 ? 26 : 27)) : (w < 1 << 29 ? (w < 1 << 28 ? 28 : 29) : (w < 1 << 30 ? 30 : 31))))); } private readonly static byte[] bitLengths = {0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8}; public int compareTo(Object o) { return compareTo((BigInteger)o); } /** * unsigned comparison on two arrays - note the arrays may * start with leading zeros. */ private int compareTo(int xIndx, int[] x, int yIndx, int[] y) { while (xIndx != x.Length && x[xIndx] == 0) { xIndx++; } while (yIndx != y.Length && y[yIndx] == 0) { yIndx++; } if ((x.Length - xIndx) < (y.Length - yIndx)) { return -1; } if ((x.Length - xIndx) > (y.Length - yIndx)) { return 1; } // lengths of magnitudes the same, test the magnitude values while (xIndx < x.Length) { long v1 = (long)(x[xIndx++]) & IMASK; long v2 = (long)(y[yIndx++]) & IMASK; if (v1 < v2) { return -1; } if (v1 > v2) { return 1; } } return 0; } public int compareTo(BigInteger val) { if (sign < val.sign) return -1; if (sign > val.sign) return 1; return compareTo(0, magnitude, 0, val.magnitude); } /** * return z = x / y - done in place (z value preserved, x contains the * remainder) */ private int[] divide(int[] x, int[] y) { int xyCmp = compareTo(0, x, 0, y); int[] count; if (xyCmp > 0) { int[] c; int shift = bitLength(0, x) - bitLength(0, y); if (shift > 1) { c = shiftLeft(y, shift - 1); count = shiftLeft(ONE.magnitude, shift - 1); if (shift % 32 == 0) { // Special case where the shift is the size of an int. int[] countSpecial = new int[shift / 32 + 1]; Array.Copy(count, 0, countSpecial, 1, countSpecial.Length - 1); countSpecial[0] = 0; count = countSpecial; } } else { c = new int[x.Length]; count = new int[1]; Array.Copy(y, 0, c, c.Length - y.Length, y.Length); count[0] = 1; } int[] iCount = new int[count.Length]; subtract(0, x, 0, c); Array.Copy(count, 0, iCount, 0, count.Length); int xStart = 0; int cStart = 0; int iCountStart = 0; for (; ; ) { int cmp = compareTo(xStart, x, cStart, c); while (cmp >= 0) { subtract(xStart, x, cStart, c); add(count, iCount); cmp = compareTo(xStart, x, cStart, c); } xyCmp = compareTo(xStart, x, 0, y); if (xyCmp > 0) { if (x[xStart] == 0) { xStart++; } shift = bitLength(cStart, c) - bitLength(xStart, x); if (shift == 0) { c = shiftRightOne(cStart, c); iCount = shiftRightOne(iCountStart, iCount); } else { c = shiftRight(cStart, c, shift); iCount = shiftRight(iCountStart, iCount, shift); } if (c[cStart] == 0) { cStart++; } if (iCount[iCountStart] == 0) { iCountStart++; } } else if (xyCmp == 0) { add(count, ONE.magnitude); for (int i = xStart; i != x.Length; i++) { x[i] = 0; } break; } else { break; } } } else if (xyCmp == 0) { count = new int[1]; count[0] = 1; } else { count = new int[1]; count[0] = 0; } return count; } public BigInteger divide(BigInteger val) //throws ArithmeticException { if (val.sign == 0) { throw new ArithmeticException("Divide by zero"); } if (sign == 0) { return BigInteger.ZERO; } if (val.compareTo(BigInteger.ONE) == 0) { return this; } int[] mag = new int[this.magnitude.Length]; Array.Copy(this.magnitude, 0, mag, 0, mag.Length); return new BigInteger(this.sign * val.sign, divide(mag, val.magnitude)); } public BigInteger[] divideAndRemainder(BigInteger val) //throws ArithmeticException { if (val.sign == 0) { throw new ArithmeticException("Divide by zero"); } BigInteger[] biggies = new BigInteger[2]; if (sign == 0) { biggies[0] = biggies[1] = BigInteger.ZERO; return biggies; } if (val.compareTo(BigInteger.ONE) == 0) { biggies[0] = this; biggies[1] = BigInteger.ZERO; return biggies; } int[] remainder = new int[this.magnitude.Length]; Array.Copy(this.magnitude, 0, remainder, 0, remainder.Length); int[] quotient = divide(remainder, val.magnitude); biggies[0] = new BigInteger(this.sign * val.sign, quotient); biggies[1] = new BigInteger(this.sign, remainder); return biggies; } public override bool Equals(Object val) { if (val == this) return true; if (!(typeof(BigInteger).IsInstanceOfType(val))) return false; BigInteger biggie = (BigInteger)val; if (biggie.sign != sign || biggie.magnitude.Length != magnitude.Length) return false; for (int i = 0; i < magnitude.Length; i++) { if (biggie.magnitude[i] != magnitude[i]) return false; } return true; } public BigInteger gcd(BigInteger val) { if (val.sign == 0) return this.abs(); else if (sign == 0) return val.abs(); BigInteger r; BigInteger u = this; BigInteger v = val; while (v.sign != 0) { r = u.mod(v); u = v; v = r; } return u; } public override int GetHashCode() { return 0; } public int intValue() { if (magnitude.Length == 0) { return 0; } if (sign < 0) { return -magnitude[magnitude.Length - 1]; } else { return magnitude[magnitude.Length - 1]; } } /** * return whether or not a BigInteger is probably prime with a * probability of 1 - (1/2)**certainty. *
* From Knuth Vol 2, pg 395.
*/
public bool isProbablePrime(int certainty)
{
if (certainty == 0)
{
return true;
}
BigInteger n = this.abs();
if (n.Equals(TWO))
{
return true;
}
if (n.Equals(ONE) || !n.testBit(0))
{
return false;
}
if ((certainty & 0x1) == 1)
{
certainty = certainty / 2 + 1;
}
else
{
certainty /= 2;
}
//
// let n = 1 + 2^kq
//
BigInteger q = n.subtract(ONE);
int k = q.getLowestSetBit();
q = q.shiftRight(k);
Random rnd = new Random();
for (int i = 0; i <= certainty; i++)
{
BigInteger x;
do
{
x = new BigInteger(n.bitLength(), rnd);
} while (x.compareTo(ONE) <= 0 || x.compareTo(n) >= 0);
int j = 0;
BigInteger y = x.modPow(q, n);
while (!((j == 0 && y.Equals(ONE)) || y.Equals(n.subtract(ONE))))
{
if (j > 0 && y.Equals(ONE))
{
return false;
}
if (++j == k)
{
return false;
}
y = y.modPow(TWO, n);
}
}
return true;
}
public long longValue()
{
long val = 0;
if (magnitude.Length == 0)
{
return 0;
}
if (magnitude.Length > 1)
{
val = ((long)magnitude[magnitude.Length - 2] << 32)
| (magnitude[magnitude.Length - 1] & IMASK);
}
else
{
val = (magnitude[magnitude.Length - 1] & IMASK);
}
if (sign < 0)
{
return -val;
}
else
{
return val;
}
}
public BigInteger max(BigInteger val)
{
return (compareTo(val) > 0) ? this : val;
}
public BigInteger min(BigInteger val)
{
return (compareTo(val) < 0) ? this : val;
}
public BigInteger mod(BigInteger m) //throws ArithmeticException
{
if (m.sign <= 0)
{
throw new ArithmeticException("BigInteger: modulus is not positive");
}
BigInteger biggie = this.remainder(m);
return (biggie.sign >= 0 ? biggie : biggie.add(m));
}
public BigInteger modInverse(BigInteger m) //throws ArithmeticException
{
if (m.sign != 1)
{
throw new ArithmeticException("Modulus must be positive");
}
BigInteger x = new BigInteger();
BigInteger y = new BigInteger();
BigInteger gcd = BigInteger.extEuclid(this, m, x, y);
if (!gcd.Equals(BigInteger.ONE))
{
throw new ArithmeticException("Numbers not relatively prime.");
}
if (x.compareTo(BigInteger.ZERO) < 0)
{
x = x.add(m);
}
return x;
}
/**
* Calculate the numbers u1, u2, and u3 such that:
*
* u1 * a + u2 * b = u3
*
* where u3 is the greatest common divider of a and b.
* a and b using the extended Euclid algorithm (refer p. 323
* of The Art of Computer Programming vol 2, 2nd ed).
* This also seems to have the side effect of calculating
* some form of multiplicative inverse.
*
* @param a First number to calculate gcd for
* @param b Second number to calculate gcd for
* @param u1Out the return object for the u1 value
* @param u2Out the return object for the u2 value
* @return The greatest common divisor of a and b
*/
private static BigInteger extEuclid(BigInteger a, BigInteger b, BigInteger u1Out,
BigInteger u2Out)
{
BigInteger res;
BigInteger u1 = BigInteger.ONE;
BigInteger u3 = a;
BigInteger v1 = BigInteger.ZERO;
BigInteger v3 = b;
while (v3.compareTo(BigInteger.ZERO) > 0)
{
BigInteger q,
tn;
//tv;
q = u3.divide(v3);
tn = u1.subtract(v1.multiply(q));
u1 = v1;
v1 = tn;
tn = u3.subtract(v3.multiply(q));
u3 = v3;
v3 = tn;
}
u1Out.sign = u1.sign;
u1Out.magnitude = u1.magnitude;
res = u3.subtract(u1.multiply(a)).divide(b);
u2Out.sign = res.sign;
u2Out.magnitude = res.magnitude;
return u3;
}
/**
* zero out the array x
*/
private void zero(int[] x)
{
for (int i = 0; i != x.Length; i++)
{
x[i] = 0;
}
}
public BigInteger modPow(
BigInteger exponent,
BigInteger m)
//throws ArithmeticException
{
int[] zVal = null;
int[] yAccum = null;
int[] yVal;
// Montgomery exponentiation is only possible if the modulus is odd,
// but AFAIK, this is always the case for crypto algo's
bool useMonty = ((m.magnitude[m.magnitude.Length - 1] & 1) == 1);
long mQ = 0;
if (useMonty)
{
mQ = m.getMQuote();
// tmp = this * R mod m
BigInteger tmp = this.shiftLeft(32 * m.magnitude.Length).mod(m);
zVal = tmp.magnitude;
useMonty = (zVal.Length == m.magnitude.Length);
if (useMonty)
{
yAccum = new int[m.magnitude.Length + 1];
}
}
if (!useMonty)
{
if (magnitude.Length <= m.magnitude.Length)
{
//zAccum = new int[m.magnitude.Length * 2];
zVal = new int[m.magnitude.Length];
Array.Copy(magnitude, 0, zVal, zVal.Length - magnitude.Length,
magnitude.Length);
}
else
{
//
// in normal practice we'll never see this...
//
BigInteger tmp = this.remainder(m);
//zAccum = new int[m.magnitude.Length * 2];
zVal = new int[m.magnitude.Length];
Array.Copy(tmp.magnitude, 0, zVal, zVal.Length - tmp.magnitude.Length,
tmp.magnitude.Length);
}
yAccum = new int[m.magnitude.Length * 2];
}
yVal = new int[m.magnitude.Length];
//
// from LSW to MSW
//
for (int i = 0; i < exponent.magnitude.Length; i++)
{
int v = exponent.magnitude[i];
int bits = 0;
if (i == 0)
{
while (v > 0)
{
v <<= 1;
bits++;
}
//
// first time in initialise y
//
Array.Copy(zVal, 0, yVal, 0, zVal.Length);
v <<= 1;
bits++;
}
while (v != 0)
{
if (useMonty)
{
// Montgomery square algo doesn't exist, and a normal
// square followed by a Montgomery reduction proved to
// be almost as heavy as a Montgomery mulitply.
multiplyMonty(yAccum, yVal, yVal, m.magnitude, mQ);
}
else
{
square(yAccum, yVal);
remainder(yAccum, m.magnitude);
Array.Copy(yAccum, yAccum.Length - yVal.Length, yVal, 0, yVal.Length);
zero(yAccum);
}
bits++;
if (v < 0)
{
if (useMonty)
{
multiplyMonty(yAccum, yVal, zVal, m.magnitude, mQ);
}
else
{
multiply(yAccum, yVal, zVal);
remainder(yAccum, m.magnitude);
Array.Copy(yAccum, yAccum.Length - yVal.Length, yVal, 0,
yVal.Length);
zero(yAccum);
}
}
v <<= 1;
}
while (bits < 32)
{
if (useMonty)
{
multiplyMonty(yAccum, yVal, yVal, m.magnitude, mQ);
}
else
{
square(yAccum, yVal);
remainder(yAccum, m.magnitude);
Array.Copy(yAccum, yAccum.Length - yVal.Length, yVal, 0, yVal.Length);
zero(yAccum);
}
bits++;
}
}
if (useMonty)
{
// Return y * R^(-1) mod m by doing y * 1 * R^(-1) mod m
zero(zVal);
zVal[zVal.Length - 1] = 1;
multiplyMonty(yAccum, yVal, zVal, m.magnitude, mQ);
}
return new BigInteger(1, yVal);
}
/**
* return w with w = x * x - w is assumed to have enough space.
*/
private int[] square(int[] w, int[] x)
{
long u1,
u2,
c;
if (w.Length != 2 * x.Length)
{
throw new ArgumentException("no I don't think so...");
}
for (int i = x.Length - 1; i != 0; i--)
{
long v = (x[i] & IMASK);
u1 = v * v;
u2 = (long)((ulong)u1 >> 32);
u1 = u1 & IMASK;
u1 += (w[2 * i + 1] & IMASK);
w[2 * i + 1] = (int)u1;
c = u2 + (u1 >> 32);
for (int j = i - 1; j >= 0; j--)
{
u1 = (x[j] & IMASK) * v;
u2 = (long)((ulong)u1 >> 31); // multiply by 2!
u1 = (u1 & 0x7fffffff) << 1; // multiply by 2!
u1 += (w[i + j + 1] & IMASK) + c;
w[i + j + 1] = (int)u1;
c = u2 + (long)((ulong)u1 >> 32);
}
c += w[i] & IMASK;
w[i] = (int)c;
w[i - 1] = (int)(c >> 32);
}
u1 = (x[0] & IMASK);
u1 = u1 * u1;
u2 = (long)((ulong)u1 >> 32);
u1 = u1 & IMASK;
u1 += (w[1] & IMASK);
w[1] = (int)u1;
w[0] = (int)(u2 + (u1 >> 32) + w[0]);
return w;
}
/**
* return x with x = y * z - x is assumed to have enough space.
*/
private int[] multiply(int[] x, int[] y, int[] z)
{
for (int i = z.Length - 1; i >= 0; i--)
{
long a = z[i] & IMASK;
long value = 0;
for (int j = y.Length - 1; j >= 0; j--)
{
value += a * (y[j] & IMASK) + (x[i + j + 1] & IMASK);
x[i + j + 1] = (int)value;
value = (long)((ulong)value >> 32);
}
x[i] = (int)value;
}
return x;
}
/**
* Calculate mQuote = -m^(-1) mod b with b = 2^32 (32 = word size)
*/
private long getMQuote()
{
if (mQuote != -1L)
{ // allready calculated
return mQuote;
}
if ((magnitude[magnitude.Length - 1] & 1) == 0)
{
return -1L; // not for even numbers
}
byte[] bytes = { 1, 0, 0, 0, 0 };
BigInteger b = new BigInteger(1, bytes); // 2^32
mQuote = this.negate().mod(b).modInverse(b).longValue();
return mQuote;
}
/**
* Montgomery multiplication: a = x * y * R^(-1) mod m
*
* Based algorithm 14.36 of Handbook of Applied Cryptography.
*
*