diff options
Diffstat (limited to 'crypto/src/math')
| -rw-r--r-- | crypto/src/math/BigInteger.cs | 3575 | ||||
| -rw-r--r-- | crypto/src/math/ec/ECAlgorithms.cs | 93 | ||||
| -rw-r--r-- | crypto/src/math/ec/ECCurve.cs | 651 | ||||
| -rw-r--r-- | crypto/src/math/ec/ECFieldElement.cs | 1253 | ||||
| -rw-r--r-- | crypto/src/math/ec/ECPoint.cs | 572 | ||||
| -rw-r--r-- | crypto/src/math/ec/IntArray.cs | 485 | ||||
| -rw-r--r-- | crypto/src/math/ec/abc/SimpleBigDecimal.cs | 241 | ||||
| -rw-r--r-- | crypto/src/math/ec/abc/Tnaf.cs | 834 | ||||
| -rw-r--r-- | crypto/src/math/ec/abc/ZTauElement.cs | 36 | ||||
| -rw-r--r-- | crypto/src/math/ec/multiplier/ECMultiplier.cs | 18 | ||||
| -rw-r--r-- | crypto/src/math/ec/multiplier/FpNafMultiplier.cs | 39 | ||||
| -rw-r--r-- | crypto/src/math/ec/multiplier/PreCompInfo.cs | 11 | ||||
| -rw-r--r-- | crypto/src/math/ec/multiplier/ReferenceMultiplier.cs | 30 | ||||
| -rw-r--r-- | crypto/src/math/ec/multiplier/WNafMultiplier.cs | 241 | ||||
| -rw-r--r-- | crypto/src/math/ec/multiplier/WNafPreCompInfo.cs | 46 | ||||
| -rw-r--r-- | crypto/src/math/ec/multiplier/WTauNafMultiplier.cs | 120 | ||||
| -rw-r--r-- | crypto/src/math/ec/multiplier/WTauNafPreCompInfo.cs | 41 |
17 files changed, 8286 insertions, 0 deletions
diff --git a/crypto/src/math/BigInteger.cs b/crypto/src/math/BigInteger.cs new file mode 100644 index 000000000..00f8f399d --- /dev/null +++ b/crypto/src/math/BigInteger.cs @@ -0,0 +1,3575 @@ +using System; +using System.Collections; +using System.Diagnostics; +using System.Globalization; +using System.Text; + +using Org.BouncyCastle.Utilities; + +namespace Org.BouncyCastle.Math +{ +#if !(NETCF_1_0 || NETCF_2_0 || SILVERLIGHT) + [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 Random RandomSource = new Random(); + + /* + * 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; + } + + 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); + + 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] |= (byte)(1 << (7 - xBits)); + + // 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)) + break; + + if (bitLength > 32) + { + for (int rep = 0; rep < 10000; ++rep) + { + int n = 33 + random.Next(bitLength - 2); + this.magnitude[this.magnitude.Length - (n >> 5)] ^= (1 << (n & 31)); + this.magnitude[this.magnitude.Length - 1] ^= ((random.Next() + 1) << 1); + this.mQuote = 0; + + if (CheckProbablePrime(certainty, random)) + 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. + // + private 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. + * <p>From Knuth Vol 2, pg 395.</p> + */ + public bool IsProbablePrime( + int certainty) + { + 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); + } + + private bool CheckProbablePrime( + int certainty, + Random random) + { + 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); +// 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) + { + Debug.Assert(certainty > 0); + Debug.Assert(BitLength > 2); + Debug.Assert(TestBit(0)); + + // 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; + } + } + + certainty -= 2; // composites pass for only 1/4 possible 'a' + } + while (certainty > 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 + * @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, + out BigInteger u1Out) + //BigInteger u2Out) + { + BigInteger u1 = One; + BigInteger u3 = a; + BigInteger v1 = Zero; + BigInteger v3 = b; + + while (v3.sign > 0) + { + BigInteger[] q = u3.DivideAndRemainder(v3); + + BigInteger tmp = v1.Multiply(q[0]); + BigInteger tn = u1.Subtract(tmp); + u1 = v1; + v1 = tn; + + u3 = v3; + v3 = q[1]; + } + + //if (u1Out != null) + //{ + // u1Out.sign = u1.sign; + // u1Out.magnitude = u1.magnitude; + //} + u1Out = u1; + + //if (u2Out != null) + //{ + // BigInteger tmp = u1.Multiply(a); + // tmp = u3.Subtract(tmp); + // BigInteger res = tmp.Divide(b); + // u2Out.sign = res.sign; + // u2Out.magnitude = res.magnitude; + //} + + 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 + * <br/> + * Based algorithm 14.36 of Handbook of Applied Cryptography. + * <br/> + * <li> m, x, y should have length n </li> + * <li> a should have length (n + 1) </li> + * <li> b = 2^32, R = b^n </li> + * <br/> + * The result is put in x + * <br/> + * NOTE: the indices of x, y, m, a different in HAC and in Java + */ + private static void MultiplyMonty(int[] a, int[] x, int[] y, int[] m, uint mDash, bool smallMontyModulus) + // mDash = -m^(-1) mod b + { + int n = m.Length; + + if (n == 1) + { + x[0] = (int)MultiplyMontyNIsOne((uint)x[0], (uint)y[0], (uint)m[0], mDash); + return; + } + + uint y0 = (uint)y[n - 1]; + + { + ulong xi = (uint)x[n - 1]; + + ulong carry = xi * y0; + ulong t = (uint)carry * mDash; + + ulong prod2 = t * (uint)m[n - 1]; + carry += (uint)prod2; + Debug.Assert((uint)carry == 0); + carry = (carry >> 32) + (prod2 >> 32); + + for (int j = n - 2; j >= 0; --j) + { + ulong prod1 = xi * (uint)y[j]; + prod2 = t * (uint)m[j]; + + carry += (prod1 & UIMASK) + (uint)prod2; + a[j + 2] = (int)carry; + carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32); + } + + a[1] = (int)carry; + a[0] = (int)(carry >> 32); + } + + for (int i = n - 2; i >= 0; --i) + { + uint a0 = (uint)a[n]; + ulong xi = (uint)x[i]; + + ulong prod1 = xi * y0; + ulong carry = (prod1 & UIMASK) + a0; + ulong t = (uint)carry * mDash; + + ulong prod2 = t * (uint)m[n - 1]; + carry += (uint)prod2; + Debug.Assert((uint)carry == 0); + carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32); + + for (int j = n - 2; j >= 0; --j) + { + prod1 = xi * (uint)y[j]; + prod2 = t * (uint)m[j]; + + carry += (prod1 & UIMASK) + (uint)prod2 + (uint)a[j + 1]; + a[j + 2] = (int)carry; + carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32); + } + + carry += (uint)a[0]; + a[1] = (int)carry; + a[0] = (int)(carry >> 32); + } + + if (!smallMontyModulus && CompareTo(0, a, 0, m) >= 0) + { + Subtract(0, a, 0, m); + } + + Array.Copy(a, 1, x, 0, n); + } + + private static void SquareMonty(int[] a, int[] x, int[] m, uint mDash, bool smallMontyModulus) + // mDash = -m^(-1) mod b + { + int n = m.Length; + + if (n == 1) + { + uint xVal = (uint)x[0]; + x[0] = (int)MultiplyMontyNIsOne(xVal, xVal, (uint)m[0], mDash); + return; + } + + ulong x0 = (uint)x[n - 1]; + + { + ulong carry = x0 * x0; + ulong t = (uint)carry * mDash; + + ulong prod2 = t * (uint)m[n - 1]; + carry += (uint)prod2; + Debug.Assert((uint)carry == 0); + carry = (carry >> 32) + (prod2 >> 32); + + for (int j = n - 2; j >= 0; --j) + { + ulong prod1 = x0 * (uint)x[j]; + prod2 = t * (uint)m[j]; + + carry += (prod2 & UIMASK) + ((uint)prod1 << 1); + a[j + 2] = (int)carry; + carry = (carry >> 32) + (prod1 >> 31) + (prod2 >> 32); + } + + a[1] = (int)carry; + a[0] = (int)(carry >> 32); + } + + for (int i = n - 2; i >= 0; --i) + { + uint a0 = (uint)a[n]; + ulong t = a0 * mDash; + + ulong carry = t * (uint)m[n - 1] + a0; + Debug.Assert((uint)carry == 0); + carry >>= 32; + + for (int j = n - 2; j > i; --j) + { + carry += t * (uint)m[j] + (uint)a[j + 1]; + a[j + 2] = (int)carry; + carry >>= 32; + } + + ulong xi = (uint)x[i]; + + { + ulong prod1 = xi * xi; + ulong prod2 = t * (uint)m[i]; + + carry += (prod1 & UIMASK) + (uint)prod2 + (uint)a[i + 1]; + a[i + 2] = (int)carry; + carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32); + } + + for (int j = i - 1; j >= 0; --j) + { + ulong prod1 = xi * (uint)x[j]; + ulong prod2 = t * (uint)m[j]; + + carry += (prod2 & UIMASK) + ((uint)prod1 << 1) + (uint)a[j + 1]; + a[j + 2] = (int)carry; + carry = (carry >> 32) + (prod1 >> 31) + (prod2 >> 32); + } + + carry += (uint)a[0]; + a[1] = (int)carry; + a[0] = (int)(carry >> 32); + } + + if (!smallMontyModulus && CompareTo(0, a, 0, m) >= 0) + { + Subtract(0, a, 0, m); + } + + Array.Copy(a, 1, x, 0, n); + } + + private static uint MultiplyMontyNIsOne(uint x, uint y, uint m, uint mDash) + { + ulong carry = (ulong)x * y; + uint t = (uint)carry * mDash; + ulong um = m; + ulong prod2 = um * t; + carry += (uint)prod2; + Debug.Assert((uint)carry == 0); + carry = (carry >> 32) + (prod2 >> 32); + if (carry > um) + { + carry -= um; + } + Debug.Assert(carry < um); + return (uint)carry; + } + + public BigInteger Multiply( + BigInteger val) + { + if (val == this) + return Square(); + + if ((sign & val.sign) == 0) + return Zero; + + if (val.QuickPow2Check()) // val is power of two + { + BigInteger result = this.ShiftLeft(val.Abs().BitLength - 1); + return val.sign > 0 ? result : result.Negate(); + } + + if (this.QuickPow2Check()) // this is power of two + { + BigInteger result = val.ShiftLeft(this.Abs().BitLength - 1); + return this.sign > 0 ? result : result.Negate(); + } + + int resLength = magnitude.Length + val.magnitude.Length; + int[] res = new int[resLength]; + + Multiply(res, this.magnitude, val.magnitude); + + int resSign = sign ^ val.sign ^ 1; + return new BigInteger(resSign, res, true); + } + + public BigInteger Square() + { + if (sign == 0) + return Zero; + if (this.QuickPow2Check()) + return ShiftLeft(Abs().BitLength - 1); + int resLength = magnitude.Length << 1; + if ((uint)magnitude[0] >> 16 == 0) + --resLength; + int[] res = new int[resLength]; + Square(res, magnitude); + return new BigInteger(1, res, false); + } + + public BigInteger Negate() + { + if (sign == 0) + return this; + + return new BigInteger(-sign, magnitude, false); + } + + public BigInteger NextProbablePrime() + { + if (sign < 0) + throw new ArithmeticException("Cannot be called on value < 0"); + + if (CompareTo(Two) < 0) + return Two; + + BigInteger n = Inc().SetBit(0); + + while (!n.CheckProbablePrime(100, RandomSource)) + { + n = n.Add(Two); + } + + return n; + } + + public BigInteger Not() + { + return Inc().Negate(); + } + + public BigInteger Pow(int exp) + { + if (exp <= 0) + { + if (exp < 0) + throw new ArithmeticException("Negative exponent"); + + return One; + } + + if (sign == 0) + { + return this; + } + + if (QuickPow2Check()) + { + long powOf2 = (long)exp * (BitLength - 1); + if (powOf2 > Int32.MaxValue) + { + throw new ArithmeticException("Result too large"); + } + return One.ShiftLeft((int)powOf2); + } + + BigInteger y = One; + BigInteger z = this; + + for (;;) + { + if ((exp & 0x1) == 1) + { + y = y.Multiply(z); + } + exp >>= 1; + if (exp == 0) break; + z = z.Multiply(z); + } + + return y; + } + + public static BigInteger ProbablePrime( + int bitLength, + Random random) + { + return new BigInteger(bitLength, 100, random); + } + + private int Remainder( + int m) + { + Debug.Assert(m > 0); + + long acc = 0; + for (int pos = 0; pos < magnitude.Length; ++pos) + { + long posVal = (uint) magnitude[pos]; + acc = (acc << 32 | posVal) % m; + } + + return (int) acc; + } + + /** + * return x = x % y - done in place (y value preserved) + */ + private static int[] Remainder( + 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); + + if (xyCmp > 0) + { + int yBitLength = CalcBitLength(1, yStart, y); + int xBitLength = CalcBitLength(1, xStart, x); + int shift = xBitLength - yBitLength; + + int[] c; + int cStart = 0; + int cBitLength = yBitLength; + if (shift > 0) + { + c = ShiftLeft(y, shift); + cBitLength += shift; + Debug.Assert(c[0] != 0); + } + else + { + int len = y.Length - yStart; + c = new int[len]; + Array.Copy(y, yStart, c, 0, len); + } + + for (;;) + { + if (cBitLength < xBitLength + || CompareNoLeadingZeroes(xStart, x, cStart, c) >= 0) + { + Subtract(xStart, x, cStart, c); + + while (x[xStart] == 0) + { + if (++xStart == x.Length) + return x; + } + + //xBitLength = CalcBitLength(xStart, x); + xBitLength = 32 * (x.Length - xStart - 1) + BitLen(x[xStart]); + + if (xBitLength <= yBitLength) + { + if (xBitLength < yBitLength) + return x; + + 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; + } + else + { + ShiftRightInPlace(cStart, c, shift); + cBitLength -= shift; + } + + //cStart = c.Length - ((cBitLength + 31) / 32); + while (c[cStart] == 0) + { + ++cStart; + } + } + } + + if (xyCmp == 0) + { + Array.Clear(x, xStart, x.Length - xStart); + } + + return x; + } + + public BigInteger Remainder( + BigInteger n) + { + if (n.sign == 0) + throw new ArithmeticException("Division by zero error"); + + if (this.sign == 0) + return Zero; + + // For small values, use fast remainder method + if (n.magnitude.Length == 1) + { + int val = n.magnitude[0]; + + if (val > 0) + { + if (val == 1) + return Zero; + + // TODO Make this func work on uint, and handle val == 1? + int rem = Remainder(val); + + return rem == 0 + ? Zero + : new BigInteger(sign, new int[]{ rem }, false); + } + } + + if (CompareNoLeadingZeroes(0, magnitude, 0, n.magnitude) < 0) + return this; + + int[] result; + if (n.QuickPow2Check()) // n is power of two + { + // TODO Move before small values branch above? + result = LastNBits(n.Abs().BitLength - 1); + } + else + { + result = (int[]) this.magnitude.Clone(); + result = Remainder(result, n.magnitude); + } + + return new BigInteger(sign, result, true); + } + + private int[] LastNBits( + int n) + { + if (n < 1) + return ZeroMagnitude; + + int numWords = (n + BitsPerInt - 1) / BitsPerInt; + numWords = System.Math.Min(numWords, this.magnitude.Length); + int[] result = new int[numWords]; + + Array.Copy(this.magnitude, this.magnitude.Length - numWords, result, 0, numWords); + + int excessBits = (numWords << 5) - n; + if (excessBits > 0) + { + result[0] &= (int)(UInt32.MaxValue >> excessBits); + } + + return result; + } + + private BigInteger DivideWords(int w) + { + Debug.Assert(w >= 0); + int n = magnitude.Length; + if (w >= n) + return Zero; + int[] mag = new int[n - w]; + Array.Copy(magnitude, 0, mag, 0, n - w); + return new BigInteger(sign, mag, false); + } + + private BigInteger RemainderWords(int w) + { + Debug.Assert(w >= 0); + int n = magnitude.Length; + if (w >= n) + return this; + int[] mag = new int[w]; + Array.Copy(magnitude, n - w, mag, 0, w); + return new BigInteger(sign, mag, false); + } + + /** + * do a left shift - this returns a new array. + */ + private static int[] ShiftLeft( + int[] mag, + int n) + { + int nInts = (int)((uint)n >> 5); + int nBits = n & 0x1f; + int magLen = mag.Length; + int[] newMag; + + if (nBits == 0) + { + newMag = new int[magLen + nInts]; + mag.CopyTo(newMag, 0); + } + else + { + int i = 0; + int nBits2 = 32 - nBits; + int highBits = (int)((uint)mag[0] >> nBits2); + + if (highBits != 0) + { + newMag = new int[magLen + nInts + 1]; + newMag[i++] = highBits; + } + else + { + newMag = new int[magLen + nInts]; + } + + int m = mag[0]; + for (int j = 0; j < magLen - 1; j++) + { + int next = mag[j + 1]; + + newMag[i++] = (m << nBits) | (int)((uint)next >> nBits2); + m = next; + } + + newMag[i] = mag[magLen - 1] << nBits; + } + + return newMag; + } + + private static int ShiftLeftOneInPlace(int[] x, int carry) + { + Debug.Assert(carry == 0 || carry == 1); + int pos = x.Length; + while (--pos >= 0) + { + uint val = (uint)x[pos]; + x[pos] = (int)(val << 1) | carry; + carry = (int)(val >> 31); + } + return carry; + } + + public BigInteger ShiftLeft( + int n) + { + if (sign == 0 || magnitude.Length == 0) + return Zero; + + if (n == 0) + return this; + + if (n < 0) + return ShiftRight(-n); + + BigInteger result = new BigInteger(sign, ShiftLeft(magnitude, n), true); + + if (this.nBits != -1) + { + result.nBits = sign > 0 + ? this.nBits + : this.nBits + n; + } + + if (this.nBitLength != -1) + { + result.nBitLength = this.nBitLength + n; + } + + return result; + } + + /** + * do a right shift - this does it in place. + */ + private static void ShiftRightInPlace( + int start, + int[] mag, + int n) + { + int nInts = (int)((uint)n >> 5) + start; + int nBits = n & 0x1f; + int magEnd = mag.Length - 1; + + if (nInts != start) + { + int delta = (nInts - start); + + for (int i = magEnd; i >= nInts; i--) + { + mag[i] = mag[i - delta]; + } + for (int i = nInts - 1; i >= start; i--) + { + mag[i] = 0; + } + } + + if (nBits != 0) + { + int nBits2 = 32 - nBits; + int m = mag[magEnd]; + + for (int i = magEnd; i > nInts; --i) + { + int next = mag[i - 1]; + + mag[i] = (int)((uint)m >> nBits) | (next << nBits2); + m = next; + } + + mag[nInts] = (int)((uint)mag[nInts] >> nBits); + } + } + + /** + * do a right shift by one - this does it in place. + */ + private static void ShiftRightOneInPlace( + int start, + int[] mag) + { + int i = mag.Length; + int m = mag[i - 1]; + + while (--i > start) + { + int next = mag[i - 1]; + mag[i] = ((int)((uint)m >> 1)) | (next << 31); + m = next; + } + + mag[start] = (int)((uint)mag[start] >> 1); + } + + public BigInteger ShiftRight( + int n) + { + if (n == 0) + return this; + + if (n < 0) + return ShiftLeft(-n); + + if (n >= BitLength) + return (this.sign < 0 ? One.Negate() : Zero); + +// int[] res = (int[]) this.magnitude.Clone(); +// +// ShiftRightInPlace(0, res, n); +// +// return new BigInteger(this.sign, res, true); + + int resultLength = (BitLength - n + 31) >> 5; + int[] res = new int[resultLength]; + + int numInts = n >> 5; + int numBits = n & 31; + + if (numBits == 0) + { + Array.Copy(this.magnitude, 0, res, 0, res.Length); + } + else + { + int numBits2 = 32 - numBits; + + int magPos = this.magnitude.Length - 1 - numInts; + for (int i = resultLength - 1; i >= 0; --i) + { + res[i] = (int)((uint) this.magnitude[magPos--] >> numBits); + + if (magPos >= 0) + { + res[i] |= this.magnitude[magPos] << numBits2; + } + } + } + + Debug.Assert(res[0] != 0); + + return new BigInteger(this.sign, res, false); + } + + public int SignValue + { + get { return sign; } + } + + /** + * returns x = x - y - we assume x is >= y + */ + private static int[] Subtract( + int xStart, + int[] x, + int yStart, + int[] y) + { + Debug.Assert(yStart < y.Length); + Debug.Assert(x.Length - xStart >= y.Length - yStart); + + int iT = x.Length; + int iV = y.Length; + long m; + int borrow = 0; + + do + { + m = (x[--iT] & IMASK) - (y[--iV] & IMASK) + borrow; + x[iT] = (int) m; + +// borrow = (m < 0) ? -1 : 0; + borrow = (int)(m >> 63); + } + while (iV > yStart); + + if (borrow != 0) + { + while (--x[--iT] == -1) + { + } + } + + return x; + } + + public BigInteger Subtract( + BigInteger n) + { + if (n.sign == 0) + return this; + + if (this.sign == 0) + return n.Negate(); + + if (this.sign != n.sign) + return Add(n.Negate()); + + int compare = CompareNoLeadingZeroes(0, magnitude, 0, n.magnitude); + if (compare == 0) + return Zero; + + BigInteger bigun, lilun; + if (compare < 0) + { + bigun = n; + lilun = this; + } + else + { + bigun = this; + lilun = n; + } + + return new BigInteger(this.sign * compare, doSubBigLil(bigun.magnitude, lilun.magnitude), true); + } + + private static int[] doSubBigLil( + int[] bigMag, + int[] lilMag) + { + int[] res = (int[]) bigMag.Clone(); + + return Subtract(0, res, 0, lilMag); + } + + public byte[] ToByteArray() + { + return ToByteArray(false); + } + + public byte[] ToByteArrayUnsigned() + { + return ToByteArray(true); + } + + private byte[] ToByteArray( + bool unsigned) + { + if (sign == 0) + return unsigned ? ZeroEncoding : new byte[1]; + + int nBits = (unsigned && sign > 0) + ? BitLength + : BitLength + 1; + + int nBytes = GetByteLength(nBits); + byte[] bytes = new byte[nBytes]; + + int magIndex = magnitude.Length; + int bytesIndex = bytes.Length; + + if (sign > 0) + { + while (magIndex > 1) + { + uint mag = (uint) magnitude[--magIndex]; + bytes[--bytesIndex] = (byte) mag; + bytes[--bytesIndex] = (byte)(mag >> 8); + bytes[--bytesIndex] = (byte)(mag >> 16); + bytes[--bytesIndex] = (byte)(mag >> 24); + } + + uint lastMag = (uint) magnitude[0]; + while (lastMag > byte.MaxValue) + { + bytes[--bytesIndex] = (byte) lastMag; + lastMag >>= 8; + } + + bytes[--bytesIndex] = (byte) lastMag; + } + else // sign < 0 + { + bool carry = true; + + while (magIndex > 1) + { + uint mag = ~((uint) magnitude[--magIndex]); + + if (carry) + { + carry = (++mag == uint.MinValue); + } + + bytes[--bytesIndex] = (byte) mag; + bytes[--bytesIndex] = (byte)(mag >> 8); + bytes[--bytesIndex] = (byte)(mag >> 16); + bytes[--bytesIndex] = (byte)(mag >> 24); + } + + uint lastMag = (uint) magnitude[0]; + + if (carry) + { + // Never wraps because magnitude[0] != 0 + --lastMag; + } + + while (lastMag > byte.MaxValue) + { + bytes[--bytesIndex] = (byte) ~lastMag; + lastMag >>= 8; + } + + bytes[--bytesIndex] = (byte) ~lastMag; + + if (bytesIndex > 0) + { + bytes[--bytesIndex] = byte.MaxValue; + } + } + + return bytes; + } + + public override string ToString() + { + return ToString(10); + } + + public string ToString(int radix) + { + // TODO Make this method work for other radices (ideally 2 <= radix <= 36 as in Java) + + switch (radix) + { + case 2: + case 8: + case 10: + case 16: + break; + default: + throw new FormatException("Only bases 2, 8, 10, 16 are allowed"); + } + + // NB: Can only happen to internally managed instances + if (magnitude == null) + return "null"; + + if (sign == 0) + return "0"; + + + // NOTE: This *should* be unnecessary, since the magnitude *should* never have leading zero digits + int firstNonZero = 0; + while (firstNonZero < magnitude.Length) + { + if (magnitude[firstNonZero] != 0) + { + break; + } + ++firstNonZero; + } + + if (firstNonZero == magnitude.Length) + { + return "0"; + } + + + StringBuilder sb = new StringBuilder(); + if (sign == -1) + { + sb.Append('-'); + } + + switch (radix) + { + case 2: + { + int pos = firstNonZero; + sb.Append(Convert.ToString(magnitude[pos], 2)); + while (++pos < magnitude.Length) + { + AppendZeroExtendedString(sb, Convert.ToString(magnitude[pos], 2), 32); + } + break; + } + case 8: + { + int mask = (1 << 30) - 1; + BigInteger u = this.Abs(); + int bits = u.BitLength; + IList S = Platform.CreateArrayList(); + while (bits > 30) + { + S.Add(Convert.ToString(u.IntValue & mask, 8)); + u = u.ShiftRight(30); + bits -= 30; + } + sb.Append(Convert.ToString(u.IntValue, 8)); + for (int i = S.Count - 1; i >= 0; --i) + { + AppendZeroExtendedString(sb, (string)S[i], 10); + } + break; + } + case 16: + { + int pos = firstNonZero; + sb.Append(Convert.ToString(magnitude[pos], 16)); + while (++pos < magnitude.Length) + { + AppendZeroExtendedString(sb, Convert.ToString(magnitude[pos], 16), 8); + } + break; + } + // TODO This could work for other radices if there is an alternative to Convert.ToString method + //default: + case 10: + { + BigInteger q = this.Abs(); + if (q.BitLength < 64) + { + sb.Append(Convert.ToString(q.LongValue, radix)); + break; + } + + // Based on algorithm 1a from chapter 4.4 in Seminumerical Algorithms (Knuth) + + // Work out the largest power of 'rdx' that is a positive 64-bit integer + // TODO possibly cache power/exponent against radix? + long limit = Int64.MaxValue / radix; + long power = radix; + int exponent = 1; + while (power <= limit) + { + power *= radix; + ++exponent; + } + + BigInteger bigPower = BigInteger.ValueOf(power); + + IList S = Platform.CreateArrayList(); + while (q.CompareTo(bigPower) >= 0) + { + BigInteger[] qr = q.DivideAndRemainder(bigPower); + S.Add(Convert.ToString(qr[1].LongValue, radix)); + q = qr[0]; + } + + sb.Append(Convert.ToString(q.LongValue, radix)); + for (int i = S.Count - 1; i >= 0; --i) + { + AppendZeroExtendedString(sb, (string)S[i], exponent); + } + break; + } + } + + return sb.ToString(); + } + + private static void AppendZeroExtendedString(StringBuilder sb, string s, int minLength) + { + for (int len = s.Length; len < minLength; ++len) + { + sb.Append('0'); + } + sb.Append(s); + } + + private static BigInteger CreateUValueOf( + ulong value) + { + int msw = (int)(value >> 32); + int lsw = (int)value; + + if (msw != 0) + return new BigInteger(1, new int[] { msw, lsw }, false); + + if (lsw != 0) + { + BigInteger n = new BigInteger(1, new int[] { lsw }, false); + // Check for a power of two + if ((lsw & -lsw) == lsw) + { + n.nBits = 1; + } + return n; + } + + return Zero; + } + + private static BigInteger CreateValueOf( + long value) + { + if (value < 0) + { + if (value == long.MinValue) + return CreateValueOf(~value).Not(); + + return CreateValueOf(-value).Negate(); + } + + return CreateUValueOf((ulong)value); + } + + public static BigInteger ValueOf( + long value) + { + if (value >= 0 && value < SMALL_CONSTANTS.Length) + { + return SMALL_CONSTANTS[value]; + } + + return CreateValueOf(value); + } + + public int GetLowestSetBit() + { + if (this.sign == 0) + return -1; + + return GetLowestSetBitMaskFirst(-1); + } + + private int GetLowestSetBitMaskFirst(int firstWordMask) + { + int w = magnitude.Length, offset = 0; + + uint word = (uint)(magnitude[--w] & firstWordMask); + Debug.Assert(magnitude[0] != 0); + + while (word == 0) + { + word = (uint)magnitude[--w]; + offset += 32; + } + + while ((word & 0xFF) == 0) + { + word >>= 8; + offset += 8; + } + + while ((word & 1) == 0) + { + word >>= 1; + ++offset; + } + + return offset; + } + + public bool TestBit( + int n) + { + if (n < 0) + throw new ArithmeticException("Bit position must not be negative"); + + if (sign < 0) + return !Not().TestBit(n); + + int wordNum = n / 32; + if (wordNum >= magnitude.Length) + return false; + + int word = magnitude[magnitude.Length - 1 - wordNum]; + return ((word >> (n % 32)) & 1) > 0; + } + + public BigInteger Or( + BigInteger value) + { + if (this.sign == 0) + return value; + + if (value.sign == 0) + return this; + + 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 Xor( + BigInteger value) + { + if (this.sign == 0) + return value; + + if (value.sign == 0) + return this; + + int[] aMag = this.sign > 0 + ? this.magnitude + : Add(One).magnitude; + + int[] bMag = value.sign > 0 + ? value.magnitude + : value.Add(One).magnitude; + + // TODO Can just replace with sign != value.sign? + bool resultNeg = (sign < 0 && value.sign >= 0) || (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 SetBit( + int n) + { + if (n < 0) + throw new ArithmeticException("Bit address less than zero"); + + if (TestBit(n)) + return this; + + // TODO Handle negative values and zero + if (sign > 0 && n < (BitLength - 1)) + return FlipExistingBit(n); + + return Or(One.ShiftLeft(n)); + } + + public BigInteger ClearBit( + int n) + { + if (n < 0) + throw new ArithmeticException("Bit address less than zero"); + + if (!TestBit(n)) + return this; + + // TODO Handle negative values + if (sign > 0 && n < (BitLength - 1)) + return FlipExistingBit(n); + + return AndNot(One.ShiftLeft(n)); + } + + public BigInteger FlipBit( + int n) + { + if (n < 0) + throw new ArithmeticException("Bit address less than zero"); + + // TODO Handle negative values and zero + if (sign > 0 && n < (BitLength - 1)) + return FlipExistingBit(n); + + return Xor(One.ShiftLeft(n)); + } + + private BigInteger FlipExistingBit( + int n) + { + Debug.Assert(sign > 0); + Debug.Assert(n >= 0); + Debug.Assert(n < BitLength - 1); + + int[] mag = (int[]) this.magnitude.Clone(); + mag[mag.Length - 1 - (n >> 5)] ^= (1 << (n & 31)); // Flip bit + //mag[mag.Length - 1 - (n / 32)] ^= (1 << (n % 32)); + return new BigInteger(this.sign, mag, false); + } + } +} diff --git a/crypto/src/math/ec/ECAlgorithms.cs b/crypto/src/math/ec/ECAlgorithms.cs new file mode 100644 index 000000000..be4fd1b14 --- /dev/null +++ b/crypto/src/math/ec/ECAlgorithms.cs @@ -0,0 +1,93 @@ +using System; + +using Org.BouncyCastle.Math; + +namespace Org.BouncyCastle.Math.EC +{ + public class ECAlgorithms + { + public static ECPoint SumOfTwoMultiplies(ECPoint P, BigInteger a, + ECPoint Q, BigInteger b) + { + ECCurve c = P.Curve; + if (!c.Equals(Q.Curve)) + throw new ArgumentException("P and Q must be on same curve"); + + // Point multiplication for Koblitz curves (using WTNAF) beats Shamir's trick + if (c is F2mCurve) + { + F2mCurve f2mCurve = (F2mCurve) c; + if (f2mCurve.IsKoblitz) + { + return P.Multiply(a).Add(Q.Multiply(b)); + } + } + + return ImplShamirsTrick(P, a, Q, b); + } + + /* + * "Shamir's Trick", originally due to E. G. Straus + * (Addition chains of vectors. American Mathematical Monthly, + * 71(7):806-808, Aug./Sept. 1964) + * + * Input: The points P, Q, scalar k = (km?, ... , k1, k0) + * and scalar l = (lm?, ... , l1, l0). + * Output: R = k * P + l * Q. + * 1: Z <- P + Q + * 2: R <- O + * 3: for i from m-1 down to 0 do + * 4: R <- R + R {point doubling} + * 5: if (ki = 1) and (li = 0) then R <- R + P end if + * 6: if (ki = 0) and (li = 1) then R <- R + Q end if + * 7: if (ki = 1) and (li = 1) then R <- R + Z end if + * 8: end for + * 9: return R + */ + public static ECPoint ShamirsTrick( + ECPoint P, + BigInteger k, + ECPoint Q, + BigInteger l) + { + if (!P.Curve.Equals(Q.Curve)) + throw new ArgumentException("P and Q must be on same curve"); + + return ImplShamirsTrick(P, k, Q, l); + } + + private static ECPoint ImplShamirsTrick(ECPoint P, BigInteger k, + ECPoint Q, BigInteger l) + { + int m = System.Math.Max(k.BitLength, l.BitLength); + ECPoint Z = P.Add(Q); + ECPoint R = P.Curve.Infinity; + + for (int i = m - 1; i >= 0; --i) + { + R = R.Twice(); + + if (k.TestBit(i)) + { + if (l.TestBit(i)) + { + R = R.Add(Z); + } + else + { + R = R.Add(P); + } + } + else + { + if (l.TestBit(i)) + { + R = R.Add(Q); + } + } + } + + return R; + } + } +} diff --git a/crypto/src/math/ec/ECCurve.cs b/crypto/src/math/ec/ECCurve.cs new file mode 100644 index 000000000..396d42f28 --- /dev/null +++ b/crypto/src/math/ec/ECCurve.cs @@ -0,0 +1,651 @@ +using System; +using System.Collections; + +using Org.BouncyCastle.Math.EC.Abc; + +namespace Org.BouncyCastle.Math.EC +{ + /// <remarks>Base class for an elliptic curve.</remarks> + public abstract class ECCurve + { + internal ECFieldElement a, b; + + public abstract int FieldSize { get; } + public abstract ECFieldElement FromBigInteger(BigInteger x); + public abstract ECPoint CreatePoint(BigInteger x, BigInteger y, bool withCompression); + public abstract ECPoint Infinity { get; } + + public ECFieldElement A + { + get { return a; } + } + + public ECFieldElement B + { + get { return b; } + } + + public override bool Equals( + object obj) + { + if (obj == this) + return true; + + ECCurve other = obj as ECCurve; + + if (other == null) + return false; + + return Equals(other); + } + + protected bool Equals( + ECCurve other) + { + return a.Equals(other.a) && b.Equals(other.b); + } + + public override int GetHashCode() + { + return a.GetHashCode() ^ b.GetHashCode(); + } + + protected abstract ECPoint DecompressPoint(int yTilde, BigInteger X1); + + /** + * Decode a point on this curve from its ASN.1 encoding. The different + * encodings are taken account of, including point compression for + * <code>F<sub>p</sub></code> (X9.62 s 4.2.1 pg 17). + * @return The decoded point. + */ + public virtual ECPoint DecodePoint(byte[] encoded) + { + ECPoint p = null; + int expectedLength = (FieldSize + 7) / 8; + + switch (encoded[0]) + { + case 0x00: // infinity + { + if (encoded.Length != 1) + throw new ArgumentException("Incorrect length for infinity encoding", "encoded"); + + p = Infinity; + break; + } + + case 0x02: // compressed + case 0x03: // compressed + { + if (encoded.Length != (expectedLength + 1)) + throw new ArgumentException("Incorrect length for compressed encoding", "encoded"); + + int yTilde = encoded[0] & 1; + BigInteger X1 = new BigInteger(1, encoded, 1, expectedLength); + + p = DecompressPoint(yTilde, X1); + break; + } + + case 0x04: // uncompressed + case 0x06: // hybrid + case 0x07: // hybrid + { + if (encoded.Length != (2 * expectedLength + 1)) + throw new ArgumentException("Incorrect length for uncompressed/hybrid encoding", "encoded"); + + BigInteger X1 = new BigInteger(1, encoded, 1, expectedLength); + BigInteger Y1 = new BigInteger(1, encoded, 1 + expectedLength, expectedLength); + + p = CreatePoint(X1, Y1, false); + break; + } + + default: + throw new FormatException("Invalid point encoding " + encoded[0]); + } + + return p; + } + } + + /** + * Elliptic curve over Fp + */ + public class FpCurve : ECCurve + { + private readonly BigInteger q; + private readonly FpPoint infinity; + + public FpCurve(BigInteger q, BigInteger a, BigInteger b) + { + this.q = q; + this.a = FromBigInteger(a); + this.b = FromBigInteger(b); + this.infinity = new FpPoint(this, null, null); + } + + public BigInteger Q + { + get { return q; } + } + + public override ECPoint Infinity + { + get { return infinity; } + } + + public override int FieldSize + { + get { return q.BitLength; } + } + + public override ECFieldElement FromBigInteger(BigInteger x) + { + return new FpFieldElement(this.q, x); + } + + public override ECPoint CreatePoint( + BigInteger X1, + BigInteger Y1, + bool withCompression) + { + // TODO Validation of X1, Y1? + return new FpPoint( + this, + FromBigInteger(X1), + FromBigInteger(Y1), + withCompression); + } + + protected override ECPoint DecompressPoint( + int yTilde, + BigInteger X1) + { + ECFieldElement x = FromBigInteger(X1); + ECFieldElement alpha = x.Multiply(x.Square().Add(a)).Add(b); + ECFieldElement beta = alpha.Sqrt(); + + // + // if we can't find a sqrt we haven't got a point on the + // curve - run! + // + if (beta == null) + throw new ArithmeticException("Invalid point compression"); + + BigInteger betaValue = beta.ToBigInteger(); + int bit0 = betaValue.TestBit(0) ? 1 : 0; + + if (bit0 != yTilde) + { + // Use the other root + beta = FromBigInteger(q.Subtract(betaValue)); + } + + return new FpPoint(this, x, beta, true); + } + + public override bool Equals( + object obj) + { + if (obj == this) + return true; + + FpCurve other = obj as FpCurve; + + if (other == null) + return false; + + return Equals(other); + } + + protected bool Equals( + FpCurve other) + { + return base.Equals(other) && q.Equals(other.q); + } + + public override int GetHashCode() + { + return base.GetHashCode() ^ q.GetHashCode(); + } + } + + /** + * Elliptic curves over F2m. The Weierstrass equation is given by + * <code>y<sup>2</sup> + xy = x<sup>3</sup> + ax<sup>2</sup> + b</code>. + */ + public class F2mCurve : ECCurve + { + /** + * The exponent <code>m</code> of <code>F<sub>2<sup>m</sup></sub></code>. + */ + private readonly int m; + + /** + * TPB: The integer <code>k</code> where <code>x<sup>m</sup> + + * x<sup>k</sup> + 1</code> represents the reduction polynomial + * <code>f(z)</code>.<br/> + * PPB: The integer <code>k1</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>.<br/> + */ + private readonly int k1; + + /** + * TPB: Always set to <code>0</code><br/> + * PPB: The integer <code>k2</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>.<br/> + */ + private readonly int k2; + + /** + * TPB: Always set to <code>0</code><br/> + * PPB: The integer <code>k3</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>.<br/> + */ + private readonly int k3; + + /** + * The order of the base point of the curve. + */ + private readonly BigInteger n; + + /** + * The cofactor of the curve. + */ + private readonly BigInteger h; + + /** + * The point at infinity on this curve. + */ + private readonly F2mPoint infinity; + + /** + * The parameter <code>μ</code> of the elliptic curve if this is + * a Koblitz curve. + */ + private sbyte mu = 0; + + /** + * The auxiliary values <code>s<sub>0</sub></code> and + * <code>s<sub>1</sub></code> used for partial modular reduction for + * Koblitz curves. + */ + private BigInteger[] si = null; + + /** + * Constructor for Trinomial Polynomial Basis (TPB). + * @param m The exponent <code>m</code> of + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param k The integer <code>k</code> where <code>x<sup>m</sup> + + * x<sup>k</sup> + 1</code> represents the reduction + * polynomial <code>f(z)</code>. + * @param a The coefficient <code>a</code> in the Weierstrass equation + * for non-supersingular elliptic curves over + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param b The coefficient <code>b</code> in the Weierstrass equation + * for non-supersingular elliptic curves over + * <code>F<sub>2<sup>m</sup></sub></code>. + */ + public F2mCurve( + int m, + int k, + BigInteger a, + BigInteger b) + : this(m, k, 0, 0, a, b, null, null) + { + } + + /** + * Constructor for Trinomial Polynomial Basis (TPB). + * @param m The exponent <code>m</code> of + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param k The integer <code>k</code> where <code>x<sup>m</sup> + + * x<sup>k</sup> + 1</code> represents the reduction + * polynomial <code>f(z)</code>. + * @param a The coefficient <code>a</code> in the Weierstrass equation + * for non-supersingular elliptic curves over + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param b The coefficient <code>b</code> in the Weierstrass equation + * for non-supersingular elliptic curves over + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param n The order of the main subgroup of the elliptic curve. + * @param h The cofactor of the elliptic curve, i.e. + * <code>#E<sub>a</sub>(F<sub>2<sup>m</sup></sub>) = h * n</code>. + */ + public F2mCurve( + int m, + int k, + BigInteger a, + BigInteger b, + BigInteger n, + BigInteger h) + : this(m, k, 0, 0, a, b, n, h) + { + } + + /** + * Constructor for Pentanomial Polynomial Basis (PPB). + * @param m The exponent <code>m</code> of + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param k1 The integer <code>k1</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>. + * @param k2 The integer <code>k2</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>. + * @param k3 The integer <code>k3</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>. + * @param a The coefficient <code>a</code> in the Weierstrass equation + * for non-supersingular elliptic curves over + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param b The coefficient <code>b</code> in the Weierstrass equation + * for non-supersingular elliptic curves over + * <code>F<sub>2<sup>m</sup></sub></code>. + */ + public F2mCurve( + int m, + int k1, + int k2, + int k3, + BigInteger a, + BigInteger b) + : this(m, k1, k2, k3, a, b, null, null) + { + } + + /** + * Constructor for Pentanomial Polynomial Basis (PPB). + * @param m The exponent <code>m</code> of + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param k1 The integer <code>k1</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>. + * @param k2 The integer <code>k2</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>. + * @param k3 The integer <code>k3</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>. + * @param a The coefficient <code>a</code> in the Weierstrass equation + * for non-supersingular elliptic curves over + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param b The coefficient <code>b</code> in the Weierstrass equation + * for non-supersingular elliptic curves over + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param n The order of the main subgroup of the elliptic curve. + * @param h The cofactor of the elliptic curve, i.e. + * <code>#E<sub>a</sub>(F<sub>2<sup>m</sup></sub>) = h * n</code>. + */ + public F2mCurve( + int m, + int k1, + int k2, + int k3, + BigInteger a, + BigInteger b, + BigInteger n, + BigInteger h) + { + this.m = m; + this.k1 = k1; + this.k2 = k2; + this.k3 = k3; + this.n = n; + this.h = h; + this.infinity = new F2mPoint(this, null, null); + + if (k1 == 0) + throw new ArgumentException("k1 must be > 0"); + + if (k2 == 0) + { + if (k3 != 0) + throw new ArgumentException("k3 must be 0 if k2 == 0"); + } + else + { + if (k2 <= k1) + throw new ArgumentException("k2 must be > k1"); + + if (k3 <= k2) + throw new ArgumentException("k3 must be > k2"); + } + + this.a = FromBigInteger(a); + this.b = FromBigInteger(b); + } + + public override ECPoint Infinity + { + get { return infinity; } + } + + public override int FieldSize + { + get { return m; } + } + + public override ECFieldElement FromBigInteger(BigInteger x) + { + return new F2mFieldElement(this.m, this.k1, this.k2, this.k3, x); + } + + /** + * Returns true if this is a Koblitz curve (ABC curve). + * @return true if this is a Koblitz curve (ABC curve), false otherwise + */ + public bool IsKoblitz + { + get + { + return n != null && h != null + && (a.ToBigInteger().Equals(BigInteger.Zero) + || a.ToBigInteger().Equals(BigInteger.One)) + && b.ToBigInteger().Equals(BigInteger.One); + } + } + + /** + * Returns the parameter <code>μ</code> of the elliptic curve. + * @return <code>μ</code> of the elliptic curve. + * @throws ArgumentException if the given ECCurve is not a + * Koblitz curve. + */ + internal sbyte GetMu() + { + if (mu == 0) + { + lock (this) + { + if (mu == 0) + { + mu = Tnaf.GetMu(this); + } + } + } + + return mu; + } + + /** + * @return the auxiliary values <code>s<sub>0</sub></code> and + * <code>s<sub>1</sub></code> used for partial modular reduction for + * Koblitz curves. + */ + internal BigInteger[] GetSi() + { + if (si == null) + { + lock (this) + { + if (si == null) + { + si = Tnaf.GetSi(this); + } + } + } + return si; + } + + public override ECPoint CreatePoint( + BigInteger X1, + BigInteger Y1, + bool withCompression) + { + // TODO Validation of X1, Y1? + return new F2mPoint( + this, + FromBigInteger(X1), + FromBigInteger(Y1), + withCompression); + } + + protected override ECPoint DecompressPoint( + int yTilde, + BigInteger X1) + { + ECFieldElement xp = FromBigInteger(X1); + ECFieldElement yp = null; + if (xp.ToBigInteger().SignValue == 0) + { + yp = (F2mFieldElement)b; + for (int i = 0; i < m - 1; i++) + { + yp = yp.Square(); + } + } + else + { + ECFieldElement beta = xp.Add(a).Add(b.Multiply(xp.Square().Invert())); + ECFieldElement z = solveQuadradicEquation(beta); + + if (z == null) + throw new ArithmeticException("Invalid point compression"); + + int zBit = z.ToBigInteger().TestBit(0) ? 1 : 0; + if (zBit != yTilde) + { + z = z.Add(FromBigInteger(BigInteger.One)); + } + + yp = xp.Multiply(z); + } + + return new F2mPoint(this, xp, yp, true); + } + + /** + * Solves a quadratic equation <code>z<sup>2</sup> + z = beta</code>(X9.62 + * D.1.6) The other solution is <code>z + 1</code>. + * + * @param beta + * The value to solve the qradratic equation for. + * @return the solution for <code>z<sup>2</sup> + z = beta</code> or + * <code>null</code> if no solution exists. + */ + private ECFieldElement solveQuadradicEquation(ECFieldElement beta) + { + if (beta.ToBigInteger().SignValue == 0) + { + return FromBigInteger(BigInteger.Zero); + } + + ECFieldElement z = null; + ECFieldElement gamma = FromBigInteger(BigInteger.Zero); + + while (gamma.ToBigInteger().SignValue == 0) + { + ECFieldElement t = FromBigInteger(new BigInteger(m, new Random())); + z = FromBigInteger(BigInteger.Zero); + + ECFieldElement w = beta; + for (int i = 1; i <= m - 1; i++) + { + ECFieldElement w2 = w.Square(); + z = z.Square().Add(w2.Multiply(t)); + w = w2.Add(beta); + } + if (w.ToBigInteger().SignValue != 0) + { + return null; + } + gamma = z.Square().Add(z); + } + return z; + } + + public override bool Equals( + object obj) + { + if (obj == this) + return true; + + F2mCurve other = obj as F2mCurve; + + if (other == null) + return false; + + return Equals(other); + } + + protected bool Equals( + F2mCurve other) + { + return m == other.m + && k1 == other.k1 + && k2 == other.k2 + && k3 == other.k3 + && base.Equals(other); + } + + public override int GetHashCode() + { + return base.GetHashCode() ^ m ^ k1 ^ k2 ^ k3; + } + + public int M + { + get { return m; } + } + + /** + * Return true if curve uses a Trinomial basis. + * + * @return true if curve Trinomial, false otherwise. + */ + public bool IsTrinomial() + { + return k2 == 0 && k3 == 0; + } + + public int K1 + { + get { return k1; } + } + + public int K2 + { + get { return k2; } + } + + public int K3 + { + get { return k3; } + } + + public BigInteger N + { + get { return n; } + } + + public BigInteger H + { + get { return h; } + } + } +} diff --git a/crypto/src/math/ec/ECFieldElement.cs b/crypto/src/math/ec/ECFieldElement.cs new file mode 100644 index 000000000..5235c6c0e --- /dev/null +++ b/crypto/src/math/ec/ECFieldElement.cs @@ -0,0 +1,1253 @@ +using System; +using System.Diagnostics; + +using Org.BouncyCastle.Utilities; + +namespace Org.BouncyCastle.Math.EC +{ + public abstract class ECFieldElement + { + public abstract BigInteger ToBigInteger(); + public abstract string FieldName { get; } + public abstract int FieldSize { get; } + public abstract ECFieldElement Add(ECFieldElement b); + public abstract ECFieldElement Subtract(ECFieldElement b); + public abstract ECFieldElement Multiply(ECFieldElement b); + public abstract ECFieldElement Divide(ECFieldElement b); + public abstract ECFieldElement Negate(); + public abstract ECFieldElement Square(); + public abstract ECFieldElement Invert(); + public abstract ECFieldElement Sqrt(); + + public override bool Equals( + object obj) + { + if (obj == this) + return true; + + ECFieldElement other = obj as ECFieldElement; + + if (other == null) + return false; + + return Equals(other); + } + + protected bool Equals( + ECFieldElement other) + { + return ToBigInteger().Equals(other.ToBigInteger()); + } + + public override int GetHashCode() + { + return ToBigInteger().GetHashCode(); + } + + public override string ToString() + { + return this.ToBigInteger().ToString(2); + } + } + + public class FpFieldElement + : ECFieldElement + { + private readonly BigInteger q, x; + + public FpFieldElement( + BigInteger q, + BigInteger x) + { + if (x.CompareTo(q) >= 0) + throw new ArgumentException("x value too large in field element"); + + this.q = q; + this.x = x; + } + + public override BigInteger ToBigInteger() + { + return x; + } + + /** + * return the field name for this field. + * + * @return the string "Fp". + */ + public override string FieldName + { + get { return "Fp"; } + } + + public override int FieldSize + { + get { return q.BitLength; } + } + + public BigInteger Q + { + get { return q; } + } + + public override ECFieldElement Add( + ECFieldElement b) + { + return new FpFieldElement(q, x.Add(b.ToBigInteger()).Mod(q)); + } + + public override ECFieldElement Subtract( + ECFieldElement b) + { + return new FpFieldElement(q, x.Subtract(b.ToBigInteger()).Mod(q)); + } + + public override ECFieldElement Multiply( + ECFieldElement b) + { + return new FpFieldElement(q, x.Multiply(b.ToBigInteger()).Mod(q)); + } + + public override ECFieldElement Divide( + ECFieldElement b) + { + return new FpFieldElement(q, x.Multiply(b.ToBigInteger().ModInverse(q)).Mod(q)); + } + + public override ECFieldElement Negate() + { + return new FpFieldElement(q, x.Negate().Mod(q)); + } + + public override ECFieldElement Square() + { + return new FpFieldElement(q, x.Multiply(x).Mod(q)); + } + + public override ECFieldElement Invert() + { + return new FpFieldElement(q, x.ModInverse(q)); + } + + // D.1.4 91 + /** + * return a sqrt root - the routine verifies that the calculation + * returns the right value - if none exists it returns null. + */ + public override ECFieldElement Sqrt() + { + if (!q.TestBit(0)) + throw Platform.CreateNotImplementedException("even value of q"); + + // p mod 4 == 3 + if (q.TestBit(1)) + { + // TODO Can this be optimised (inline the Square?) + // z = g^(u+1) + p, p = 4u + 3 + ECFieldElement z = new FpFieldElement(q, x.ModPow(q.ShiftRight(2).Add(BigInteger.One), q)); + + return this.Equals(z.Square()) ? z : null; + } + + // p mod 4 == 1 + BigInteger qMinusOne = q.Subtract(BigInteger.One); + + BigInteger legendreExponent = qMinusOne.ShiftRight(1); + if (!(x.ModPow(legendreExponent, q).Equals(BigInteger.One))) + return null; + + BigInteger u = qMinusOne.ShiftRight(2); + BigInteger k = u.ShiftLeft(1).Add(BigInteger.One); + + BigInteger Q = this.x; + BigInteger fourQ = Q.ShiftLeft(2).Mod(q); + + BigInteger U, V; + do + { + Random rand = new Random(); + BigInteger P; + do + { + P = new BigInteger(q.BitLength, rand); + } + while (P.CompareTo(q) >= 0 + || !(P.Multiply(P).Subtract(fourQ).ModPow(legendreExponent, q).Equals(qMinusOne))); + + BigInteger[] result = fastLucasSequence(q, P, Q, k); + U = result[0]; + V = result[1]; + + if (V.Multiply(V).Mod(q).Equals(fourQ)) + { + // Integer division by 2, mod q + if (V.TestBit(0)) + { + V = V.Add(q); + } + + V = V.ShiftRight(1); + + Debug.Assert(V.Multiply(V).Mod(q).Equals(x)); + + return new FpFieldElement(q, V); + } + } + while (U.Equals(BigInteger.One) || U.Equals(qMinusOne)); + + return null; + + +// BigInteger qMinusOne = q.Subtract(BigInteger.One); +// +// BigInteger legendreExponent = qMinusOne.ShiftRight(1); +// if (!(x.ModPow(legendreExponent, q).Equals(BigInteger.One))) +// return null; +// +// Random rand = new Random(); +// BigInteger fourX = x.ShiftLeft(2); +// +// BigInteger r; +// do +// { +// r = new BigInteger(q.BitLength, rand); +// } +// while (r.CompareTo(q) >= 0 +// || !(r.Multiply(r).Subtract(fourX).ModPow(legendreExponent, q).Equals(qMinusOne))); +// +// BigInteger n1 = qMinusOne.ShiftRight(2); +// BigInteger n2 = n1.Add(BigInteger.One); +// +// BigInteger wOne = WOne(r, x, q); +// BigInteger wSum = W(n1, wOne, q).Add(W(n2, wOne, q)).Mod(q); +// BigInteger twoR = r.ShiftLeft(1); +// +// BigInteger root = twoR.ModPow(q.Subtract(BigInteger.Two), q) +// .Multiply(x).Mod(q) +// .Multiply(wSum).Mod(q); +// +// return new FpFieldElement(q, root); + } + +// private static BigInteger W(BigInteger n, BigInteger wOne, BigInteger p) +// { +// if (n.Equals(BigInteger.One)) +// return wOne; +// +// bool isEven = !n.TestBit(0); +// n = n.ShiftRight(1); +// if (isEven) +// { +// BigInteger w = W(n, wOne, p); +// return w.Multiply(w).Subtract(BigInteger.Two).Mod(p); +// } +// BigInteger w1 = W(n.Add(BigInteger.One), wOne, p); +// BigInteger w2 = W(n, wOne, p); +// return w1.Multiply(w2).Subtract(wOne).Mod(p); +// } +// +// private BigInteger WOne(BigInteger r, BigInteger x, BigInteger p) +// { +// return r.Multiply(r).Multiply(x.ModPow(q.Subtract(BigInteger.Two), q)).Subtract(BigInteger.Two).Mod(p); +// } + + private static BigInteger[] fastLucasSequence( + BigInteger p, + BigInteger P, + BigInteger Q, + BigInteger k) + { + // TODO Research and apply "common-multiplicand multiplication here" + + int n = k.BitLength; + int s = k.GetLowestSetBit(); + + Debug.Assert(k.TestBit(s)); + + BigInteger Uh = BigInteger.One; + BigInteger Vl = BigInteger.Two; + BigInteger Vh = P; + BigInteger Ql = BigInteger.One; + BigInteger Qh = BigInteger.One; + + for (int j = n - 1; j >= s + 1; --j) + { + Ql = Ql.Multiply(Qh).Mod(p); + + if (k.TestBit(j)) + { + Qh = Ql.Multiply(Q).Mod(p); + Uh = Uh.Multiply(Vh).Mod(p); + Vl = Vh.Multiply(Vl).Subtract(P.Multiply(Ql)).Mod(p); + Vh = Vh.Multiply(Vh).Subtract(Qh.ShiftLeft(1)).Mod(p); + } + else + { + Qh = Ql; + Uh = Uh.Multiply(Vl).Subtract(Ql).Mod(p); + Vh = Vh.Multiply(Vl).Subtract(P.Multiply(Ql)).Mod(p); + Vl = Vl.Multiply(Vl).Subtract(Ql.ShiftLeft(1)).Mod(p); + } + } + + Ql = Ql.Multiply(Qh).Mod(p); + Qh = Ql.Multiply(Q).Mod(p); + Uh = Uh.Multiply(Vl).Subtract(Ql).Mod(p); + Vl = Vh.Multiply(Vl).Subtract(P.Multiply(Ql)).Mod(p); + Ql = Ql.Multiply(Qh).Mod(p); + + for (int j = 1; j <= s; ++j) + { + Uh = Uh.Multiply(Vl).Mod(p); + Vl = Vl.Multiply(Vl).Subtract(Ql.ShiftLeft(1)).Mod(p); + Ql = Ql.Multiply(Ql).Mod(p); + } + + return new BigInteger[]{ Uh, Vl }; + } + +// private static BigInteger[] verifyLucasSequence( +// BigInteger p, +// BigInteger P, +// BigInteger Q, +// BigInteger k) +// { +// BigInteger[] actual = fastLucasSequence(p, P, Q, k); +// BigInteger[] plus1 = fastLucasSequence(p, P, Q, k.Add(BigInteger.One)); +// BigInteger[] plus2 = fastLucasSequence(p, P, Q, k.Add(BigInteger.Two)); +// +// BigInteger[] check = stepLucasSequence(p, P, Q, actual, plus1); +// +// Debug.Assert(check[0].Equals(plus2[0])); +// Debug.Assert(check[1].Equals(plus2[1])); +// +// return actual; +// } +// +// private static BigInteger[] stepLucasSequence( +// BigInteger p, +// BigInteger P, +// BigInteger Q, +// BigInteger[] backTwo, +// BigInteger[] backOne) +// { +// return new BigInteger[] +// { +// P.Multiply(backOne[0]).Subtract(Q.Multiply(backTwo[0])).Mod(p), +// P.Multiply(backOne[1]).Subtract(Q.Multiply(backTwo[1])).Mod(p) +// }; +// } + + public override bool Equals( + object obj) + { + if (obj == this) + return true; + + FpFieldElement other = obj as FpFieldElement; + + if (other == null) + return false; + + return Equals(other); + } + + protected bool Equals( + FpFieldElement other) + { + return q.Equals(other.q) && base.Equals(other); + } + + public override int GetHashCode() + { + return q.GetHashCode() ^ base.GetHashCode(); + } + } + +// /** +// * Class representing the Elements of the finite field +// * <code>F<sub>2<sup>m</sup></sub></code> in polynomial basis (PB) +// * representation. Both trinomial (Tpb) and pentanomial (Ppb) polynomial +// * basis representations are supported. Gaussian normal basis (GNB) +// * representation is not supported. +// */ +// public class F2mFieldElement +// : ECFieldElement +// { +// /** +// * Indicates gaussian normal basis representation (GNB). Number chosen +// * according to X9.62. GNB is not implemented at present. +// */ +// public const int Gnb = 1; +// +// /** +// * Indicates trinomial basis representation (Tpb). Number chosen +// * according to X9.62. +// */ +// public const int Tpb = 2; +// +// /** +// * Indicates pentanomial basis representation (Ppb). Number chosen +// * according to X9.62. +// */ +// public const int Ppb = 3; +// +// /** +// * Tpb or Ppb. +// */ +// private int representation; +// +// /** +// * The exponent <code>m</code> of <code>F<sub>2<sup>m</sup></sub></code>. +// */ +// private int m; +// +// /** +// * Tpb: The integer <code>k</code> where <code>x<sup>m</sup> + +// * x<sup>k</sup> + 1</code> represents the reduction polynomial +// * <code>f(z)</code>.<br/> +// * Ppb: The integer <code>k1</code> where <code>x<sup>m</sup> + +// * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> +// * represents the reduction polynomial <code>f(z)</code>.<br/> +// */ +// private int k1; +// +// /** +// * Tpb: Always set to <code>0</code><br/> +// * Ppb: The integer <code>k2</code> where <code>x<sup>m</sup> + +// * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> +// * represents the reduction polynomial <code>f(z)</code>.<br/> +// */ +// private int k2; +// +// /** +// * Tpb: Always set to <code>0</code><br/> +// * Ppb: The integer <code>k3</code> where <code>x<sup>m</sup> + +// * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> +// * represents the reduction polynomial <code>f(z)</code>.<br/> +// */ +// private int k3; +// +// /** +// * Constructor for Ppb. +// * @param m The exponent <code>m</code> of +// * <code>F<sub>2<sup>m</sup></sub></code>. +// * @param k1 The integer <code>k1</code> where <code>x<sup>m</sup> + +// * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> +// * represents the reduction polynomial <code>f(z)</code>. +// * @param k2 The integer <code>k2</code> where <code>x<sup>m</sup> + +// * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> +// * represents the reduction polynomial <code>f(z)</code>. +// * @param k3 The integer <code>k3</code> where <code>x<sup>m</sup> + +// * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> +// * represents the reduction polynomial <code>f(z)</code>. +// * @param x The BigInteger representing the value of the field element. +// */ +// public F2mFieldElement( +// int m, +// int k1, +// int k2, +// int k3, +// BigInteger x) +// : base(x) +// { +// if ((k2 == 0) && (k3 == 0)) +// { +// this.representation = Tpb; +// } +// else +// { +// if (k2 >= k3) +// throw new ArgumentException("k2 must be smaller than k3"); +// if (k2 <= 0) +// throw new ArgumentException("k2 must be larger than 0"); +// +// this.representation = Ppb; +// } +// +// if (x.SignValue < 0) +// throw new ArgumentException("x value cannot be negative"); +// +// this.m = m; +// this.k1 = k1; +// this.k2 = k2; +// this.k3 = k3; +// } +// +// /** +// * Constructor for Tpb. +// * @param m The exponent <code>m</code> of +// * <code>F<sub>2<sup>m</sup></sub></code>. +// * @param k The integer <code>k</code> where <code>x<sup>m</sup> + +// * x<sup>k</sup> + 1</code> represents the reduction +// * polynomial <code>f(z)</code>. +// * @param x The BigInteger representing the value of the field element. +// */ +// public F2mFieldElement( +// int m, +// int k, +// BigInteger x) +// : this(m, k, 0, 0, x) +// { +// // Set k1 to k, and set k2 and k3 to 0 +// } +// +// public override string FieldName +// { +// get { return "F2m"; } +// } +// +// /** +// * Checks, if the ECFieldElements <code>a</code> and <code>b</code> +// * are elements of the same field <code>F<sub>2<sup>m</sup></sub></code> +// * (having the same representation). +// * @param a field element. +// * @param b field element to be compared. +// * @throws ArgumentException if <code>a</code> and <code>b</code> +// * are not elements of the same field +// * <code>F<sub>2<sup>m</sup></sub></code> (having the same +// * representation). +// */ +// public static void CheckFieldElements( +// ECFieldElement a, +// ECFieldElement b) +// { +// if (!(a is F2mFieldElement) || !(b is F2mFieldElement)) +// { +// throw new ArgumentException("Field elements are not " +// + "both instances of F2mFieldElement"); +// } +// +// if ((a.x.SignValue < 0) || (b.x.SignValue < 0)) +// { +// throw new ArgumentException( +// "x value may not be negative"); +// } +// +// F2mFieldElement aF2m = (F2mFieldElement)a; +// F2mFieldElement bF2m = (F2mFieldElement)b; +// +// if ((aF2m.m != bF2m.m) || (aF2m.k1 != bF2m.k1) +// || (aF2m.k2 != bF2m.k2) || (aF2m.k3 != bF2m.k3)) +// { +// throw new ArgumentException("Field elements are not " +// + "elements of the same field F2m"); +// } +// +// if (aF2m.representation != bF2m.representation) +// { +// // Should never occur +// throw new ArgumentException( +// "One of the field " +// + "elements are not elements has incorrect representation"); +// } +// } +// +// /** +// * Computes <code>z * a(z) mod f(z)</code>, where <code>f(z)</code> is +// * the reduction polynomial of <code>this</code>. +// * @param a The polynomial <code>a(z)</code> to be multiplied by +// * <code>z mod f(z)</code>. +// * @return <code>z * a(z) mod f(z)</code> +// */ +// private BigInteger multZModF( +// BigInteger a) +// { +// // Left-shift of a(z) +// BigInteger az = a.ShiftLeft(1); +// if (az.TestBit(this.m)) +// { +// // If the coefficient of z^m in a(z) Equals 1, reduction +// // modulo f(z) is performed: Add f(z) to to a(z): +// // Step 1: Unset mth coeffient of a(z) +// az = az.ClearBit(this.m); +// +// // Step 2: Add r(z) to a(z), where r(z) is defined as +// // f(z) = z^m + r(z), and k1, k2, k3 are the positions of +// // the non-zero coefficients in r(z) +// az = az.FlipBit(0); +// az = az.FlipBit(this.k1); +// if (this.representation == Ppb) +// { +// az = az.FlipBit(this.k2); +// az = az.FlipBit(this.k3); +// } +// } +// return az; +// } +// +// public override ECFieldElement Add( +// ECFieldElement b) +// { +// // No check performed here for performance reasons. Instead the +// // elements involved are checked in ECPoint.F2m +// // checkFieldElements(this, b); +// if (b.x.SignValue == 0) +// return this; +// +// return new F2mFieldElement(this.m, this.k1, this.k2, this.k3, this.x.Xor(b.x)); +// } +// +// public override ECFieldElement Subtract( +// ECFieldElement b) +// { +// // Addition and subtraction are the same in F2m +// return Add(b); +// } +// +// public override ECFieldElement Multiply( +// ECFieldElement b) +// { +// // Left-to-right shift-and-add field multiplication in F2m +// // Input: Binary polynomials a(z) and b(z) of degree at most m-1 +// // Output: c(z) = a(z) * b(z) mod f(z) +// +// // No check performed here for performance reasons. Instead the +// // elements involved are checked in ECPoint.F2m +// // checkFieldElements(this, b); +// BigInteger az = this.x; +// BigInteger bz = b.x; +// BigInteger cz; +// +// // Compute c(z) = a(z) * b(z) mod f(z) +// if (az.TestBit(0)) +// { +// cz = bz; +// } +// else +// { +// cz = BigInteger.Zero; +// } +// +// for (int i = 1; i < this.m; i++) +// { +// // b(z) := z * b(z) mod f(z) +// bz = multZModF(bz); +// +// if (az.TestBit(i)) +// { +// // If the coefficient of x^i in a(z) Equals 1, b(z) is added +// // to c(z) +// cz = cz.Xor(bz); +// } +// } +// return new F2mFieldElement(m, this.k1, this.k2, this.k3, cz); +// } +// +// +// public override ECFieldElement Divide( +// ECFieldElement b) +// { +// // There may be more efficient implementations +// ECFieldElement bInv = b.Invert(); +// return Multiply(bInv); +// } +// +// public override ECFieldElement Negate() +// { +// // -x == x holds for all x in F2m +// return this; +// } +// +// public override ECFieldElement Square() +// { +// // Naive implementation, can probably be speeded up using modular +// // reduction +// return Multiply(this); +// } +// +// public override ECFieldElement Invert() +// { +// // Inversion in F2m using the extended Euclidean algorithm +// // Input: A nonzero polynomial a(z) of degree at most m-1 +// // Output: a(z)^(-1) mod f(z) +// +// // u(z) := a(z) +// BigInteger uz = this.x; +// if (uz.SignValue <= 0) +// { +// throw new ArithmeticException("x is zero or negative, " + +// "inversion is impossible"); +// } +// +// // v(z) := f(z) +// BigInteger vz = BigInteger.One.ShiftLeft(m); +// vz = vz.SetBit(0); +// vz = vz.SetBit(this.k1); +// if (this.representation == Ppb) +// { +// vz = vz.SetBit(this.k2); +// vz = vz.SetBit(this.k3); +// } +// +// // g1(z) := 1, g2(z) := 0 +// BigInteger g1z = BigInteger.One; +// BigInteger g2z = BigInteger.Zero; +// +// // while u != 1 +// while (uz.SignValue != 0) +// { +// // j := deg(u(z)) - deg(v(z)) +// int j = uz.BitLength - vz.BitLength; +// +// // If j < 0 then: u(z) <-> v(z), g1(z) <-> g2(z), j := -j +// if (j < 0) +// { +// BigInteger uzCopy = uz; +// uz = vz; +// vz = uzCopy; +// +// BigInteger g1zCopy = g1z; +// g1z = g2z; +// g2z = g1zCopy; +// +// j = -j; +// } +// +// // u(z) := u(z) + z^j * v(z) +// // Note, that no reduction modulo f(z) is required, because +// // deg(u(z) + z^j * v(z)) <= max(deg(u(z)), j + deg(v(z))) +// // = max(deg(u(z)), deg(u(z)) - deg(v(z)) + deg(v(z)) +// // = deg(u(z)) +// uz = uz.Xor(vz.ShiftLeft(j)); +// +// // g1(z) := g1(z) + z^j * g2(z) +// g1z = g1z.Xor(g2z.ShiftLeft(j)); +// // if (g1z.BitLength() > this.m) { +// // throw new ArithmeticException( +// // "deg(g1z) >= m, g1z = " + g1z.ToString(2)); +// // } +// } +// return new F2mFieldElement(this.m, this.k1, this.k2, this.k3, g2z); +// } +// +// public override ECFieldElement Sqrt() +// { +// throw new ArithmeticException("Not implemented"); +// } +// +// /** +// * @return the representation of the field +// * <code>F<sub>2<sup>m</sup></sub></code>, either of +// * {@link F2mFieldElement.Tpb} (trinomial +// * basis representation) or +// * {@link F2mFieldElement.Ppb} (pentanomial +// * basis representation). +// */ +// public int Representation +// { +// get { return this.representation; } +// } +// +// /** +// * @return the degree <code>m</code> of the reduction polynomial +// * <code>f(z)</code>. +// */ +// public int M +// { +// get { return this.m; } +// } +// +// /** +// * @return Tpb: The integer <code>k</code> where <code>x<sup>m</sup> + +// * x<sup>k</sup> + 1</code> represents the reduction polynomial +// * <code>f(z)</code>.<br/> +// * Ppb: The integer <code>k1</code> where <code>x<sup>m</sup> + +// * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> +// * represents the reduction polynomial <code>f(z)</code>.<br/> +// */ +// public int K1 +// { +// get { return this.k1; } +// } +// +// /** +// * @return Tpb: Always returns <code>0</code><br/> +// * Ppb: The integer <code>k2</code> where <code>x<sup>m</sup> + +// * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> +// * represents the reduction polynomial <code>f(z)</code>.<br/> +// */ +// public int K2 +// { +// get { return this.k2; } +// } +// +// /** +// * @return Tpb: Always set to <code>0</code><br/> +// * Ppb: The integer <code>k3</code> where <code>x<sup>m</sup> + +// * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> +// * represents the reduction polynomial <code>f(z)</code>.<br/> +// */ +// public int K3 +// { +// get { return this.k3; } +// } +// +// public override bool Equals( +// object obj) +// { +// if (obj == this) +// return true; +// +// F2mFieldElement other = obj as F2mFieldElement; +// +// if (other == null) +// return false; +// +// return Equals(other); +// } +// +// protected bool Equals( +// F2mFieldElement other) +// { +// return m == other.m +// && k1 == other.k1 +// && k2 == other.k2 +// && k3 == other.k3 +// && representation == other.representation +// && base.Equals(other); +// } +// +// public override int GetHashCode() +// { +// return base.GetHashCode() ^ m ^ k1 ^ k2 ^ k3; +// } +// } + + /** + * Class representing the Elements of the finite field + * <code>F<sub>2<sup>m</sup></sub></code> in polynomial basis (PB) + * representation. Both trinomial (Tpb) and pentanomial (Ppb) polynomial + * basis representations are supported. Gaussian normal basis (GNB) + * representation is not supported. + */ + public class F2mFieldElement + : ECFieldElement + { + /** + * Indicates gaussian normal basis representation (GNB). Number chosen + * according to X9.62. GNB is not implemented at present. + */ + public const int Gnb = 1; + + /** + * Indicates trinomial basis representation (Tpb). Number chosen + * according to X9.62. + */ + public const int Tpb = 2; + + /** + * Indicates pentanomial basis representation (Ppb). Number chosen + * according to X9.62. + */ + public const int Ppb = 3; + + /** + * Tpb or Ppb. + */ + private int representation; + + /** + * The exponent <code>m</code> of <code>F<sub>2<sup>m</sup></sub></code>. + */ + private int m; + + /** + * Tpb: The integer <code>k</code> where <code>x<sup>m</sup> + + * x<sup>k</sup> + 1</code> represents the reduction polynomial + * <code>f(z)</code>.<br/> + * Ppb: The integer <code>k1</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>.<br/> + */ + private int k1; + + /** + * Tpb: Always set to <code>0</code><br/> + * Ppb: The integer <code>k2</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>.<br/> + */ + private int k2; + + /** + * Tpb: Always set to <code>0</code><br/> + * Ppb: The integer <code>k3</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>.<br/> + */ + private int k3; + + /** + * The <code>IntArray</code> holding the bits. + */ + private IntArray x; + + /** + * The number of <code>int</code>s required to hold <code>m</code> bits. + */ + private readonly int t; + + /** + * Constructor for Ppb. + * @param m The exponent <code>m</code> of + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param k1 The integer <code>k1</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>. + * @param k2 The integer <code>k2</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>. + * @param k3 The integer <code>k3</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>. + * @param x The BigInteger representing the value of the field element. + */ + public F2mFieldElement( + int m, + int k1, + int k2, + int k3, + BigInteger x) + { + // t = m / 32 rounded up to the next integer + this.t = (m + 31) >> 5; + this.x = new IntArray(x, t); + + if ((k2 == 0) && (k3 == 0)) + { + this.representation = Tpb; + } + else + { + if (k2 >= k3) + throw new ArgumentException("k2 must be smaller than k3"); + if (k2 <= 0) + throw new ArgumentException("k2 must be larger than 0"); + + this.representation = Ppb; + } + + if (x.SignValue < 0) + throw new ArgumentException("x value cannot be negative"); + + this.m = m; + this.k1 = k1; + this.k2 = k2; + this.k3 = k3; + } + + /** + * Constructor for Tpb. + * @param m The exponent <code>m</code> of + * <code>F<sub>2<sup>m</sup></sub></code>. + * @param k The integer <code>k</code> where <code>x<sup>m</sup> + + * x<sup>k</sup> + 1</code> represents the reduction + * polynomial <code>f(z)</code>. + * @param x The BigInteger representing the value of the field element. + */ + public F2mFieldElement( + int m, + int k, + BigInteger x) + : this(m, k, 0, 0, x) + { + // Set k1 to k, and set k2 and k3 to 0 + } + + private F2mFieldElement(int m, int k1, int k2, int k3, IntArray x) + { + t = (m + 31) >> 5; + this.x = x; + this.m = m; + this.k1 = k1; + this.k2 = k2; + this.k3 = k3; + + if ((k2 == 0) && (k3 == 0)) + { + this.representation = Tpb; + } + else + { + this.representation = Ppb; + } + } + + public override BigInteger ToBigInteger() + { + return x.ToBigInteger(); + } + + public override string FieldName + { + get { return "F2m"; } + } + + public override int FieldSize + { + get { return m; } + } + + /** + * Checks, if the ECFieldElements <code>a</code> and <code>b</code> + * are elements of the same field <code>F<sub>2<sup>m</sup></sub></code> + * (having the same representation). + * @param a field element. + * @param b field element to be compared. + * @throws ArgumentException if <code>a</code> and <code>b</code> + * are not elements of the same field + * <code>F<sub>2<sup>m</sup></sub></code> (having the same + * representation). + */ + public static void CheckFieldElements( + ECFieldElement a, + ECFieldElement b) + { + if (!(a is F2mFieldElement) || !(b is F2mFieldElement)) + { + throw new ArgumentException("Field elements are not " + + "both instances of F2mFieldElement"); + } + + F2mFieldElement aF2m = (F2mFieldElement)a; + F2mFieldElement bF2m = (F2mFieldElement)b; + + if ((aF2m.m != bF2m.m) || (aF2m.k1 != bF2m.k1) + || (aF2m.k2 != bF2m.k2) || (aF2m.k3 != bF2m.k3)) + { + throw new ArgumentException("Field elements are not " + + "elements of the same field F2m"); + } + + if (aF2m.representation != bF2m.representation) + { + // Should never occur + throw new ArgumentException( + "One of the field " + + "elements are not elements has incorrect representation"); + } + } + + public override ECFieldElement Add( + ECFieldElement b) + { + // No check performed here for performance reasons. Instead the + // elements involved are checked in ECPoint.F2m + // checkFieldElements(this, b); + IntArray iarrClone = (IntArray) this.x.Copy(); + F2mFieldElement bF2m = (F2mFieldElement) b; + iarrClone.AddShifted(bF2m.x, 0); + return new F2mFieldElement(m, k1, k2, k3, iarrClone); + } + + public override ECFieldElement Subtract( + ECFieldElement b) + { + // Addition and subtraction are the same in F2m + return Add(b); + } + + public override ECFieldElement Multiply( + ECFieldElement b) + { + // Right-to-left comb multiplication in the IntArray + // Input: Binary polynomials a(z) and b(z) of degree at most m-1 + // Output: c(z) = a(z) * b(z) mod f(z) + + // No check performed here for performance reasons. Instead the + // elements involved are checked in ECPoint.F2m + // checkFieldElements(this, b); + F2mFieldElement bF2m = (F2mFieldElement) b; + IntArray mult = x.Multiply(bF2m.x, m); + mult.Reduce(m, new int[]{k1, k2, k3}); + return new F2mFieldElement(m, k1, k2, k3, mult); + } + + public override ECFieldElement Divide( + ECFieldElement b) + { + // There may be more efficient implementations + ECFieldElement bInv = b.Invert(); + return Multiply(bInv); + } + + public override ECFieldElement Negate() + { + // -x == x holds for all x in F2m + return this; + } + + public override ECFieldElement Square() + { + IntArray squared = x.Square(m); + squared.Reduce(m, new int[]{k1, k2, k3}); + return new F2mFieldElement(m, k1, k2, k3, squared); + } + + public override ECFieldElement Invert() + { + // Inversion in F2m using the extended Euclidean algorithm + // Input: A nonzero polynomial a(z) of degree at most m-1 + // Output: a(z)^(-1) mod f(z) + + // u(z) := a(z) + IntArray uz = (IntArray)this.x.Copy(); + + // v(z) := f(z) + IntArray vz = new IntArray(t); + vz.SetBit(m); + vz.SetBit(0); + vz.SetBit(this.k1); + if (this.representation == Ppb) + { + vz.SetBit(this.k2); + vz.SetBit(this.k3); + } + + // g1(z) := 1, g2(z) := 0 + IntArray g1z = new IntArray(t); + g1z.SetBit(0); + IntArray g2z = new IntArray(t); + + // while u != 0 + while (uz.GetUsedLength() > 0) +// while (uz.bitLength() > 1) + { + // j := deg(u(z)) - deg(v(z)) + int j = uz.BitLength - vz.BitLength; + + // If j < 0 then: u(z) <-> v(z), g1(z) <-> g2(z), j := -j + if (j < 0) + { + IntArray uzCopy = uz; + uz = vz; + vz = uzCopy; + + IntArray g1zCopy = g1z; + g1z = g2z; + g2z = g1zCopy; + + j = -j; + } + + // u(z) := u(z) + z^j * v(z) + // Note, that no reduction modulo f(z) is required, because + // deg(u(z) + z^j * v(z)) <= max(deg(u(z)), j + deg(v(z))) + // = max(deg(u(z)), deg(u(z)) - deg(v(z)) + deg(v(z)) + // = deg(u(z)) + // uz = uz.xor(vz.ShiftLeft(j)); + // jInt = n / 32 + int jInt = j >> 5; + // jInt = n % 32 + int jBit = j & 0x1F; + IntArray vzShift = vz.ShiftLeft(jBit); + uz.AddShifted(vzShift, jInt); + + // g1(z) := g1(z) + z^j * g2(z) +// g1z = g1z.xor(g2z.ShiftLeft(j)); + IntArray g2zShift = g2z.ShiftLeft(jBit); + g1z.AddShifted(g2zShift, jInt); + } + return new F2mFieldElement(this.m, this.k1, this.k2, this.k3, g2z); + } + + public override ECFieldElement Sqrt() + { + throw new ArithmeticException("Not implemented"); + } + + /** + * @return the representation of the field + * <code>F<sub>2<sup>m</sup></sub></code>, either of + * {@link F2mFieldElement.Tpb} (trinomial + * basis representation) or + * {@link F2mFieldElement.Ppb} (pentanomial + * basis representation). + */ + public int Representation + { + get { return this.representation; } + } + + /** + * @return the degree <code>m</code> of the reduction polynomial + * <code>f(z)</code>. + */ + public int M + { + get { return this.m; } + } + + /** + * @return Tpb: The integer <code>k</code> where <code>x<sup>m</sup> + + * x<sup>k</sup> + 1</code> represents the reduction polynomial + * <code>f(z)</code>.<br/> + * Ppb: The integer <code>k1</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>.<br/> + */ + public int K1 + { + get { return this.k1; } + } + + /** + * @return Tpb: Always returns <code>0</code><br/> + * Ppb: The integer <code>k2</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>.<br/> + */ + public int K2 + { + get { return this.k2; } + } + + /** + * @return Tpb: Always set to <code>0</code><br/> + * Ppb: The integer <code>k3</code> where <code>x<sup>m</sup> + + * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code> + * represents the reduction polynomial <code>f(z)</code>.<br/> + */ + public int K3 + { + get { return this.k3; } + } + + public override bool Equals( + object obj) + { + if (obj == this) + return true; + + F2mFieldElement other = obj as F2mFieldElement; + + if (other == null) + return false; + + return Equals(other); + } + + protected bool Equals( + F2mFieldElement other) + { + return m == other.m + && k1 == other.k1 + && k2 == other.k2 + && k3 == other.k3 + && representation == other.representation + && base.Equals(other); + } + + public override int GetHashCode() + { + return m.GetHashCode() + ^ k1.GetHashCode() + ^ k2.GetHashCode() + ^ k3.GetHashCode() + ^ representation.GetHashCode() + ^ base.GetHashCode(); + } + } +} diff --git a/crypto/src/math/ec/ECPoint.cs b/crypto/src/math/ec/ECPoint.cs new file mode 100644 index 000000000..6a06be4d1 --- /dev/null +++ b/crypto/src/math/ec/ECPoint.cs @@ -0,0 +1,572 @@ +using System; +using System.Collections; +using System.Diagnostics; + +using Org.BouncyCastle.Asn1.X9; + +using Org.BouncyCastle.Math.EC.Multiplier; + +namespace Org.BouncyCastle.Math.EC +{ + /** + * base class for points on elliptic curves. + */ + public abstract class ECPoint + { + internal readonly ECCurve curve; + internal readonly ECFieldElement x, y; + internal readonly bool withCompression; + internal ECMultiplier multiplier = null; + internal PreCompInfo preCompInfo = null; + + protected internal ECPoint( + ECCurve curve, + ECFieldElement x, + ECFieldElement y, + bool withCompression) + { + if (curve == null) + throw new ArgumentNullException("curve"); + + this.curve = curve; + this.x = x; + this.y = y; + this.withCompression = withCompression; + } + + public ECCurve Curve + { + get { return curve; } + } + + public ECFieldElement X + { + get { return x; } + } + + public ECFieldElement Y + { + get { return y; } + } + + public bool IsInfinity + { + get { return x == null && y == null; } + } + + public bool IsCompressed + { + get { return withCompression; } + } + + public override bool Equals( + object obj) + { + if (obj == this) + return true; + + ECPoint o = obj as ECPoint; + + if (o == null) + return false; + + if (this.IsInfinity) + return o.IsInfinity; + + return x.Equals(o.x) && y.Equals(o.y); + } + + public override int GetHashCode() + { + if (this.IsInfinity) + return 0; + + return x.GetHashCode() ^ y.GetHashCode(); + } + +// /** +// * Mainly for testing. Explicitly set the <code>ECMultiplier</code>. +// * @param multiplier The <code>ECMultiplier</code> to be used to multiply +// * this <code>ECPoint</code>. +// */ +// internal void SetECMultiplier( +// ECMultiplier multiplier) +// { +// this.multiplier = multiplier; +// } + + /** + * Sets the <code>PreCompInfo</code>. Used by <code>ECMultiplier</code>s + * to save the precomputation for this <code>ECPoint</code> to store the + * precomputation result for use by subsequent multiplication. + * @param preCompInfo The values precomputed by the + * <code>ECMultiplier</code>. + */ + internal void SetPreCompInfo( + PreCompInfo preCompInfo) + { + this.preCompInfo = preCompInfo; + } + + public virtual byte[] GetEncoded() + { + return GetEncoded(withCompression); + } + + public abstract byte[] GetEncoded(bool compressed); + + public abstract ECPoint Add(ECPoint b); + public abstract ECPoint Subtract(ECPoint b); + public abstract ECPoint Negate(); + public abstract ECPoint Twice(); + public abstract ECPoint Multiply(BigInteger b); + + /** + * Sets the appropriate <code>ECMultiplier</code>, unless already set. + */ + internal virtual void AssertECMultiplier() + { + if (this.multiplier == null) + { + lock (this) + { + if (this.multiplier == null) + { + this.multiplier = new FpNafMultiplier(); + } + } + } + } + } + + public abstract class ECPointBase + : ECPoint + { + protected internal ECPointBase( + ECCurve curve, + ECFieldElement x, + ECFieldElement y, + bool withCompression) + : base(curve, x, y, withCompression) + { + } + + protected internal abstract bool YTilde { get; } + + /** + * return the field element encoded with point compression. (S 4.3.6) + */ + public override byte[] GetEncoded(bool compressed) + { + if (this.IsInfinity) + return new byte[1]; + + // Note: some of the tests rely on calculating byte length from the field element + // (since the test cases use mismatching fields for curve/elements) + int byteLength = X9IntegerConverter.GetByteLength(x); + byte[] X = X9IntegerConverter.IntegerToBytes(this.X.ToBigInteger(), byteLength); + byte[] PO; + + if (compressed) + { + PO = new byte[1 + X.Length]; + + PO[0] = (byte)(YTilde ? 0x03 : 0x02); + } + else + { + byte[] Y = X9IntegerConverter.IntegerToBytes(this.Y.ToBigInteger(), byteLength); + PO = new byte[1 + X.Length + Y.Length]; + + PO[0] = 0x04; + + Y.CopyTo(PO, 1 + X.Length); + } + + X.CopyTo(PO, 1); + + return PO; + } + + /** + * Multiplies this <code>ECPoint</code> by the given number. + * @param k The multiplicator. + * @return <code>k * this</code>. + */ + public override ECPoint Multiply( + BigInteger k) + { + if (k.SignValue < 0) + throw new ArgumentException("The multiplicator cannot be negative", "k"); + + if (this.IsInfinity) + return this; + + if (k.SignValue == 0) + return this.curve.Infinity; + + AssertECMultiplier(); + return this.multiplier.Multiply(this, k, preCompInfo); + } + } + + /** + * Elliptic curve points over Fp + */ + public class FpPoint + : ECPointBase + { + /** + * Create a point which encodes with point compression. + * + * @param curve the curve to use + * @param x affine x co-ordinate + * @param y affine y co-ordinate + */ + public FpPoint( + ECCurve curve, + ECFieldElement x, + ECFieldElement y) + : this(curve, x, y, false) + { + } + + /** + * Create a point that encodes with or without point compresion. + * + * @param curve the curve to use + * @param x affine x co-ordinate + * @param y affine y co-ordinate + * @param withCompression if true encode with point compression + */ + public FpPoint( + ECCurve curve, + ECFieldElement x, + ECFieldElement y, + bool withCompression) + : base(curve, x, y, withCompression) + { + if ((x == null) != (y == null)) + throw new ArgumentException("Exactly one of the field elements is null"); + } + + protected internal override bool YTilde + { + get + { + return this.Y.ToBigInteger().TestBit(0); + } + } + + // B.3 pg 62 + public override ECPoint Add( + ECPoint b) + { + if (this.IsInfinity) + return b; + + if (b.IsInfinity) + return this; + + // Check if b = this or b = -this + if (this.x.Equals(b.x)) + { + if (this.y.Equals(b.y)) + { + // this = b, i.e. this must be doubled + return this.Twice(); + } + + Debug.Assert(this.y.Equals(b.y.Negate())); + + // this = -b, i.e. the result is the point at infinity + return this.curve.Infinity; + } + + ECFieldElement gamma = b.y.Subtract(this.y).Divide(b.x.Subtract(this.x)); + + ECFieldElement x3 = gamma.Square().Subtract(this.x).Subtract(b.x); + ECFieldElement y3 = gamma.Multiply(this.x.Subtract(x3)).Subtract(this.y); + + return new FpPoint(curve, x3, y3, withCompression); + } + + // B.3 pg 62 + public override ECPoint Twice() + { + // Twice identity element (point at infinity) is identity + if (this.IsInfinity) + return this; + + // if y1 == 0, then (x1, y1) == (x1, -y1) + // and hence this = -this and thus 2(x1, y1) == infinity + if (this.y.ToBigInteger().SignValue == 0) + return this.curve.Infinity; + + ECFieldElement TWO = this.curve.FromBigInteger(BigInteger.Two); + ECFieldElement THREE = this.curve.FromBigInteger(BigInteger.Three); + ECFieldElement gamma = this.x.Square().Multiply(THREE).Add(curve.a).Divide(y.Multiply(TWO)); + + ECFieldElement x3 = gamma.Square().Subtract(this.x.Multiply(TWO)); + ECFieldElement y3 = gamma.Multiply(this.x.Subtract(x3)).Subtract(this.y); + + return new FpPoint(curve, x3, y3, this.withCompression); + } + + // D.3.2 pg 102 (see Note:) + public override ECPoint Subtract( + ECPoint b) + { + if (b.IsInfinity) + return this; + + // Add -b + return Add(b.Negate()); + } + + public override ECPoint Negate() + { + return new FpPoint(this.curve, this.x, this.y.Negate(), this.withCompression); + } + + /** + * Sets the default <code>ECMultiplier</code>, unless already set. + */ + internal override void AssertECMultiplier() + { + if (this.multiplier == null) + { + lock (this) + { + if (this.multiplier == null) + { + this.multiplier = new WNafMultiplier(); + } + } + } + } + } + + /** + * Elliptic curve points over F2m + */ + public class F2mPoint + : ECPointBase + { + /** + * @param curve base curve + * @param x x point + * @param y y point + */ + public F2mPoint( + ECCurve curve, + ECFieldElement x, + ECFieldElement y) + : this(curve, x, y, false) + { + } + + /** + * @param curve base curve + * @param x x point + * @param y y point + * @param withCompression true if encode with point compression. + */ + public F2mPoint( + ECCurve curve, + ECFieldElement x, + ECFieldElement y, + bool withCompression) + : base(curve, x, y, withCompression) + { + if ((x != null && y == null) || (x == null && y != null)) + { + throw new ArgumentException("Exactly one of the field elements is null"); + } + + if (x != null) + { + // Check if x and y are elements of the same field + F2mFieldElement.CheckFieldElements(this.x, this.y); + + // Check if x and a are elements of the same field + F2mFieldElement.CheckFieldElements(this.x, this.curve.A); + } + } + + /** + * Constructor for point at infinity + */ + [Obsolete("Use ECCurve.Infinity property")] + public F2mPoint( + ECCurve curve) + : this(curve, null, null) + { + } + + protected internal override bool YTilde + { + get + { + // X9.62 4.2.2 and 4.3.6: + // if x = 0 then ypTilde := 0, else ypTilde is the rightmost + // bit of y * x^(-1) + return this.X.ToBigInteger().SignValue != 0 + && this.Y.Multiply(this.X.Invert()).ToBigInteger().TestBit(0); + } + } + + /** + * Check, if two <code>ECPoint</code>s can be added or subtracted. + * @param a The first <code>ECPoint</code> to check. + * @param b The second <code>ECPoint</code> to check. + * @throws IllegalArgumentException if <code>a</code> and <code>b</code> + * cannot be added. + */ + private static void CheckPoints( + ECPoint a, + ECPoint b) + { + // Check, if points are on the same curve + if (!a.curve.Equals(b.curve)) + throw new ArgumentException("Only points on the same curve can be added or subtracted"); + +// F2mFieldElement.CheckFieldElements(a.x, b.x); + } + + /* (non-Javadoc) + * @see org.bouncycastle.math.ec.ECPoint#add(org.bouncycastle.math.ec.ECPoint) + */ + public override ECPoint Add(ECPoint b) + { + CheckPoints(this, b); + return AddSimple((F2mPoint) b); + } + + /** + * Adds another <code>ECPoints.F2m</code> to <code>this</code> without + * checking if both points are on the same curve. Used by multiplication + * algorithms, because there all points are a multiple of the same point + * and hence the checks can be omitted. + * @param b The other <code>ECPoints.F2m</code> to add to + * <code>this</code>. + * @return <code>this + b</code> + */ + internal F2mPoint AddSimple(F2mPoint b) + { + if (this.IsInfinity) + return b; + + if (b.IsInfinity) + return this; + + F2mFieldElement x2 = (F2mFieldElement) b.X; + F2mFieldElement y2 = (F2mFieldElement) b.Y; + + // Check if b == this or b == -this + if (this.x.Equals(x2)) + { + // this == b, i.e. this must be doubled + if (this.y.Equals(y2)) + return (F2mPoint) this.Twice(); + + // this = -other, i.e. the result is the point at infinity + return (F2mPoint) this.curve.Infinity; + } + + ECFieldElement xSum = this.x.Add(x2); + + F2mFieldElement lambda + = (F2mFieldElement)(this.y.Add(y2)).Divide(xSum); + + F2mFieldElement x3 + = (F2mFieldElement)lambda.Square().Add(lambda).Add(xSum).Add(this.curve.A); + + F2mFieldElement y3 + = (F2mFieldElement)lambda.Multiply(this.x.Add(x3)).Add(x3).Add(this.y); + + return new F2mPoint(curve, x3, y3, withCompression); + } + + /* (non-Javadoc) + * @see org.bouncycastle.math.ec.ECPoint#subtract(org.bouncycastle.math.ec.ECPoint) + */ + public override ECPoint Subtract( + ECPoint b) + { + CheckPoints(this, b); + return SubtractSimple((F2mPoint) b); + } + + /** + * Subtracts another <code>ECPoints.F2m</code> from <code>this</code> + * without checking if both points are on the same curve. Used by + * multiplication algorithms, because there all points are a multiple + * of the same point and hence the checks can be omitted. + * @param b The other <code>ECPoints.F2m</code> to subtract from + * <code>this</code>. + * @return <code>this - b</code> + */ + internal F2mPoint SubtractSimple( + F2mPoint b) + { + if (b.IsInfinity) + return this; + + // Add -b + return AddSimple((F2mPoint) b.Negate()); + } + + /* (non-Javadoc) + * @see Org.BouncyCastle.Math.EC.ECPoint#twice() + */ + public override ECPoint Twice() + { + // Twice identity element (point at infinity) is identity + if (this.IsInfinity) + return this; + + // if x1 == 0, then (x1, y1) == (x1, x1 + y1) + // and hence this = -this and thus 2(x1, y1) == infinity + if (this.x.ToBigInteger().SignValue == 0) + return this.curve.Infinity; + + F2mFieldElement lambda = (F2mFieldElement) this.x.Add(this.y.Divide(this.x)); + F2mFieldElement x2 = (F2mFieldElement)lambda.Square().Add(lambda).Add(this.curve.A); + ECFieldElement ONE = this.curve.FromBigInteger(BigInteger.One); + F2mFieldElement y2 = (F2mFieldElement)this.x.Square().Add( + x2.Multiply(lambda.Add(ONE))); + + return new F2mPoint(this.curve, x2, y2, withCompression); + } + + public override ECPoint Negate() + { + return new F2mPoint(curve, this.x, this.x.Add(this.y), withCompression); + } + + /** + * Sets the appropriate <code>ECMultiplier</code>, unless already set. + */ + internal override void AssertECMultiplier() + { + if (this.multiplier == null) + { + lock (this) + { + if (this.multiplier == null) + { + if (((F2mCurve) this.curve).IsKoblitz) + { + this.multiplier = new WTauNafMultiplier(); + } + else + { + this.multiplier = new WNafMultiplier(); + } + } + } + } + } + } +} diff --git a/crypto/src/math/ec/IntArray.cs b/crypto/src/math/ec/IntArray.cs new file mode 100644 index 000000000..1089966a8 --- /dev/null +++ b/crypto/src/math/ec/IntArray.cs @@ -0,0 +1,485 @@ +using System; +using System.Text; + +namespace Org.BouncyCastle.Math.EC +{ + internal class IntArray + { + // TODO make m fixed for the IntArray, and hence compute T once and for all + + // TODO Use uint's internally? + private int[] m_ints; + + public IntArray(int intLen) + { + m_ints = new int[intLen]; + } + + private IntArray(int[] ints) + { + m_ints = ints; + } + + public IntArray(BigInteger bigInt) + : this(bigInt, 0) + { + } + + public IntArray(BigInteger bigInt, int minIntLen) + { + if (bigInt.SignValue == -1) + throw new ArgumentException("Only positive Integers allowed", "bigint"); + + if (bigInt.SignValue == 0) + { + m_ints = new int[] { 0 }; + return; + } + + byte[] barr = bigInt.ToByteArrayUnsigned(); + int barrLen = barr.Length; + + int intLen = (barrLen + 3) / 4; + m_ints = new int[System.Math.Max(intLen, minIntLen)]; + + int rem = barrLen % 4; + int barrI = 0; + + if (0 < rem) + { + int temp = (int) barr[barrI++]; + while (barrI < rem) + { + temp = temp << 8 | (int) barr[barrI++]; + } + m_ints[--intLen] = temp; + } + + while (intLen > 0) + { + int temp = (int) barr[barrI++]; + for (int i = 1; i < 4; i++) + { + temp = temp << 8 | (int) barr[barrI++]; + } + m_ints[--intLen] = temp; + } + } + + public int GetUsedLength() + { + int highestIntPos = m_ints.Length; + + if (highestIntPos < 1) + return 0; + + // Check if first element will act as sentinel + if (m_ints[0] != 0) + { + while (m_ints[--highestIntPos] == 0) + { + } + return highestIntPos + 1; + } + + do + { + if (m_ints[--highestIntPos] != 0) + { + return highestIntPos + 1; + } + } + while (highestIntPos > 0); + + return 0; + } + + public int BitLength + { + get + { + // JDK 1.5: see Integer.numberOfLeadingZeros() + int intLen = GetUsedLength(); + if (intLen == 0) + return 0; + + int last = intLen - 1; + uint highest = (uint) m_ints[last]; + int bits = (last << 5) + 1; + + // A couple of binary search steps + if (highest > 0x0000ffff) + { + if (highest > 0x00ffffff) + { + bits += 24; + highest >>= 24; + } + else + { + bits += 16; + highest >>= 16; + } + } + else if (highest > 0x000000ff) + { + bits += 8; + highest >>= 8; + } + + while (highest > 1) + { + ++bits; + highest >>= 1; + } + + return bits; + } + } + + private int[] resizedInts(int newLen) + { + int[] newInts = new int[newLen]; + int oldLen = m_ints.Length; + int copyLen = oldLen < newLen ? oldLen : newLen; + Array.Copy(m_ints, 0, newInts, 0, copyLen); + return newInts; + } + + public BigInteger ToBigInteger() + { + int usedLen = GetUsedLength(); + if (usedLen == 0) + { + return BigInteger.Zero; + } + + int highestInt = m_ints[usedLen - 1]; + byte[] temp = new byte[4]; + int barrI = 0; + bool trailingZeroBytesDone = false; + for (int j = 3; j >= 0; j--) + { + byte thisByte = (byte)((int)((uint) highestInt >> (8 * j))); + if (trailingZeroBytesDone || (thisByte != 0)) + { + trailingZeroBytesDone = true; + temp[barrI++] = thisByte; + } + } + + int barrLen = 4 * (usedLen - 1) + barrI; + byte[] barr = new byte[barrLen]; + for (int j = 0; j < barrI; j++) + { + barr[j] = temp[j]; + } + // Highest value int is done now + + for (int iarrJ = usedLen - 2; iarrJ >= 0; iarrJ--) + { + for (int j = 3; j >= 0; j--) + { + barr[barrI++] = (byte)((int)((uint)m_ints[iarrJ] >> (8 * j))); + } + } + return new BigInteger(1, barr); + } + + public void ShiftLeft() + { + int usedLen = GetUsedLength(); + if (usedLen == 0) + { + return; + } + if (m_ints[usedLen - 1] < 0) + { + // highest bit of highest used byte is set, so shifting left will + // make the IntArray one byte longer + usedLen++; + if (usedLen > m_ints.Length) + { + // make the m_ints one byte longer, because we need one more + // byte which is not available in m_ints + m_ints = resizedInts(m_ints.Length + 1); + } + } + + bool carry = false; + for (int i = 0; i < usedLen; i++) + { + // nextCarry is true if highest bit is set + bool nextCarry = m_ints[i] < 0; + m_ints[i] <<= 1; + if (carry) + { + // set lowest bit + m_ints[i] |= 1; + } + carry = nextCarry; + } + } + + public IntArray ShiftLeft(int n) + { + int usedLen = GetUsedLength(); + if (usedLen == 0) + { + return this; + } + + if (n == 0) + { + return this; + } + + if (n > 31) + { + throw new ArgumentException("shiftLeft() for max 31 bits " + + ", " + n + "bit shift is not possible", "n"); + } + + int[] newInts = new int[usedLen + 1]; + + int nm32 = 32 - n; + newInts[0] = m_ints[0] << n; + for (int i = 1; i < usedLen; i++) + { + newInts[i] = (m_ints[i] << n) | (int)((uint)m_ints[i - 1] >> nm32); + } + newInts[usedLen] = (int)((uint)m_ints[usedLen - 1] >> nm32); + + return new IntArray(newInts); + } + + public void AddShifted(IntArray other, int shift) + { + int usedLenOther = other.GetUsedLength(); + int newMinUsedLen = usedLenOther + shift; + if (newMinUsedLen > m_ints.Length) + { + m_ints = resizedInts(newMinUsedLen); + //Console.WriteLine("Resize required"); + } + + for (int i = 0; i < usedLenOther; i++) + { + m_ints[i + shift] ^= other.m_ints[i]; + } + } + + public int Length + { + get { return m_ints.Length; } + } + + public bool TestBit(int n) + { + // theInt = n / 32 + int theInt = n >> 5; + // theBit = n % 32 + int theBit = n & 0x1F; + int tester = 1 << theBit; + return ((m_ints[theInt] & tester) != 0); + } + + public void FlipBit(int n) + { + // theInt = n / 32 + int theInt = n >> 5; + // theBit = n % 32 + int theBit = n & 0x1F; + int flipper = 1 << theBit; + m_ints[theInt] ^= flipper; + } + + public void SetBit(int n) + { + // theInt = n / 32 + int theInt = n >> 5; + // theBit = n % 32 + int theBit = n & 0x1F; + int setter = 1 << theBit; + m_ints[theInt] |= setter; + } + + public IntArray Multiply(IntArray other, int m) + { + // Lenght of c is 2m bits rounded up to the next int (32 bit) + int t = (m + 31) >> 5; + if (m_ints.Length < t) + { + m_ints = resizedInts(t); + } + + IntArray b = new IntArray(other.resizedInts(other.Length + 1)); + IntArray c = new IntArray((m + m + 31) >> 5); + // IntArray c = new IntArray(t + t); + int testBit = 1; + for (int k = 0; k < 32; k++) + { + for (int j = 0; j < t; j++) + { + if ((m_ints[j] & testBit) != 0) + { + // The kth bit of m_ints[j] is set + c.AddShifted(b, j); + } + } + testBit <<= 1; + b.ShiftLeft(); + } + return c; + } + + // public IntArray multiplyLeftToRight(IntArray other, int m) { + // // Lenght of c is 2m bits rounded up to the next int (32 bit) + // int t = (m + 31) / 32; + // if (m_ints.Length < t) { + // m_ints = resizedInts(t); + // } + // + // IntArray b = new IntArray(other.resizedInts(other.getLength() + 1)); + // IntArray c = new IntArray((m + m + 31) / 32); + // // IntArray c = new IntArray(t + t); + // int testBit = 1 << 31; + // for (int k = 31; k >= 0; k--) { + // for (int j = 0; j < t; j++) { + // if ((m_ints[j] & testBit) != 0) { + // // The kth bit of m_ints[j] is set + // c.addShifted(b, j); + // } + // } + // testBit >>>= 1; + // if (k > 0) { + // c.shiftLeft(); + // } + // } + // return c; + // } + + // TODO note, redPol.Length must be 3 for TPB and 5 for PPB + public void Reduce(int m, int[] redPol) + { + for (int i = m + m - 2; i >= m; i--) + { + if (TestBit(i)) + { + int bit = i - m; + FlipBit(bit); + FlipBit(i); + int l = redPol.Length; + while (--l >= 0) + { + FlipBit(redPol[l] + bit); + } + } + } + m_ints = resizedInts((m + 31) >> 5); + } + + public IntArray Square(int m) + { + // TODO make the table static readonly + int[] table = { 0x0, 0x1, 0x4, 0x5, 0x10, 0x11, 0x14, 0x15, 0x40, + 0x41, 0x44, 0x45, 0x50, 0x51, 0x54, 0x55 }; + + int t = (m + 31) >> 5; + if (m_ints.Length < t) + { + m_ints = resizedInts(t); + } + + IntArray c = new IntArray(t + t); + + // TODO twice the same code, put in separate private method + for (int i = 0; i < t; i++) + { + int v0 = 0; + for (int j = 0; j < 4; j++) + { + v0 = (int)((uint) v0 >> 8); + int u = (int)((uint)m_ints[i] >> (j * 4)) & 0xF; + int w = table[u] << 24; + v0 |= w; + } + c.m_ints[i + i] = v0; + + v0 = 0; + int upper = (int)((uint) m_ints[i] >> 16); + for (int j = 0; j < 4; j++) + { + v0 = (int)((uint) v0 >> 8); + int u = (int)((uint)upper >> (j * 4)) & 0xF; + int w = table[u] << 24; + v0 |= w; + } + c.m_ints[i + i + 1] = v0; + } + return c; + } + + public override bool Equals(object o) + { + if (!(o is IntArray)) + { + return false; + } + IntArray other = (IntArray) o; + int usedLen = GetUsedLength(); + if (other.GetUsedLength() != usedLen) + { + return false; + } + for (int i = 0; i < usedLen; i++) + { + if (m_ints[i] != other.m_ints[i]) + { + return false; + } + } + return true; + } + + public override int GetHashCode() + { + int i = GetUsedLength(); + int hc = i; + while (--i >= 0) + { + hc *= 17; + hc ^= m_ints[i]; + } + return hc; + } + + internal IntArray Copy() + { + return new IntArray((int[]) m_ints.Clone()); + } + + public override string ToString() + { + int usedLen = GetUsedLength(); + if (usedLen == 0) + { + return "0"; + } + + StringBuilder sb = new StringBuilder(Convert.ToString(m_ints[usedLen - 1], 2)); + for (int iarrJ = usedLen - 2; iarrJ >= 0; iarrJ--) + { + string hexString = Convert.ToString(m_ints[iarrJ], 2); + + // Add leading zeroes, except for highest significant int + for (int i = hexString.Length; i < 8; i++) + { + hexString = "0" + hexString; + } + sb.Append(hexString); + } + return sb.ToString(); + } + } +} diff --git a/crypto/src/math/ec/abc/SimpleBigDecimal.cs b/crypto/src/math/ec/abc/SimpleBigDecimal.cs new file mode 100644 index 000000000..d5664dbfd --- /dev/null +++ b/crypto/src/math/ec/abc/SimpleBigDecimal.cs @@ -0,0 +1,241 @@ +using System; +using System.Text; + +namespace Org.BouncyCastle.Math.EC.Abc +{ + /** + * Class representing a simple version of a big decimal. A + * <code>SimpleBigDecimal</code> is basically a + * {@link java.math.BigInteger BigInteger} with a few digits on the right of + * the decimal point. The number of (binary) digits on the right of the decimal + * point is called the <code>scale</code> of the <code>SimpleBigDecimal</code>. + * Unlike in {@link java.math.BigDecimal BigDecimal}, the scale is not adjusted + * automatically, but must be set manually. All <code>SimpleBigDecimal</code>s + * taking part in the same arithmetic operation must have equal scale. The + * result of a multiplication of two <code>SimpleBigDecimal</code>s returns a + * <code>SimpleBigDecimal</code> with double scale. + */ + internal class SimpleBigDecimal + // : Number + { + // private static final long serialVersionUID = 1L; + + private readonly BigInteger bigInt; + private readonly int scale; + + /** + * Returns a <code>SimpleBigDecimal</code> representing the same numerical + * value as <code>value</code>. + * @param value The value of the <code>SimpleBigDecimal</code> to be + * created. + * @param scale The scale of the <code>SimpleBigDecimal</code> to be + * created. + * @return The such created <code>SimpleBigDecimal</code>. + */ + public static SimpleBigDecimal GetInstance(BigInteger val, int scale) + { + return new SimpleBigDecimal(val.ShiftLeft(scale), scale); + } + + /** + * Constructor for <code>SimpleBigDecimal</code>. The value of the + * constructed <code>SimpleBigDecimal</code> Equals <code>bigInt / + * 2<sup>scale</sup></code>. + * @param bigInt The <code>bigInt</code> value parameter. + * @param scale The scale of the constructed <code>SimpleBigDecimal</code>. + */ + public SimpleBigDecimal(BigInteger bigInt, int scale) + { + if (scale < 0) + throw new ArgumentException("scale may not be negative"); + + this.bigInt = bigInt; + this.scale = scale; + } + + private SimpleBigDecimal(SimpleBigDecimal limBigDec) + { + bigInt = limBigDec.bigInt; + scale = limBigDec.scale; + } + + private void CheckScale(SimpleBigDecimal b) + { + if (scale != b.scale) + throw new ArgumentException("Only SimpleBigDecimal of same scale allowed in arithmetic operations"); + } + + public SimpleBigDecimal AdjustScale(int newScale) + { + if (newScale < 0) + throw new ArgumentException("scale may not be negative"); + + if (newScale == scale) + return this; + + return new SimpleBigDecimal(bigInt.ShiftLeft(newScale - scale), newScale); + } + + public SimpleBigDecimal Add(SimpleBigDecimal b) + { + CheckScale(b); + return new SimpleBigDecimal(bigInt.Add(b.bigInt), scale); + } + + public SimpleBigDecimal Add(BigInteger b) + { + return new SimpleBigDecimal(bigInt.Add(b.ShiftLeft(scale)), scale); + } + + public SimpleBigDecimal Negate() + { + return new SimpleBigDecimal(bigInt.Negate(), scale); + } + + public SimpleBigDecimal Subtract(SimpleBigDecimal b) + { + return Add(b.Negate()); + } + + public SimpleBigDecimal Subtract(BigInteger b) + { + return new SimpleBigDecimal(bigInt.Subtract(b.ShiftLeft(scale)), scale); + } + + public SimpleBigDecimal Multiply(SimpleBigDecimal b) + { + CheckScale(b); + return new SimpleBigDecimal(bigInt.Multiply(b.bigInt), scale + scale); + } + + public SimpleBigDecimal Multiply(BigInteger b) + { + return new SimpleBigDecimal(bigInt.Multiply(b), scale); + } + + public SimpleBigDecimal Divide(SimpleBigDecimal b) + { + CheckScale(b); + BigInteger dividend = bigInt.ShiftLeft(scale); + return new SimpleBigDecimal(dividend.Divide(b.bigInt), scale); + } + + public SimpleBigDecimal Divide(BigInteger b) + { + return new SimpleBigDecimal(bigInt.Divide(b), scale); + } + + public SimpleBigDecimal ShiftLeft(int n) + { + return new SimpleBigDecimal(bigInt.ShiftLeft(n), scale); + } + + public int CompareTo(SimpleBigDecimal val) + { + CheckScale(val); + return bigInt.CompareTo(val.bigInt); + } + + public int CompareTo(BigInteger val) + { + return bigInt.CompareTo(val.ShiftLeft(scale)); + } + + public BigInteger Floor() + { + return bigInt.ShiftRight(scale); + } + + public BigInteger Round() + { + SimpleBigDecimal oneHalf = new SimpleBigDecimal(BigInteger.One, 1); + return Add(oneHalf.AdjustScale(scale)).Floor(); + } + + public int IntValue + { + get { return Floor().IntValue; } + } + + public long LongValue + { + get { return Floor().LongValue; } + } + +// public double doubleValue() +// { +// return new Double(ToString()).doubleValue(); +// } +// +// public float floatValue() +// { +// return new Float(ToString()).floatValue(); +// } + + public int Scale + { + get { return scale; } + } + + public override string ToString() + { + if (scale == 0) + return bigInt.ToString(); + + BigInteger floorBigInt = Floor(); + + BigInteger fract = bigInt.Subtract(floorBigInt.ShiftLeft(scale)); + if (bigInt.SignValue < 0) + { + fract = BigInteger.One.ShiftLeft(scale).Subtract(fract); + } + + if ((floorBigInt.SignValue == -1) && (!(fract.Equals(BigInteger.Zero)))) + { + floorBigInt = floorBigInt.Add(BigInteger.One); + } + string leftOfPoint = floorBigInt.ToString(); + + char[] fractCharArr = new char[scale]; + string fractStr = fract.ToString(2); + int fractLen = fractStr.Length; + int zeroes = scale - fractLen; + for (int i = 0; i < zeroes; i++) + { + fractCharArr[i] = '0'; + } + for (int j = 0; j < fractLen; j++) + { + fractCharArr[zeroes + j] = fractStr[j]; + } + string rightOfPoint = new string(fractCharArr); + + StringBuilder sb = new StringBuilder(leftOfPoint); + sb.Append("."); + sb.Append(rightOfPoint); + + return sb.ToString(); + } + + public override bool Equals( + object obj) + { + if (this == obj) + return true; + + SimpleBigDecimal other = obj as SimpleBigDecimal; + + if (other == null) + return false; + + return bigInt.Equals(other.bigInt) + && scale == other.scale; + } + + public override int GetHashCode() + { + return bigInt.GetHashCode() ^ scale; + } + + } +} diff --git a/crypto/src/math/ec/abc/Tnaf.cs b/crypto/src/math/ec/abc/Tnaf.cs new file mode 100644 index 000000000..225fc3075 --- /dev/null +++ b/crypto/src/math/ec/abc/Tnaf.cs @@ -0,0 +1,834 @@ +using System; + +namespace Org.BouncyCastle.Math.EC.Abc +{ + /** + * Class holding methods for point multiplication based on the window + * τ-adic nonadjacent form (WTNAF). The algorithms are based on the + * paper "Improved Algorithms for Arithmetic on Anomalous Binary Curves" + * by Jerome A. Solinas. The paper first appeared in the Proceedings of + * Crypto 1997. + */ + internal class Tnaf + { + private static readonly BigInteger MinusOne = BigInteger.One.Negate(); + private static readonly BigInteger MinusTwo = BigInteger.Two.Negate(); + private static readonly BigInteger MinusThree = BigInteger.Three.Negate(); + private static readonly BigInteger Four = BigInteger.ValueOf(4); + + /** + * The window width of WTNAF. The standard value of 4 is slightly less + * than optimal for running time, but keeps space requirements for + * precomputation low. For typical curves, a value of 5 or 6 results in + * a better running time. When changing this value, the + * <code>α<sub>u</sub></code>'s must be computed differently, see + * e.g. "Guide to Elliptic Curve Cryptography", Darrel Hankerson, + * Alfred Menezes, Scott Vanstone, Springer-Verlag New York Inc., 2004, + * p. 121-122 + */ + public const sbyte Width = 4; + + /** + * 2<sup>4</sup> + */ + public const sbyte Pow2Width = 16; + + /** + * The <code>α<sub>u</sub></code>'s for <code>a=0</code> as an array + * of <code>ZTauElement</code>s. + */ + public static readonly ZTauElement[] Alpha0 = + { + null, + new ZTauElement(BigInteger.One, BigInteger.Zero), null, + new ZTauElement(MinusThree, MinusOne), null, + new ZTauElement(MinusOne, MinusOne), null, + new ZTauElement(BigInteger.One, MinusOne), null + }; + + /** + * The <code>α<sub>u</sub></code>'s for <code>a=0</code> as an array + * of TNAFs. + */ + public static readonly sbyte[][] Alpha0Tnaf = + { + null, new sbyte[]{1}, null, new sbyte[]{-1, 0, 1}, null, new sbyte[]{1, 0, 1}, null, new sbyte[]{-1, 0, 0, 1} + }; + + /** + * The <code>α<sub>u</sub></code>'s for <code>a=1</code> as an array + * of <code>ZTauElement</code>s. + */ + public static readonly ZTauElement[] Alpha1 = + { + null, + new ZTauElement(BigInteger.One, BigInteger.Zero), null, + new ZTauElement(MinusThree, BigInteger.One), null, + new ZTauElement(MinusOne, BigInteger.One), null, + new ZTauElement(BigInteger.One, BigInteger.One), null + }; + + /** + * The <code>α<sub>u</sub></code>'s for <code>a=1</code> as an array + * of TNAFs. + */ + public static readonly sbyte[][] Alpha1Tnaf = + { + null, new sbyte[]{1}, null, new sbyte[]{-1, 0, 1}, null, new sbyte[]{1, 0, 1}, null, new sbyte[]{-1, 0, 0, -1} + }; + + /** + * Computes the norm of an element <code>λ</code> of + * <code><b>Z</b>[τ]</code>. + * @param mu The parameter <code>μ</code> of the elliptic curve. + * @param lambda The element <code>λ</code> of + * <code><b>Z</b>[τ]</code>. + * @return The norm of <code>λ</code>. + */ + public static BigInteger Norm(sbyte mu, ZTauElement lambda) + { + BigInteger norm; + + // s1 = u^2 + BigInteger s1 = lambda.u.Multiply(lambda.u); + + // s2 = u * v + BigInteger s2 = lambda.u.Multiply(lambda.v); + + // s3 = 2 * v^2 + BigInteger s3 = lambda.v.Multiply(lambda.v).ShiftLeft(1); + + if (mu == 1) + { + norm = s1.Add(s2).Add(s3); + } + else if (mu == -1) + { + norm = s1.Subtract(s2).Add(s3); + } + else + { + throw new ArgumentException("mu must be 1 or -1"); + } + + return norm; + } + + /** + * Computes the norm of an element <code>λ</code> of + * <code><b>R</b>[τ]</code>, where <code>λ = u + vτ</code> + * and <code>u</code> and <code>u</code> are real numbers (elements of + * <code><b>R</b></code>). + * @param mu The parameter <code>μ</code> of the elliptic curve. + * @param u The real part of the element <code>λ</code> of + * <code><b>R</b>[τ]</code>. + * @param v The <code>τ</code>-adic part of the element + * <code>λ</code> of <code><b>R</b>[τ]</code>. + * @return The norm of <code>λ</code>. + */ + public static SimpleBigDecimal Norm(sbyte mu, SimpleBigDecimal u, SimpleBigDecimal v) + { + SimpleBigDecimal norm; + + // s1 = u^2 + SimpleBigDecimal s1 = u.Multiply(u); + + // s2 = u * v + SimpleBigDecimal s2 = u.Multiply(v); + + // s3 = 2 * v^2 + SimpleBigDecimal s3 = v.Multiply(v).ShiftLeft(1); + + if (mu == 1) + { + norm = s1.Add(s2).Add(s3); + } + else if (mu == -1) + { + norm = s1.Subtract(s2).Add(s3); + } + else + { + throw new ArgumentException("mu must be 1 or -1"); + } + + return norm; + } + + /** + * Rounds an element <code>λ</code> of <code><b>R</b>[τ]</code> + * to an element of <code><b>Z</b>[τ]</code>, such that their difference + * has minimal norm. <code>λ</code> is given as + * <code>λ = λ<sub>0</sub> + λ<sub>1</sub>τ</code>. + * @param lambda0 The component <code>λ<sub>0</sub></code>. + * @param lambda1 The component <code>λ<sub>1</sub></code>. + * @param mu The parameter <code>μ</code> of the elliptic curve. Must + * equal 1 or -1. + * @return The rounded element of <code><b>Z</b>[τ]</code>. + * @throws ArgumentException if <code>lambda0</code> and + * <code>lambda1</code> do not have same scale. + */ + public static ZTauElement Round(SimpleBigDecimal lambda0, + SimpleBigDecimal lambda1, sbyte mu) + { + int scale = lambda0.Scale; + if (lambda1.Scale != scale) + throw new ArgumentException("lambda0 and lambda1 do not have same scale"); + + if (!((mu == 1) || (mu == -1))) + throw new ArgumentException("mu must be 1 or -1"); + + BigInteger f0 = lambda0.Round(); + BigInteger f1 = lambda1.Round(); + + SimpleBigDecimal eta0 = lambda0.Subtract(f0); + SimpleBigDecimal eta1 = lambda1.Subtract(f1); + + // eta = 2*eta0 + mu*eta1 + SimpleBigDecimal eta = eta0.Add(eta0); + if (mu == 1) + { + eta = eta.Add(eta1); + } + else + { + // mu == -1 + eta = eta.Subtract(eta1); + } + + // check1 = eta0 - 3*mu*eta1 + // check2 = eta0 + 4*mu*eta1 + SimpleBigDecimal threeEta1 = eta1.Add(eta1).Add(eta1); + SimpleBigDecimal fourEta1 = threeEta1.Add(eta1); + SimpleBigDecimal check1; + SimpleBigDecimal check2; + if (mu == 1) + { + check1 = eta0.Subtract(threeEta1); + check2 = eta0.Add(fourEta1); + } + else + { + // mu == -1 + check1 = eta0.Add(threeEta1); + check2 = eta0.Subtract(fourEta1); + } + + sbyte h0 = 0; + sbyte h1 = 0; + + // if eta >= 1 + if (eta.CompareTo(BigInteger.One) >= 0) + { + if (check1.CompareTo(MinusOne) < 0) + { + h1 = mu; + } + else + { + h0 = 1; + } + } + else + { + // eta < 1 + if (check2.CompareTo(BigInteger.Two) >= 0) + { + h1 = mu; + } + } + + // if eta < -1 + if (eta.CompareTo(MinusOne) < 0) + { + if (check1.CompareTo(BigInteger.One) >= 0) + { + h1 = (sbyte)-mu; + } + else + { + h0 = -1; + } + } + else + { + // eta >= -1 + if (check2.CompareTo(MinusTwo) < 0) + { + h1 = (sbyte)-mu; + } + } + + BigInteger q0 = f0.Add(BigInteger.ValueOf(h0)); + BigInteger q1 = f1.Add(BigInteger.ValueOf(h1)); + return new ZTauElement(q0, q1); + } + + /** + * Approximate division by <code>n</code>. For an integer + * <code>k</code>, the value <code>λ = s k / n</code> is + * computed to <code>c</code> bits of accuracy. + * @param k The parameter <code>k</code>. + * @param s The curve parameter <code>s<sub>0</sub></code> or + * <code>s<sub>1</sub></code>. + * @param vm The Lucas Sequence element <code>V<sub>m</sub></code>. + * @param a The parameter <code>a</code> of the elliptic curve. + * @param m The bit length of the finite field + * <code><b>F</b><sub>m</sub></code>. + * @param c The number of bits of accuracy, i.e. the scale of the returned + * <code>SimpleBigDecimal</code>. + * @return The value <code>λ = s k / n</code> computed to + * <code>c</code> bits of accuracy. + */ + public static SimpleBigDecimal ApproximateDivisionByN(BigInteger k, + BigInteger s, BigInteger vm, sbyte a, int m, int c) + { + int _k = (m + 5)/2 + c; + BigInteger ns = k.ShiftRight(m - _k - 2 + a); + + BigInteger gs = s.Multiply(ns); + + BigInteger hs = gs.ShiftRight(m); + + BigInteger js = vm.Multiply(hs); + + BigInteger gsPlusJs = gs.Add(js); + BigInteger ls = gsPlusJs.ShiftRight(_k-c); + if (gsPlusJs.TestBit(_k-c-1)) + { + // round up + ls = ls.Add(BigInteger.One); + } + + return new SimpleBigDecimal(ls, c); + } + + /** + * Computes the <code>τ</code>-adic NAF (non-adjacent form) of an + * element <code>λ</code> of <code><b>Z</b>[τ]</code>. + * @param mu The parameter <code>μ</code> of the elliptic curve. + * @param lambda The element <code>λ</code> of + * <code><b>Z</b>[τ]</code>. + * @return The <code>τ</code>-adic NAF of <code>λ</code>. + */ + public static sbyte[] TauAdicNaf(sbyte mu, ZTauElement lambda) + { + if (!((mu == 1) || (mu == -1))) + throw new ArgumentException("mu must be 1 or -1"); + + BigInteger norm = Norm(mu, lambda); + + // Ceiling of log2 of the norm + int log2Norm = norm.BitLength; + + // If length(TNAF) > 30, then length(TNAF) < log2Norm + 3.52 + int maxLength = log2Norm > 30 ? log2Norm + 4 : 34; + + // The array holding the TNAF + sbyte[] u = new sbyte[maxLength]; + int i = 0; + + // The actual length of the TNAF + int length = 0; + + BigInteger r0 = lambda.u; + BigInteger r1 = lambda.v; + + while(!((r0.Equals(BigInteger.Zero)) && (r1.Equals(BigInteger.Zero)))) + { + // If r0 is odd + if (r0.TestBit(0)) + { + u[i] = (sbyte) BigInteger.Two.Subtract((r0.Subtract(r1.ShiftLeft(1))).Mod(Four)).IntValue; + + // r0 = r0 - u[i] + if (u[i] == 1) + { + r0 = r0.ClearBit(0); + } + else + { + // u[i] == -1 + r0 = r0.Add(BigInteger.One); + } + length = i; + } + else + { + u[i] = 0; + } + + BigInteger t = r0; + BigInteger s = r0.ShiftRight(1); + if (mu == 1) + { + r0 = r1.Add(s); + } + else + { + // mu == -1 + r0 = r1.Subtract(s); + } + + r1 = t.ShiftRight(1).Negate(); + i++; + } + + length++; + + // Reduce the TNAF array to its actual length + sbyte[] tnaf = new sbyte[length]; + Array.Copy(u, 0, tnaf, 0, length); + return tnaf; + } + + /** + * Applies the operation <code>τ()</code> to an + * <code>F2mPoint</code>. + * @param p The F2mPoint to which <code>τ()</code> is applied. + * @return <code>τ(p)</code> + */ + public static F2mPoint Tau(F2mPoint p) + { + if (p.IsInfinity) + return p; + + ECFieldElement x = p.X; + ECFieldElement y = p.Y; + + return new F2mPoint(p.Curve, x.Square(), y.Square(), p.IsCompressed); + } + + /** + * Returns the parameter <code>μ</code> of the elliptic curve. + * @param curve The elliptic curve from which to obtain <code>μ</code>. + * The curve must be a Koblitz curve, i.e. <code>a</code> Equals + * <code>0</code> or <code>1</code> and <code>b</code> Equals + * <code>1</code>. + * @return <code>μ</code> of the elliptic curve. + * @throws ArgumentException if the given ECCurve is not a Koblitz + * curve. + */ + public static sbyte GetMu(F2mCurve curve) + { + BigInteger a = curve.A.ToBigInteger(); + + sbyte mu; + if (a.SignValue == 0) + { + mu = -1; + } + else if (a.Equals(BigInteger.One)) + { + mu = 1; + } + else + { + throw new ArgumentException("No Koblitz curve (ABC), TNAF multiplication not possible"); + } + return mu; + } + + /** + * Calculates the Lucas Sequence elements <code>U<sub>k-1</sub></code> and + * <code>U<sub>k</sub></code> or <code>V<sub>k-1</sub></code> and + * <code>V<sub>k</sub></code>. + * @param mu The parameter <code>μ</code> of the elliptic curve. + * @param k The index of the second element of the Lucas Sequence to be + * returned. + * @param doV If set to true, computes <code>V<sub>k-1</sub></code> and + * <code>V<sub>k</sub></code>, otherwise <code>U<sub>k-1</sub></code> and + * <code>U<sub>k</sub></code>. + * @return An array with 2 elements, containing <code>U<sub>k-1</sub></code> + * and <code>U<sub>k</sub></code> or <code>V<sub>k-1</sub></code> + * and <code>V<sub>k</sub></code>. + */ + public static BigInteger[] GetLucas(sbyte mu, int k, bool doV) + { + if (!(mu == 1 || mu == -1)) + throw new ArgumentException("mu must be 1 or -1"); + + BigInteger u0; + BigInteger u1; + BigInteger u2; + + if (doV) + { + u0 = BigInteger.Two; + u1 = BigInteger.ValueOf(mu); + } + else + { + u0 = BigInteger.Zero; + u1 = BigInteger.One; + } + + for (int i = 1; i < k; i++) + { + // u2 = mu*u1 - 2*u0; + BigInteger s = null; + if (mu == 1) + { + s = u1; + } + else + { + // mu == -1 + s = u1.Negate(); + } + + u2 = s.Subtract(u0.ShiftLeft(1)); + u0 = u1; + u1 = u2; + // System.out.println(i + ": " + u2); + // System.out.println(); + } + + BigInteger[] retVal = {u0, u1}; + return retVal; + } + + /** + * Computes the auxiliary value <code>t<sub>w</sub></code>. If the width is + * 4, then for <code>mu = 1</code>, <code>t<sub>w</sub> = 6</code> and for + * <code>mu = -1</code>, <code>t<sub>w</sub> = 10</code> + * @param mu The parameter <code>μ</code> of the elliptic curve. + * @param w The window width of the WTNAF. + * @return the auxiliary value <code>t<sub>w</sub></code> + */ + public static BigInteger GetTw(sbyte mu, int w) + { + if (w == 4) + { + if (mu == 1) + { + return BigInteger.ValueOf(6); + } + else + { + // mu == -1 + return BigInteger.ValueOf(10); + } + } + else + { + // For w <> 4, the values must be computed + BigInteger[] us = GetLucas(mu, w, false); + BigInteger twoToW = BigInteger.Zero.SetBit(w); + BigInteger u1invert = us[1].ModInverse(twoToW); + BigInteger tw; + tw = BigInteger.Two.Multiply(us[0]).Multiply(u1invert).Mod(twoToW); + //System.out.println("mu = " + mu); + //System.out.println("tw = " + tw); + return tw; + } + } + + /** + * Computes the auxiliary values <code>s<sub>0</sub></code> and + * <code>s<sub>1</sub></code> used for partial modular reduction. + * @param curve The elliptic curve for which to compute + * <code>s<sub>0</sub></code> and <code>s<sub>1</sub></code>. + * @throws ArgumentException if <code>curve</code> is not a + * Koblitz curve (Anomalous Binary Curve, ABC). + */ + public static BigInteger[] GetSi(F2mCurve curve) + { + if (!curve.IsKoblitz) + throw new ArgumentException("si is defined for Koblitz curves only"); + + int m = curve.M; + int a = curve.A.ToBigInteger().IntValue; + sbyte mu = curve.GetMu(); + int h = curve.H.IntValue; + int index = m + 3 - a; + BigInteger[] ui = GetLucas(mu, index, false); + + BigInteger dividend0; + BigInteger dividend1; + if (mu == 1) + { + dividend0 = BigInteger.One.Subtract(ui[1]); + dividend1 = BigInteger.One.Subtract(ui[0]); + } + else if (mu == -1) + { + dividend0 = BigInteger.One.Add(ui[1]); + dividend1 = BigInteger.One.Add(ui[0]); + } + else + { + throw new ArgumentException("mu must be 1 or -1"); + } + + BigInteger[] si = new BigInteger[2]; + + if (h == 2) + { + si[0] = dividend0.ShiftRight(1); + si[1] = dividend1.ShiftRight(1).Negate(); + } + else if (h == 4) + { + si[0] = dividend0.ShiftRight(2); + si[1] = dividend1.ShiftRight(2).Negate(); + } + else + { + throw new ArgumentException("h (Cofactor) must be 2 or 4"); + } + + return si; + } + + /** + * Partial modular reduction modulo + * <code>(τ<sup>m</sup> - 1)/(τ - 1)</code>. + * @param k The integer to be reduced. + * @param m The bitlength of the underlying finite field. + * @param a The parameter <code>a</code> of the elliptic curve. + * @param s The auxiliary values <code>s<sub>0</sub></code> and + * <code>s<sub>1</sub></code>. + * @param mu The parameter μ of the elliptic curve. + * @param c The precision (number of bits of accuracy) of the partial + * modular reduction. + * @return <code>ρ := k partmod (τ<sup>m</sup> - 1)/(τ - 1)</code> + */ + public static ZTauElement PartModReduction(BigInteger k, int m, sbyte a, + BigInteger[] s, sbyte mu, sbyte c) + { + // d0 = s[0] + mu*s[1]; mu is either 1 or -1 + BigInteger d0; + if (mu == 1) + { + d0 = s[0].Add(s[1]); + } + else + { + d0 = s[0].Subtract(s[1]); + } + + BigInteger[] v = GetLucas(mu, m, true); + BigInteger vm = v[1]; + + SimpleBigDecimal lambda0 = ApproximateDivisionByN( + k, s[0], vm, a, m, c); + + SimpleBigDecimal lambda1 = ApproximateDivisionByN( + k, s[1], vm, a, m, c); + + ZTauElement q = Round(lambda0, lambda1, mu); + + // r0 = n - d0*q0 - 2*s1*q1 + BigInteger r0 = k.Subtract(d0.Multiply(q.u)).Subtract( + BigInteger.ValueOf(2).Multiply(s[1]).Multiply(q.v)); + + // r1 = s1*q0 - s0*q1 + BigInteger r1 = s[1].Multiply(q.u).Subtract(s[0].Multiply(q.v)); + + return new ZTauElement(r0, r1); + } + + /** + * Multiplies a {@link org.bouncycastle.math.ec.F2mPoint F2mPoint} + * by a <code>BigInteger</code> using the reduced <code>τ</code>-adic + * NAF (RTNAF) method. + * @param p The F2mPoint to Multiply. + * @param k The <code>BigInteger</code> by which to Multiply <code>p</code>. + * @return <code>k * p</code> + */ + public static F2mPoint MultiplyRTnaf(F2mPoint p, BigInteger k) + { + F2mCurve curve = (F2mCurve) p.Curve; + int m = curve.M; + sbyte a = (sbyte) curve.A.ToBigInteger().IntValue; + sbyte mu = curve.GetMu(); + BigInteger[] s = curve.GetSi(); + ZTauElement rho = PartModReduction(k, m, a, s, mu, (sbyte)10); + + return MultiplyTnaf(p, rho); + } + + /** + * Multiplies a {@link org.bouncycastle.math.ec.F2mPoint F2mPoint} + * by an element <code>λ</code> of <code><b>Z</b>[τ]</code> + * using the <code>τ</code>-adic NAF (TNAF) method. + * @param p The F2mPoint to Multiply. + * @param lambda The element <code>λ</code> of + * <code><b>Z</b>[τ]</code>. + * @return <code>λ * p</code> + */ + public static F2mPoint MultiplyTnaf(F2mPoint p, ZTauElement lambda) + { + F2mCurve curve = (F2mCurve)p.Curve; + sbyte mu = curve.GetMu(); + sbyte[] u = TauAdicNaf(mu, lambda); + + F2mPoint q = MultiplyFromTnaf(p, u); + + return q; + } + + /** + * Multiplies a {@link org.bouncycastle.math.ec.F2mPoint F2mPoint} + * by an element <code>λ</code> of <code><b>Z</b>[τ]</code> + * using the <code>τ</code>-adic NAF (TNAF) method, given the TNAF + * of <code>λ</code>. + * @param p The F2mPoint to Multiply. + * @param u The the TNAF of <code>λ</code>.. + * @return <code>λ * p</code> + */ + public static F2mPoint MultiplyFromTnaf(F2mPoint p, sbyte[] u) + { + F2mCurve curve = (F2mCurve)p.Curve; + F2mPoint q = (F2mPoint) curve.Infinity; + for (int i = u.Length - 1; i >= 0; i--) + { + q = Tau(q); + if (u[i] == 1) + { + q = (F2mPoint)q.AddSimple(p); + } + else if (u[i] == -1) + { + q = (F2mPoint)q.SubtractSimple(p); + } + } + return q; + } + + /** + * Computes the <code>[τ]</code>-adic window NAF of an element + * <code>λ</code> of <code><b>Z</b>[τ]</code>. + * @param mu The parameter μ of the elliptic curve. + * @param lambda The element <code>λ</code> of + * <code><b>Z</b>[τ]</code> of which to compute the + * <code>[τ]</code>-adic NAF. + * @param width The window width of the resulting WNAF. + * @param pow2w 2<sup>width</sup>. + * @param tw The auxiliary value <code>t<sub>w</sub></code>. + * @param alpha The <code>α<sub>u</sub></code>'s for the window width. + * @return The <code>[τ]</code>-adic window NAF of + * <code>λ</code>. + */ + public static sbyte[] TauAdicWNaf(sbyte mu, ZTauElement lambda, + sbyte width, BigInteger pow2w, BigInteger tw, ZTauElement[] alpha) + { + if (!((mu == 1) || (mu == -1))) + throw new ArgumentException("mu must be 1 or -1"); + + BigInteger norm = Norm(mu, lambda); + + // Ceiling of log2 of the norm + int log2Norm = norm.BitLength; + + // If length(TNAF) > 30, then length(TNAF) < log2Norm + 3.52 + int maxLength = log2Norm > 30 ? log2Norm + 4 + width : 34 + width; + + // The array holding the TNAF + sbyte[] u = new sbyte[maxLength]; + + // 2^(width - 1) + BigInteger pow2wMin1 = pow2w.ShiftRight(1); + + // Split lambda into two BigIntegers to simplify calculations + BigInteger r0 = lambda.u; + BigInteger r1 = lambda.v; + int i = 0; + + // while lambda <> (0, 0) + while (!((r0.Equals(BigInteger.Zero))&&(r1.Equals(BigInteger.Zero)))) + { + // if r0 is odd + if (r0.TestBit(0)) + { + // uUnMod = r0 + r1*tw Mod 2^width + BigInteger uUnMod + = r0.Add(r1.Multiply(tw)).Mod(pow2w); + + sbyte uLocal; + // if uUnMod >= 2^(width - 1) + if (uUnMod.CompareTo(pow2wMin1) >= 0) + { + uLocal = (sbyte) uUnMod.Subtract(pow2w).IntValue; + } + else + { + uLocal = (sbyte) uUnMod.IntValue; + } + // uLocal is now in [-2^(width-1), 2^(width-1)-1] + + u[i] = uLocal; + bool s = true; + if (uLocal < 0) + { + s = false; + uLocal = (sbyte)-uLocal; + } + // uLocal is now >= 0 + + if (s) + { + r0 = r0.Subtract(alpha[uLocal].u); + r1 = r1.Subtract(alpha[uLocal].v); + } + else + { + r0 = r0.Add(alpha[uLocal].u); + r1 = r1.Add(alpha[uLocal].v); + } + } + else + { + u[i] = 0; + } + + BigInteger t = r0; + + if (mu == 1) + { + r0 = r1.Add(r0.ShiftRight(1)); + } + else + { + // mu == -1 + r0 = r1.Subtract(r0.ShiftRight(1)); + } + r1 = t.ShiftRight(1).Negate(); + i++; + } + return u; + } + + /** + * Does the precomputation for WTNAF multiplication. + * @param p The <code>ECPoint</code> for which to do the precomputation. + * @param a The parameter <code>a</code> of the elliptic curve. + * @return The precomputation array for <code>p</code>. + */ + public static F2mPoint[] GetPreComp(F2mPoint p, sbyte a) + { + F2mPoint[] pu; + pu = new F2mPoint[16]; + pu[1] = p; + sbyte[][] alphaTnaf; + if (a == 0) + { + alphaTnaf = Tnaf.Alpha0Tnaf; + } + else + { + // a == 1 + alphaTnaf = Tnaf.Alpha1Tnaf; + } + + int precompLen = alphaTnaf.Length; + for (int i = 3; i < precompLen; i = i + 2) + { + pu[i] = Tnaf.MultiplyFromTnaf(p, alphaTnaf[i]); + } + + return pu; + } + } +} diff --git a/crypto/src/math/ec/abc/ZTauElement.cs b/crypto/src/math/ec/abc/ZTauElement.cs new file mode 100644 index 000000000..4fcbf1bdf --- /dev/null +++ b/crypto/src/math/ec/abc/ZTauElement.cs @@ -0,0 +1,36 @@ +namespace Org.BouncyCastle.Math.EC.Abc +{ + /** + * Class representing an element of <code><b>Z</b>[τ]</code>. Let + * <code>λ</code> be an element of <code><b>Z</b>[τ]</code>. Then + * <code>λ</code> is given as <code>λ = u + vτ</code>. The + * components <code>u</code> and <code>v</code> may be used directly, there + * are no accessor methods. + * Immutable class. + */ + internal class ZTauElement + { + /** + * The "real" part of <code>λ</code>. + */ + public readonly BigInteger u; + + /** + * The "<code>τ</code>-adic" part of <code>λ</code>. + */ + public readonly BigInteger v; + + /** + * Constructor for an element <code>λ</code> of + * <code><b>Z</b>[τ]</code>. + * @param u The "real" part of <code>λ</code>. + * @param v The "<code>τ</code>-adic" part of + * <code>λ</code>. + */ + public ZTauElement(BigInteger u, BigInteger v) + { + this.u = u; + this.v = v; + } + } +} diff --git a/crypto/src/math/ec/multiplier/ECMultiplier.cs b/crypto/src/math/ec/multiplier/ECMultiplier.cs new file mode 100644 index 000000000..c6d768ea8 --- /dev/null +++ b/crypto/src/math/ec/multiplier/ECMultiplier.cs @@ -0,0 +1,18 @@ +namespace Org.BouncyCastle.Math.EC.Multiplier +{ + /** + * Interface for classes encapsulating a point multiplication algorithm + * for <code>ECPoint</code>s. + */ + internal interface ECMultiplier + { + /** + * Multiplies the <code>ECPoint p</code> by <code>k</code>, i.e. + * <code>p</code> is added <code>k</code> times to itself. + * @param p The <code>ECPoint</code> to be multiplied. + * @param k The factor by which <code>p</code> i multiplied. + * @return <code>p</code> multiplied by <code>k</code>. + */ + ECPoint Multiply(ECPoint p, BigInteger k, PreCompInfo preCompInfo); + } +} diff --git a/crypto/src/math/ec/multiplier/FpNafMultiplier.cs b/crypto/src/math/ec/multiplier/FpNafMultiplier.cs new file mode 100644 index 000000000..f5a98501a --- /dev/null +++ b/crypto/src/math/ec/multiplier/FpNafMultiplier.cs @@ -0,0 +1,39 @@ +namespace Org.BouncyCastle.Math.EC.Multiplier +{ + /** + * Class implementing the NAF (Non-Adjacent Form) multiplication algorithm. + */ + internal class FpNafMultiplier + : ECMultiplier + { + /** + * D.3.2 pg 101 + * @see org.bouncycastle.math.ec.multiplier.ECMultiplier#multiply(org.bouncycastle.math.ec.ECPoint, java.math.BigInteger) + */ + public ECPoint Multiply(ECPoint p, BigInteger k, PreCompInfo preCompInfo) + { + // TODO Probably should try to add this + // BigInteger e = k.Mod(n); // n == order of p + BigInteger e = k; + BigInteger h = e.Multiply(BigInteger.Three); + + ECPoint neg = p.Negate(); + ECPoint R = p; + + for (int i = h.BitLength - 2; i > 0; --i) + { + R = R.Twice(); + + bool hBit = h.TestBit(i); + bool eBit = e.TestBit(i); + + if (hBit != eBit) + { + R = R.Add(hBit ? p : neg); + } + } + + return R; + } + } +} diff --git a/crypto/src/math/ec/multiplier/PreCompInfo.cs b/crypto/src/math/ec/multiplier/PreCompInfo.cs new file mode 100644 index 000000000..d379508c8 --- /dev/null +++ b/crypto/src/math/ec/multiplier/PreCompInfo.cs @@ -0,0 +1,11 @@ +namespace Org.BouncyCastle.Math.EC.Multiplier +{ + /** + * Interface for classes storing precomputation data for multiplication + * algorithms. Used as a Memento (see GOF patterns) for + * <code>WNafMultiplier</code>. + */ + internal interface PreCompInfo + { + } +} diff --git a/crypto/src/math/ec/multiplier/ReferenceMultiplier.cs b/crypto/src/math/ec/multiplier/ReferenceMultiplier.cs new file mode 100644 index 000000000..cdccffc2d --- /dev/null +++ b/crypto/src/math/ec/multiplier/ReferenceMultiplier.cs @@ -0,0 +1,30 @@ +namespace Org.BouncyCastle.Math.EC.Multiplier +{ + internal class ReferenceMultiplier + : ECMultiplier + { + /** + * Simple shift-and-add multiplication. Serves as reference implementation + * to verify (possibly faster) implementations in + * {@link org.bouncycastle.math.ec.ECPoint ECPoint}. + * + * @param p The point to multiply. + * @param k The factor by which to multiply. + * @return The result of the point multiplication <code>k * p</code>. + */ + public ECPoint Multiply(ECPoint p, BigInteger k, PreCompInfo preCompInfo) + { + ECPoint q = p.Curve.Infinity; + int t = k.BitLength; + for (int i = 0; i < t; i++) + { + if (k.TestBit(i)) + { + q = q.Add(p); + } + p = p.Twice(); + } + return q; + } + } +} diff --git a/crypto/src/math/ec/multiplier/WNafMultiplier.cs b/crypto/src/math/ec/multiplier/WNafMultiplier.cs new file mode 100644 index 000000000..b5cf34ba8 --- /dev/null +++ b/crypto/src/math/ec/multiplier/WNafMultiplier.cs @@ -0,0 +1,241 @@ +using System; + +namespace Org.BouncyCastle.Math.EC.Multiplier +{ + /** + * Class implementing the WNAF (Window Non-Adjacent Form) multiplication + * algorithm. + */ + internal class WNafMultiplier + : ECMultiplier + { + /** + * Computes the Window NAF (non-adjacent Form) of an integer. + * @param width The width <code>w</code> of the Window NAF. The width is + * defined as the minimal number <code>w</code>, such that for any + * <code>w</code> consecutive digits in the resulting representation, at + * most one is non-zero. + * @param k The integer of which the Window NAF is computed. + * @return The Window NAF of the given width, such that the following holds: + * <code>k = −<sub>i=0</sub><sup>l-1</sup> k<sub>i</sub>2<sup>i</sup> + * </code>, where the <code>k<sub>i</sub></code> denote the elements of the + * returned <code>sbyte[]</code>. + */ + public sbyte[] WindowNaf(sbyte width, BigInteger k) + { + // The window NAF is at most 1 element longer than the binary + // representation of the integer k. sbyte can be used instead of short or + // int unless the window width is larger than 8. For larger width use + // short or int. However, a width of more than 8 is not efficient for + // m = log2(q) smaller than 2305 Bits. Note: Values for m larger than + // 1000 Bits are currently not used in practice. + sbyte[] wnaf = new sbyte[k.BitLength + 1]; + + // 2^width as short and BigInteger + short pow2wB = (short)(1 << width); + BigInteger pow2wBI = BigInteger.ValueOf(pow2wB); + + int i = 0; + + // The actual length of the WNAF + int length = 0; + + // while k >= 1 + while (k.SignValue > 0) + { + // if k is odd + if (k.TestBit(0)) + { + // k Mod 2^width + BigInteger remainder = k.Mod(pow2wBI); + + // if remainder > 2^(width - 1) - 1 + if (remainder.TestBit(width - 1)) + { + wnaf[i] = (sbyte)(remainder.IntValue - pow2wB); + } + else + { + wnaf[i] = (sbyte)remainder.IntValue; + } + // wnaf[i] is now in [-2^(width-1), 2^(width-1)-1] + + k = k.Subtract(BigInteger.ValueOf(wnaf[i])); + length = i; + } + else + { + wnaf[i] = 0; + } + + // k = k/2 + k = k.ShiftRight(1); + i++; + } + + length++; + + // Reduce the WNAF array to its actual length + sbyte[] wnafShort = new sbyte[length]; + Array.Copy(wnaf, 0, wnafShort, 0, length); + return wnafShort; + } + + /** + * Multiplies <code>this</code> by an integer <code>k</code> using the + * Window NAF method. + * @param k The integer by which <code>this</code> is multiplied. + * @return A new <code>ECPoint</code> which equals <code>this</code> + * multiplied by <code>k</code>. + */ + public ECPoint Multiply(ECPoint p, BigInteger k, PreCompInfo preCompInfo) + { + WNafPreCompInfo wnafPreCompInfo; + + if ((preCompInfo != null) && (preCompInfo is WNafPreCompInfo)) + { + wnafPreCompInfo = (WNafPreCompInfo)preCompInfo; + } + else + { + // Ignore empty PreCompInfo or PreCompInfo of incorrect type + wnafPreCompInfo = new WNafPreCompInfo(); + } + + // floor(log2(k)) + int m = k.BitLength; + + // width of the Window NAF + sbyte width; + + // Required length of precomputation array + int reqPreCompLen; + + // Determine optimal width and corresponding length of precomputation + // array based on literature values + if (m < 13) + { + width = 2; + reqPreCompLen = 1; + } + else + { + if (m < 41) + { + width = 3; + reqPreCompLen = 2; + } + else + { + if (m < 121) + { + width = 4; + reqPreCompLen = 4; + } + else + { + if (m < 337) + { + width = 5; + reqPreCompLen = 8; + } + else + { + if (m < 897) + { + width = 6; + reqPreCompLen = 16; + } + else + { + if (m < 2305) + { + width = 7; + reqPreCompLen = 32; + } + else + { + width = 8; + reqPreCompLen = 127; + } + } + } + } + } + } + + // The length of the precomputation array + int preCompLen = 1; + + ECPoint[] preComp = wnafPreCompInfo.GetPreComp(); + ECPoint twiceP = wnafPreCompInfo.GetTwiceP(); + + // Check if the precomputed ECPoints already exist + if (preComp == null) + { + // Precomputation must be performed from scratch, create an empty + // precomputation array of desired length + preComp = new ECPoint[]{ p }; + } + else + { + // Take the already precomputed ECPoints to start with + preCompLen = preComp.Length; + } + + if (twiceP == null) + { + // Compute twice(p) + twiceP = p.Twice(); + } + + if (preCompLen < reqPreCompLen) + { + // Precomputation array must be made bigger, copy existing preComp + // array into the larger new preComp array + ECPoint[] oldPreComp = preComp; + preComp = new ECPoint[reqPreCompLen]; + Array.Copy(oldPreComp, 0, preComp, 0, preCompLen); + + for (int i = preCompLen; i < reqPreCompLen; i++) + { + // Compute the new ECPoints for the precomputation array. + // The values 1, 3, 5, ..., 2^(width-1)-1 times p are + // computed + preComp[i] = twiceP.Add(preComp[i - 1]); + } + } + + // Compute the Window NAF of the desired width + sbyte[] wnaf = WindowNaf(width, k); + int l = wnaf.Length; + + // Apply the Window NAF to p using the precomputed ECPoint values. + ECPoint q = p.Curve.Infinity; + for (int i = l - 1; i >= 0; i--) + { + q = q.Twice(); + + if (wnaf[i] != 0) + { + if (wnaf[i] > 0) + { + q = q.Add(preComp[(wnaf[i] - 1)/2]); + } + else + { + // wnaf[i] < 0 + q = q.Subtract(preComp[(-wnaf[i] - 1)/2]); + } + } + } + + // Set PreCompInfo in ECPoint, such that it is available for next + // multiplication. + wnafPreCompInfo.SetPreComp(preComp); + wnafPreCompInfo.SetTwiceP(twiceP); + p.SetPreCompInfo(wnafPreCompInfo); + return q; + } + } +} diff --git a/crypto/src/math/ec/multiplier/WNafPreCompInfo.cs b/crypto/src/math/ec/multiplier/WNafPreCompInfo.cs new file mode 100644 index 000000000..d9305dace --- /dev/null +++ b/crypto/src/math/ec/multiplier/WNafPreCompInfo.cs @@ -0,0 +1,46 @@ +namespace Org.BouncyCastle.Math.EC.Multiplier +{ + /** + * Class holding precomputation data for the WNAF (Window Non-Adjacent Form) + * algorithm. + */ + internal class WNafPreCompInfo + : PreCompInfo + { + /** + * Array holding the precomputed <code>ECPoint</code>s used for the Window + * NAF multiplication in <code> + * {@link org.bouncycastle.math.ec.multiplier.WNafMultiplier.multiply() + * WNafMultiplier.multiply()}</code>. + */ + private ECPoint[] preComp = null; + + /** + * Holds an <code>ECPoint</code> representing twice(this). Used for the + * Window NAF multiplication in <code> + * {@link org.bouncycastle.math.ec.multiplier.WNafMultiplier.multiply() + * WNafMultiplier.multiply()}</code>. + */ + private ECPoint twiceP = null; + + internal ECPoint[] GetPreComp() + { + return preComp; + } + + internal void SetPreComp(ECPoint[] preComp) + { + this.preComp = preComp; + } + + internal ECPoint GetTwiceP() + { + return twiceP; + } + + internal void SetTwiceP(ECPoint twiceThis) + { + this.twiceP = twiceThis; + } + } +} diff --git a/crypto/src/math/ec/multiplier/WTauNafMultiplier.cs b/crypto/src/math/ec/multiplier/WTauNafMultiplier.cs new file mode 100644 index 000000000..f1a605770 --- /dev/null +++ b/crypto/src/math/ec/multiplier/WTauNafMultiplier.cs @@ -0,0 +1,120 @@ +using System; + +using Org.BouncyCastle.Math.EC.Abc; + +namespace Org.BouncyCastle.Math.EC.Multiplier +{ + /** + * Class implementing the WTNAF (Window + * <code>τ</code>-adic Non-Adjacent Form) algorithm. + */ + internal class WTauNafMultiplier + : ECMultiplier + { + /** + * Multiplies a {@link org.bouncycastle.math.ec.F2mPoint F2mPoint} + * by <code>k</code> using the reduced <code>τ</code>-adic NAF (RTNAF) + * method. + * @param p The F2mPoint to multiply. + * @param k The integer by which to multiply <code>k</code>. + * @return <code>p</code> multiplied by <code>k</code>. + */ + public ECPoint Multiply(ECPoint point, BigInteger k, PreCompInfo preCompInfo) + { + if (!(point is F2mPoint)) + throw new ArgumentException("Only F2mPoint can be used in WTauNafMultiplier"); + + F2mPoint p = (F2mPoint)point; + + F2mCurve curve = (F2mCurve) p.Curve; + int m = curve.M; + sbyte a = (sbyte) curve.A.ToBigInteger().IntValue; + sbyte mu = curve.GetMu(); + BigInteger[] s = curve.GetSi(); + + ZTauElement rho = Tnaf.PartModReduction(k, m, a, s, mu, (sbyte)10); + + return MultiplyWTnaf(p, rho, preCompInfo, a, mu); + } + + /** + * Multiplies a {@link org.bouncycastle.math.ec.F2mPoint F2mPoint} + * by an element <code>λ</code> of <code><b>Z</b>[τ]</code> using + * the <code>τ</code>-adic NAF (TNAF) method. + * @param p The F2mPoint to multiply. + * @param lambda The element <code>λ</code> of + * <code><b>Z</b>[τ]</code> of which to compute the + * <code>[τ]</code>-adic NAF. + * @return <code>p</code> multiplied by <code>λ</code>. + */ + private F2mPoint MultiplyWTnaf(F2mPoint p, ZTauElement lambda, + PreCompInfo preCompInfo, sbyte a, sbyte mu) + { + ZTauElement[] alpha; + if (a == 0) + { + alpha = Tnaf.Alpha0; + } + else + { + // a == 1 + alpha = Tnaf.Alpha1; + } + + BigInteger tw = Tnaf.GetTw(mu, Tnaf.Width); + + sbyte[]u = Tnaf.TauAdicWNaf(mu, lambda, Tnaf.Width, + BigInteger.ValueOf(Tnaf.Pow2Width), tw, alpha); + + return MultiplyFromWTnaf(p, u, preCompInfo); + } + + /** + * Multiplies a {@link org.bouncycastle.math.ec.F2mPoint F2mPoint} + * by an element <code>λ</code> of <code><b>Z</b>[τ]</code> + * using the window <code>τ</code>-adic NAF (TNAF) method, given the + * WTNAF of <code>λ</code>. + * @param p The F2mPoint to multiply. + * @param u The the WTNAF of <code>λ</code>.. + * @return <code>λ * p</code> + */ + private static F2mPoint MultiplyFromWTnaf(F2mPoint p, sbyte[] u, + PreCompInfo preCompInfo) + { + F2mCurve curve = (F2mCurve)p.Curve; + sbyte a = (sbyte) curve.A.ToBigInteger().IntValue; + + F2mPoint[] pu; + if ((preCompInfo == null) || !(preCompInfo is WTauNafPreCompInfo)) + { + pu = Tnaf.GetPreComp(p, a); + p.SetPreCompInfo(new WTauNafPreCompInfo(pu)); + } + else + { + pu = ((WTauNafPreCompInfo)preCompInfo).GetPreComp(); + } + + // q = infinity + F2mPoint q = (F2mPoint) p.Curve.Infinity; + for (int i = u.Length - 1; i >= 0; i--) + { + q = Tnaf.Tau(q); + if (u[i] != 0) + { + if (u[i] > 0) + { + q = q.AddSimple(pu[u[i]]); + } + else + { + // u[i] < 0 + q = q.SubtractSimple(pu[-u[i]]); + } + } + } + + return q; + } + } +} diff --git a/crypto/src/math/ec/multiplier/WTauNafPreCompInfo.cs b/crypto/src/math/ec/multiplier/WTauNafPreCompInfo.cs new file mode 100644 index 000000000..cede4a05d --- /dev/null +++ b/crypto/src/math/ec/multiplier/WTauNafPreCompInfo.cs @@ -0,0 +1,41 @@ +namespace Org.BouncyCastle.Math.EC.Multiplier +{ + /** + * Class holding precomputation data for the WTNAF (Window + * <code>τ</code>-adic Non-Adjacent Form) algorithm. + */ + internal class WTauNafPreCompInfo + : PreCompInfo + { + /** + * Array holding the precomputed <code>F2mPoint</code>s used for the + * WTNAF multiplication in <code> + * {@link org.bouncycastle.math.ec.multiplier.WTauNafMultiplier.multiply() + * WTauNafMultiplier.multiply()}</code>. + */ + private readonly F2mPoint[] preComp; + + /** + * Constructor for <code>WTauNafPreCompInfo</code> + * @param preComp Array holding the precomputed <code>F2mPoint</code>s + * used for the WTNAF multiplication in <code> + * {@link org.bouncycastle.math.ec.multiplier.WTauNafMultiplier.multiply() + * WTauNafMultiplier.multiply()}</code>. + */ + internal WTauNafPreCompInfo(F2mPoint[] preComp) + { + this.preComp = preComp; + } + + /** + * @return the array holding the precomputed <code>F2mPoint</code>s + * used for the WTNAF multiplication in <code> + * {@link org.bouncycastle.math.ec.multiplier.WTauNafMultiplier.multiply() + * WTauNafMultiplier.multiply()}</code>. + */ + internal F2mPoint[] GetPreComp() + { + return preComp; + } + } +} |