summary refs log tree commit diff
path: root/Crypto/src/math
diff options
context:
space:
mode:
Diffstat (limited to 'Crypto/src/math')
-rw-r--r--Crypto/src/math/BigInteger.cs3141
-rw-r--r--Crypto/src/math/ec/ECAlgorithms.cs93
-rw-r--r--Crypto/src/math/ec/ECCurve.cs661
-rw-r--r--Crypto/src/math/ec/ECFieldElement.cs1253
-rw-r--r--Crypto/src/math/ec/ECPoint.cs567
-rw-r--r--Crypto/src/math/ec/IntArray.cs485
-rw-r--r--Crypto/src/math/ec/abc/SimpleBigDecimal.cs241
-rw-r--r--Crypto/src/math/ec/abc/Tnaf.cs834
-rw-r--r--Crypto/src/math/ec/abc/ZTauElement.cs36
-rw-r--r--Crypto/src/math/ec/multiplier/ECMultiplier.cs18
-rw-r--r--Crypto/src/math/ec/multiplier/FpNafMultiplier.cs39
-rw-r--r--Crypto/src/math/ec/multiplier/PreCompInfo.cs11
-rw-r--r--Crypto/src/math/ec/multiplier/ReferenceMultiplier.cs30
-rw-r--r--Crypto/src/math/ec/multiplier/WNafMultiplier.cs241
-rw-r--r--Crypto/src/math/ec/multiplier/WNafPreCompInfo.cs46
-rw-r--r--Crypto/src/math/ec/multiplier/WTauNafMultiplier.cs120
-rw-r--r--Crypto/src/math/ec/multiplier/WTauNafPreCompInfo.cs41
17 files changed, 7857 insertions, 0 deletions
diff --git a/Crypto/src/math/BigInteger.cs b/Crypto/src/math/BigInteger.cs
new file mode 100644
index 000000000..d52c0f83c
--- /dev/null
+++ b/Crypto/src/math/BigInteger.cs
@@ -0,0 +1,3141 @@
+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 || PORTABLE)
+	[Serializable]
+#endif
+    public class BigInteger
+	{
+		// The primes b/w 2 and ~2^10
+		/*
+				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
+		*/
+
+		// Each list has a product < 2^31
+		private 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 },
+		};
+
+		private static readonly int[] primeProducts;
+
+		private const long IMASK = 0xffffffffL;
+		private static readonly ulong UIMASK = (ulong)IMASK;
+
+		private static readonly int[] ZeroMagnitude = new int[0];
+		private static readonly byte[] ZeroEncoding = new byte[0];
+
+		public static readonly BigInteger Zero = new BigInteger(0, ZeroMagnitude, false);
+		public static readonly BigInteger One = createUValueOf(1);
+		public static readonly BigInteger Two = createUValueOf(2);
+		public static readonly BigInteger Three = createUValueOf(3);
+		public static readonly BigInteger Ten = createUValueOf(10);
+
+		private static readonly int chunk2 = 1; // TODO Parse 64 bits at a time
+		private static readonly BigInteger radix2 = ValueOf(2);
+		private static readonly BigInteger radix2E = radix2.Pow(chunk2);
+
+		private static readonly int chunk10 = 19;
+		private static readonly BigInteger radix10 = ValueOf(10);
+		private static readonly BigInteger radix10E = radix10.Pow(chunk10);
+
+		private static readonly int chunk16 = 16;
+		private static readonly BigInteger radix16 = ValueOf(16);
+		private static readonly BigInteger radix16E = radix16.Pow(chunk16);
+
+		private static readonly Random RandomSource = new Random();
+
+		private const int BitsPerByte = 8;
+		private const int BitsPerInt = 32;
+		private const int BytesPerInt = 4;
+
+		static BigInteger()
+		{
+			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 sign; // -1 means -ve; +1 means +ve; 0 means 0;
+		private int[] magnitude; // array of ints with [0] being the most significant
+		private int nBits = -1; // cache BitCount() value
+		private int nBitLength = -1; // cache calcBitLength() value
+		private long mQuote = -1L; // -m^(-1) mod b, b = 2^32 (see Montgomery mult.)
+
+		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 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, 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 > 1)
+								throw new FormatException("Bad character in radix 2 string: " + s);
+
+							// TODO Parse 64 bits at a time
+							b = b.ShiftLeft(1);
+							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 == 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
+			b[0] &= rndMask[BitsPerByte * nBytes - sizeInBits];
+
+			this.magnitude = MakeMagnitude(b, 0, b.Length);
+			this.sign = this.magnitude.Length < 1 ? 0 : 1;
+		}
+
+		private static readonly byte[] rndMask = { 255, 127, 63, 31, 15, 7, 3, 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 = rndMask[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 = -1L;
+
+				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 = -1L;
+
+						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 += bitCounts[(byte) magnitude[i]];
+							sum += bitCounts[(byte)(magnitude[i] >> 8)];
+							sum += bitCounts[(byte)(magnitude[i] >> 16)];
+							sum += bitCounts[(byte)(magnitude[i] >> 24)];
+						}
+						nBits = sum;
+					}
+				}
+
+				return nBits;
+			}
+		}
+
+		private readonly static byte[] bitCounts =
+		{
+			0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1,
+			2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4,
+			4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3,
+			4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 2, 3, 3, 4, 3, 4, 4, 5,
+			3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 1, 2, 2, 3, 2,
+			3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 2, 3,
+			3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6,
+			7, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6,
+			5, 6, 6, 7, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 4, 5, 5, 6, 5, 6, 6, 7, 5,
+			6, 6, 7, 6, 7, 7, 8
+		};
+
+		private int calcBitLength(
+			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(0, magnitude);
+				}
+
+				return nBitLength;
+			}
+		}
+
+		//
+		// BitLen(value) is the number of bits in value.
+		//
+		private static int BitLen(
+			int w)
+		{
+			// Binary search - decision tree (5 tests, rarely 6)
+			return (w < 1 << 15 ? (w < 1 << 7
+				? (w < 1 << 3 ? (w < 1 << 1
+				? (w < 1 << 0 ? (w < 0 ? 32 : 0) : 1)
+				: (w < 1 << 2 ? 2 : 3)) : (w < 1 << 5
+				? (w < 1 << 4 ? 4 : 5)
+				: (w < 1 << 6 ? 6 : 7)))
+				: (w < 1 << 11
+				? (w < 1 << 9 ? (w < 1 << 8 ? 8 : 9) : (w < 1 << 10 ? 10 : 11))
+				: (w < 1 << 13 ? (w < 1 << 12 ? 12 : 13) : (w < 1 << 14 ? 14 : 15)))) : (w < 1 << 23 ? (w < 1 << 19
+				? (w < 1 << 17 ? (w < 1 << 16 ? 16 : 17) : (w < 1 << 18 ? 18 : 19))
+				: (w < 1 << 21 ? (w < 1 << 20 ? 20 : 21) : (w < 1 << 22 ? 22 : 23))) : (w < 1 << 27
+				? (w < 1 << 25 ? (w < 1 << 24 ? 24 : 25) : (w < 1 << 26 ? 26 : 27))
+				: (w < 1 << 29 ? (w < 1 << 28 ? 28 : 29) : (w < 1 << 30 ? 30 : 31)))));
+		}
+
+//		private readonly static byte[] bitLengths =
+//		{
+//			0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,
+//			5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
+//			6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+//			7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+//			7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8,
+//			8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+//			8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+//			8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+//			8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+//			8, 8, 8, 8, 8, 8, 8, 8
+//		};
+
+		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(yStart, y);
+				int xBitLength = calcBitLength(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;
+
+			if (biggie.sign != sign || biggie.magnitude.Length != magnitude.Length)
+				return false;
+
+			for (int i = 0; i < magnitude.Length; i++)
+			{
+				if (biggie.magnitude[i] != magnitude[i])
+				{
+					return false;
+				}
+			}
+
+			return true;
+		}
+
+		public BigInteger Gcd(
+			BigInteger 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
+			{
+				return sign == 0 ? 0
+					: sign > 0 ? magnitude[magnitude.Length - 1]
+					: -magnitude[magnitude.Length - 1];
+			}
+		}
+
+		/**
+		 * 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;
+		}
+
+		internal 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;
+			BigInteger nMinusOne = n.Subtract(One);
+			int s = nMinusOne.GetLowestSetBit();
+			BigInteger r = nMinusOne.ShiftRight(s);
+
+			Debug.Assert(s >= 1);
+
+			do
+			{
+				// TODO Make a method for random BigIntegers in range 0 < x < n)
+				// - Method can be optimized by only replacing examined bits at each trial
+				BigInteger a;
+				do
+				{
+					a = new BigInteger(n.BitLength, random);
+				}
+				while (a.CompareTo(One) <= 0 || a.CompareTo(nMinusOne) >= 0);
+
+				BigInteger y = a.ModPow(r, n);
+
+				if (!y.Equals(One))
+				{
+					int j = 0;
+					while (!y.Equals(nMinusOne))
+					{
+						if (++j == s)
+							return false;
+
+						y = y.ModPow(Two, n);
+
+						if (y.Equals(One))
+							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;
+
+				long v;
+				if (magnitude.Length > 1)
+				{
+					v = ((long)magnitude[magnitude.Length - 2] << 32)
+						| (magnitude[magnitude.Length - 1] & IMASK);
+				}
+				else
+				{
+					v = (magnitude[magnitude.Length - 1] & IMASK);
+				}
+
+				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);
+//				}
+//			}
+
+            BigInteger x;
+            BigInteger gcd = ExtEuclid(this.Mod(m), m, out x);
+
+			if (!gcd.Equals(One))
+				throw new ArithmeticException("Numbers not relatively prime.");
+
+			if (x.sign < 0)
+			{
+                x = x.Add(m);
+			}
+
+			return x;
+		}
+
+		/**
+		 * Calculate the numbers u1, u2, and u3 such that:
+		 *
+		 * u1 * a + u2 * b = u3
+		 *
+		 * where u3 is the greatest common divider of a and b.
+		 * a and b using the extended Euclid algorithm (refer p. 323
+		 * of The Art of Computer Programming vol 2, 2nd ed).
+		 * This also seems to have the side effect of calculating
+		 * some form of multiplicative inverse.
+		 *
+		 * @param a    First number to calculate gcd for
+		 * @param b    Second number to calculate gcd for
+		 * @param u1Out      the return object for the u1 value
+		 * @param u2Out      the return object for the u2 value
+		 * @return     The greatest common divisor of a and b
+		 */
+		private static BigInteger ExtEuclid(
+			BigInteger	    a,
+			BigInteger	    b,
+			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 exponent,
+			BigInteger m)
+		{
+			if (m.sign < 1)
+				throw new ArithmeticException("Modulus must be positive");
+
+			if (m.Equals(One))
+				return Zero;
+
+			if (exponent.sign == 0)
+				return One;
+
+			if (sign == 0)
+				return Zero;
+
+			int[] zVal = null;
+			int[] yAccum = null;
+			int[] yVal;
+
+			// Montgomery exponentiation is only possible if the modulus is odd,
+			// but AFAIK, this is always the case for crypto algo's
+			bool useMonty = ((m.magnitude[m.magnitude.Length - 1] & 1) == 1);
+			long mQ = 0;
+			if (useMonty)
+			{
+				mQ = m.GetMQuote();
+
+				// tmp = this * R mod m
+				BigInteger tmp = ShiftLeft(32 * m.magnitude.Length).Mod(m);
+				zVal = tmp.magnitude;
+
+				useMonty = (zVal.Length <= m.magnitude.Length);
+
+				if (useMonty)
+				{
+					yAccum = new int[m.magnitude.Length + 1];
+					if (zVal.Length < m.magnitude.Length)
+					{
+						int[] longZ = new int[m.magnitude.Length];
+						zVal.CopyTo(longZ, longZ.Length - zVal.Length);
+						zVal = longZ;
+					}
+				}
+			}
+
+			if (!useMonty)
+			{
+				if (magnitude.Length <= m.magnitude.Length)
+				{
+					//zAccum = new int[m.magnitude.Length * 2];
+					zVal = new int[m.magnitude.Length];
+					magnitude.CopyTo(zVal, zVal.Length - magnitude.Length);
+				}
+				else
+				{
+					//
+					// in normal practice we'll never see this...
+					//
+					BigInteger tmp = Remainder(m);
+
+					//zAccum = new int[m.magnitude.Length * 2];
+					zVal = new int[m.magnitude.Length];
+					tmp.magnitude.CopyTo(zVal, zVal.Length - tmp.magnitude.Length);
+				}
+
+				yAccum = new int[m.magnitude.Length * 2];
+			}
+
+			yVal = new int[m.magnitude.Length];
+
+			//
+			// from LSW to MSW
+			//
+			for (int i = 0; i < exponent.magnitude.Length; i++)
+			{
+				int v = exponent.magnitude[i];
+				int bits = 0;
+
+				if (i == 0)
+				{
+					while (v > 0)
+					{
+						v <<= 1;
+						bits++;
+					}
+
+					//
+					// first time in initialise y
+					//
+					zVal.CopyTo(yVal, 0);
+
+					v <<= 1;
+					bits++;
+				}
+
+				while (v != 0)
+				{
+					if (useMonty)
+					{
+						// Montgomery square algo doesn't exist, and a normal
+						// square followed by a Montgomery reduction proved to
+						// be almost as heavy as a Montgomery mulitply.
+						MultiplyMonty(yAccum, yVal, yVal, m.magnitude, mQ);
+					}
+					else
+					{
+						Square(yAccum, yVal);
+						Remainder(yAccum, m.magnitude);
+						Array.Copy(yAccum, yAccum.Length - yVal.Length, yVal, 0, yVal.Length);
+						ZeroOut(yAccum);
+					}
+					bits++;
+
+					if (v < 0)
+					{
+						if (useMonty)
+						{
+							MultiplyMonty(yAccum, yVal, zVal, m.magnitude, mQ);
+						}
+						else
+						{
+							Multiply(yAccum, yVal, zVal);
+							Remainder(yAccum, m.magnitude);
+							Array.Copy(yAccum, yAccum.Length - yVal.Length, yVal, 0,
+								yVal.Length);
+							ZeroOut(yAccum);
+						}
+					}
+
+					v <<= 1;
+				}
+
+				while (bits < 32)
+				{
+					if (useMonty)
+					{
+						MultiplyMonty(yAccum, yVal, yVal, m.magnitude, mQ);
+					}
+					else
+					{
+						Square(yAccum, yVal);
+						Remainder(yAccum, m.magnitude);
+						Array.Copy(yAccum, yAccum.Length - yVal.Length, yVal, 0, yVal.Length);
+						ZeroOut(yAccum);
+					}
+					bits++;
+				}
+			}
+
+			if (useMonty)
+			{
+				// Return y * R^(-1) mod m by doing y * 1 * R^(-1) mod m
+				ZeroOut(zVal);
+				zVal[zVal.Length - 1] = 1;
+				MultiplyMonty(yAccum, yVal, zVal, m.magnitude, mQ);
+			}
+
+			BigInteger result = new BigInteger(1, yVal, true);
+
+			return exponent.sign > 0
+				?	result
+				:	result.ModInverse(m);
+		}
+
+		/**
+		 * 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 u1, u2, c;
+
+			int wBase = w.Length - 1;
+
+			for (int i = x.Length - 1; i != 0; i--)
+			{
+				ulong v = (ulong)(uint) x[i];
+
+				u1 = v * v;
+				u2 = u1 >> 32;
+				u1 = (uint) u1;
+
+				u1 += (ulong)(uint) w[wBase];
+
+				w[wBase] = (int)(uint) u1;
+				c = u2 + (u1 >> 32);
+
+				for (int j = i - 1; j >= 0; j--)
+				{
+					--wBase;
+					u1 = v * (ulong)(uint) x[j];
+					u2 = u1 >> 31; // multiply by 2!
+					u1 = (uint)(u1 << 1); // multiply by 2!
+					u1 += c + (ulong)(uint) w[wBase];
+
+					w[wBase] = (int)(uint) u1;
+					c = u2 + (u1 >> 32);
+				}
+
+				c += (ulong)(uint) w[--wBase];
+				w[wBase] = (int)(uint) c;
+
+				if (--wBase >= 0)
+				{
+					w[wBase] = (int)(uint)(c >> 32);
+				}
+				else
+				{
+					Debug.Assert((uint)(c >> 32) == 0);
+				}
+				wBase += i;
+			}
+
+			u1 = (ulong)(uint) x[0];
+			u1 = u1 * u1;
+			u2 = u1 >> 32;
+			u1 = u1 & IMASK;
+
+			u1 += (ulong)(uint) w[wBase];
+
+			w[wBase] = (int)(uint) u1;
+			if (--wBase >= 0)
+			{
+				w[wBase] = (int)(uint)(u2 + (u1 >> 32) + (ulong)(uint) w[wBase]);
+			}
+			else
+			{
+				Debug.Assert((uint)(u2 + (u1 >> 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;
+		}
+
+		private static long FastExtEuclid(
+			long	a,
+			long	b,
+			long[]	uOut)
+		{
+			long u1 = 1;
+			long u3 = a;
+			long v1 = 0;
+			long v3 = b;
+
+			while (v3 > 0)
+			{
+				long q, tn;
+
+				q = u3 / v3;
+
+				tn = u1 - (v1 * q);
+				u1 = v1;
+				v1 = tn;
+
+				tn = u3 - (v3 * q);
+				u3 = v3;
+				v3 = tn;
+			}
+
+			uOut[0] = u1;
+			uOut[1] = (u3 - (u1 * a)) / b;
+
+			return u3;
+		}
+
+		private static long FastModInverse(
+			long	v,
+			long	m)
+		{
+			if (m < 1)
+				throw new ArithmeticException("Modulus must be positive");
+
+			long[] x = new long[2];
+			long gcd = FastExtEuclid(v, m, x);
+
+			if (gcd != 1)
+				throw new ArithmeticException("Numbers not relatively prime.");
+
+			if (x[0] < 0)
+			{
+				x[0] += m;
+			}
+
+			return x[0];
+		}
+
+//		private static BigInteger MQuoteB = One.ShiftLeft(32);
+//		private static BigInteger MQuoteBSub1 = MQuoteB.Subtract(One);
+
+		/**
+		 * Calculate mQuote = -m^(-1) mod b with b = 2^32 (32 = word size)
+		 */
+		private long GetMQuote()
+		{
+			Debug.Assert(this.sign > 0);
+
+			if (mQuote != -1)
+			{
+				return mQuote; // already calculated
+			}
+
+			if (magnitude.Length == 0 || (magnitude[magnitude.Length - 1] & 1) == 0)
+			{
+				return -1; // not for even numbers
+			}
+
+			long v = (((~this.magnitude[this.magnitude.Length - 1]) | 1) & 0xffffffffL);
+			mQuote = FastModInverse(v, 0x100000000L);
+
+			return mQuote;
+		}
+
+		/**
+		 * 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,
+			long	mQuote)
+			// mQuote = -m^(-1) mod b
+		{
+			if (m.Length == 1)
+			{
+				x[0] = (int)MultiplyMontyNIsOne((uint)x[0], (uint)y[0], (uint)m[0], (ulong)mQuote);
+				return;
+			}
+
+			int n = m.Length;
+			int nMinus1 = n - 1;
+			long y_0 = y[nMinus1] & IMASK;
+
+			// 1. a = 0 (Notation: a = (a_{n} a_{n-1} ... a_{0})_{b} )
+			Array.Clear(a, 0, n + 1);
+
+			// 2. for i from 0 to (n - 1) do the following:
+			for (int i = n; i > 0; i--)
+			{
+				long x_i = x[i - 1] & IMASK;
+
+				// 2.1 u = ((a[0] + (x[i] * y[0]) * mQuote) mod b
+				long u = ((((a[n] & IMASK) + ((x_i * y_0) & IMASK)) & IMASK) * mQuote) & IMASK;
+
+				// 2.2 a = (a + x_i * y + u * m) / b
+				long prod1 = x_i * y_0;
+				long prod2 = u * (m[nMinus1] & IMASK);
+				long tmp = (a[n] & IMASK) + (prod1 & IMASK) + (prod2 & IMASK);
+				long carry = (long)((ulong)prod1 >> 32) + (long)((ulong)prod2 >> 32) + (long)((ulong)tmp >> 32);
+				for (int j = nMinus1; j > 0; j--)
+				{
+					prod1 = x_i * (y[j - 1] & IMASK);
+					prod2 = u * (m[j - 1] & IMASK);
+					tmp = (a[j] & IMASK) + (prod1 & IMASK) + (prod2 & IMASK) + (carry & IMASK);
+					carry = (long)((ulong)carry >> 32) + (long)((ulong)prod1 >> 32) +
+						(long)((ulong)prod2 >> 32) + (long)((ulong)tmp >> 32);
+					a[j + 1] = (int)tmp; // division by b
+				}
+				carry += (a[0] & IMASK);
+				a[1] = (int)carry;
+				a[0] = (int)((ulong)carry >> 32); // OJO!!!!!
+			}
+
+			// 3. if x >= m the x = x - m
+			if (CompareTo(0, a, 0, m) >= 0)
+			{
+				Subtract(0, a, 0, m);
+			}
+
+			// put the result in x
+			Array.Copy(a, 1, x, 0, n);
+		}
+
+		private static uint MultiplyMontyNIsOne(
+			uint	x,
+			uint	y,
+			uint	m,
+			ulong	mQuote)
+		{
+			ulong um = m;
+			ulong prod1 = (ulong)x * (ulong)y;
+			ulong u = (prod1 * mQuote) & UIMASK;
+			ulong prod2 = u * um;
+			ulong tmp = (prod1 & UIMASK) + (prod2 & UIMASK);
+			ulong carry = (prod1 >> 32) + (prod2 >> 32) + (tmp >> 32);
+
+			if (carry > um)
+			{
+				carry -= um;
+			}
+
+			return (uint)(carry & UIMASK);
+		}
+
+		public BigInteger Multiply(
+			BigInteger val)
+		{
+			if (sign == 0 || 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 = (this.BitLength + val.BitLength) / BitsPerInt + 1;
+			int[] res = new int[resLength];
+
+			if (val == this)
+			{
+				Square(res, this.magnitude);
+			}
+			else
+			{
+				Multiply(res, this.magnitude, val.magnitude);
+			}
+
+			return new BigInteger(sign * val.sign, res, true);
+		}
+
+		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)
+			{
+				throw new ArithmeticException("Negative exponent");
+			}
+
+			if (exp == 0)
+			{
+				return One;
+			}
+
+			if (sign == 0 || Equals(One))
+			{
+				return this;
+			}
+
+			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 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(yStart, y);
+				int xBitLength = calcBitLength(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 hiBits = n % 32;
+			if (hiBits != 0)
+			{
+				result[0] &= ~(-1 << hiBits);
+			}
+
+			return result;
+		}
+
+		/**
+		 * 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;
+		}
+
+		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 <= 16)
+
+			switch (radix)
+			{
+				case 2:
+				case 10:
+				case 16:
+					break;
+				default:
+					throw new FormatException("Only bases 2, 10, 16 are allowed");
+			}
+
+			// NB: Can only happen to internally managed instances
+			if (magnitude == null)
+				return "null";
+
+			if (sign == 0)
+				return "0";
+
+			Debug.Assert(magnitude.Length > 0);
+
+			StringBuilder sb = new StringBuilder();
+
+			if (radix == 16)
+			{
+				sb.Append(magnitude[0].ToString("x"));
+
+				for (int i = 1; i < magnitude.Length; i++)
+				{
+					sb.Append(magnitude[i].ToString("x8"));
+				}
+			}
+			else if (radix == 2)
+			{
+				sb.Append('1');
+
+				for (int i = BitLength - 2; i >= 0; --i)
+				{
+					sb.Append(TestBit(i) ? '1' : '0');
+				}
+			}
+			else
+			{
+				// This is algorithm 1a from chapter 4.4 in Seminumerical Algorithms, slow but it works
+				IList S = Platform.CreateArrayList();
+				BigInteger bs = ValueOf(radix);
+
+				// The sign is handled separatly.
+				// Notice however that for this to work, radix 16 _MUST_ be a special case,
+				// unless we want to enter a recursion well. In their infinite wisdom, why did not
+				// the Sun engineers made a c'tor for BigIntegers taking a BigInteger as parameter?
+				// (Answer: Becuase Sun's BigIntger is clonable, something bouncycastle's isn't.)
+//				BigInteger u = new BigInteger(Abs().ToString(16), 16);
+				BigInteger u = this.Abs();
+				BigInteger b;
+
+				while (u.sign != 0)
+				{
+					b = u.Mod(bs);
+					if (b.sign == 0)
+					{
+						S.Add("0");
+					}
+					else
+					{
+						// see how to interact with different bases
+						S.Add(b.magnitude[0].ToString("d"));
+					}
+					u = u.Divide(bs);
+				}
+
+				// Then pop the stack
+                for (int i = S.Count - 1; i >= 0; --i)
+                {
+                    sb.Append((string)S[i]);
+                }
+			}
+
+			string s = sb.ToString();
+
+			Debug.Assert(s.Length > 0);
+
+			// Strip leading zeros. (We know this number is not all zeroes though)
+			if (s[0] == '0')
+			{
+				int nonZeroPos = 0;
+				while (s[++nonZeroPos] == '0') {}
+
+				s = s.Substring(nonZeroPos);
+			}
+
+			if (sign == -1)
+			{
+				s = "-" + s;
+			}
+
+			return 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);
+
+//			// store value into a byte array
+//			byte[] b = new byte[8];
+//			for (int i = 0; i < 8; i++)
+//			{
+//				b[7 - i] = (byte)value;
+//				value >>= 8;
+//			}
+//
+//			return new BigInteger(b);
+		}
+
+		public static BigInteger ValueOf(
+			long value)
+		{
+			switch (value)
+			{
+				case 0:
+					return Zero;
+				case 1:
+					return One;
+				case 2:
+					return Two;
+				case 3:
+					return Three;
+				case 10:
+					return Ten;
+			}
+
+			return createValueOf(value);
+		}
+
+		public int GetLowestSetBit()
+		{
+			if (this.sign == 0)
+				return -1;
+
+			int w = magnitude.Length;
+
+			while (--w > 0)
+			{
+				if (magnitude[w] != 0)
+					break;
+			}
+
+			int word = (int) magnitude[w];
+			Debug.Assert(word != 0);
+
+			int b = (word & 0x0000FFFF) == 0
+				?	(word & 0x00FF0000) == 0
+					?	7
+					:	15
+				:	(word & 0x000000FF) == 0
+					?	23
+					:	31;
+
+			while (b > 0)
+			{
+				if ((word << b) == int.MinValue)
+					break;
+
+				b--;
+			}
+
+			return ((magnitude.Length - w) * 32 - (b + 1));
+		}
+
+		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..4dd5e74e2
--- /dev/null
+++ b/Crypto/src/math/ec/ECCurve.cs
@@ -0,0 +1,661 @@
+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 DecodePoint(byte[] encoded);
+		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();
+		}
+	}
+
+	public abstract class ECCurveBase : ECCurve
+	{
+		protected internal ECCurveBase()
+		{
+		}
+
+		protected internal 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 override 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, encoded.Length - 1);
+
+					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 : ECCurveBase
+    {
+        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 internal 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 : ECCurveBase
+    {
+        /**
+         * 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>&#956;</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>&#956;</code> of the elliptic curve.
+		 * @return <code>&#956;</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 internal 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..f95f09974
--- /dev/null
+++ b/Crypto/src/math/ec/ECPoint.cs
@@ -0,0 +1,567 @@
+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 abstract byte[] GetEncoded();
+
+		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()
+		{
+			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 (withCompression)
+			{
+				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) || (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);
+		}
+
+		// 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
+	* &#964;-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>&#945;<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>&#945;<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>&#945;<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>&#945;<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>&#945;<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>&#955;</code> of
+		* <code><b>Z</b>[&#964;]</code>.
+		* @param mu The parameter <code>&#956;</code> of the elliptic curve.
+		* @param lambda The element <code>&#955;</code> of
+		* <code><b>Z</b>[&#964;]</code>.
+		* @return The norm of <code>&#955;</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>&#955;</code> of
+		* <code><b>R</b>[&#964;]</code>, where <code>&#955; = u + v&#964;</code>
+		* and <code>u</code> and <code>u</code> are real numbers (elements of
+		* <code><b>R</b></code>). 
+		* @param mu The parameter <code>&#956;</code> of the elliptic curve.
+		* @param u The real part of the element <code>&#955;</code> of
+		* <code><b>R</b>[&#964;]</code>.
+		* @param v The <code>&#964;</code>-adic part of the element
+		* <code>&#955;</code> of <code><b>R</b>[&#964;]</code>.
+		* @return The norm of <code>&#955;</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>&#955;</code> of <code><b>R</b>[&#964;]</code>
+		* to an element of <code><b>Z</b>[&#964;]</code>, such that their difference
+		* has minimal norm. <code>&#955;</code> is given as
+		* <code>&#955; = &#955;<sub>0</sub> + &#955;<sub>1</sub>&#964;</code>.
+		* @param lambda0 The component <code>&#955;<sub>0</sub></code>.
+		* @param lambda1 The component <code>&#955;<sub>1</sub></code>.
+		* @param mu The parameter <code>&#956;</code> of the elliptic curve. Must
+		* equal 1 or -1.
+		* @return The rounded element of <code><b>Z</b>[&#964;]</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>&#955; = 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>&#955; = 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>&#964;</code>-adic NAF (non-adjacent form) of an
+		* element <code>&#955;</code> of <code><b>Z</b>[&#964;]</code>.
+		* @param mu The parameter <code>&#956;</code> of the elliptic curve.
+		* @param lambda The element <code>&#955;</code> of
+		* <code><b>Z</b>[&#964;]</code>.
+		* @return The <code>&#964;</code>-adic NAF of <code>&#955;</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>&#964;()</code> to an
+		* <code>F2mPoint</code>. 
+		* @param p The F2mPoint to which <code>&#964;()</code> is applied.
+		* @return <code>&#964;(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>&#956;</code> of the elliptic curve.
+		* @param curve The elliptic curve from which to obtain <code>&#956;</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>&#956;</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>&#956;</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>&#956;</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>(&#964;<sup>m</sup> - 1)/(&#964; - 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 &#956; of the elliptic curve.
+		* @param c The precision (number of bits of accuracy) of the partial
+		* modular reduction.
+		* @return <code>&#961; := k partmod (&#964;<sup>m</sup> - 1)/(&#964; - 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>&#964;</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>&#955;</code> of <code><b>Z</b>[&#964;]</code>
+		* using the <code>&#964;</code>-adic NAF (TNAF) method.
+		* @param p The F2mPoint to Multiply.
+		* @param lambda The element <code>&#955;</code> of
+		* <code><b>Z</b>[&#964;]</code>.
+		* @return <code>&#955; * 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>&#955;</code> of <code><b>Z</b>[&#964;]</code>
+		* using the <code>&#964;</code>-adic NAF (TNAF) method, given the TNAF
+		* of <code>&#955;</code>.
+		* @param p The F2mPoint to Multiply.
+		* @param u The the TNAF of <code>&#955;</code>..
+		* @return <code>&#955; * 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>[&#964;]</code>-adic window NAF of an element
+		* <code>&#955;</code> of <code><b>Z</b>[&#964;]</code>.
+		* @param mu The parameter &#956; of the elliptic curve.
+		* @param lambda The element <code>&#955;</code> of
+		* <code><b>Z</b>[&#964;]</code> of which to compute the
+		* <code>[&#964;]</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>&#945;<sub>u</sub></code>'s for the window width.
+		* @return The <code>[&#964;]</code>-adic window NAF of
+		* <code>&#955;</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>[&#964;]</code>. Let
+	* <code>&#955;</code> be an element of <code><b>Z</b>[&#964;]</code>. Then
+	* <code>&#955;</code> is given as <code>&#955; = u + v&#964;</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 &quot;real&quot; part of <code>&#955;</code>.
+		*/
+		public readonly BigInteger u;
+
+		/**
+		* The &quot;<code>&#964;</code>-adic&quot; part of <code>&#955;</code>.
+		*/
+		public readonly BigInteger v;
+
+		/**
+		* Constructor for an element <code>&#955;</code> of
+		* <code><b>Z</b>[&#964;]</code>.
+		* @param u The &quot;real&quot; part of <code>&#955;</code>.
+		* @param v The &quot;<code>&#964;</code>-adic&quot; part of
+		* <code>&#955;</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 = &#8722;<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>&#964;</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>&#964;</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>&#955;</code> of <code><b>Z</b>[&#964;]</code> using
+		* the <code>&#964;</code>-adic NAF (TNAF) method.
+		* @param p The F2mPoint to multiply.
+		* @param lambda The element <code>&#955;</code> of
+		* <code><b>Z</b>[&#964;]</code> of which to compute the
+		* <code>[&#964;]</code>-adic NAF.
+		* @return <code>p</code> multiplied by <code>&#955;</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>&#955;</code> of <code><b>Z</b>[&#964;]</code>
+		* using the window <code>&#964;</code>-adic NAF (TNAF) method, given the
+		* WTNAF of <code>&#955;</code>.
+		* @param p The F2mPoint to multiply.
+		* @param u The the WTNAF of <code>&#955;</code>..
+		* @return <code>&#955; * 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>&#964;</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;
+		}
+	}
+}
generated by cgit-pink 1.4.1 (git 2.36.1) at 2024-11-10 03:12:59 +0000