#if !BESTHTTP_DISABLE_ALTERNATE_SSL && (!UNITY_WEBGL || UNITY_EDITOR) using System; using System.Collections; using System.Diagnostics; using System.Globalization; using System.Text; using Org.BouncyCastle.Security; using Org.BouncyCastle.Utilities; namespace Org.BouncyCastle.Math { #if !(NETCF_1_0 || NETCF_2_0 || SILVERLIGHT || NETFX_CORE || PORTABLE) [Serializable] #endif public class BigInteger { // The first few odd primes /* 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069 1087 1091 1093 1097 1103 1109 1117 1123 1129 1151 1153 1163 1171 1181 1187 1193 1201 1213 1217 1223 1229 1231 1237 1249 1259 1277 1279 1283 1289 */ // Each list has a product < 2^31 internal static readonly int[][] primeLists = new int[][] { new int[]{ 3, 5, 7, 11, 13, 17, 19, 23 }, new int[]{ 29, 31, 37, 41, 43 }, new int[]{ 47, 53, 59, 61, 67 }, new int[]{ 71, 73, 79, 83 }, new int[]{ 89, 97, 101, 103 }, new int[]{ 107, 109, 113, 127 }, new int[]{ 131, 137, 139, 149 }, new int[]{ 151, 157, 163, 167 }, new int[]{ 173, 179, 181, 191 }, new int[]{ 193, 197, 199, 211 }, new int[]{ 223, 227, 229 }, new int[]{ 233, 239, 241 }, new int[]{ 251, 257, 263 }, new int[]{ 269, 271, 277 }, new int[]{ 281, 283, 293 }, new int[]{ 307, 311, 313 }, new int[]{ 317, 331, 337 }, new int[]{ 347, 349, 353 }, new int[]{ 359, 367, 373 }, new int[]{ 379, 383, 389 }, new int[]{ 397, 401, 409 }, new int[]{ 419, 421, 431 }, new int[]{ 433, 439, 443 }, new int[]{ 449, 457, 461 }, new int[]{ 463, 467, 479 }, new int[]{ 487, 491, 499 }, new int[]{ 503, 509, 521 }, new int[]{ 523, 541, 547 }, new int[]{ 557, 563, 569 }, new int[]{ 571, 577, 587 }, new int[]{ 593, 599, 601 }, new int[]{ 607, 613, 617 }, new int[]{ 619, 631, 641 }, new int[]{ 643, 647, 653 }, new int[]{ 659, 661, 673 }, new int[]{ 677, 683, 691 }, new int[]{ 701, 709, 719 }, new int[]{ 727, 733, 739 }, new int[]{ 743, 751, 757 }, new int[]{ 761, 769, 773 }, new int[]{ 787, 797, 809 }, new int[]{ 811, 821, 823 }, new int[]{ 827, 829, 839 }, new int[]{ 853, 857, 859 }, new int[]{ 863, 877, 881 }, new int[]{ 883, 887, 907 }, new int[]{ 911, 919, 929 }, new int[]{ 937, 941, 947 }, new int[]{ 953, 967, 971 }, new int[]{ 977, 983, 991 }, new int[]{ 997, 1009, 1013 }, new int[]{ 1019, 1021, 1031 }, new int[]{ 1033, 1039, 1049 }, new int[]{ 1051, 1061, 1063 }, new int[]{ 1069, 1087, 1091 }, new int[]{ 1093, 1097, 1103 }, new int[]{ 1109, 1117, 1123 }, new int[]{ 1129, 1151, 1153 }, new int[]{ 1163, 1171, 1181 }, new int[]{ 1187, 1193, 1201 }, new int[]{ 1213, 1217, 1223 }, new int[]{ 1229, 1231, 1237 }, new int[]{ 1249, 1259, 1277 }, new int[]{ 1279, 1283, 1289 }, }; internal static readonly int[] primeProducts; private const long IMASK = 0xFFFFFFFFL; private const ulong UIMASK = 0xFFFFFFFFUL; private static readonly int[] ZeroMagnitude = new int[0]; private static readonly byte[] ZeroEncoding = new byte[0]; private static readonly BigInteger[] SMALL_CONSTANTS = new BigInteger[17]; public static readonly BigInteger Zero; public static readonly BigInteger One; public static readonly BigInteger Two; public static readonly BigInteger Three; public static readonly BigInteger Ten; //private readonly static byte[] BitCountTable = //{ // 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 readonly static byte[] BitLengthTable = { 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 }; // TODO Parse radix-2 64 bits at a time and radix-8 63 bits at a time private const int chunk2 = 1, chunk8 = 1, chunk10 = 19, chunk16 = 16; private static readonly BigInteger radix2, radix2E, radix8, radix8E, radix10, radix10E, radix16, radix16E; private static readonly SecureRandom RandomSource = new SecureRandom(); /* * These are the threshold bit-lengths (of an exponent) where we increase the window size. * They are calculated according to the expected savings in multiplications. * Some squares will also be saved on average, but we offset these against the extra storage costs. */ private static readonly int[] ExpWindowThresholds = { 7, 25, 81, 241, 673, 1793, 4609, Int32.MaxValue }; private const int BitsPerByte = 8; private const int BitsPerInt = 32; private const int BytesPerInt = 4; static BigInteger() { Zero = new BigInteger(0, ZeroMagnitude, false); Zero.nBits = 0; Zero.nBitLength = 0; SMALL_CONSTANTS[0] = Zero; for (uint i = 1; i < SMALL_CONSTANTS.Length; ++i) { SMALL_CONSTANTS[i] = CreateUValueOf(i); } One = SMALL_CONSTANTS[1]; Two = SMALL_CONSTANTS[2]; Three = SMALL_CONSTANTS[3]; Ten = SMALL_CONSTANTS[10]; radix2 = ValueOf(2); radix2E = radix2.Pow(chunk2); radix8 = ValueOf(8); radix8E = radix8.Pow(chunk8); radix10 = ValueOf(10); radix10E = radix10.Pow(chunk10); radix16 = ValueOf(16); radix16E = radix16.Pow(chunk16); primeProducts = new int[primeLists.Length]; for (int i = 0; i < primeLists.Length; ++i) { int[] primeList = primeLists[i]; int product = primeList[0]; for (int j = 1; j < primeList.Length; ++j) { product *= primeList[j]; } primeProducts[i] = product; } } private int[] magnitude; // array of ints with [0] being the most significant private int sign; // -1 means -ve; +1 means +ve; 0 means 0; private int nBits = -1; // cache BitCount() value private int nBitLength = -1; // cache BitLength() value private int mQuote = 0; // -m^(-1) mod b, b = 2^32 (see Montgomery mult.), 0 when uninitialised private static int GetByteLength( int nBits) { return (nBits + BitsPerByte - 1) / BitsPerByte; } internal static BigInteger Arbitrary(int sizeInBits) { return new BigInteger(sizeInBits, RandomSource); } private BigInteger( int signum, int[] mag, bool checkMag) { if (checkMag) { int i = 0; while (i < mag.Length && mag[i] == 0) { ++i; } if (i == mag.Length) { this.sign = 0; this.magnitude = ZeroMagnitude; } else { this.sign = signum; if (i == 0) { this.magnitude = mag; } else { // strip leading 0 words this.magnitude = new int[mag.Length - i]; Array.Copy(mag, i, this.magnitude, 0, this.magnitude.Length); } } } else { this.sign = signum; this.magnitude = mag; } } public BigInteger( string value) : this(value, 10) { } public BigInteger( string str, int radix) { if (str.Length == 0) throw new FormatException("Zero length BigInteger"); NumberStyles style; int chunk; BigInteger r; BigInteger rE; switch (radix) { case 2: // Is there anyway to restrict to binary digits? style = NumberStyles.Integer; chunk = chunk2; r = radix2; rE = radix2E; break; case 8: // Is there anyway to restrict to octal digits? style = NumberStyles.Integer; chunk = chunk8; r = radix8; rE = radix8E; break; case 10: // This style seems to handle spaces and minus sign already (our processing redundant?) style = NumberStyles.Integer; chunk = chunk10; r = radix10; rE = radix10E; break; case 16: // TODO Should this be HexNumber? style = NumberStyles.AllowHexSpecifier; chunk = chunk16; r = radix16; rE = radix16E; break; default: throw new FormatException("Only bases 2, 8, 10, or 16 allowed"); } int index = 0; sign = 1; if (str[0] == '-') { if (str.Length == 1) throw new FormatException("Zero length BigInteger"); sign = -1; index = 1; } // strip leading zeros from the string str while (index < str.Length && Int32.Parse(str[index].ToString(), style) == 0) { index++; } if (index >= str.Length) { // zero value - we're done sign = 0; magnitude = ZeroMagnitude; return; } ////// // could we work out the max number of ints required to store // str.Length digits in the given base, then allocate that // storage in one hit?, then Generate the magnitude in one hit too? ////// BigInteger b = Zero; int next = index + chunk; if (next <= str.Length) { do { string s = str.Substring(index, chunk); ulong i = ulong.Parse(s, style); BigInteger bi = CreateUValueOf(i); switch (radix) { case 2: // TODO Need this because we are parsing in radix 10 above if (i >= 2) throw new FormatException("Bad character in radix 2 string: " + s); // TODO Parse 64 bits at a time b = b.ShiftLeft(1); break; case 8: // TODO Need this because we are parsing in radix 10 above if (i >= 8) throw new FormatException("Bad character in radix 8 string: " + s); // TODO Parse 63 bits at a time b = b.ShiftLeft(3); break; case 16: b = b.ShiftLeft(64); break; default: b = b.Multiply(rE); break; } b = b.Add(bi); index = next; next += chunk; } while (next <= str.Length); } if (index < str.Length) { string s = str.Substring(index); ulong i = ulong.Parse(s, style); BigInteger bi = CreateUValueOf(i); if (b.sign > 0) { if (radix == 2) { // NB: Can't reach here since we are parsing one char at a time Debug.Assert(false); // TODO Parse all bits at once // b = b.ShiftLeft(s.Length); } else if (radix == 8) { // NB: Can't reach here since we are parsing one char at a time Debug.Assert(false); // TODO Parse all bits at once // b = b.ShiftLeft(s.Length * 3); } else if (radix == 16) { b = b.ShiftLeft(s.Length << 2); } else { b = b.Multiply(r.Pow(s.Length)); } b = b.Add(bi); } else { b = bi; } } // Note: This is the previous (slower) algorithm // while (index < value.Length) // { // char c = value[index]; // string s = c.ToString(); // int i = Int32.Parse(s, style); // // b = b.Multiply(r).Add(ValueOf(i)); // index++; // } magnitude = b.magnitude; } public BigInteger( byte[] bytes) : this(bytes, 0, bytes.Length) { } public BigInteger( byte[] bytes, int offset, int length) { if (length == 0) throw new FormatException("Zero length BigInteger"); // TODO Move this processing into MakeMagnitude (provide sign argument) if ((sbyte)bytes[offset] < 0) { this.sign = -1; int end = offset + length; int iBval; // strip leading sign bytes for (iBval = offset; iBval < end && ((sbyte)bytes[iBval] == -1); iBval++) { } if (iBval >= end) { this.magnitude = One.magnitude; } else { int numBytes = end - iBval; byte[] inverse = new byte[numBytes]; int index = 0; while (index < numBytes) { inverse[index++] = (byte)~bytes[iBval++]; } Debug.Assert(iBval == end); while (inverse[--index] == byte.MaxValue) { inverse[index] = byte.MinValue; } inverse[index]++; this.magnitude = MakeMagnitude(inverse, 0, inverse.Length); } } else { // strip leading zero bytes and return magnitude bytes this.magnitude = MakeMagnitude(bytes, offset, length); this.sign = this.magnitude.Length > 0 ? 1 : 0; } } private static int[] MakeMagnitude( byte[] bytes, int offset, int length) { int end = offset + length; // strip leading zeros int firstSignificant; for (firstSignificant = offset; firstSignificant < end && bytes[firstSignificant] == 0; firstSignificant++) { } if (firstSignificant >= end) { return ZeroMagnitude; } int nInts = (end - firstSignificant + 3) / BytesPerInt; int bCount = (end - firstSignificant) % BytesPerInt; if (bCount == 0) { bCount = BytesPerInt; } if (nInts < 1) { return ZeroMagnitude; } int[] mag = new int[nInts]; int v = 0; int magnitudeIndex = 0; for (int i = firstSignificant; i < end; ++i) { v <<= 8; v |= bytes[i] & 0xff; bCount--; if (bCount <= 0) { mag[magnitudeIndex] = v; magnitudeIndex++; bCount = BytesPerInt; v = 0; } } if (magnitudeIndex < mag.Length) { mag[magnitudeIndex] = v; } return mag; } public BigInteger( int sign, byte[] bytes) : this(sign, bytes, 0, bytes.Length) { } public BigInteger( int sign, byte[] bytes, int offset, int length) { if (sign < -1 || sign > 1) throw new FormatException("Invalid sign value"); if (sign == 0) { this.sign = 0; this.magnitude = ZeroMagnitude; } else { // copy bytes this.magnitude = MakeMagnitude(bytes, offset, length); this.sign = this.magnitude.Length < 1 ? 0 : sign; } } public BigInteger( int sizeInBits, Random random) { if (sizeInBits < 0) throw new ArgumentException("sizeInBits must be non-negative"); this.nBits = -1; this.nBitLength = -1; if (sizeInBits == 0) { this.sign = 0; this.magnitude = ZeroMagnitude; return; } int nBytes = GetByteLength(sizeInBits); byte[] b = new byte[nBytes]; random.NextBytes(b); // strip off any excess bits in the MSB int xBits = BitsPerByte * nBytes - sizeInBits; b[0] &= (byte)(255U >> xBits); this.magnitude = MakeMagnitude(b, 0, b.Length); this.sign = this.magnitude.Length < 1 ? 0 : 1; } public BigInteger( int bitLength, int certainty, Random random) { if (bitLength < 2) throw new ArithmeticException("bitLength < 2"); this.sign = 1; this.nBitLength = bitLength; if (bitLength == 2) { this.magnitude = random.Next(2) == 0 ? Two.magnitude : Three.magnitude; return; } int nBytes = GetByteLength(bitLength); byte[] b = new byte[nBytes]; int xBits = BitsPerByte * nBytes - bitLength; byte mask = (byte)(255U >> xBits); byte lead = (byte)(1 << (7 - xBits)); for (;;) { random.NextBytes(b); // strip off any excess bits in the MSB b[0] &= mask; // ensure the leading bit is 1 (to meet the strength requirement) b[0] |= lead; // ensure the trailing bit is 1 (i.e. must be odd) b[nBytes - 1] |= 1; this.magnitude = MakeMagnitude(b, 0, b.Length); this.nBits = -1; this.mQuote = 0; if (certainty < 1) break; if (CheckProbablePrime(certainty, random, true)) break; for (int j = 1; j < (magnitude.Length - 1); ++j) { this.magnitude[j] ^= random.Next(); if (CheckProbablePrime(certainty, random, true)) return; } } } public BigInteger Abs() { return sign >= 0 ? this : Negate(); } /** * return a = a + b - b preserved. */ private static int[] AddMagnitudes( int[] a, int[] b) { int tI = a.Length - 1; int vI = b.Length - 1; long m = 0; while (vI >= 0) { m += ((long)(uint)a[tI] + (long)(uint)b[vI--]); a[tI--] = (int)m; m = (long)((ulong)m >> 32); } if (m != 0) { while (tI >= 0 && ++a[tI--] == 0) { } } return a; } public BigInteger Add( BigInteger value) { if (this.sign == 0) return value; if (this.sign != value.sign) { if (value.sign == 0) return this; if (value.sign < 0) return Subtract(value.Negate()); return value.Subtract(Negate()); } return AddToMagnitude(value.magnitude); } private BigInteger AddToMagnitude( int[] magToAdd) { int[] big, small; if (this.magnitude.Length < magToAdd.Length) { big = magToAdd; small = this.magnitude; } else { big = this.magnitude; small = magToAdd; } // Conservatively avoid over-allocation when no overflow possible uint limit = uint.MaxValue; if (big.Length == small.Length) limit -= (uint) small[0]; bool possibleOverflow = (uint) big[0] >= limit; int[] bigCopy; if (possibleOverflow) { bigCopy = new int[big.Length + 1]; big.CopyTo(bigCopy, 1); } else { bigCopy = (int[]) big.Clone(); } bigCopy = AddMagnitudes(bigCopy, small); return new BigInteger(this.sign, bigCopy, possibleOverflow); } public BigInteger And( BigInteger value) { if (this.sign == 0 || value.sign == 0) { return Zero; } int[] aMag = this.sign > 0 ? this.magnitude : Add(One).magnitude; int[] bMag = value.sign > 0 ? value.magnitude : value.Add(One).magnitude; bool resultNeg = sign < 0 && value.sign < 0; int resultLength = System.Math.Max(aMag.Length, bMag.Length); int[] resultMag = new int[resultLength]; int aStart = resultMag.Length - aMag.Length; int bStart = resultMag.Length - bMag.Length; for (int i = 0; i < resultMag.Length; ++i) { int aWord = i >= aStart ? aMag[i - aStart] : 0; int bWord = i >= bStart ? bMag[i - bStart] : 0; if (this.sign < 0) { aWord = ~aWord; } if (value.sign < 0) { bWord = ~bWord; } resultMag[i] = aWord & bWord; if (resultNeg) { resultMag[i] = ~resultMag[i]; } } BigInteger result = new BigInteger(1, resultMag, true); // TODO Optimise this case if (resultNeg) { result = result.Not(); } return result; } public BigInteger AndNot( BigInteger val) { return And(val.Not()); } public int BitCount { get { if (nBits == -1) { if (sign < 0) { // TODO Optimise this case nBits = Not().BitCount; } else { int sum = 0; for (int i = 0; i < magnitude.Length; ++i) { sum += BitCnt(magnitude[i]); } nBits = sum; } } return nBits; } } public static int BitCnt(int i) { uint u = (uint)i; u = u - ((u >> 1) & 0x55555555); u = (u & 0x33333333) + ((u >> 2) & 0x33333333); u = (u + (u >> 4)) & 0x0f0f0f0f; u += (u >> 8); u += (u >> 16); u &= 0x3f; return (int)u; } private static int CalcBitLength(int sign, int indx, int[] mag) { for (;;) { if (indx >= mag.Length) return 0; if (mag[indx] != 0) break; ++indx; } // bit length for everything after the first int int bitLength = 32 * ((mag.Length - indx) - 1); // and determine bitlength of first int int firstMag = mag[indx]; bitLength += BitLen(firstMag); // Check for negative powers of two if (sign < 0 && ((firstMag & -firstMag) == firstMag)) { do { if (++indx >= mag.Length) { --bitLength; break; } } while (mag[indx] == 0); } return bitLength; } public int BitLength { get { if (nBitLength == -1) { nBitLength = sign == 0 ? 0 : CalcBitLength(sign, 0, magnitude); } return nBitLength; } } // // BitLen(value) is the number of bits in value. // internal static int BitLen(int w) { uint v = (uint)w; uint t = v >> 24; if (t != 0) return 24 + BitLengthTable[t]; t = v >> 16; if (t != 0) return 16 + BitLengthTable[t]; t = v >> 8; if (t != 0) return 8 + BitLengthTable[t]; return BitLengthTable[v]; } private bool QuickPow2Check() { return sign > 0 && nBits == 1; } public int CompareTo( object obj) { return CompareTo((BigInteger)obj); } /** * unsigned comparison on two arrays - note the arrays may * start with leading zeros. */ private static 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++; } return CompareNoLeadingZeroes(xIndx, x, yIndx, y); } private static int CompareNoLeadingZeroes( int xIndx, int[] x, int yIndx, int[] y) { int diff = (x.Length - y.Length) - (xIndx - yIndx); if (diff != 0) { return diff < 0 ? -1 : 1; } // lengths of magnitudes the same, test the magnitude values while (xIndx < x.Length) { uint v1 = (uint)x[xIndx++]; uint v2 = (uint)y[yIndx++]; if (v1 != v2) return v1 < v2 ? -1 : 1; } return 0; } public int CompareTo( BigInteger value) { return sign < value.sign ? -1 : sign > value.sign ? 1 : sign == 0 ? 0 : sign * CompareNoLeadingZeroes(0, magnitude, 0, value.magnitude); } /** * return z = x / y - done in place (z value preserved, x contains the * remainder) */ private int[] Divide( int[] x, int[] y) { int xStart = 0; while (xStart < x.Length && x[xStart] == 0) { ++xStart; } int yStart = 0; while (yStart < y.Length && y[yStart] == 0) { ++yStart; } Debug.Assert(yStart < y.Length); int xyCmp = CompareNoLeadingZeroes(xStart, x, yStart, y); int[] count; if (xyCmp > 0) { int yBitLength = CalcBitLength(1, yStart, y); int xBitLength = CalcBitLength(1, xStart, x); int shift = xBitLength - yBitLength; int[] iCount; int iCountStart = 0; int[] c; int cStart = 0; int cBitLength = yBitLength; if (shift > 0) { // iCount = ShiftLeft(One.magnitude, shift); iCount = new int[(shift >> 5) + 1]; iCount[0] = 1 << (shift % 32); c = ShiftLeft(y, shift); cBitLength += shift; } else { iCount = new int[] { 1 }; int len = y.Length - yStart; c = new int[len]; Array.Copy(y, yStart, c, 0, len); } count = new int[iCount.Length]; for (;;) { if (cBitLength < xBitLength || CompareNoLeadingZeroes(xStart, x, cStart, c) >= 0) { Subtract(xStart, x, cStart, c); AddMagnitudes(count, iCount); while (x[xStart] == 0) { if (++xStart == x.Length) return count; } //xBitLength = CalcBitLength(xStart, x); xBitLength = 32 * (x.Length - xStart - 1) + BitLen(x[xStart]); if (xBitLength <= yBitLength) { if (xBitLength < yBitLength) return count; xyCmp = CompareNoLeadingZeroes(xStart, x, yStart, y); if (xyCmp <= 0) break; } } shift = cBitLength - xBitLength; // NB: The case where c[cStart] is 1-bit is harmless if (shift == 1) { uint firstC = (uint) c[cStart] >> 1; uint firstX = (uint) x[xStart]; if (firstC > firstX) ++shift; } if (shift < 2) { ShiftRightOneInPlace(cStart, c); --cBitLength; ShiftRightOneInPlace(iCountStart, iCount); } else { ShiftRightInPlace(cStart, c, shift); cBitLength -= shift; ShiftRightInPlace(iCountStart, iCount, shift); } //cStart = c.Length - ((cBitLength + 31) / 32); while (c[cStart] == 0) { ++cStart; } while (iCount[iCountStart] == 0) { ++iCountStart; } } } else { count = new int[1]; } if (xyCmp == 0) { AddMagnitudes(count, One.magnitude); Array.Clear(x, xStart, x.Length - xStart); } return count; } public BigInteger Divide( BigInteger val) { if (val.sign == 0) throw new ArithmeticException("Division by zero error"); if (sign == 0) return Zero; if (val.QuickPow2Check()) // val is power of two { BigInteger result = this.Abs().ShiftRight(val.Abs().BitLength - 1); return val.sign == this.sign ? result : result.Negate(); } int[] mag = (int[]) this.magnitude.Clone(); return new BigInteger(this.sign * val.sign, Divide(mag, val.magnitude), true); } public BigInteger[] DivideAndRemainder( BigInteger val) { if (val.sign == 0) throw new ArithmeticException("Division by zero error"); BigInteger[] biggies = new BigInteger[2]; if (sign == 0) { biggies[0] = Zero; biggies[1] = Zero; } else if (val.QuickPow2Check()) // val is power of two { int e = val.Abs().BitLength - 1; BigInteger quotient = this.Abs().ShiftRight(e); int[] remainder = this.LastNBits(e); biggies[0] = val.sign == this.sign ? quotient : quotient.Negate(); biggies[1] = new BigInteger(this.sign, remainder, true); } else { int[] remainder = (int[]) this.magnitude.Clone(); int[] quotient = Divide(remainder, val.magnitude); biggies[0] = new BigInteger(this.sign * val.sign, quotient, true); biggies[1] = new BigInteger(this.sign, remainder, true); } return biggies; } public override bool Equals( object obj) { if (obj == this) return true; BigInteger biggie = obj as BigInteger; if (biggie == null) return false; return sign == biggie.sign && IsEqualMagnitude(biggie); } private bool IsEqualMagnitude(BigInteger x) { //int[] xMag = x.magnitude; if (magnitude.Length != x.magnitude.Length) return false; for (int i = 0; i < magnitude.Length; i++) { if (magnitude[i] != x.magnitude[i]) return false; } return true; } public BigInteger Gcd( BigInteger value) { if (value.sign == 0) return Abs(); if (sign == 0) return value.Abs(); BigInteger r; BigInteger u = this; BigInteger v = value; while (v.sign != 0) { r = u.Mod(v); u = v; v = r; } return u; } public override int GetHashCode() { int hc = magnitude.Length; if (magnitude.Length > 0) { hc ^= magnitude[0]; if (magnitude.Length > 1) { hc ^= magnitude[magnitude.Length - 1]; } } return sign < 0 ? ~hc : hc; } // TODO Make public? private BigInteger Inc() { if (this.sign == 0) return One; if (this.sign < 0) return new BigInteger(-1, doSubBigLil(this.magnitude, One.magnitude), true); return AddToMagnitude(One.magnitude); } public int IntValue { get { if (sign == 0) return 0; int n = magnitude.Length; int v = magnitude[n - 1]; return sign < 0 ? -v : v; } } /** * 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) { return IsProbablePrime(certainty, false); } internal bool IsProbablePrime(int certainty, bool randomlySelected) { if (certainty <= 0) return true; BigInteger n = Abs(); if (!n.TestBit(0)) return n.Equals(Two); if (n.Equals(One)) return false; return n.CheckProbablePrime(certainty, RandomSource, randomlySelected); } private bool CheckProbablePrime(int certainty, Random random, bool randomlySelected) { Debug.Assert(certainty > 0); Debug.Assert(CompareTo(Two) > 0); Debug.Assert(TestBit(0)); // Try to reduce the penalty for really small numbers int numLists = System.Math.Min(BitLength - 1, primeLists.Length); for (int i = 0; i < numLists; ++i) { int test = Remainder(primeProducts[i]); int[] primeList = primeLists[i]; for (int j = 0; j < primeList.Length; ++j) { int prime = primeList[j]; int qRem = test % prime; if (qRem == 0) { // We may find small numbers in the list return BitLength < 16 && IntValue == prime; } } } // TODO Special case for < 10^16 (RabinMiller fixed list) // if (BitLength < 30) // { // RabinMiller against 2, 3, 5, 7, 11, 13, 23 is sufficient // } // TODO Is it worth trying to create a hybrid of these two? return RabinMillerTest(certainty, random, randomlySelected); // return SolovayStrassenTest(certainty, random); // bool rbTest = RabinMillerTest(certainty, random); // bool ssTest = SolovayStrassenTest(certainty, random); // // Debug.Assert(rbTest == ssTest); // // return rbTest; } public bool RabinMillerTest(int certainty, Random random) { return RabinMillerTest(certainty, random, false); } internal bool RabinMillerTest(int certainty, Random random, bool randomlySelected) { int bits = BitLength; Debug.Assert(certainty > 0); Debug.Assert(bits > 2); Debug.Assert(TestBit(0)); int iterations = ((certainty - 1) / 2) + 1; if (randomlySelected) { int itersFor100Cert = bits >= 1024 ? 4 : bits >= 512 ? 8 : bits >= 256 ? 16 : 50; if (certainty < 100) { iterations = System.Math.Min(itersFor100Cert, iterations); } else { iterations -= 50; iterations += itersFor100Cert; } } // let n = 1 + d . 2^s BigInteger n = this; int s = n.GetLowestSetBitMaskFirst(-1 << 1); Debug.Assert(s >= 1); BigInteger r = n.ShiftRight(s); // NOTE: Avoid conversion to/from Montgomery form and check for R/-R as result instead BigInteger montRadix = One.ShiftLeft(32 * n.magnitude.Length).Remainder(n); BigInteger minusMontRadix = n.Subtract(montRadix); do { BigInteger a; do { a = new BigInteger(n.BitLength, random); } while (a.sign == 0 || a.CompareTo(n) >= 0 || a.IsEqualMagnitude(montRadix) || a.IsEqualMagnitude(minusMontRadix)); BigInteger y = ModPowMonty(a, r, n, false); if (!y.Equals(montRadix)) { int j = 0; while (!y.Equals(minusMontRadix)) { if (++j == s) return false; y = ModPowMonty(y, Two, n, false); if (y.Equals(montRadix)) return false; } } } while (--iterations > 0); return true; } // private bool SolovayStrassenTest( // int certainty, // Random random) // { // Debug.Assert(certainty > 0); // Debug.Assert(CompareTo(Two) > 0); // Debug.Assert(TestBit(0)); // // BigInteger n = this; // BigInteger nMinusOne = n.Subtract(One); // BigInteger e = nMinusOne.ShiftRight(1); // // do // { // BigInteger a; // do // { // a = new BigInteger(nBitLength, random); // } // // NB: Spec says 0 < x < n, but 1 is trivial // while (a.CompareTo(One) <= 0 || a.CompareTo(n) >= 0); // // // // TODO Check this is redundant given the way Jacobi() works? //// if (!a.Gcd(n).Equals(One)) //// return false; // // int x = Jacobi(a, n); // // if (x == 0) // return false; // // BigInteger check = a.ModPow(e, n); // // if (x == 1 && !check.Equals(One)) // return false; // // if (x == -1 && !check.Equals(nMinusOne)) // return false; // // --certainty; // } // while (certainty > 0); // // return true; // } // // private static int Jacobi( // BigInteger a, // BigInteger b) // { // Debug.Assert(a.sign >= 0); // Debug.Assert(b.sign > 0); // Debug.Assert(b.TestBit(0)); // Debug.Assert(a.CompareTo(b) < 0); // // int totalS = 1; // for (;;) // { // if (a.sign == 0) // return 0; // // if (a.Equals(One)) // break; // // int e = a.GetLowestSetBit(); // // int bLsw = b.magnitude[b.magnitude.Length - 1]; // if ((e & 1) != 0 && ((bLsw & 7) == 3 || (bLsw & 7) == 5)) // totalS = -totalS; // // // TODO Confirm this is faster than later a1.Equals(One) test // if (a.BitLength == e + 1) // break; // BigInteger a1 = a.ShiftRight(e); //// if (a1.Equals(One)) //// break; // // int a1Lsw = a1.magnitude[a1.magnitude.Length - 1]; // if ((bLsw & 3) == 3 && (a1Lsw & 3) == 3) // totalS = -totalS; // //// a = b.Mod(a1); // a = b.Remainder(a1); // b = a1; // } // return totalS; // } public long LongValue { get { if (sign == 0) return 0; int n = magnitude.Length; long v = magnitude[n - 1] & IMASK; if (n > 1) { v |= (magnitude[n - 2] & IMASK) << 32; } return sign < 0 ? -v : v; } } public BigInteger Max( BigInteger value) { return CompareTo(value) > 0 ? this : value; } public BigInteger Min( BigInteger value) { return CompareTo(value) < 0 ? this : value; } public BigInteger Mod( BigInteger m) { if (m.sign < 1) throw new ArithmeticException("Modulus must be positive"); BigInteger biggie = Remainder(m); return (biggie.sign >= 0 ? biggie : biggie.Add(m)); } public BigInteger ModInverse( BigInteger m) { if (m.sign < 1) throw new ArithmeticException("Modulus must be positive"); // TODO Too slow at the moment // // "Fast Key Exchange with Elliptic Curve Systems" R.Schoeppel // if (m.TestBit(0)) // { // //The Almost Inverse Algorithm // int k = 0; // BigInteger B = One, C = Zero, F = this, G = m, tmp; // // for (;;) // { // // While F is even, do F=F/u, C=C*u, k=k+1. // int zeroes = F.GetLowestSetBit(); // if (zeroes > 0) // { // F = F.ShiftRight(zeroes); // C = C.ShiftLeft(zeroes); // k += zeroes; // } // // // If F = 1, then return B,k. // if (F.Equals(One)) // { // BigInteger half = m.Add(One).ShiftRight(1); // BigInteger halfK = half.ModPow(BigInteger.ValueOf(k), m); // return B.Multiply(halfK).Mod(m); // } // // if (F.CompareTo(G) < 0) // { // tmp = G; G = F; F = tmp; // tmp = B; B = C; C = tmp; // } // // F = F.Add(G); // B = B.Add(C); // } // } if (m.QuickPow2Check()) { return ModInversePow2(m); } BigInteger d = this.Remainder(m); BigInteger x; BigInteger gcd = ExtEuclid(d, m, out x); if (!gcd.Equals(One)) throw new ArithmeticException("Numbers not relatively prime."); if (x.sign < 0) { x = x.Add(m); } return x; } private BigInteger ModInversePow2(BigInteger m) { Debug.Assert(m.SignValue > 0); Debug.Assert(m.BitCount == 1); if (!TestBit(0)) { throw new ArithmeticException("Numbers not relatively prime."); } int pow = m.BitLength - 1; long inv64 = ModInverse64(LongValue); if (pow < 64) { inv64 &= ((1L << pow) - 1); } BigInteger x = BigInteger.ValueOf(inv64); if (pow > 64) { BigInteger d = this.Remainder(m); int bitsCorrect = 64; do { BigInteger t = x.Multiply(d).Remainder(m); x = x.Multiply(Two.Subtract(t)).Remainder(m); bitsCorrect <<= 1; } while (bitsCorrect < pow); } if (x.sign < 0) { x = x.Add(m); } return x; } private static int ModInverse32(int d) { // Newton's method with initial estimate "correct to 4 bits" Debug.Assert((d & 1) != 0); int x = d + (((d + 1) & 4) << 1); // d.x == 1 mod 2**4 Debug.Assert(((d * x) & 15) == 1); x *= 2 - d * x; // d.x == 1 mod 2**8 x *= 2 - d * x; // d.x == 1 mod 2**16 x *= 2 - d * x; // d.x == 1 mod 2**32 Debug.Assert(d * x == 1); return x; } private static long ModInverse64(long d) { // Newton's method with initial estimate "correct to 4 bits" Debug.Assert((d & 1L) != 0); long x = d + (((d + 1L) & 4L) << 1); // d.x == 1 mod 2**4 Debug.Assert(((d * x) & 15L) == 1L); x *= 2 - d * x; // d.x == 1 mod 2**8 x *= 2 - d * x; // d.x == 1 mod 2**16 x *= 2 - d * x; // d.x == 1 mod 2**32 x *= 2 - d * x; // d.x == 1 mod 2**64 Debug.Assert(d * x == 1L); 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 * @return The greatest common divisor of a and b */ private static BigInteger ExtEuclid(BigInteger a, BigInteger b, out BigInteger u1Out) { BigInteger u1 = One, v1 = Zero; BigInteger u3 = a, v3 = b; if (v3.sign > 0) { for (;;) { BigInteger[] q = u3.DivideAndRemainder(v3); u3 = v3; v3 = q[1]; BigInteger oldU1 = u1; u1 = v1; if (v3.sign <= 0) break; v1 = oldU1.Subtract(v1.Multiply(q[0])); } } u1Out = u1; return u3; } private static void ZeroOut( int[] x) { Array.Clear(x, 0, x.Length); } public BigInteger ModPow(BigInteger e, BigInteger m) { if (m.sign < 1) throw new ArithmeticException("Modulus must be positive"); if (m.Equals(One)) return Zero; if (e.sign == 0) return One; if (sign == 0) return Zero; bool negExp = e.sign < 0; if (negExp) e = e.Negate(); BigInteger result = this.Mod(m); if (!e.Equals(One)) { if ((m.magnitude[m.magnitude.Length - 1] & 1) == 0) { result = ModPowBarrett(result, e, m); } else { result = ModPowMonty(result, e, m, true); } } if (negExp) result = result.ModInverse(m); return result; } private static BigInteger ModPowBarrett(BigInteger b, BigInteger e, BigInteger m) { int k = m.magnitude.Length; BigInteger mr = One.ShiftLeft((k + 1) << 5); BigInteger yu = One.ShiftLeft(k << 6).Divide(m); // Sliding window from MSW to LSW int extraBits = 0, expLength = e.BitLength; while (expLength > ExpWindowThresholds[extraBits]) { ++extraBits; } int numPowers = 1 << extraBits; BigInteger[] oddPowers = new BigInteger[numPowers]; oddPowers[0] = b; BigInteger b2 = ReduceBarrett(b.Square(), m, mr, yu); for (int i = 1; i < numPowers; ++i) { oddPowers[i] = ReduceBarrett(oddPowers[i - 1].Multiply(b2), m, mr, yu); } int[] windowList = GetWindowList(e.magnitude, extraBits); Debug.Assert(windowList.Length > 0); int window = windowList[0]; int mult = window & 0xFF, lastZeroes = window >> 8; BigInteger y; if (mult == 1) { y = b2; --lastZeroes; } else { y = oddPowers[mult >> 1]; } int windowPos = 1; while ((window = windowList[windowPos++]) != -1) { mult = window & 0xFF; int bits = lastZeroes + BitLengthTable[mult]; for (int j = 0; j < bits; ++j) { y = ReduceBarrett(y.Square(), m, mr, yu); } y = ReduceBarrett(y.Multiply(oddPowers[mult >> 1]), m, mr, yu); lastZeroes = window >> 8; } for (int i = 0; i < lastZeroes; ++i) { y = ReduceBarrett(y.Square(), m, mr, yu); } return y; } private static BigInteger ReduceBarrett(BigInteger x, BigInteger m, BigInteger mr, BigInteger yu) { int xLen = x.BitLength, mLen = m.BitLength; if (xLen < mLen) return x; if (xLen - mLen > 1) { int k = m.magnitude.Length; BigInteger q1 = x.DivideWords(k - 1); BigInteger q2 = q1.Multiply(yu); // TODO Only need partial multiplication here BigInteger q3 = q2.DivideWords(k + 1); BigInteger r1 = x.RemainderWords(k + 1); BigInteger r2 = q3.Multiply(m); // TODO Only need partial multiplication here BigInteger r3 = r2.RemainderWords(k + 1); x = r1.Subtract(r3); if (x.sign < 0) { x = x.Add(mr); } } while (x.CompareTo(m) >= 0) { x = x.Subtract(m); } return x; } private static BigInteger ModPowMonty(BigInteger b, BigInteger e, BigInteger m, bool convert) { int n = m.magnitude.Length; int powR = 32 * n; bool smallMontyModulus = m.BitLength + 2 <= powR; uint mDash = (uint)m.GetMQuote(); // tmp = this * R mod m if (convert) { b = b.ShiftLeft(powR).Remainder(m); } int[] yAccum = new int[n + 1]; int[] zVal = b.magnitude; Debug.Assert(zVal.Length <= n); if (zVal.Length < n) { int[] tmp = new int[n]; zVal.CopyTo(tmp, n - zVal.Length); zVal = tmp; } // Sliding window from MSW to LSW int extraBits = 0; // Filter the common case of small RSA exponents with few bits set if (e.magnitude.Length > 1 || e.BitCount > 2) { int expLength = e.BitLength; while (expLength > ExpWindowThresholds[extraBits]) { ++extraBits; } } int numPowers = 1 << extraBits; int[][] oddPowers = new int[numPowers][]; oddPowers[0] = zVal; int[] zSquared = Arrays.Clone(zVal); SquareMonty(yAccum, zSquared, m.magnitude, mDash, smallMontyModulus); for (int i = 1; i < numPowers; ++i) { oddPowers[i] = Arrays.Clone(oddPowers[i - 1]); MultiplyMonty(yAccum, oddPowers[i], zSquared, m.magnitude, mDash, smallMontyModulus); } int[] windowList = GetWindowList(e.magnitude, extraBits); Debug.Assert(windowList.Length > 1); int window = windowList[0]; int mult = window & 0xFF, lastZeroes = window >> 8; int[] yVal; if (mult == 1) { yVal = zSquared; --lastZeroes; } else { yVal = Arrays.Clone(oddPowers[mult >> 1]); } int windowPos = 1; while ((window = windowList[windowPos++]) != -1) { mult = window & 0xFF; int bits = lastZeroes + BitLengthTable[mult]; for (int j = 0; j < bits; ++j) { SquareMonty(yAccum, yVal, m.magnitude, mDash, smallMontyModulus); } MultiplyMonty(yAccum, yVal, oddPowers[mult >> 1], m.magnitude, mDash, smallMontyModulus); lastZeroes = window >> 8; } for (int i = 0; i < lastZeroes; ++i) { SquareMonty(yAccum, yVal, m.magnitude, mDash, smallMontyModulus); } if (convert) { // Return y * R^(-1) mod m MontgomeryReduce(yVal, m.magnitude, mDash); } else if (smallMontyModulus && CompareTo(0, yVal, 0, m.magnitude) >= 0) { Subtract(0, yVal, 0, m.magnitude); } return new BigInteger(1, yVal, true); } private static int[] GetWindowList(int[] mag, int extraBits) { int v = mag[0]; Debug.Assert(v != 0); int leadingBits = BitLen(v); int resultSize = (((mag.Length - 1) << 5) + leadingBits) / (1 + extraBits) + 2; int[] result = new int[resultSize]; int resultPos = 0; int bitPos = 33 - leadingBits; v <<= bitPos; int mult = 1, multLimit = 1 << extraBits; int zeroes = 0; int i = 0; for (; ; ) { for (; bitPos < 32; ++bitPos) { if (mult < multLimit) { mult = (mult << 1) | (int)((uint)v >> 31); } else if (v < 0) { result[resultPos++] = CreateWindowEntry(mult, zeroes); mult = 1; zeroes = 0; } else { ++zeroes; } v <<= 1; } if (++i == mag.Length) { result[resultPos++] = CreateWindowEntry(mult, zeroes); break; } v = mag[i]; bitPos = 0; } result[resultPos] = -1; return result; } private static int CreateWindowEntry(int mult, int zeroes) { while ((mult & 1) == 0) { mult >>= 1; ++zeroes; } return mult | (zeroes << 8); } /** * return w with w = x * x - w is assumed to have enough space. */ private static int[] Square( int[] w, int[] x) { // Note: this method allows w to be only (2 * x.Length - 1) words if result will fit // if (w.Length != 2 * x.Length) // throw new ArgumentException("no I don't think so..."); ulong c; int wBase = w.Length - 1; for (int i = x.Length - 1; i > 0; --i) { ulong v = (uint)x[i]; c = v * v + (uint)w[wBase]; w[wBase] = (int)c; c >>= 32; for (int j = i - 1; j >= 0; --j) { ulong prod = v * (uint)x[j]; c += ((uint)w[--wBase] & UIMASK) + ((uint)prod << 1); w[wBase] = (int)c; c = (c >> 32) + (prod >> 31); } c += (uint)w[--wBase]; w[wBase] = (int)c; if (--wBase >= 0) { w[wBase] = (int)(c >> 32); } else { Debug.Assert((c >> 32) == 0); } wBase += i; } c = (uint)x[0]; c = c * c + (uint)w[wBase]; w[wBase] = (int)c; if (--wBase >= 0) { w[wBase] += (int)(c >> 32); } else { Debug.Assert((c >> 32) == 0); } return w; } /** * return x with x = y * z - x is assumed to have enough space. */ private static int[] Multiply(int[] x, int[] y, int[] z) { int i = z.Length; if (i < 1) return x; int xBase = x.Length - y.Length; do { long a = z[--i] & IMASK; long val = 0; if (a != 0) { for (int j = y.Length - 1; j >= 0; j--) { val += a * (y[j] & IMASK) + (x[xBase + j] & IMASK); x[xBase + j] = (int)val; val = (long)((ulong)val >> 32); } } --xBase; if (xBase >= 0) { x[xBase] = (int)val; } else { Debug.Assert(val == 0); } } while (i > 0); return x; } /** * Calculate mQuote = -m^(-1) mod b with b = 2^32 (32 = word size) */ private int GetMQuote() { if (mQuote != 0) { return mQuote; // already calculated } Debug.Assert(this.sign > 0); int d = -magnitude[magnitude.Length - 1]; Debug.Assert((d & 1) != 0); return mQuote = ModInverse32(d); } private static void MontgomeryReduce(int[] x, int[] m, uint mDash) // mDash = -m^(-1) mod b { // NOTE: Not a general purpose reduction (which would allow x up to twice the bitlength of m) Debug.Assert(x.Length == m.Length); int n = m.Length; for (int i = n - 1; i >= 0; --i) { uint x0 = (uint)x[n - 1]; ulong t = x0 * mDash; ulong carry = t * (uint)m[n - 1] + x0; Debug.Assert((uint)carry == 0); carry >>= 32; for (int j = n - 2; j >= 0; --j) { carry += t * (uint)m[j] + (uint)x[j]; x[j + 1] = (int)carry; carry >>= 32; } x[0] = (int)carry; Debug.Assert(carry >> 32 == 0); } if (CompareTo(0, x, 0, m) >= 0) { Subtract(0, x, 0, m); } } /** * Montgomery multiplication: a = x * y * R^(-1) mod m *