Initial import of old CVS repository

17 files changed, 8286 insertions, 0 deletions

diff --git a/crypto/src/math/BigInteger.cs b/crypto/src/math/BigInteger.cs
new file mode 100644
index 000000000..00f8f399d
--- /dev/null
+++ b/crypto/src/math/BigInteger.cs
@@ -0,0 +1,3575 @@
+using System;
+using System.Collections;
+using System.Diagnostics;
+using System.Globalization;
+using System.Text;
+
+using Org.BouncyCastle.Utilities;
+
+namespace Org.BouncyCastle.Math
+{
+#if !(NETCF_1_0 || NETCF_2_0 || SILVERLIGHT)
+	[Serializable]
+#endif
+	public class BigInteger
+	{
+		// The first few odd primes
+        /*
+                3   5   7   11  13  17  19  23  29
+            31  37  41  43  47  53  59  61  67  71
+            73  79  83  89  97  101 103 107 109 113
+            127 131 137 139 149 151 157 163 167 173
+            179 181 191 193 197 199 211 223 227 229
+            233 239 241 251 257 263 269 271 277 281
+            283 293 307 311 313 317 331 337 347 349
+            353 359 367 373 379 383 389 397 401 409
+            419 421 431 433 439 443 449 457 461 463
+            467 479 487 491 499 503 509 521 523 541
+            547 557 563 569 571 577 587 593 599 601
+            607 613 617 619 631 641 643 647 653 659
+            661 673 677 683 691 701 709 719 727 733
+            739 743 751 757 761 769 773 787 797 809
+            811 821 823 827 829 839 853 857 859 863
+            877 881 883 887 907 911 919 929 937 941
+            947 953 967 971 977 983 991 997 1009
+            1013 1019 1021 1031 1033 1039 1049 1051
+            1061 1063 1069 1087 1091 1093 1097 1103
+            1109 1117 1123 1129 1151 1153 1163 1171
+            1181 1187 1193 1201 1213 1217 1223 1229
+            1231 1237 1249 1259 1277 1279 1283 1289
+        */
+
+        // Each list has a product < 2^31
+		internal static readonly int[][] primeLists = new int[][]
+		{
+			new int[]{ 3, 5, 7, 11, 13, 17, 19, 23 },
+			new int[]{ 29, 31, 37, 41, 43 },
+			new int[]{ 47, 53, 59, 61, 67 },
+			new int[]{ 71, 73, 79, 83 },
+			new int[]{ 89, 97, 101, 103 },
+
+			new int[]{ 107, 109, 113, 127 },
+			new int[]{ 131, 137, 139, 149 },
+			new int[]{ 151, 157, 163, 167 },
+			new int[]{ 173, 179, 181, 191 },
+			new int[]{ 193, 197, 199, 211 },
+
+			new int[]{ 223, 227, 229 },
+			new int[]{ 233, 239, 241 },
+			new int[]{ 251, 257, 263 },
+			new int[]{ 269, 271, 277 },
+			new int[]{ 281, 283, 293 },
+
+			new int[]{ 307, 311, 313 },
+			new int[]{ 317, 331, 337 },
+			new int[]{ 347, 349, 353 },
+			new int[]{ 359, 367, 373 },
+			new int[]{ 379, 383, 389 },
+
+			new int[]{ 397, 401, 409 },
+			new int[]{ 419, 421, 431 },
+			new int[]{ 433, 439, 443 },
+			new int[]{ 449, 457, 461 },
+			new int[]{ 463, 467, 479 },
+
+			new int[]{ 487, 491, 499 },
+			new int[]{ 503, 509, 521 },
+			new int[]{ 523, 541, 547 },
+			new int[]{ 557, 563, 569 },
+			new int[]{ 571, 577, 587 },
+
+			new int[]{ 593, 599, 601 },
+			new int[]{ 607, 613, 617 },
+			new int[]{ 619, 631, 641 },
+			new int[]{ 643, 647, 653 },
+			new int[]{ 659, 661, 673 },
+
+			new int[]{ 677, 683, 691 },
+			new int[]{ 701, 709, 719 },
+			new int[]{ 727, 733, 739 },
+			new int[]{ 743, 751, 757 },
+			new int[]{ 761, 769, 773 },
+
+			new int[]{ 787, 797, 809 },
+			new int[]{ 811, 821, 823 },
+			new int[]{ 827, 829, 839 },
+			new int[]{ 853, 857, 859 },
+			new int[]{ 863, 877, 881 },
+
+			new int[]{ 883, 887, 907 },
+			new int[]{ 911, 919, 929 },
+			new int[]{ 937, 941, 947 },
+			new int[]{ 953, 967, 971 },
+			new int[]{ 977, 983, 991 },
+
+            new int[]{ 997, 1009, 1013 },
+            new int[]{ 1019, 1021, 1031 },
+            new int[]{ 1033, 1039, 1049 },
+            new int[]{ 1051, 1061, 1063 },
+            new int[]{ 1069, 1087, 1091 },
+
+            new int[]{ 1093, 1097, 1103 },
+            new int[]{ 1109, 1117, 1123 },
+            new int[]{ 1129, 1151, 1153 },
+            new int[]{ 1163, 1171, 1181 },
+            new int[]{ 1187, 1193, 1201 },
+
+            new int[]{ 1213, 1217, 1223 },
+            new int[]{ 1229, 1231, 1237 },
+            new int[]{ 1249, 1259, 1277 },
+            new int[]{ 1279, 1283, 1289 },
+        };
+
+        internal static readonly int[] primeProducts;
+
+		private const long IMASK = 0xFFFFFFFFL;
+        private const ulong UIMASK = 0xFFFFFFFFUL;
+
+        private static readonly int[] ZeroMagnitude = new int[0];
+		private static readonly byte[] ZeroEncoding = new byte[0];
+
+        private static readonly BigInteger[] SMALL_CONSTANTS = new BigInteger[17];
+		public static readonly BigInteger Zero;
+		public static readonly BigInteger One;
+        public static readonly BigInteger Two;
+        public static readonly BigInteger Three;
+        public static readonly BigInteger Ten;
+
+        //private readonly static byte[] BitCountTable =
+        //{
+        //    0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
+        //    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
+        //    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
+        //    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
+        //    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
+        //    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
+        //    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
+        //    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
+        //    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
+        //    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
+        //    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
+        //    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
+        //    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
+        //    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
+        //    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
+        //    4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
+        //};
+
+        private readonly static byte[] BitLengthTable =
+		{
+			0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,
+            5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
+            6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
+            6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
+            7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+            7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+            7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+            7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+            8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+            8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+            8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+            8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+            8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+            8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+            8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+            8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8
+		};
+
+        // TODO Parse radix-2 64 bits at a time and radix-8 63 bits at a time
+        private const int chunk2 = 1, chunk8 = 1, chunk10 = 19, chunk16 = 16;
+		private static readonly BigInteger radix2, radix2E, radix8, radix8E, radix10, radix10E, radix16, radix16E;
+
+        private static readonly Random RandomSource = new Random();
+
+        /*
+         * These are the threshold bit-lengths (of an exponent) where we increase the window size.
+         * They are calculated according to the expected savings in multiplications.
+         * Some squares will also be saved on average, but we offset these against the extra storage costs.
+         */
+        private static readonly int[] ExpWindowThresholds = { 7, 25, 81, 241, 673, 1793, 4609, Int32.MaxValue };
+
+        private const int BitsPerByte = 8;
+		private const int BitsPerInt = 32;
+		private const int BytesPerInt = 4;
+
+		static BigInteger()
+		{
+            Zero = new BigInteger(0, ZeroMagnitude, false);
+            Zero.nBits = 0; Zero.nBitLength = 0;
+
+            SMALL_CONSTANTS[0] = Zero;
+            for (uint i = 1; i < SMALL_CONSTANTS.Length; ++i)
+            {
+                SMALL_CONSTANTS[i] = CreateUValueOf(i);
+            }
+
+            One = SMALL_CONSTANTS[1];
+            Two = SMALL_CONSTANTS[2];
+            Three = SMALL_CONSTANTS[3];
+            Ten = SMALL_CONSTANTS[10];
+
+            radix2 = ValueOf(2);
+            radix2E = radix2.Pow(chunk2);
+
+            radix8 = ValueOf(8);
+            radix8E = radix8.Pow(chunk8);
+
+            radix10 = ValueOf(10);
+		    radix10E = radix10.Pow(chunk10);
+
+            radix16 = ValueOf(16);
+            radix16E = radix16.Pow(chunk16);
+
+            primeProducts = new int[primeLists.Length];
+
+			for (int i = 0; i < primeLists.Length; ++i)
+			{
+				int[] primeList = primeLists[i];
+				int product = primeList[0];
+				for (int j = 1; j < primeList.Length; ++j)
+				{
+					product *= primeList[j];
+				}
+				primeProducts[i] = product;
+			}
+		}
+
+		private int[] magnitude; // array of ints with [0] being the most significant
+        private int sign; // -1 means -ve; +1 means +ve; 0 means 0;
+		private int nBits = -1; // cache BitCount() value
+		private int nBitLength = -1; // cache BitLength() value
+        private int mQuote = 0; // -m^(-1) mod b, b = 2^32 (see Montgomery mult.), 0 when uninitialised
+
+        private static int GetByteLength(
+			int nBits)
+		{
+			return (nBits + BitsPerByte - 1) / BitsPerByte;
+		}
+
+		private BigInteger(
+			int		signum,
+			int[]	mag,
+			bool	checkMag)
+		{
+			if (checkMag)
+			{
+				int i = 0;
+				while (i < mag.Length && mag[i] == 0)
+				{
+					++i;
+				}
+
+				if (i == mag.Length)
+				{
+					this.sign = 0;
+					this.magnitude = ZeroMagnitude;
+				}
+				else
+				{
+					this.sign = signum;
+
+					if (i == 0)
+					{
+						this.magnitude = mag;
+					}
+					else
+					{
+						// strip leading 0 words
+						this.magnitude = new int[mag.Length - i];
+						Array.Copy(mag, i, this.magnitude, 0, this.magnitude.Length);
+					}
+				}
+			}
+			else
+			{
+				this.sign = signum;
+				this.magnitude = mag;
+			}
+		}
+
+		public BigInteger(
+			string value)
+			: this(value, 10)
+		{
+		}
+
+		public BigInteger(
+			string	str,
+			int		radix)
+		{
+			if (str.Length == 0)
+				throw new FormatException("Zero length BigInteger");
+
+			NumberStyles style;
+			int chunk;
+			BigInteger r;
+			BigInteger rE;
+
+			switch (radix)
+			{
+				case 2:
+					// Is there anyway to restrict to binary digits?
+					style = NumberStyles.Integer;
+					chunk = chunk2;
+					r = radix2;
+					rE = radix2E;
+					break;
+                case 8:
+                    // Is there anyway to restrict to octal digits?
+                    style = NumberStyles.Integer;
+                    chunk = chunk8;
+                    r = radix8;
+                    rE = radix8E;
+                    break;
+                case 10:
+					// This style seems to handle spaces and minus sign already (our processing redundant?)
+					style = NumberStyles.Integer;
+					chunk = chunk10;
+					r = radix10;
+					rE = radix10E;
+					break;
+				case 16:
+					// TODO Should this be HexNumber?
+					style = NumberStyles.AllowHexSpecifier;
+					chunk = chunk16;
+					r = radix16;
+					rE = radix16E;
+					break;
+				default:
+					throw new FormatException("Only bases 2, 8, 10, or 16 allowed");
+			}
+
+
+			int index = 0;
+			sign = 1;
+
+			if (str[0] == '-')
+			{
+				if (str.Length == 1)
+					throw new FormatException("Zero length BigInteger");
+
+				sign = -1;
+				index = 1;
+			}
+
+			// strip leading zeros from the string str
+			while (index < str.Length && Int32.Parse(str[index].ToString(), style) == 0)
+			{
+				index++;
+			}
+
+			if (index >= str.Length)
+			{
+				// zero value - we're done
+				sign = 0;
+				magnitude = ZeroMagnitude;
+				return;
+			}
+
+			//////
+			// could we work out the max number of ints required to store
+			// str.Length digits in the given base, then allocate that
+			// storage in one hit?, then Generate the magnitude in one hit too?
+			//////
+
+			BigInteger b = Zero;
+
+
+			int next = index + chunk;
+
+			if (next <= str.Length)
+			{
+				do
+				{
+					string s = str.Substring(index, chunk);
+					ulong i = ulong.Parse(s, style);
+					BigInteger bi = CreateUValueOf(i);
+
+					switch (radix)
+					{
+						case 2:
+							// TODO Need this because we are parsing in radix 10 above
+							if (i >= 2)
+								throw new FormatException("Bad character in radix 2 string: " + s);
+
+							// TODO Parse 64 bits at a time
+							b = b.ShiftLeft(1);
+							break;
+                        case 8:
+                            // TODO Need this because we are parsing in radix 10 above
+                            if (i >= 8)
+                                throw new FormatException("Bad character in radix 8 string: " + s);
+
+                            // TODO Parse 63 bits at a time
+                            b = b.ShiftLeft(3);
+                            break;
+                        case 16:
+							b = b.ShiftLeft(64);
+							break;
+						default:
+							b = b.Multiply(rE);
+							break;
+					}
+
+					b = b.Add(bi);
+
+					index = next;
+					next += chunk;
+				}
+				while (next <= str.Length);
+			}
+
+			if (index < str.Length)
+			{
+				string s = str.Substring(index);
+				ulong i = ulong.Parse(s, style);
+				BigInteger bi = CreateUValueOf(i);
+
+				if (b.sign > 0)
+				{
+					if (radix == 2)
+					{
+						// NB: Can't reach here since we are parsing one char at a time
+						Debug.Assert(false);
+
+						// TODO Parse all bits at once
+//						b = b.ShiftLeft(s.Length);
+					}
+                    else if (radix == 8)
+                    {
+                        // NB: Can't reach here since we are parsing one char at a time
+                        Debug.Assert(false);
+
+                        // TODO Parse all bits at once
+//						b = b.ShiftLeft(s.Length * 3);
+                    }
+                    else if (radix == 16)
+					{
+						b = b.ShiftLeft(s.Length << 2);
+					}
+					else
+					{
+						b = b.Multiply(r.Pow(s.Length));
+					}
+
+					b = b.Add(bi);
+				}
+				else
+				{
+					b = bi;
+				}
+			}
+
+            // Note: This is the previous (slower) algorithm
+//			while (index < value.Length)
+//            {
+//				char c = value[index];
+//				string s = c.ToString();
+//				int i = Int32.Parse(s, style);
+//
+//                b = b.Multiply(r).Add(ValueOf(i));
+//                index++;
+//            }
+
+			magnitude = b.magnitude;
+		}
+
+		public BigInteger(
+			byte[] bytes)
+			: this(bytes, 0, bytes.Length)
+		{
+		}
+
+		public BigInteger(
+			byte[]	bytes,
+			int		offset,
+			int		length)
+		{
+			if (length == 0)
+				throw new FormatException("Zero length BigInteger");
+
+			// TODO Move this processing into MakeMagnitude (provide sign argument)
+			if ((sbyte)bytes[offset] < 0)
+			{
+				this.sign = -1;
+
+				int end = offset + length;
+
+				int iBval;
+				// strip leading sign bytes
+				for (iBval = offset; iBval < end && ((sbyte)bytes[iBval] == -1); iBval++)
+				{
+				}
+
+				if (iBval >= end)
+				{
+					this.magnitude = One.magnitude;
+				}
+				else
+				{
+					int numBytes = end - iBval;
+					byte[] inverse = new byte[numBytes];
+
+					int index = 0;
+					while (index < numBytes)
+					{
+						inverse[index++] = (byte)~bytes[iBval++];
+					}
+
+					Debug.Assert(iBval == end);
+
+					while (inverse[--index] == byte.MaxValue)
+					{
+						inverse[index] = byte.MinValue;
+					}
+
+					inverse[index]++;
+
+					this.magnitude = MakeMagnitude(inverse, 0, inverse.Length);
+				}
+			}
+			else
+			{
+				// strip leading zero bytes and return magnitude bytes
+				this.magnitude = MakeMagnitude(bytes, offset, length);
+				this.sign = this.magnitude.Length > 0 ? 1 : 0;
+			}
+		}
+
+		private static int[] MakeMagnitude(
+			byte[]	bytes,
+			int		offset,
+			int		length)
+		{
+			int end = offset + length;
+
+			// strip leading zeros
+			int firstSignificant;
+			for (firstSignificant = offset; firstSignificant < end
+				&& bytes[firstSignificant] == 0; firstSignificant++)
+			{
+			}
+
+			if (firstSignificant >= end)
+			{
+				return ZeroMagnitude;
+			}
+
+			int nInts = (end - firstSignificant + 3) / BytesPerInt;
+			int bCount = (end - firstSignificant) % BytesPerInt;
+			if (bCount == 0)
+			{
+				bCount = BytesPerInt;
+			}
+
+			if (nInts < 1)
+			{
+				return ZeroMagnitude;
+			}
+
+			int[] mag = new int[nInts];
+
+			int v = 0;
+			int magnitudeIndex = 0;
+			for (int i = firstSignificant; i < end; ++i)
+			{
+				v <<= 8;
+				v |= bytes[i] & 0xff;
+				bCount--;
+				if (bCount <= 0)
+				{
+					mag[magnitudeIndex] = v;
+					magnitudeIndex++;
+					bCount = BytesPerInt;
+					v = 0;
+				}
+			}
+
+			if (magnitudeIndex < mag.Length)
+			{
+				mag[magnitudeIndex] = v;
+			}
+
+			return mag;
+		}
+
+		public BigInteger(
+			int		sign,
+			byte[]	bytes)
+			: this(sign, bytes, 0, bytes.Length)
+		{
+		}
+
+		public BigInteger(
+			int		sign,
+			byte[]	bytes,
+			int		offset,
+			int		length)
+		{
+			if (sign < -1 || sign > 1)
+				throw new FormatException("Invalid sign value");
+
+			if (sign == 0)
+			{
+				this.sign = 0;
+				this.magnitude = ZeroMagnitude;
+			}
+			else
+			{
+				// copy bytes
+				this.magnitude = MakeMagnitude(bytes, offset, length);
+				this.sign = this.magnitude.Length < 1 ? 0 : sign;
+			}
+		}
+
+		public BigInteger(
+			int		sizeInBits,
+			Random	random)
+		{
+			if (sizeInBits < 0)
+				throw new ArgumentException("sizeInBits must be non-negative");
+
+			this.nBits = -1;
+			this.nBitLength = -1;
+
+			if (sizeInBits == 0)
+			{
+				this.sign = 0;
+				this.magnitude = ZeroMagnitude;
+				return;
+			}
+
+			int nBytes = GetByteLength(sizeInBits);
+			byte[] b = new byte[nBytes];
+			random.NextBytes(b);
+
+			// strip off any excess bits in the MSB
+            int xBits = BitsPerByte * nBytes - sizeInBits;
+            b[0] &= (byte)(255U >> xBits);
+
+            this.magnitude = MakeMagnitude(b, 0, b.Length);
+			this.sign = this.magnitude.Length < 1 ? 0 : 1;
+		}
+
+        public BigInteger(
+			int		bitLength,
+			int		certainty,
+			Random	random)
+		{
+			if (bitLength < 2)
+				throw new ArithmeticException("bitLength < 2");
+
+			this.sign = 1;
+			this.nBitLength = bitLength;
+
+			if (bitLength == 2)
+			{
+				this.magnitude = random.Next(2) == 0
+					?	Two.magnitude
+					:	Three.magnitude;
+				return;
+			}
+             
+			int nBytes = GetByteLength(bitLength);
+			byte[] b = new byte[nBytes];
+
+			int xBits = BitsPerByte * nBytes - bitLength;
+            byte mask = (byte)(255U >> xBits);
+
+            for (;;)
+			{
+				random.NextBytes(b);
+
+				// strip off any excess bits in the MSB
+				b[0] &= mask;
+
+				// ensure the leading bit is 1 (to meet the strength requirement)
+				b[0] |= (byte)(1 << (7 - xBits));
+
+				// ensure the trailing bit is 1 (i.e. must be odd)
+				b[nBytes - 1] |= 1;
+
+				this.magnitude = MakeMagnitude(b, 0, b.Length);
+				this.nBits = -1;
+				this.mQuote = 0;
+
+				if (certainty < 1)
+					break;
+
+				if (CheckProbablePrime(certainty, random))
+					break;
+
+				if (bitLength > 32)
+				{
+					for (int rep = 0; rep < 10000; ++rep)
+					{
+						int n = 33 + random.Next(bitLength - 2);
+						this.magnitude[this.magnitude.Length - (n >> 5)] ^= (1 << (n & 31));
+						this.magnitude[this.magnitude.Length - 1] ^= ((random.Next() + 1) << 1);
+						this.mQuote = 0;
+
+						if (CheckProbablePrime(certainty, random))
+							return;
+					}
+				}
+			}
+		}
+
+		public BigInteger Abs()
+		{
+			return sign >= 0 ? this : Negate();
+		}
+
+		/**
+		 * return a = a + b - b preserved.
+		 */
+		private static int[] AddMagnitudes(
+			int[] a,
+			int[] b)
+		{
+			int tI = a.Length - 1;
+			int vI = b.Length - 1;
+			long m = 0;
+
+			while (vI >= 0)
+			{
+				m += ((long)(uint)a[tI] + (long)(uint)b[vI--]);
+				a[tI--] = (int)m;
+				m = (long)((ulong)m >> 32);
+			}
+
+			if (m != 0)
+			{
+				while (tI >= 0 && ++a[tI--] == 0)
+				{
+				}
+			}
+
+			return a;
+		}
+
+		public BigInteger Add(
+			BigInteger value)
+		{
+			if (this.sign == 0)
+				return value;
+
+			if (this.sign != value.sign)
+			{
+				if (value.sign == 0)
+					return this;
+
+				if (value.sign < 0)
+					return Subtract(value.Negate());
+
+				return value.Subtract(Negate());
+			}
+
+			return AddToMagnitude(value.magnitude);
+		}
+
+		private BigInteger AddToMagnitude(
+			int[] magToAdd)
+		{
+			int[] big, small;
+			if (this.magnitude.Length < magToAdd.Length)
+			{
+				big = magToAdd;
+				small = this.magnitude;
+			}
+			else
+			{
+				big = this.magnitude;
+				small = magToAdd;
+			}
+
+			// Conservatively avoid over-allocation when no overflow possible
+			uint limit = uint.MaxValue;
+			if (big.Length == small.Length)
+				limit -= (uint) small[0];
+
+			bool possibleOverflow = (uint) big[0] >= limit;
+
+			int[] bigCopy;
+			if (possibleOverflow)
+			{
+				bigCopy = new int[big.Length + 1];
+				big.CopyTo(bigCopy, 1);
+			}
+			else
+			{
+				bigCopy = (int[]) big.Clone();
+			}
+
+			bigCopy = AddMagnitudes(bigCopy, small);
+
+			return new BigInteger(this.sign, bigCopy, possibleOverflow);
+		}
+
+		public BigInteger And(
+			BigInteger value)
+		{
+			if (this.sign == 0 || value.sign == 0)
+			{
+				return Zero;
+			}
+
+			int[] aMag = this.sign > 0
+				? this.magnitude
+				: Add(One).magnitude;
+
+			int[] bMag = value.sign > 0
+				? value.magnitude
+				: value.Add(One).magnitude;
+
+			bool resultNeg = sign < 0 && value.sign < 0;
+			int resultLength = System.Math.Max(aMag.Length, bMag.Length);
+			int[] resultMag = new int[resultLength];
+
+			int aStart = resultMag.Length - aMag.Length;
+			int bStart = resultMag.Length - bMag.Length;
+
+			for (int i = 0; i < resultMag.Length; ++i)
+			{
+				int aWord = i >= aStart ? aMag[i - aStart] : 0;
+				int bWord = i >= bStart ? bMag[i - bStart] : 0;
+
+				if (this.sign < 0)
+				{
+					aWord = ~aWord;
+				}
+
+				if (value.sign < 0)
+				{
+					bWord = ~bWord;
+				}
+
+				resultMag[i] = aWord & bWord;
+
+				if (resultNeg)
+				{
+					resultMag[i] = ~resultMag[i];
+				}
+			}
+
+			BigInteger result = new BigInteger(1, resultMag, true);
+
+			// TODO Optimise this case
+			if (resultNeg)
+			{
+				result = result.Not();
+			}
+
+			return result;
+		}
+
+		public BigInteger AndNot(
+			BigInteger val)
+		{
+			return And(val.Not());
+		}
+
+		public int BitCount
+		{
+			get
+			{
+				if (nBits == -1)
+				{
+					if (sign < 0)
+					{
+						// TODO Optimise this case
+						nBits = Not().BitCount;
+					}
+					else
+					{
+						int sum = 0;
+						for (int i = 0; i < magnitude.Length; ++i)
+						{
+                            sum += BitCnt(magnitude[i]);
+						}
+						nBits = sum;
+					}
+				}
+
+				return nBits;
+			}
+		}
+
+        public static int BitCnt(int i)
+        {
+            uint u = (uint)i;
+            u = u - ((u >> 1) & 0x55555555);
+            u = (u & 0x33333333) + ((u >> 2) & 0x33333333);
+            u = (u + (u >> 4)) & 0x0f0f0f0f;
+            u += (u >> 8);
+            u += (u >> 16);
+            u &= 0x3f;
+            return (int)u;
+        }
+
+        private static int CalcBitLength(int sign, int indx, int[]	mag)
+		{
+			for (;;)
+			{
+				if (indx >= mag.Length)
+					return 0;
+
+                if (mag[indx] != 0)
+					break;
+
+				++indx;
+			}
+
+			// bit length for everything after the first int
+			int bitLength = 32 * ((mag.Length - indx) - 1);
+
+			// and determine bitlength of first int
+			int firstMag = mag[indx];
+			bitLength += BitLen(firstMag);
+
+			// Check for negative powers of two
+			if (sign < 0 && ((firstMag & -firstMag) == firstMag))
+			{
+				do
+				{
+					if (++indx >= mag.Length)
+					{
+						--bitLength;
+						break;
+					}
+				}
+				while (mag[indx] == 0);
+			}
+
+			return bitLength;
+		}
+
+		public int BitLength
+		{
+			get
+			{
+				if (nBitLength == -1)
+				{
+					nBitLength = sign == 0
+						? 0
+						: CalcBitLength(sign, 0, magnitude);
+				}
+
+				return nBitLength;
+			}
+		}
+
+		//
+		// BitLen(value) is the number of bits in value.
+		//
+		private static int BitLen(int w)
+		{
+            uint v = (uint)w;
+            uint t = v >> 24;
+            if (t != 0)
+                return 24 + BitLengthTable[t];
+            t = v >> 16;
+            if (t != 0)
+                return 16 + BitLengthTable[t];
+            t = v >> 8;
+            if (t != 0)
+                return 8 + BitLengthTable[t];
+            return BitLengthTable[v];
+		}
+
+        private bool QuickPow2Check()
+		{
+			return sign > 0 && nBits == 1;
+		}
+
+		public int CompareTo(
+			object obj)
+		{
+			return CompareTo((BigInteger)obj);
+		}
+
+		/**
+		 * unsigned comparison on two arrays - note the arrays may
+		 * start with leading zeros.
+		 */
+		private static int CompareTo(
+			int		xIndx,
+			int[]	x,
+			int		yIndx,
+			int[]	y)
+		{
+			while (xIndx != x.Length && x[xIndx] == 0)
+			{
+				xIndx++;
+			}
+
+			while (yIndx != y.Length && y[yIndx] == 0)
+			{
+				yIndx++;
+			}
+
+			return CompareNoLeadingZeroes(xIndx, x, yIndx, y);
+		}
+
+		private static int CompareNoLeadingZeroes(
+			int		xIndx,
+			int[]	x,
+			int		yIndx,
+			int[]	y)
+		{
+			int diff = (x.Length - y.Length) - (xIndx - yIndx);
+
+			if (diff != 0)
+			{
+				return diff < 0 ? -1 : 1;
+			}
+
+			// lengths of magnitudes the same, test the magnitude values
+
+			while (xIndx < x.Length)
+			{
+				uint v1 = (uint)x[xIndx++];
+				uint v2 = (uint)y[yIndx++];
+
+				if (v1 != v2)
+					return v1 < v2 ? -1 : 1;
+			}
+
+			return 0;
+		}
+
+		public int CompareTo(
+			BigInteger value)
+		{
+			return sign < value.sign ? -1
+				: sign > value.sign ? 1
+				: sign == 0 ? 0
+				: sign * CompareNoLeadingZeroes(0, magnitude, 0, value.magnitude);
+		}
+
+		/**
+		 * return z = x / y - done in place (z value preserved, x contains the
+		 * remainder)
+		 */
+		private int[] Divide(
+			int[]	x,
+			int[]	y)
+		{
+			int xStart = 0;
+			while (xStart < x.Length && x[xStart] == 0)
+			{
+				++xStart;
+			}
+
+			int yStart = 0;
+			while (yStart < y.Length && y[yStart] == 0)
+			{
+				++yStart;
+			}
+
+			Debug.Assert(yStart < y.Length);
+
+			int xyCmp = CompareNoLeadingZeroes(xStart, x, yStart, y);
+			int[] count;
+
+			if (xyCmp > 0)
+			{
+				int yBitLength = CalcBitLength(1, yStart, y);
+				int xBitLength = CalcBitLength(1, xStart, x);
+				int shift = xBitLength - yBitLength;
+
+				int[] iCount;
+				int iCountStart = 0;
+
+				int[] c;
+				int cStart = 0;
+				int cBitLength = yBitLength;
+				if (shift > 0)
+				{
+//					iCount = ShiftLeft(One.magnitude, shift);
+					iCount = new int[(shift >> 5) + 1];
+					iCount[0] = 1 << (shift % 32);
+
+					c = ShiftLeft(y, shift);
+					cBitLength += shift;
+				}
+				else
+				{
+					iCount = new int[] { 1 };
+
+					int len = y.Length - yStart;
+					c = new int[len];
+					Array.Copy(y, yStart, c, 0, len);
+				}
+
+				count = new int[iCount.Length];
+
+				for (;;)
+				{
+					if (cBitLength < xBitLength
+						|| CompareNoLeadingZeroes(xStart, x, cStart, c) >= 0)
+					{
+						Subtract(xStart, x, cStart, c);
+						AddMagnitudes(count, iCount);
+
+						while (x[xStart] == 0)
+						{
+							if (++xStart == x.Length)
+								return count;
+						}
+
+						//xBitLength = CalcBitLength(xStart, x);
+						xBitLength = 32 * (x.Length - xStart - 1) + BitLen(x[xStart]);
+
+						if (xBitLength <= yBitLength)
+						{
+							if (xBitLength < yBitLength)
+								return count;
+
+							xyCmp = CompareNoLeadingZeroes(xStart, x, yStart, y);
+
+							if (xyCmp <= 0)
+								break;
+						}
+					}
+
+					shift = cBitLength - xBitLength;
+
+					// NB: The case where c[cStart] is 1-bit is harmless
+					if (shift == 1)
+					{
+						uint firstC = (uint) c[cStart] >> 1;
+						uint firstX = (uint) x[xStart];
+						if (firstC > firstX)
+							++shift;
+					}
+
+					if (shift < 2)
+					{
+						ShiftRightOneInPlace(cStart, c);
+						--cBitLength;
+						ShiftRightOneInPlace(iCountStart, iCount);
+					}
+					else
+					{
+						ShiftRightInPlace(cStart, c, shift);
+						cBitLength -= shift;
+						ShiftRightInPlace(iCountStart, iCount, shift);
+					}
+
+					//cStart = c.Length - ((cBitLength + 31) / 32);
+					while (c[cStart] == 0)
+					{
+						++cStart;
+					}
+
+					while (iCount[iCountStart] == 0)
+					{
+						++iCountStart;
+					}
+				}
+			}
+			else
+			{
+				count = new int[1];
+			}
+
+			if (xyCmp == 0)
+			{
+				AddMagnitudes(count, One.magnitude);
+				Array.Clear(x, xStart, x.Length - xStart);
+			}
+
+			return count;
+		}
+
+		public BigInteger Divide(
+			BigInteger val)
+		{
+			if (val.sign == 0)
+				throw new ArithmeticException("Division by zero error");
+
+			if (sign == 0)
+				return Zero;
+
+			if (val.QuickPow2Check()) // val is power of two
+			{
+				BigInteger result = this.Abs().ShiftRight(val.Abs().BitLength - 1);
+				return val.sign == this.sign ? result : result.Negate();
+			}
+
+			int[] mag = (int[]) this.magnitude.Clone();
+
+			return new BigInteger(this.sign * val.sign, Divide(mag, val.magnitude), true);
+		}
+
+		public BigInteger[] DivideAndRemainder(
+			BigInteger val)
+		{
+			if (val.sign == 0)
+				throw new ArithmeticException("Division by zero error");
+
+			BigInteger[] biggies = new BigInteger[2];
+
+			if (sign == 0)
+			{
+				biggies[0] = Zero;
+				biggies[1] = Zero;
+			}
+			else if (val.QuickPow2Check()) // val is power of two
+			{
+				int e = val.Abs().BitLength - 1;
+				BigInteger quotient = this.Abs().ShiftRight(e);
+				int[] remainder = this.LastNBits(e);
+
+				biggies[0] = val.sign == this.sign ? quotient : quotient.Negate();
+				biggies[1] = new BigInteger(this.sign, remainder, true);
+			}
+			else
+			{
+				int[] remainder = (int[]) this.magnitude.Clone();
+				int[] quotient = Divide(remainder, val.magnitude);
+
+				biggies[0] = new BigInteger(this.sign * val.sign, quotient, true);
+				biggies[1] = new BigInteger(this.sign, remainder, true);
+			}
+
+			return biggies;
+		}
+
+		public override bool Equals(
+			object obj)
+		{
+			if (obj == this)
+				return true;
+
+			BigInteger biggie = obj as BigInteger;
+			if (biggie == null)
+				return false;
+
+            return sign == biggie.sign && IsEqualMagnitude(biggie);
+		}
+
+        private bool IsEqualMagnitude(BigInteger x)
+        {
+            int[] xMag = x.magnitude;
+            if (magnitude.Length != x.magnitude.Length)
+                return false;
+            for (int i = 0; i < magnitude.Length; i++)
+            {
+                if (magnitude[i] != x.magnitude[i])
+                    return false;
+            }
+            return true;
+        }
+
+        public BigInteger Gcd(
+			BigInteger value)
+		{
+			if (value.sign == 0)
+				return Abs();
+
+			if (sign == 0)
+				return value.Abs();
+
+			BigInteger r;
+			BigInteger u = this;
+			BigInteger v = value;
+
+			while (v.sign != 0)
+			{
+				r = u.Mod(v);
+				u = v;
+				v = r;
+			}
+
+			return u;
+		}
+
+		public override int GetHashCode()
+		{
+			int hc = magnitude.Length;
+			if (magnitude.Length > 0)
+			{
+				hc ^= magnitude[0];
+
+				if (magnitude.Length > 1)
+				{
+					hc ^= magnitude[magnitude.Length - 1];
+				}
+			}
+
+			return sign < 0 ? ~hc : hc;
+		}
+
+		// TODO Make public?
+		private BigInteger Inc()
+		{
+			if (this.sign == 0)
+				return One;
+
+			if (this.sign < 0)
+				return new BigInteger(-1, doSubBigLil(this.magnitude, One.magnitude), true);
+
+			return AddToMagnitude(One.magnitude);
+		}
+
+		public int IntValue
+		{
+            get
+            {
+                if (sign == 0)
+                    return 0;
+
+                int n = magnitude.Length;
+
+                int v = magnitude[n - 1];
+
+                return sign < 0 ? -v : v;
+            }
+		}
+
+        /**
+		 * return whether or not a BigInteger is probably prime with a
+		 * probability of 1 - (1/2)**certainty.
+		 * <p>From Knuth Vol 2, pg 395.</p>
+		 */
+		public bool IsProbablePrime(
+			int certainty)
+		{
+			if (certainty <= 0)
+				return true;
+
+			BigInteger n = Abs();
+
+			if (!n.TestBit(0))
+				return n.Equals(Two);
+
+			if (n.Equals(One))
+				return false;
+
+			return n.CheckProbablePrime(certainty, RandomSource);
+		}
+
+		private bool CheckProbablePrime(
+			int		certainty,
+			Random	random)
+		{
+			Debug.Assert(certainty > 0);
+			Debug.Assert(CompareTo(Two) > 0);
+			Debug.Assert(TestBit(0));
+
+
+			// Try to reduce the penalty for really small numbers
+			int numLists = System.Math.Min(BitLength - 1, primeLists.Length);
+
+			for (int i = 0; i < numLists; ++i)
+			{
+				int test = Remainder(primeProducts[i]);
+
+				int[] primeList = primeLists[i];
+				for (int j = 0; j < primeList.Length; ++j)
+				{
+					int prime = primeList[j];
+					int qRem = test % prime;
+					if (qRem == 0)
+					{
+						// We may find small numbers in the list
+						return BitLength < 16 && IntValue == prime;
+					}
+				}
+			}
+
+
+			// TODO Special case for < 10^16 (RabinMiller fixed list)
+//			if (BitLength < 30)
+//			{
+//				RabinMiller against 2, 3, 5, 7, 11, 13, 23 is sufficient
+//			}
+
+
+			// TODO Is it worth trying to create a hybrid of these two?
+			return RabinMillerTest(certainty, random);
+//			return SolovayStrassenTest(certainty, random);
+
+//			bool rbTest = RabinMillerTest(certainty, random);
+//			bool ssTest = SolovayStrassenTest(certainty, random);
+//
+//			Debug.Assert(rbTest == ssTest);
+//
+//			return rbTest;
+		}
+
+		public bool RabinMillerTest(int certainty, Random random)
+		{
+			Debug.Assert(certainty > 0);
+			Debug.Assert(BitLength > 2);
+			Debug.Assert(TestBit(0));
+
+			// let n = 1 + d . 2^s
+            BigInteger n = this;
+            int s = n.GetLowestSetBitMaskFirst(-1 << 1);
+            Debug.Assert(s >= 1);
+            BigInteger r = n.ShiftRight(s);
+
+            // NOTE: Avoid conversion to/from Montgomery form and check for R/-R as result instead
+
+            BigInteger montRadix = One.ShiftLeft(32 * n.magnitude.Length).Remainder(n);
+            BigInteger minusMontRadix = n.Subtract(montRadix);
+
+            do
+			{
+				BigInteger a;
+                do
+                {
+                    a = new BigInteger(n.BitLength, random);
+                }
+                while (a.sign == 0 || a.CompareTo(n) >= 0
+                    || a.IsEqualMagnitude(montRadix) || a.IsEqualMagnitude(minusMontRadix));
+
+				BigInteger y = ModPowMonty(a, r, n, false);
+
+                if (!y.Equals(montRadix))
+                {
+                    int j = 0;
+                    while (!y.Equals(minusMontRadix))
+					{
+						if (++j == s)
+							return false;
+
+						y = ModPowMonty(y, Two, n, false);
+
+                        if (y.Equals(montRadix))
+							return false;
+					}
+				}
+
+				certainty -= 2; // composites pass for only 1/4 possible 'a'
+			}
+			while (certainty > 0);
+
+			return true;
+		}
+
+//		private bool SolovayStrassenTest(
+//			int		certainty,
+//			Random	random)
+//		{
+//			Debug.Assert(certainty > 0);
+//			Debug.Assert(CompareTo(Two) > 0);
+//			Debug.Assert(TestBit(0));
+//
+//			BigInteger n = this;
+//			BigInteger nMinusOne = n.Subtract(One);
+//			BigInteger e = nMinusOne.ShiftRight(1);
+//
+//			do
+//			{
+//				BigInteger a;
+//				do
+//				{
+//					a = new BigInteger(nBitLength, random);
+//				}
+//				// NB: Spec says 0 < x < n, but 1 is trivial
+//				while (a.CompareTo(One) <= 0 || a.CompareTo(n) >= 0);
+//
+//
+//				// TODO Check this is redundant given the way Jacobi() works?
+////				if (!a.Gcd(n).Equals(One))
+////					return false;
+//
+//				int x = Jacobi(a, n);
+//
+//				if (x == 0)
+//					return false;
+//
+//				BigInteger check = a.ModPow(e, n);
+//
+//				if (x == 1 && !check.Equals(One))
+//					return false;
+//
+//				if (x == -1 && !check.Equals(nMinusOne))
+//					return false;
+//
+//				--certainty;
+//			}
+//			while (certainty > 0);
+//
+//			return true;
+//		}
+//
+//		private static int Jacobi(
+//			BigInteger	a,
+//			BigInteger	b)
+//		{
+//			Debug.Assert(a.sign >= 0);
+//			Debug.Assert(b.sign > 0);
+//			Debug.Assert(b.TestBit(0));
+//			Debug.Assert(a.CompareTo(b) < 0);
+//
+//			int totalS = 1;
+//			for (;;)
+//			{
+//				if (a.sign == 0)
+//					return 0;
+//
+//				if (a.Equals(One))
+//					break;
+//
+//				int e = a.GetLowestSetBit();
+//
+//				int bLsw = b.magnitude[b.magnitude.Length - 1];
+//				if ((e & 1) != 0 && ((bLsw & 7) == 3 || (bLsw & 7) == 5))
+//					totalS = -totalS;
+//
+//				// TODO Confirm this is faster than later a1.Equals(One) test
+//				if (a.BitLength == e + 1)
+//					break;
+//				BigInteger a1 = a.ShiftRight(e);
+////				if (a1.Equals(One))
+////					break;
+//
+//				int a1Lsw = a1.magnitude[a1.magnitude.Length - 1];
+//				if ((bLsw & 3) == 3 && (a1Lsw & 3) == 3)
+//					totalS = -totalS;
+//
+////				a = b.Mod(a1);
+//				a = b.Remainder(a1);
+//				b = a1;
+//			}
+//			return totalS;
+//		}
+
+		public long LongValue
+		{
+            get
+            {
+                if (sign == 0)
+                    return 0;
+
+                int n = magnitude.Length;
+
+                long v = magnitude[n - 1] & IMASK;
+                if (n > 1)
+                {
+                    v |= (magnitude[n - 2] & IMASK) << 32;
+                }
+
+                return sign < 0 ? -v : v;
+            }
+        }
+
+        public BigInteger Max(
+			BigInteger value)
+		{
+			return CompareTo(value) > 0 ? this : value;
+		}
+
+		public BigInteger Min(
+			BigInteger value)
+		{
+			return CompareTo(value) < 0 ? this : value;
+		}
+
+		public BigInteger Mod(
+			BigInteger m)
+		{
+			if (m.sign < 1)
+				throw new ArithmeticException("Modulus must be positive");
+
+			BigInteger biggie = Remainder(m);
+
+			return (biggie.sign >= 0 ? biggie : biggie.Add(m));
+		}
+
+		public BigInteger ModInverse(
+			BigInteger m)
+		{
+			if (m.sign < 1)
+				throw new ArithmeticException("Modulus must be positive");
+
+			// TODO Too slow at the moment
+//			// "Fast Key Exchange with Elliptic Curve Systems" R.Schoeppel
+//			if (m.TestBit(0))
+//			{
+//				//The Almost Inverse Algorithm
+//				int k = 0;
+//				BigInteger B = One, C = Zero, F = this, G = m, tmp;
+//
+//				for (;;)
+//				{
+//					// While F is even, do F=F/u, C=C*u, k=k+1.
+//					int zeroes = F.GetLowestSetBit();
+//					if (zeroes > 0)
+//					{
+//						F = F.ShiftRight(zeroes);
+//						C = C.ShiftLeft(zeroes);
+//						k += zeroes;
+//					}
+//
+//					// If F = 1, then return B,k.
+//					if (F.Equals(One))
+//					{
+//						BigInteger half = m.Add(One).ShiftRight(1);
+//						BigInteger halfK = half.ModPow(BigInteger.ValueOf(k), m);
+//						return B.Multiply(halfK).Mod(m);
+//					}
+//
+//					if (F.CompareTo(G) < 0)
+//					{
+//						tmp = G; G = F; F = tmp;
+//						tmp = B; B = C; C = tmp;
+//					}
+//
+//					F = F.Add(G);
+//					B = B.Add(C);
+//				}
+//			}
+
+            if (m.QuickPow2Check())
+            {
+                return ModInversePow2(m);
+            }
+
+            BigInteger d = this.Remainder(m);
+            BigInteger x;
+            BigInteger gcd = ExtEuclid(d, m, out x);
+
+			if (!gcd.Equals(One))
+				throw new ArithmeticException("Numbers not relatively prime.");
+
+			if (x.sign < 0)
+			{
+                x = x.Add(m);
+			}
+
+			return x;
+		}
+
+        private BigInteger ModInversePow2(BigInteger m)
+        {
+            Debug.Assert(m.SignValue > 0);
+            Debug.Assert(m.BitCount == 1);
+
+            if (!TestBit(0))
+            {
+                throw new ArithmeticException("Numbers not relatively prime.");
+            }
+
+            int pow = m.BitLength - 1;
+
+            long inv64 = ModInverse64(LongValue);
+            if (pow < 64)
+            {
+                inv64 &= ((1L << pow) - 1);
+            }
+
+            BigInteger x = BigInteger.ValueOf(inv64);
+
+            if (pow > 64)
+            {
+                BigInteger d = this.Remainder(m);
+                int bitsCorrect = 64;
+
+                do
+                {
+                    BigInteger t = x.Multiply(d).Remainder(m);
+                    x = x.Multiply(Two.Subtract(t)).Remainder(m);
+                    bitsCorrect <<= 1;
+                }
+                while (bitsCorrect < pow);
+
+                if (x.sign < 0)
+                {
+                    x = x.Add(m);
+                }
+            }
+
+            return x;
+        }
+
+        private static int ModInverse32(int d)
+        {
+            // Newton's method with initial estimate "correct to 4 bits"
+            Debug.Assert((d & 1) != 0);
+            int x = d + (((d + 1) & 4) << 1);   // d.x == 1 mod 2**4
+            Debug.Assert(((d * x) & 15) == 1);
+            x *= 2 - d * x;                     // d.x == 1 mod 2**8
+            x *= 2 - d * x;                     // d.x == 1 mod 2**16
+            x *= 2 - d * x;                     // d.x == 1 mod 2**32
+            Debug.Assert(d * x == 1);
+            return x;
+        }
+
+        private static long ModInverse64(long d)
+        {
+            // Newton's method with initial estimate "correct to 4 bits"
+            Debug.Assert((d & 1L) != 0);
+            long x = d + (((d + 1L) & 4L) << 1);    // d.x == 1 mod 2**4
+            Debug.Assert(((d * x) & 15L) == 1L);
+            x *= 2 - d * x;                         // d.x == 1 mod 2**8
+            x *= 2 - d * x;                         // d.x == 1 mod 2**16
+            x *= 2 - d * x;                         // d.x == 1 mod 2**32
+            x *= 2 - d * x;                         // d.x == 1 mod 2**64
+            Debug.Assert(d * x == 1L);
+            return x;
+        }
+
+		/**
+		 * Calculate the numbers u1, u2, and u3 such that:
+		 *
+		 * u1 * a + u2 * b = u3
+		 *
+		 * where u3 is the greatest common divider of a and b.
+		 * a and b using the extended Euclid algorithm (refer p. 323
+		 * of The Art of Computer Programming vol 2, 2nd ed).
+		 * This also seems to have the side effect of calculating
+		 * some form of multiplicative inverse.
+		 *
+		 * @param a    First number to calculate gcd for
+		 * @param b    Second number to calculate gcd for
+		 * @param u1Out      the return object for the u1 value
+		 * @param u2Out      the return object for the u2 value
+		 * @return     The greatest common divisor of a and b
+		 */
+		private static BigInteger ExtEuclid(
+			BigInteger	    a,
+			BigInteger	    b,
+			out BigInteger  u1Out)
+            //BigInteger	    u2Out)
+		{
+			BigInteger u1 = One;
+			BigInteger u3 = a;
+			BigInteger v1 = Zero;
+			BigInteger v3 = b;
+
+			while (v3.sign > 0)
+			{
+				BigInteger[] q = u3.DivideAndRemainder(v3);
+
+				BigInteger tmp = v1.Multiply(q[0]);
+				BigInteger tn = u1.Subtract(tmp);
+				u1 = v1;
+				v1 = tn;
+
+				u3 = v3;
+				v3 = q[1];
+			}
+
+            //if (u1Out != null)
+            //{
+            //    u1Out.sign = u1.sign;
+            //    u1Out.magnitude = u1.magnitude;
+            //}
+            u1Out = u1;
+
+            //if (u2Out != null)
+            //{
+            //    BigInteger tmp = u1.Multiply(a);
+            //    tmp = u3.Subtract(tmp);
+            //    BigInteger res = tmp.Divide(b);
+            //    u2Out.sign = res.sign;
+            //    u2Out.magnitude = res.magnitude;
+            //}
+
+			return u3;
+		}
+
+		private static void ZeroOut(
+			int[] x)
+		{
+			Array.Clear(x, 0, x.Length);
+		}
+
+        public BigInteger ModPow(BigInteger e, BigInteger m)
+        {
+            if (m.sign < 1)
+                throw new ArithmeticException("Modulus must be positive");
+
+            if (m.Equals(One))
+                return Zero;
+
+            if (e.sign == 0)
+                return One;
+
+            if (sign == 0)
+                return Zero;
+
+            bool negExp = e.sign < 0;
+            if (negExp)
+                e = e.Negate();
+
+            BigInteger result = this.Mod(m);
+            if (!e.Equals(One))
+            {
+                if ((m.magnitude[m.magnitude.Length - 1] & 1) == 0)
+                {
+                    result = ModPowBarrett(result, e, m);
+                }
+                else
+                {
+                    result = ModPowMonty(result, e, m, true);
+                }
+            }
+
+            if (negExp)
+                result = result.ModInverse(m);
+
+            return result;
+        }
+
+        private static BigInteger ModPowBarrett(BigInteger b, BigInteger e, BigInteger m)
+        {
+            int k = m.magnitude.Length;
+            BigInteger mr = One.ShiftLeft((k + 1) << 5);
+            BigInteger yu = One.ShiftLeft(k << 6).Divide(m);
+
+            // Sliding window from MSW to LSW
+            int extraBits = 0, expLength = e.BitLength;
+            while (expLength > ExpWindowThresholds[extraBits])
+            {
+                ++extraBits;
+            }
+
+            int numPowers = 1 << extraBits;
+            BigInteger[] oddPowers = new BigInteger[numPowers];
+            oddPowers[0] = b;
+
+            BigInteger b2 = ReduceBarrett(b.Square(), m, mr, yu);
+
+            for (int i = 1; i < numPowers; ++i)
+            {
+                oddPowers[i] = ReduceBarrett(oddPowers[i - 1].Multiply(b2), m, mr, yu);
+            }
+
+            int[] windowList = GetWindowList(e.magnitude, extraBits);
+            Debug.Assert(windowList.Length > 0);
+
+            int window = windowList[0];
+            int mult = window & 0xFF, lastZeroes = window >> 8;
+
+            BigInteger y;
+            if (mult == 1)
+            {
+                y = b2;
+                --lastZeroes;
+            }
+            else
+            {
+                y = oddPowers[mult >> 1];
+            }
+
+            int windowPos = 1;
+            while ((window = windowList[windowPos++]) != -1)
+            {
+                mult = window & 0xFF;
+
+                int bits = lastZeroes + BitLengthTable[mult];
+                for (int j = 0; j < bits; ++j)
+                {
+                    y = ReduceBarrett(y.Square(), m, mr, yu);
+                }
+
+                y = ReduceBarrett(y.Multiply(oddPowers[mult >> 1]), m, mr, yu);
+
+                lastZeroes = window >> 8;
+            }
+
+            for (int i = 0; i < lastZeroes; ++i)
+            {
+                y = ReduceBarrett(y.Square(), m, mr, yu);
+            }
+
+            return y;
+        }
+
+        private static BigInteger ReduceBarrett(BigInteger x, BigInteger m, BigInteger mr, BigInteger yu)
+        {
+            int xLen = x.BitLength, mLen = m.BitLength;
+            if (xLen < mLen)
+                return x;
+
+            if (xLen - mLen > 1)
+            {
+                int k = m.magnitude.Length;
+
+                BigInteger q1 = x.DivideWords(k - 1);
+                BigInteger q2 = q1.Multiply(yu); // TODO Only need partial multiplication here
+                BigInteger q3 = q2.DivideWords(k + 1);
+
+                BigInteger r1 = x.RemainderWords(k + 1);
+                BigInteger r2 = q3.Multiply(m); // TODO Only need partial multiplication here
+                BigInteger r3 = r2.RemainderWords(k + 1);
+
+                x = r1.Subtract(r3);
+                if (x.sign < 0)
+                {
+                    x = x.Add(mr);
+                }
+            }
+
+            while (x.CompareTo(m) >= 0)
+            {
+                x = x.Subtract(m);
+            }
+
+            return x;
+        }
+
+        private static BigInteger ModPowMonty(BigInteger b, BigInteger e, BigInteger m, bool convert)
+        {
+            int n = m.magnitude.Length;
+            int powR = 32 * n;
+            bool smallMontyModulus = m.BitLength + 2 <= powR;
+            uint mDash = (uint)m.GetMQuote();
+
+            // tmp = this * R mod m
+            if (convert)
+            {
+                b = b.ShiftLeft(powR).Remainder(m);
+            }
+
+            int[] yAccum = new int[n + 1];
+
+            int[] zVal = b.magnitude;
+            Debug.Assert(zVal.Length <= n);
+            if (zVal.Length < n)
+            {
+                int[] tmp = new int[n];
+                zVal.CopyTo(tmp, n - zVal.Length);
+                zVal = tmp;
+            }
+
+            // Sliding window from MSW to LSW
+
+            int extraBits = 0;
+
+            // Filter the common case of small RSA exponents with few bits set
+            if (e.magnitude.Length > 1 || e.BitCount > 2)
+            {
+                int expLength = e.BitLength;
+                while (expLength > ExpWindowThresholds[extraBits])
+                {
+                    ++extraBits;
+                }
+            }
+
+            int numPowers = 1 << extraBits;
+            int[][] oddPowers = new int[numPowers][];
+            oddPowers[0] = zVal;
+
+            int[] zSquared = Arrays.Clone(zVal);
+            SquareMonty(yAccum, zSquared, m.magnitude, mDash, smallMontyModulus);
+
+            for (int i = 1; i < numPowers; ++i)
+            {
+                oddPowers[i] = Arrays.Clone(oddPowers[i - 1]);
+                MultiplyMonty(yAccum, oddPowers[i], zSquared, m.magnitude, mDash, smallMontyModulus);
+            }
+
+            int[] windowList = GetWindowList(e.magnitude, extraBits);
+            Debug.Assert(windowList.Length > 1);
+
+            int window = windowList[0];
+            int mult = window & 0xFF, lastZeroes = window >> 8;
+
+            int[] yVal;
+            if (mult == 1)
+            {
+                yVal = zSquared;
+                --lastZeroes;
+            }
+            else
+            {
+                yVal = Arrays.Clone(oddPowers[mult >> 1]);
+            }
+
+            int windowPos = 1;
+            while ((window = windowList[windowPos++]) != -1)
+            {
+                mult = window & 0xFF;
+
+                int bits = lastZeroes + BitLengthTable[mult];
+                for (int j = 0; j < bits; ++j)
+                {
+                    SquareMonty(yAccum, yVal, m.magnitude, mDash, smallMontyModulus);
+                }
+
+                MultiplyMonty(yAccum, yVal, oddPowers[mult >> 1], m.magnitude, mDash, smallMontyModulus);
+
+                lastZeroes = window >> 8;
+            }
+
+            for (int i = 0; i < lastZeroes; ++i)
+            {
+                SquareMonty(yAccum, yVal, m.magnitude, mDash, smallMontyModulus);
+            }
+
+            if (convert)
+            {
+                // Return y * R^(-1) mod m
+                MontgomeryReduce(yVal, m.magnitude, mDash);
+            }
+            else if (smallMontyModulus && CompareTo(0, yVal, 0, m.magnitude) >= 0)
+            {
+                Subtract(0, yVal, 0, m.magnitude);
+            }
+
+            return new BigInteger(1, yVal, true);
+        }
+
+        private static int[] GetWindowList(int[] mag, int extraBits)
+        {
+            int v = mag[0];
+            Debug.Assert(v != 0);
+
+            int leadingBits = BitLen(v);
+
+            int resultSize = (((mag.Length - 1) << 5) + leadingBits) / (1 + extraBits) + 2;
+            int[] result = new int[resultSize];
+            int resultPos = 0;
+
+            int bitPos = 33 - leadingBits;
+            v <<= bitPos;
+
+            int mult = 1, multLimit = 1 << extraBits;
+            int zeroes = 0;
+
+            int i = 0;
+            for (; ; )
+            {
+                for (; bitPos < 32; ++bitPos)
+                {
+                    if (mult < multLimit)
+                    {
+                        mult = (mult << 1) | (int)((uint)v >> 31);
+                    }
+                    else if (v < 0)
+                    {
+                        result[resultPos++] = CreateWindowEntry(mult, zeroes);
+                        mult = 1;
+                        zeroes = 0;
+                    }
+                    else
+                    {
+                        ++zeroes;
+                    }
+
+                    v <<= 1;
+                }
+
+                if (++i == mag.Length)
+                {
+                    result[resultPos++] = CreateWindowEntry(mult, zeroes);
+                    break;
+                }
+
+                v = mag[i];
+                bitPos = 0;
+            }
+
+            result[resultPos] = -1;
+            return result;
+        }
+
+        private static int CreateWindowEntry(int mult, int zeroes)
+        {
+            while ((mult & 1) == 0)
+            {
+                mult >>= 1;
+                ++zeroes;
+            }
+
+            return mult | (zeroes << 8);
+        }
+
+        /**
+		 * return w with w = x * x - w is assumed to have enough space.
+		 */
+		private static int[] Square(
+			int[]	w,
+			int[]	x)
+		{
+            // Note: this method allows w to be only (2 * x.Length - 1) words if result will fit
+//			if (w.Length != 2 * x.Length)
+//				throw new ArgumentException("no I don't think so...");
+
+            ulong c;
+
+            int wBase = w.Length - 1;
+
+            for (int i = x.Length - 1; i > 0; --i)
+            {
+                ulong v = (uint)x[i];
+
+                c = v * v + (uint)w[wBase];
+                w[wBase] = (int)c;
+                c >>= 32;
+
+                for (int j = i - 1; j >= 0; --j)
+                {
+                    ulong prod = v * (uint)x[j];
+
+                    c += ((uint)w[--wBase] & UIMASK) + ((uint)prod << 1);
+                    w[wBase] = (int)c;
+                    c = (c >> 32) + (prod >> 31);
+                }
+
+                c += (uint)w[--wBase];
+                w[wBase] = (int)c;
+
+                if (--wBase >= 0)
+                {
+                    w[wBase] = (int)(c >> 32);
+                }
+                else
+                {
+                    Debug.Assert((c >> 32) == 0);
+                }
+
+                wBase += i;
+            }
+
+            c = (uint)x[0];
+
+            c = c * c + (uint)w[wBase];
+            w[wBase] = (int)c;
+
+            if (--wBase >= 0)
+            {
+                w[wBase] += (int)(c >> 32);
+            }
+            else
+            {
+                Debug.Assert((c >> 32) == 0);
+            }
+
+            return w;
+		}
+
+        /**
+		 * return x with x = y * z - x is assumed to have enough space.
+		 */
+		private static int[] Multiply(int[]	x, int[] y, int[] z)
+		{
+			int i = z.Length;
+
+			if (i < 1)
+				return x;
+
+			int xBase = x.Length - y.Length;
+
+			do
+			{
+				long a = z[--i] & IMASK;
+				long val = 0;
+
+				if (a != 0)
+				{
+					for (int j = y.Length - 1; j >= 0; j--)
+					{
+						val += a * (y[j] & IMASK) + (x[xBase + j] & IMASK);
+	
+						x[xBase + j] = (int)val;
+	
+						val = (long)((ulong)val >> 32);
+					}
+				}
+
+				--xBase;
+
+				if (xBase >= 0)
+				{
+					x[xBase] = (int)val;
+				}
+				else
+				{
+					Debug.Assert(val == 0);
+				}
+			}
+			while (i > 0);
+
+			return x;
+		}
+
+        /**
+		 * Calculate mQuote = -m^(-1) mod b with b = 2^32 (32 = word size)
+		 */
+        private int GetMQuote()
+		{
+			if (mQuote != 0)
+			{
+				return mQuote; // already calculated
+			}
+
+            Debug.Assert(this.sign > 0);
+
+            int d = -magnitude[magnitude.Length - 1];
+
+            Debug.Assert((d & 1) != 0);
+
+            return mQuote = ModInverse32(d);
+		}
+
+        private static void MontgomeryReduce(int[] x, int[] m, uint mDash) // mDash = -m^(-1) mod b
+        {
+            // NOTE: Not a general purpose reduction (which would allow x up to twice the bitlength of m)
+            Debug.Assert(x.Length == m.Length);
+
+            int n = m.Length;
+
+            for (int i = n - 1; i >= 0; --i)
+            {
+                uint x0 = (uint)x[n - 1];
+                ulong t = x0 * mDash;
+
+                ulong carry = t * (uint)m[n - 1] + x0;
+                Debug.Assert((uint)carry == 0);
+                carry >>= 32;
+
+                for (int j = n - 2; j >= 0; --j)
+                {
+                    carry += t * (uint)m[j] + (uint)x[j];
+                    x[j + 1] = (int)carry;
+                    carry >>= 32;
+                }
+
+                x[0] = (int)carry;
+                Debug.Assert(carry >> 32 == 0);
+            }
+
+            if (CompareTo(0, x, 0, m) >= 0)
+            {
+                Subtract(0, x, 0, m);
+            }
+        }
+
+        /**
+		 * Montgomery multiplication: a = x * y * R^(-1) mod m
+		 * <br/>
+		 * Based algorithm 14.36 of Handbook of Applied Cryptography.
+		 * <br/>
+		 * <li> m, x, y should have length n </li>
+		 * <li> a should have length (n + 1) </li>
+		 * <li> b = 2^32, R = b^n </li>
+		 * <br/>
+		 * The result is put in x
+		 * <br/>
+		 * NOTE: the indices of x, y, m, a different in HAC and in Java
+		 */
+		private static void MultiplyMonty(int[]	a, int[] x, int[] y, int[] m, uint mDash, bool smallMontyModulus)
+			// mDash = -m^(-1) mod b
+		{
+            int n = m.Length;
+
+            if (n == 1)
+            {
+                x[0] = (int)MultiplyMontyNIsOne((uint)x[0], (uint)y[0], (uint)m[0], mDash);
+                return;
+            }
+
+            uint y0 = (uint)y[n - 1];
+
+            {
+                ulong xi = (uint)x[n - 1];
+
+                ulong carry = xi * y0;
+                ulong t = (uint)carry * mDash;
+
+                ulong prod2 = t * (uint)m[n - 1];
+                carry += (uint)prod2;
+                Debug.Assert((uint)carry == 0);
+                carry = (carry >> 32) + (prod2 >> 32);
+
+                for (int j = n - 2; j >= 0; --j)
+                {
+                    ulong prod1 = xi * (uint)y[j];
+                    prod2 = t * (uint)m[j];
+
+                    carry += (prod1 & UIMASK) + (uint)prod2;
+                    a[j + 2] = (int)carry;
+                    carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32);
+                }
+
+                a[1] = (int)carry;
+                a[0] = (int)(carry >> 32);
+            }
+
+            for (int i = n - 2; i >= 0; --i)
+            {
+                uint a0 = (uint)a[n];
+                ulong xi = (uint)x[i];
+
+                ulong prod1 = xi * y0;
+                ulong carry = (prod1 & UIMASK) + a0;
+                ulong t = (uint)carry * mDash;
+
+                ulong prod2 = t * (uint)m[n - 1];
+                carry += (uint)prod2;
+                Debug.Assert((uint)carry == 0);
+                carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32);
+
+                for (int j = n - 2; j >= 0; --j)
+                {
+                    prod1 = xi * (uint)y[j];
+                    prod2 = t * (uint)m[j];
+
+                    carry += (prod1 & UIMASK) + (uint)prod2 + (uint)a[j + 1];
+                    a[j + 2] = (int)carry;
+                    carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32);
+                }
+
+                carry += (uint)a[0];
+                a[1] = (int)carry;
+                a[0] = (int)(carry >> 32);
+            }
+
+            if (!smallMontyModulus && CompareTo(0, a, 0, m) >= 0)
+            {
+                Subtract(0, a, 0, m);
+            }
+
+            Array.Copy(a, 1, x, 0, n);
+        }
+
+        private static void SquareMonty(int[] a, int[] x, int[] m, uint mDash, bool smallMontyModulus)
+            // mDash = -m^(-1) mod b
+        {
+            int n = m.Length;
+
+            if (n == 1)
+            {
+                uint xVal = (uint)x[0];
+                x[0] = (int)MultiplyMontyNIsOne(xVal, xVal, (uint)m[0], mDash);
+                return;
+            }
+
+            ulong x0 = (uint)x[n - 1];
+
+            {
+                ulong carry = x0 * x0;
+                ulong t = (uint)carry * mDash;
+
+                ulong prod2 = t * (uint)m[n - 1];
+                carry += (uint)prod2;
+                Debug.Assert((uint)carry == 0);
+                carry = (carry >> 32) + (prod2 >> 32);
+
+                for (int j = n - 2; j >= 0; --j)
+                {
+                    ulong prod1 = x0 * (uint)x[j];
+                    prod2 = t * (uint)m[j];
+
+                    carry += (prod2 & UIMASK) + ((uint)prod1 << 1);
+                    a[j + 2] = (int)carry;
+                    carry = (carry >> 32) + (prod1 >> 31) + (prod2 >> 32);
+                }
+
+                a[1] = (int)carry;
+                a[0] = (int)(carry >> 32);
+            }
+
+            for (int i = n - 2; i >= 0; --i)
+            {
+                uint a0 = (uint)a[n];
+                ulong t = a0 * mDash;
+
+                ulong carry = t * (uint)m[n - 1] + a0;
+                Debug.Assert((uint)carry == 0);
+                carry >>= 32;
+
+                for (int j = n - 2; j > i; --j)
+                {
+                    carry += t * (uint)m[j] + (uint)a[j + 1];
+                    a[j + 2] = (int)carry;
+                    carry >>= 32;
+                }
+
+                ulong xi = (uint)x[i];
+
+                {
+                    ulong prod1 = xi * xi;
+                    ulong prod2 = t * (uint)m[i];
+
+                    carry += (prod1 & UIMASK) + (uint)prod2 + (uint)a[i + 1];
+                    a[i + 2] = (int)carry;
+                    carry = (carry >> 32) + (prod1 >> 32) + (prod2 >> 32);
+                }
+
+                for (int j = i - 1; j >= 0; --j)
+                {
+                    ulong prod1 = xi * (uint)x[j];
+                    ulong prod2 = t * (uint)m[j];
+
+                    carry += (prod2 & UIMASK) + ((uint)prod1 << 1) + (uint)a[j + 1];
+                    a[j + 2] = (int)carry;
+                    carry = (carry >> 32) + (prod1 >> 31) + (prod2 >> 32);
+                }
+
+                carry += (uint)a[0];
+                a[1] = (int)carry;
+                a[0] = (int)(carry >> 32);
+            }
+
+            if (!smallMontyModulus && CompareTo(0, a, 0, m) >= 0)
+            {
+                Subtract(0, a, 0, m);
+            }
+
+            Array.Copy(a, 1, x, 0, n);
+        }
+
+        private static uint MultiplyMontyNIsOne(uint x, uint y, uint m, uint mDash)
+		{
+            ulong carry = (ulong)x * y;
+            uint t = (uint)carry * mDash;
+            ulong um = m;
+            ulong prod2 = um * t;
+            carry += (uint)prod2;
+            Debug.Assert((uint)carry == 0);
+            carry = (carry >> 32) + (prod2 >> 32);
+            if (carry > um)
+            {
+                carry -= um;
+            }
+            Debug.Assert(carry < um);
+            return (uint)carry;
+		}
+
+        public BigInteger Multiply(
+			BigInteger val)
+		{
+            if (val == this)
+                return Square();
+
+            if ((sign & val.sign) == 0)
+				return Zero;
+
+            if (val.QuickPow2Check()) // val is power of two
+			{
+				BigInteger result = this.ShiftLeft(val.Abs().BitLength - 1);
+				return val.sign > 0 ? result : result.Negate();
+			}
+
+			if (this.QuickPow2Check()) // this is power of two
+			{
+				BigInteger result = val.ShiftLeft(this.Abs().BitLength - 1);
+				return this.sign > 0 ? result : result.Negate();
+			}
+
+            int resLength = magnitude.Length + val.magnitude.Length;
+            int[] res = new int[resLength];
+
+            Multiply(res, this.magnitude, val.magnitude);
+
+            int resSign = sign ^ val.sign ^ 1;
+            return new BigInteger(resSign, res, true);
+		}
+
+        public BigInteger Square()
+        {
+            if (sign == 0)
+                return Zero;
+            if (this.QuickPow2Check())
+                return ShiftLeft(Abs().BitLength - 1);
+            int resLength = magnitude.Length << 1;
+            if ((uint)magnitude[0] >> 16 == 0)
+                --resLength;
+            int[] res = new int[resLength];
+            Square(res, magnitude);
+            return new BigInteger(1, res, false);
+        }
+
+        public BigInteger Negate()
+		{
+			if (sign == 0)
+				return this;
+
+			return new BigInteger(-sign, magnitude, false);
+		}
+
+		public BigInteger NextProbablePrime()
+		{
+			if (sign < 0)
+				throw new ArithmeticException("Cannot be called on value < 0");
+
+			if (CompareTo(Two) < 0)
+				return Two;
+
+			BigInteger n = Inc().SetBit(0);
+
+			while (!n.CheckProbablePrime(100, RandomSource))
+			{
+				n = n.Add(Two);
+			}
+
+			return n;
+		}
+
+		public BigInteger Not()
+		{
+			return Inc().Negate();
+		}
+
+		public BigInteger Pow(int exp)
+		{
+			if (exp <= 0)
+			{
+                if (exp < 0)
+    				throw new ArithmeticException("Negative exponent");
+
+                return One;
+			}
+
+            if (sign == 0)
+			{
+				return this;
+			}
+
+            if (QuickPow2Check())
+            {
+                long powOf2 = (long)exp * (BitLength - 1);
+                if (powOf2 > Int32.MaxValue)
+                {
+                    throw new ArithmeticException("Result too large");
+                }
+                return One.ShiftLeft((int)powOf2); 
+            }
+
+            BigInteger y = One;
+			BigInteger z = this;
+
+			for (;;)
+			{
+				if ((exp & 0x1) == 1)
+				{
+					y = y.Multiply(z);
+				}
+				exp >>= 1;
+				if (exp == 0) break;
+				z = z.Multiply(z);
+			}
+
+			return y;
+		}
+
+		public static BigInteger ProbablePrime(
+			int bitLength,
+			Random random)
+		{
+			return new BigInteger(bitLength, 100, random);
+		}
+
+		private int Remainder(
+			int m)
+		{
+			Debug.Assert(m > 0);
+
+			long acc = 0;
+			for (int pos = 0; pos < magnitude.Length; ++pos)
+			{
+				long posVal = (uint) magnitude[pos];
+				acc = (acc << 32 | posVal) % m;
+			}
+
+			return (int) acc;
+		}
+
+		/**
+		 * return x = x % y - done in place (y value preserved)
+		 */
+		private static int[] Remainder(
+			int[] x,
+			int[] y)
+		{
+			int xStart = 0;
+			while (xStart < x.Length && x[xStart] == 0)
+			{
+				++xStart;
+			}
+
+			int yStart = 0;
+			while (yStart < y.Length && y[yStart] == 0)
+			{
+				++yStart;
+			}
+
+			Debug.Assert(yStart < y.Length);
+
+			int xyCmp = CompareNoLeadingZeroes(xStart, x, yStart, y);
+
+			if (xyCmp > 0)
+			{
+				int yBitLength = CalcBitLength(1, yStart, y);
+				int xBitLength = CalcBitLength(1, xStart, x);
+				int shift = xBitLength - yBitLength;
+
+				int[] c;
+				int cStart = 0;
+				int cBitLength = yBitLength;
+				if (shift > 0)
+				{
+					c = ShiftLeft(y, shift);
+					cBitLength += shift;
+					Debug.Assert(c[0] != 0);
+				}
+				else
+				{
+					int len = y.Length - yStart;
+					c = new int[len];
+					Array.Copy(y, yStart, c, 0, len);
+				}
+
+				for (;;)
+				{
+					if (cBitLength < xBitLength
+						|| CompareNoLeadingZeroes(xStart, x, cStart, c) >= 0)
+					{
+						Subtract(xStart, x, cStart, c);
+
+						while (x[xStart] == 0)
+						{
+							if (++xStart == x.Length)
+								return x;
+						}
+
+						//xBitLength = CalcBitLength(xStart, x);
+						xBitLength = 32 * (x.Length - xStart - 1) + BitLen(x[xStart]);
+
+						if (xBitLength <= yBitLength)
+						{
+							if (xBitLength < yBitLength)
+								return x;
+
+							xyCmp = CompareNoLeadingZeroes(xStart, x, yStart, y);
+
+							if (xyCmp <= 0)
+								break;
+						}
+					}
+
+					shift = cBitLength - xBitLength;
+
+					// NB: The case where c[cStart] is 1-bit is harmless
+					if (shift == 1)
+					{
+						uint firstC = (uint) c[cStart] >> 1;
+						uint firstX = (uint) x[xStart];
+						if (firstC > firstX)
+							++shift;
+					}
+
+					if (shift < 2)
+					{
+						ShiftRightOneInPlace(cStart, c);
+						--cBitLength;
+					}
+					else
+					{
+						ShiftRightInPlace(cStart, c, shift);
+						cBitLength -= shift;
+					}
+
+					//cStart = c.Length - ((cBitLength + 31) / 32);
+					while (c[cStart] == 0)
+					{
+						++cStart;
+					}
+				}
+			}
+
+			if (xyCmp == 0)
+			{
+				Array.Clear(x, xStart, x.Length - xStart);
+			}
+
+			return x;
+		}
+
+		public BigInteger Remainder(
+			BigInteger n)
+		{
+			if (n.sign == 0)
+				throw new ArithmeticException("Division by zero error");
+
+			if (this.sign == 0)
+				return Zero;
+
+			// For small values, use fast remainder method
+			if (n.magnitude.Length == 1)
+			{
+				int val = n.magnitude[0];
+
+				if (val > 0)
+				{
+					if (val == 1)
+						return Zero;
+
+					// TODO Make this func work on uint, and handle val == 1?
+					int rem = Remainder(val);
+
+					return rem == 0
+						?	Zero
+						:	new BigInteger(sign, new int[]{ rem }, false);
+				}
+			}
+
+			if (CompareNoLeadingZeroes(0, magnitude, 0, n.magnitude) < 0)
+				return this;
+
+			int[] result;
+			if (n.QuickPow2Check())  // n is power of two
+			{
+				// TODO Move before small values branch above?
+				result = LastNBits(n.Abs().BitLength - 1);
+			}
+			else
+			{
+				result = (int[]) this.magnitude.Clone();
+				result = Remainder(result, n.magnitude);
+			}
+
+			return new BigInteger(sign, result, true);
+		}
+
+		private int[] LastNBits(
+			int n)
+		{
+			if (n < 1)
+				return ZeroMagnitude;
+
+			int numWords = (n + BitsPerInt - 1) / BitsPerInt;
+			numWords = System.Math.Min(numWords, this.magnitude.Length);
+			int[] result = new int[numWords];
+
+			Array.Copy(this.magnitude, this.magnitude.Length - numWords, result, 0, numWords);
+
+            int excessBits = (numWords << 5) - n;
+            if (excessBits > 0)
+            {
+                result[0] &= (int)(UInt32.MaxValue >> excessBits);
+            }
+
+            return result;
+		}
+
+        private BigInteger DivideWords(int w)
+        {
+            Debug.Assert(w >= 0);
+            int n = magnitude.Length;
+            if (w >= n)
+                return Zero;
+            int[] mag = new int[n - w];
+            Array.Copy(magnitude, 0, mag, 0, n - w);
+            return new BigInteger(sign, mag, false);
+        }
+
+        private BigInteger RemainderWords(int w)
+        {
+            Debug.Assert(w >= 0);
+            int n = magnitude.Length;
+            if (w >= n)
+                return this;
+            int[] mag = new int[w];
+            Array.Copy(magnitude, n - w, mag, 0, w);
+            return new BigInteger(sign, mag, false);
+        }
+
+        /**
+		 * do a left shift - this returns a new array.
+		 */
+		private static int[] ShiftLeft(
+			int[]	mag,
+			int		n)
+		{
+			int nInts = (int)((uint)n >> 5);
+			int nBits = n & 0x1f;
+			int magLen = mag.Length;
+			int[] newMag;
+
+			if (nBits == 0)
+			{
+				newMag = new int[magLen + nInts];
+				mag.CopyTo(newMag, 0);
+			}
+			else
+			{
+				int i = 0;
+				int nBits2 = 32 - nBits;
+				int highBits = (int)((uint)mag[0] >> nBits2);
+
+				if (highBits != 0)
+				{
+					newMag = new int[magLen + nInts + 1];
+					newMag[i++] = highBits;
+				}
+				else
+				{
+					newMag = new int[magLen + nInts];
+				}
+
+				int m = mag[0];
+				for (int j = 0; j < magLen - 1; j++)
+				{
+					int next = mag[j + 1];
+
+					newMag[i++] = (m << nBits) | (int)((uint)next >> nBits2);
+					m = next;
+				}
+
+				newMag[i] = mag[magLen - 1] << nBits;
+			}
+
+			return newMag;
+		}
+
+        private static int ShiftLeftOneInPlace(int[] x, int carry)
+        {
+            Debug.Assert(carry == 0 || carry == 1);
+            int pos = x.Length;
+            while (--pos >= 0)
+            {
+                uint val = (uint)x[pos];
+                x[pos] = (int)(val << 1) | carry;
+                carry = (int)(val >> 31);
+            }
+            return carry;
+        }
+
+		public BigInteger ShiftLeft(
+			int n)
+		{
+			if (sign == 0 || magnitude.Length == 0)
+				return Zero;
+
+			if (n == 0)
+				return this;
+
+			if (n < 0)
+				return ShiftRight(-n);
+
+			BigInteger result = new BigInteger(sign, ShiftLeft(magnitude, n), true);
+
+			if (this.nBits != -1)
+			{
+				result.nBits = sign > 0
+					?	this.nBits
+					:	this.nBits + n;
+			}
+
+			if (this.nBitLength != -1)
+			{
+				result.nBitLength = this.nBitLength + n;
+			}
+
+			return result;
+		}
+
+		/**
+		 * do a right shift - this does it in place.
+		 */
+		private static void ShiftRightInPlace(
+			int		start,
+			int[]	mag,
+			int		n)
+		{
+			int nInts = (int)((uint)n >> 5) + start;
+			int nBits = n & 0x1f;
+			int magEnd = mag.Length - 1;
+
+			if (nInts != start)
+			{
+				int delta = (nInts - start);
+
+				for (int i = magEnd; i >= nInts; i--)
+				{
+					mag[i] = mag[i - delta];
+				}
+				for (int i = nInts - 1; i >= start; i--)
+				{
+					mag[i] = 0;
+				}
+			}
+
+			if (nBits != 0)
+			{
+				int nBits2 = 32 - nBits;
+				int m = mag[magEnd];
+
+				for (int i = magEnd; i > nInts; --i)
+				{
+					int next = mag[i - 1];
+
+					mag[i] = (int)((uint)m >> nBits) | (next << nBits2);
+					m = next;
+				}
+
+				mag[nInts] = (int)((uint)mag[nInts] >> nBits);
+			}
+		}
+
+		/**
+		 * do a right shift by one - this does it in place.
+		 */
+		private static void ShiftRightOneInPlace(
+			int		start,
+			int[]	mag)
+		{
+			int i = mag.Length;
+			int m = mag[i - 1];
+
+			while (--i > start)
+			{
+				int next = mag[i - 1];
+				mag[i] = ((int)((uint)m >> 1)) | (next << 31);
+				m = next;
+			}
+
+			mag[start] = (int)((uint)mag[start] >> 1);
+		}
+
+        public BigInteger ShiftRight(
+			int n)
+		{
+			if (n == 0)
+				return this;
+
+			if (n < 0)
+				return ShiftLeft(-n);
+
+			if (n >= BitLength)
+				return (this.sign < 0 ? One.Negate() : Zero);
+
+//			int[] res = (int[]) this.magnitude.Clone();
+//
+//			ShiftRightInPlace(0, res, n);
+//
+//			return new BigInteger(this.sign, res, true);
+
+			int resultLength = (BitLength - n + 31) >> 5;
+			int[] res = new int[resultLength];
+
+			int numInts = n >> 5;
+			int numBits = n & 31;
+
+			if (numBits == 0)
+			{
+				Array.Copy(this.magnitude, 0, res, 0, res.Length);
+			}
+			else
+			{
+				int numBits2 = 32 - numBits;
+
+				int magPos = this.magnitude.Length - 1 - numInts;
+				for (int i = resultLength - 1; i >= 0; --i)
+				{
+					res[i] = (int)((uint) this.magnitude[magPos--] >> numBits);
+
+					if (magPos >= 0)
+					{
+						res[i] |= this.magnitude[magPos] << numBits2;
+					}
+				}
+			}
+
+			Debug.Assert(res[0] != 0);
+
+			return new BigInteger(this.sign, res, false);
+		}
+
+        public int SignValue
+		{
+			get { return sign; }
+		}
+
+		/**
+		 * returns x = x - y - we assume x is >= y
+		 */
+		private static int[] Subtract(
+			int		xStart,
+			int[]	x,
+			int		yStart,
+			int[]	y)
+		{
+			Debug.Assert(yStart < y.Length);
+			Debug.Assert(x.Length - xStart >= y.Length - yStart);
+
+			int iT = x.Length;
+			int iV = y.Length;
+			long m;
+			int borrow = 0;
+
+			do
+			{
+				m = (x[--iT] & IMASK) - (y[--iV] & IMASK) + borrow;
+				x[iT] = (int) m;
+
+//				borrow = (m < 0) ? -1 : 0;
+				borrow = (int)(m >> 63);
+			}
+			while (iV > yStart);
+
+			if (borrow != 0)
+			{
+				while (--x[--iT] == -1)
+				{
+				}
+			}
+
+			return x;
+		}
+
+		public BigInteger Subtract(
+			BigInteger n)
+		{
+			if (n.sign == 0)
+				return this;
+
+			if (this.sign == 0)
+				return n.Negate();
+
+			if (this.sign != n.sign)
+				return Add(n.Negate());
+
+			int compare = CompareNoLeadingZeroes(0, magnitude, 0, n.magnitude);
+			if (compare == 0)
+				return Zero;
+
+			BigInteger bigun, lilun;
+			if (compare < 0)
+			{
+				bigun = n;
+				lilun = this;
+			}
+			else
+			{
+				bigun = this;
+				lilun = n;
+			}
+
+			return new BigInteger(this.sign * compare, doSubBigLil(bigun.magnitude, lilun.magnitude), true);
+		}
+
+		private static int[] doSubBigLil(
+			int[]	bigMag,
+			int[]	lilMag)
+		{
+			int[] res = (int[]) bigMag.Clone();
+
+			return Subtract(0, res, 0, lilMag);
+		}
+
+		public byte[] ToByteArray()
+		{
+			return ToByteArray(false);
+		}
+
+		public byte[] ToByteArrayUnsigned()
+		{
+			return ToByteArray(true);
+		}
+
+		private byte[] ToByteArray(
+			bool unsigned)
+		{
+			if (sign == 0)
+				return unsigned ? ZeroEncoding : new byte[1];
+
+			int nBits = (unsigned && sign > 0)
+				?	BitLength
+				:	BitLength + 1;
+
+			int nBytes = GetByteLength(nBits);
+			byte[] bytes = new byte[nBytes];
+
+			int magIndex = magnitude.Length;
+			int bytesIndex = bytes.Length;
+
+			if (sign > 0)
+			{
+				while (magIndex > 1)
+				{
+					uint mag = (uint) magnitude[--magIndex];
+					bytes[--bytesIndex] = (byte) mag;
+					bytes[--bytesIndex] = (byte)(mag >> 8);
+					bytes[--bytesIndex] = (byte)(mag >> 16);
+					bytes[--bytesIndex] = (byte)(mag >> 24);
+				}
+
+				uint lastMag = (uint) magnitude[0];
+				while (lastMag > byte.MaxValue)
+				{
+					bytes[--bytesIndex] = (byte) lastMag;
+					lastMag >>= 8;
+				}
+
+				bytes[--bytesIndex] = (byte) lastMag;
+			}
+			else // sign < 0
+			{
+				bool carry = true;
+
+				while (magIndex > 1)
+				{
+					uint mag = ~((uint) magnitude[--magIndex]);
+
+					if (carry)
+					{
+						carry = (++mag == uint.MinValue);
+					}
+
+					bytes[--bytesIndex] = (byte) mag;
+					bytes[--bytesIndex] = (byte)(mag >> 8);
+					bytes[--bytesIndex] = (byte)(mag >> 16);
+					bytes[--bytesIndex] = (byte)(mag >> 24);
+				}
+
+				uint lastMag = (uint) magnitude[0];
+
+				if (carry)
+				{
+					// Never wraps because magnitude[0] != 0
+					--lastMag;
+				}
+
+				while (lastMag > byte.MaxValue)
+				{
+					bytes[--bytesIndex] = (byte) ~lastMag;
+					lastMag >>= 8;
+				}
+
+				bytes[--bytesIndex] = (byte) ~lastMag;
+
+				if (bytesIndex > 0)
+				{
+					bytes[--bytesIndex] = byte.MaxValue;
+				}
+			}
+
+			return bytes;
+		}
+
+		public override string ToString()
+		{
+			return ToString(10);
+		}
+
+        public string ToString(int radix)
+        {
+            // TODO Make this method work for other radices (ideally 2 <= radix <= 36 as in Java)
+
+            switch (radix)
+            {
+                case 2:
+                case 8:
+                case 10:
+                case 16:
+                    break;
+                default:
+                    throw new FormatException("Only bases 2, 8, 10, 16 are allowed");
+            }
+
+            // NB: Can only happen to internally managed instances
+            if (magnitude == null)
+                return "null";
+
+            if (sign == 0)
+                return "0";
+
+
+            // NOTE: This *should* be unnecessary, since the magnitude *should* never have leading zero digits
+            int firstNonZero = 0;
+            while (firstNonZero < magnitude.Length)
+            {
+                if (magnitude[firstNonZero] != 0)
+                {
+                    break;
+                }
+                ++firstNonZero;
+            }
+
+            if (firstNonZero == magnitude.Length)
+            {
+                return "0";
+            }
+
+
+            StringBuilder sb = new StringBuilder();
+            if (sign == -1)
+            {
+                sb.Append('-');
+            }
+
+            switch (radix)
+            {
+            case 2:
+            {
+                int pos = firstNonZero;
+                sb.Append(Convert.ToString(magnitude[pos], 2));
+                while (++pos < magnitude.Length)
+                {
+                    AppendZeroExtendedString(sb, Convert.ToString(magnitude[pos], 2), 32);
+                }
+                break;
+            }
+            case 8:
+            {
+                int mask = (1 << 30) - 1;
+                BigInteger u = this.Abs();
+                int bits = u.BitLength;
+                IList S = Platform.CreateArrayList();
+                while (bits > 30)
+                {
+                    S.Add(Convert.ToString(u.IntValue & mask, 8));
+                    u = u.ShiftRight(30);
+                    bits -= 30;
+                }
+                sb.Append(Convert.ToString(u.IntValue, 8));
+                for (int i = S.Count - 1; i >= 0; --i)
+                {
+                    AppendZeroExtendedString(sb, (string)S[i], 10);
+                }
+                break;
+            }
+            case 16:
+            {
+                int pos = firstNonZero;
+                sb.Append(Convert.ToString(magnitude[pos], 16));
+                while (++pos < magnitude.Length)
+                {
+                    AppendZeroExtendedString(sb, Convert.ToString(magnitude[pos], 16), 8);
+                }
+                break;
+            }
+            // TODO This could work for other radices if there is an alternative to Convert.ToString method
+            //default:
+            case 10:
+            {
+                BigInteger q = this.Abs();
+                if (q.BitLength < 64)
+                {
+                    sb.Append(Convert.ToString(q.LongValue, radix));
+                    break;
+                }
+
+                // Based on algorithm 1a from chapter 4.4 in Seminumerical Algorithms (Knuth)
+
+                // Work out the largest power of 'rdx' that is a positive 64-bit integer
+                // TODO possibly cache power/exponent against radix?
+                long limit = Int64.MaxValue / radix;
+                long power = radix;
+                int exponent = 1;
+                while (power <= limit)
+                {
+                    power *= radix;
+                    ++exponent;
+                }
+
+                BigInteger bigPower = BigInteger.ValueOf(power);
+
+                IList S = Platform.CreateArrayList();
+                while (q.CompareTo(bigPower) >= 0)
+                {
+                    BigInteger[] qr = q.DivideAndRemainder(bigPower);
+                    S.Add(Convert.ToString(qr[1].LongValue, radix));
+                    q = qr[0];
+                }
+
+                sb.Append(Convert.ToString(q.LongValue, radix));
+                for (int i = S.Count - 1; i >= 0; --i)
+                {
+                    AppendZeroExtendedString(sb, (string)S[i], exponent);
+                }
+                break;
+            }
+            }
+
+            return sb.ToString();
+        }
+
+        private static void AppendZeroExtendedString(StringBuilder sb, string s, int minLength)
+        {
+            for (int len = s.Length; len < minLength; ++len)
+            {
+                sb.Append('0');
+            }
+            sb.Append(s);
+        }
+
+        private static BigInteger CreateUValueOf(
+			ulong value)
+		{
+			int msw = (int)(value >> 32);
+			int lsw = (int)value;
+
+			if (msw != 0)
+				return new BigInteger(1, new int[] { msw, lsw }, false);
+
+			if (lsw != 0)
+			{
+				BigInteger n = new BigInteger(1, new int[] { lsw }, false);
+				// Check for a power of two
+				if ((lsw & -lsw) == lsw)
+				{
+					n.nBits = 1;
+				}
+				return n;
+			}
+
+			return Zero;
+		}
+
+		private static BigInteger CreateValueOf(
+			long value)
+		{
+			if (value < 0)
+			{
+				if (value == long.MinValue)
+					return CreateValueOf(~value).Not();
+
+				return CreateValueOf(-value).Negate();
+			}
+
+			return CreateUValueOf((ulong)value);
+		}
+
+        public static BigInteger ValueOf(
+			long value)
+		{
+            if (value >= 0 && value < SMALL_CONSTANTS.Length)
+            {
+                return SMALL_CONSTANTS[value];
+            }
+
+            return CreateValueOf(value);
+		}
+
+        public int GetLowestSetBit()
+		{
+			if (this.sign == 0)
+				return -1;
+
+            return GetLowestSetBitMaskFirst(-1);
+		}
+
+        private int GetLowestSetBitMaskFirst(int firstWordMask)
+        {
+            int w = magnitude.Length, offset = 0;
+
+            uint word = (uint)(magnitude[--w] & firstWordMask);
+            Debug.Assert(magnitude[0] != 0);
+
+            while (word == 0)
+            {
+                word = (uint)magnitude[--w];
+                offset += 32;
+            }
+
+            while ((word & 0xFF) == 0)
+            {
+                word >>= 8;
+                offset += 8;
+            }
+
+            while ((word & 1) == 0)
+            {
+                word >>= 1;
+                ++offset;
+            }
+
+            return offset;
+        }
+
+        public bool TestBit(
+			int n)
+		{
+			if (n < 0)
+				throw new ArithmeticException("Bit position must not be negative");
+
+			if (sign < 0)
+				return !Not().TestBit(n);
+
+			int wordNum = n / 32;
+			if (wordNum >= magnitude.Length)
+				return false;
+
+			int word = magnitude[magnitude.Length - 1 - wordNum];
+			return ((word >> (n % 32)) & 1) > 0;
+		}
+
+		public BigInteger Or(
+			BigInteger value)
+		{
+			if (this.sign == 0)
+				return value;
+
+			if (value.sign == 0)
+				return this;
+
+			int[] aMag = this.sign > 0
+				? this.magnitude
+				: Add(One).magnitude;
+
+			int[] bMag = value.sign > 0
+				? value.magnitude
+				: value.Add(One).magnitude;
+
+			bool resultNeg = sign < 0 || value.sign < 0;
+			int resultLength = System.Math.Max(aMag.Length, bMag.Length);
+			int[] resultMag = new int[resultLength];
+
+			int aStart = resultMag.Length - aMag.Length;
+			int bStart = resultMag.Length - bMag.Length;
+
+			for (int i = 0; i < resultMag.Length; ++i)
+			{
+				int aWord = i >= aStart ? aMag[i - aStart] : 0;
+				int bWord = i >= bStart ? bMag[i - bStart] : 0;
+
+				if (this.sign < 0)
+				{
+					aWord = ~aWord;
+				}
+
+				if (value.sign < 0)
+				{
+					bWord = ~bWord;
+				}
+
+				resultMag[i] = aWord | bWord;
+
+				if (resultNeg)
+				{
+					resultMag[i] = ~resultMag[i];
+				}
+			}
+
+			BigInteger result = new BigInteger(1, resultMag, true);
+
+			// TODO Optimise this case
+			if (resultNeg)
+			{
+				result = result.Not();
+			}
+
+			return result;
+		}
+
+		public BigInteger Xor(
+			BigInteger value)
+		{
+			if (this.sign == 0)
+				return value;
+
+			if (value.sign == 0)
+				return this;
+
+			int[] aMag = this.sign > 0
+				? this.magnitude
+				: Add(One).magnitude;
+
+			int[] bMag = value.sign > 0
+				? value.magnitude
+				: value.Add(One).magnitude;
+
+			// TODO Can just replace with sign != value.sign?
+			bool resultNeg = (sign < 0 && value.sign >= 0) || (sign >= 0 && value.sign < 0);
+			int resultLength = System.Math.Max(aMag.Length, bMag.Length);
+			int[] resultMag = new int[resultLength];
+
+			int aStart = resultMag.Length - aMag.Length;
+			int bStart = resultMag.Length - bMag.Length;
+
+			for (int i = 0; i < resultMag.Length; ++i)
+			{
+				int aWord = i >= aStart ? aMag[i - aStart] : 0;
+				int bWord = i >= bStart ? bMag[i - bStart] : 0;
+
+				if (this.sign < 0)
+				{
+					aWord = ~aWord;
+				}
+
+				if (value.sign < 0)
+				{
+					bWord = ~bWord;
+				}
+
+				resultMag[i] = aWord ^ bWord;
+
+				if (resultNeg)
+				{
+					resultMag[i] = ~resultMag[i];
+				}
+			}
+
+			BigInteger result = new BigInteger(1, resultMag, true);
+
+			// TODO Optimise this case
+			if (resultNeg)
+			{
+				result = result.Not();
+			}
+
+			return result;
+		}
+
+		public BigInteger SetBit(
+			int n)
+		{
+			if (n < 0)
+				throw new ArithmeticException("Bit address less than zero");
+
+			if (TestBit(n))
+				return this;
+
+			// TODO Handle negative values and zero
+			if (sign > 0 && n < (BitLength - 1))
+				return FlipExistingBit(n);
+
+			return Or(One.ShiftLeft(n));
+		}
+
+		public BigInteger ClearBit(
+			int n)
+		{
+			if (n < 0)
+				throw new ArithmeticException("Bit address less than zero");
+
+			if (!TestBit(n))
+				return this;
+
+			// TODO Handle negative values
+			if (sign > 0 && n < (BitLength - 1))
+				return FlipExistingBit(n);
+
+			return AndNot(One.ShiftLeft(n));
+		}
+
+		public BigInteger FlipBit(
+			int n)
+		{
+			if (n < 0)
+				throw new ArithmeticException("Bit address less than zero");
+
+			// TODO Handle negative values and zero
+			if (sign > 0 && n < (BitLength - 1))
+				return FlipExistingBit(n);
+
+			return Xor(One.ShiftLeft(n));
+		}
+
+		private BigInteger FlipExistingBit(
+			int n)
+		{
+			Debug.Assert(sign > 0);
+			Debug.Assert(n >= 0);
+			Debug.Assert(n < BitLength - 1);
+
+			int[] mag = (int[]) this.magnitude.Clone();
+			mag[mag.Length - 1 - (n >> 5)] ^= (1 << (n & 31)); // Flip bit
+			//mag[mag.Length - 1 - (n / 32)] ^= (1 << (n % 32));
+			return new BigInteger(this.sign, mag, false);
+		}
+	}
+}
diff --git a/crypto/src/math/ec/ECAlgorithms.cs b/crypto/src/math/ec/ECAlgorithms.cs
new file mode 100644
index 000000000..be4fd1b14
--- /dev/null
+++ b/crypto/src/math/ec/ECAlgorithms.cs
@@ -0,0 +1,93 @@
+using System;
+
+using Org.BouncyCastle.Math;
+
+namespace Org.BouncyCastle.Math.EC
+{
+	public class ECAlgorithms
+	{
+		public static ECPoint SumOfTwoMultiplies(ECPoint P, BigInteger a,
+			ECPoint Q, BigInteger b)
+		{
+			ECCurve c = P.Curve;
+			if (!c.Equals(Q.Curve))
+				throw new ArgumentException("P and Q must be on same curve");
+
+			// Point multiplication for Koblitz curves (using WTNAF) beats Shamir's trick
+			if (c is F2mCurve)
+			{
+				F2mCurve f2mCurve = (F2mCurve) c;
+				if (f2mCurve.IsKoblitz)
+				{
+					return P.Multiply(a).Add(Q.Multiply(b));
+				}
+			}
+
+			return ImplShamirsTrick(P, a, Q, b);
+		}
+
+		/*
+		* "Shamir's Trick", originally due to E. G. Straus
+		* (Addition chains of vectors. American Mathematical Monthly,
+		* 71(7):806-808, Aug./Sept. 1964)
+		*  
+		* Input: The points P, Q, scalar k = (km?, ... , k1, k0)
+		* and scalar l = (lm?, ... , l1, l0).
+		* Output: R = k * P + l * Q.
+		* 1: Z <- P + Q
+		* 2: R <- O
+		* 3: for i from m-1 down to 0 do
+		* 4:        R <- R + R        {point doubling}
+		* 5:        if (ki = 1) and (li = 0) then R <- R + P end if
+		* 6:        if (ki = 0) and (li = 1) then R <- R + Q end if
+		* 7:        if (ki = 1) and (li = 1) then R <- R + Z end if
+		* 8: end for
+		* 9: return R
+		*/
+		public static ECPoint ShamirsTrick(
+			ECPoint		P,
+			BigInteger	k,
+			ECPoint		Q,
+			BigInteger	l)
+		{
+			if (!P.Curve.Equals(Q.Curve))
+				throw new ArgumentException("P and Q must be on same curve");
+
+			return ImplShamirsTrick(P, k, Q, l);
+		}
+
+		private static ECPoint ImplShamirsTrick(ECPoint P, BigInteger k,
+			ECPoint Q, BigInteger l)
+		{
+			int m = System.Math.Max(k.BitLength, l.BitLength);
+			ECPoint Z = P.Add(Q);
+			ECPoint R = P.Curve.Infinity;
+
+			for (int i = m - 1; i >= 0; --i)
+			{
+				R = R.Twice();
+
+				if (k.TestBit(i))
+				{
+					if (l.TestBit(i))
+					{
+						R = R.Add(Z);
+					}
+					else
+					{
+						R = R.Add(P);
+					}
+				}
+				else
+				{
+					if (l.TestBit(i))
+					{
+						R = R.Add(Q);
+					}
+				}
+			}
+
+			return R;
+		}
+	}
+}
diff --git a/crypto/src/math/ec/ECCurve.cs b/crypto/src/math/ec/ECCurve.cs
new file mode 100644
index 000000000..396d42f28
--- /dev/null
+++ b/crypto/src/math/ec/ECCurve.cs
@@ -0,0 +1,651 @@
+using System;
+using System.Collections;
+
+using Org.BouncyCastle.Math.EC.Abc;
+
+namespace Org.BouncyCastle.Math.EC
+{
+    /// <remarks>Base class for an elliptic curve.</remarks>
+    public abstract class ECCurve
+    {
+        internal ECFieldElement a, b;
+
+        public abstract int FieldSize { get; }
+        public abstract ECFieldElement FromBigInteger(BigInteger x);
+        public abstract ECPoint CreatePoint(BigInteger x, BigInteger y, bool withCompression);
+        public abstract ECPoint Infinity { get; }
+
+        public ECFieldElement A
+        {
+            get { return a; }
+        }
+
+        public ECFieldElement B
+        {
+            get { return b; }
+        }
+
+        public override bool Equals(
+            object obj)
+        {
+            if (obj == this)
+                return true;
+
+            ECCurve other = obj as ECCurve;
+
+            if (other == null)
+                return false;
+
+            return Equals(other);
+        }
+
+        protected bool Equals(
+            ECCurve other)
+        {
+            return a.Equals(other.a) && b.Equals(other.b);
+        }
+
+        public override int GetHashCode()
+        {
+            return a.GetHashCode() ^ b.GetHashCode();
+        }
+
+        protected abstract ECPoint DecompressPoint(int yTilde, BigInteger X1);
+
+        /**
+         * Decode a point on this curve from its ASN.1 encoding. The different
+         * encodings are taken account of, including point compression for
+         * <code>F<sub>p</sub></code> (X9.62 s 4.2.1 pg 17).
+         * @return The decoded point.
+         */
+        public virtual ECPoint DecodePoint(byte[] encoded)
+        {
+            ECPoint p = null;
+            int expectedLength = (FieldSize + 7) / 8;
+
+            switch (encoded[0])
+            {
+                case 0x00: // infinity
+                {
+                    if (encoded.Length != 1)
+                        throw new ArgumentException("Incorrect length for infinity encoding", "encoded");
+
+                    p = Infinity;
+                    break;
+                }
+
+                case 0x02: // compressed
+                case 0x03: // compressed
+                {
+                    if (encoded.Length != (expectedLength + 1))
+                        throw new ArgumentException("Incorrect length for compressed encoding", "encoded");
+
+                    int yTilde = encoded[0] & 1;
+                    BigInteger X1 = new BigInteger(1, encoded, 1, expectedLength);
+
+                    p = DecompressPoint(yTilde, X1);
+                    break;
+                }
+
+                case 0x04: // uncompressed
+                case 0x06: // hybrid
+                case 0x07: // hybrid
+                {
+                    if (encoded.Length != (2 * expectedLength + 1))
+                        throw new ArgumentException("Incorrect length for uncompressed/hybrid encoding", "encoded");
+
+                    BigInteger X1 = new BigInteger(1, encoded, 1, expectedLength);
+                    BigInteger Y1 = new BigInteger(1, encoded, 1 + expectedLength, expectedLength);
+
+                    p = CreatePoint(X1, Y1, false);
+                    break;
+                }
+
+                default:
+                    throw new FormatException("Invalid point encoding " + encoded[0]);
+            }
+
+            return p;
+        }
+    }
+
+    /**
+     * Elliptic curve over Fp
+     */
+    public class FpCurve : ECCurve
+    {
+        private readonly BigInteger q;
+        private readonly FpPoint infinity;
+
+        public FpCurve(BigInteger q, BigInteger a, BigInteger b)
+        {
+            this.q = q;
+            this.a = FromBigInteger(a);
+            this.b = FromBigInteger(b);
+            this.infinity = new FpPoint(this, null, null);
+        }
+
+        public BigInteger Q
+        {
+            get { return q; }
+        }
+
+        public override ECPoint Infinity
+        {
+            get { return infinity; }
+        }
+
+        public override int FieldSize
+        {
+            get { return q.BitLength; }
+        }
+
+        public override ECFieldElement FromBigInteger(BigInteger x)
+        {
+            return new FpFieldElement(this.q, x);
+        }
+
+        public override ECPoint CreatePoint(
+            BigInteger	X1,
+            BigInteger	Y1,
+            bool		withCompression)
+        {
+            // TODO Validation of X1, Y1?
+            return new FpPoint(
+                this,
+                FromBigInteger(X1),
+                FromBigInteger(Y1),
+                withCompression);
+        }
+
+        protected override ECPoint DecompressPoint(
+            int			yTilde,
+            BigInteger	X1)
+        {
+            ECFieldElement x = FromBigInteger(X1);
+            ECFieldElement alpha = x.Multiply(x.Square().Add(a)).Add(b);
+            ECFieldElement beta = alpha.Sqrt();
+
+            //
+            // if we can't find a sqrt we haven't got a point on the
+            // curve - run!
+            //
+            if (beta == null)
+                throw new ArithmeticException("Invalid point compression");
+
+            BigInteger betaValue = beta.ToBigInteger();
+            int bit0 = betaValue.TestBit(0) ? 1 : 0;
+
+            if (bit0 != yTilde)
+            {
+                // Use the other root
+                beta = FromBigInteger(q.Subtract(betaValue));
+            }
+
+            return new FpPoint(this, x, beta, true);
+        }
+
+        public override bool Equals(
+            object obj)
+        {
+            if (obj == this)
+                return true;
+
+            FpCurve other = obj as FpCurve;
+
+            if (other == null)
+                return false;
+
+            return Equals(other);
+        }
+
+        protected bool Equals(
+            FpCurve other)
+        {
+            return base.Equals(other) && q.Equals(other.q);
+        }
+
+        public override int GetHashCode()
+        {
+            return base.GetHashCode() ^ q.GetHashCode();
+        }
+    }
+
+    /**
+     * Elliptic curves over F2m. The Weierstrass equation is given by
+     * <code>y<sup>2</sup> + xy = x<sup>3</sup> + ax<sup>2</sup> + b</code>.
+     */
+    public class F2mCurve : ECCurve
+    {
+        /**
+         * The exponent <code>m</code> of <code>F<sub>2<sup>m</sup></sub></code>.
+         */
+        private readonly int m;
+
+        /**
+         * TPB: The integer <code>k</code> where <code>x<sup>m</sup> +
+         * x<sup>k</sup> + 1</code> represents the reduction polynomial
+         * <code>f(z)</code>.<br/>
+         * PPB: The integer <code>k1</code> where <code>x<sup>m</sup> +
+         * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+         * represents the reduction polynomial <code>f(z)</code>.<br/>
+         */
+        private readonly int k1;
+
+        /**
+         * TPB: Always set to <code>0</code><br/>
+         * PPB: The integer <code>k2</code> where <code>x<sup>m</sup> +
+         * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+         * represents the reduction polynomial <code>f(z)</code>.<br/>
+         */
+        private readonly int k2;
+
+        /**
+         * TPB: Always set to <code>0</code><br/>
+         * PPB: The integer <code>k3</code> where <code>x<sup>m</sup> +
+         * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+         * represents the reduction polynomial <code>f(z)</code>.<br/>
+         */
+        private readonly int k3;
+
+        /**
+         * The order of the base point of the curve.
+         */
+        private readonly BigInteger n;
+
+        /**
+         * The cofactor of the curve.
+         */
+        private readonly BigInteger h;
+
+        /**
+         * The point at infinity on this curve.
+         */
+        private readonly F2mPoint infinity;
+
+        /**
+         * The parameter <code>&#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 override ECPoint DecompressPoint(
+            int			yTilde,
+            BigInteger	X1)
+        {
+            ECFieldElement xp = FromBigInteger(X1);
+            ECFieldElement yp = null;
+            if (xp.ToBigInteger().SignValue == 0)
+            {
+                yp = (F2mFieldElement)b;
+                for (int i = 0; i < m - 1; i++)
+                {
+                    yp = yp.Square();
+                }
+            }
+            else
+            {
+                ECFieldElement beta = xp.Add(a).Add(b.Multiply(xp.Square().Invert()));
+                ECFieldElement z = solveQuadradicEquation(beta);
+
+                if (z == null)
+                    throw new ArithmeticException("Invalid point compression");
+
+                int zBit = z.ToBigInteger().TestBit(0) ? 1 : 0;
+                if (zBit != yTilde)
+                {
+                    z = z.Add(FromBigInteger(BigInteger.One));
+                }
+
+                yp = xp.Multiply(z);
+            }
+
+            return new F2mPoint(this, xp, yp, true);
+        }
+
+        /**
+         * Solves a quadratic equation <code>z<sup>2</sup> + z = beta</code>(X9.62
+         * D.1.6) The other solution is <code>z + 1</code>.
+         *
+         * @param beta
+         *            The value to solve the qradratic equation for.
+         * @return the solution for <code>z<sup>2</sup> + z = beta</code> or
+         *         <code>null</code> if no solution exists.
+         */
+        private ECFieldElement solveQuadradicEquation(ECFieldElement beta)
+        {
+            if (beta.ToBigInteger().SignValue == 0)
+            {
+                return FromBigInteger(BigInteger.Zero);
+            }
+
+            ECFieldElement z = null;
+            ECFieldElement gamma = FromBigInteger(BigInteger.Zero);
+
+            while (gamma.ToBigInteger().SignValue == 0)
+            {
+                ECFieldElement t = FromBigInteger(new BigInteger(m, new Random()));
+                z = FromBigInteger(BigInteger.Zero);
+
+                ECFieldElement w = beta;
+                for (int i = 1; i <= m - 1; i++)
+                {
+                    ECFieldElement w2 = w.Square();
+                    z = z.Square().Add(w2.Multiply(t));
+                    w = w2.Add(beta);
+                }
+                if (w.ToBigInteger().SignValue != 0)
+                {
+                    return null;
+                }
+                gamma = z.Square().Add(z);
+            }
+            return z;
+        }
+
+        public override bool Equals(
+            object obj)
+        {
+            if (obj == this)
+                return true;
+
+            F2mCurve other = obj as F2mCurve;
+
+            if (other == null)
+                return false;
+
+            return Equals(other);
+        }
+
+        protected bool Equals(
+            F2mCurve other)
+        {
+            return m == other.m
+                && k1 == other.k1
+                && k2 == other.k2
+                && k3 == other.k3
+                && base.Equals(other);
+        }
+
+        public override int GetHashCode()
+        {
+            return base.GetHashCode() ^ m ^ k1 ^ k2 ^ k3;
+        }
+
+        public int M
+        {
+            get { return m; }
+        }
+
+        /**
+         * Return true if curve uses a Trinomial basis.
+         *
+         * @return true if curve Trinomial, false otherwise.
+         */
+        public bool IsTrinomial()
+        {
+            return k2 == 0 && k3 == 0;
+        }
+
+        public int K1
+        {
+            get { return k1; }
+        }
+
+        public int K2
+        {
+            get { return k2; }
+        }
+
+        public int K3
+        {
+            get { return k3; }
+        }
+
+        public BigInteger N
+        {
+            get { return n; }
+        }
+
+        public BigInteger H
+        {
+            get { return h; }
+        }
+    }
+}
diff --git a/crypto/src/math/ec/ECFieldElement.cs b/crypto/src/math/ec/ECFieldElement.cs
new file mode 100644
index 000000000..5235c6c0e
--- /dev/null
+++ b/crypto/src/math/ec/ECFieldElement.cs
@@ -0,0 +1,1253 @@
+using System;
+using System.Diagnostics;
+
+using Org.BouncyCastle.Utilities;
+
+namespace Org.BouncyCastle.Math.EC
+{
+	public abstract class ECFieldElement
+	{
+		public abstract BigInteger ToBigInteger();
+		public abstract string FieldName { get; }
+		public abstract int FieldSize { get; }
+		public abstract ECFieldElement Add(ECFieldElement b);
+		public abstract ECFieldElement Subtract(ECFieldElement b);
+		public abstract ECFieldElement Multiply(ECFieldElement b);
+		public abstract ECFieldElement Divide(ECFieldElement b);
+		public abstract ECFieldElement Negate();
+		public abstract ECFieldElement Square();
+		public abstract ECFieldElement Invert();
+		public abstract ECFieldElement Sqrt();
+
+		public override bool Equals(
+			object obj)
+		{
+			if (obj == this)
+				return true;
+
+			ECFieldElement other = obj as ECFieldElement;
+
+			if (other == null)
+				return false;
+
+			return Equals(other);
+		}
+
+		protected bool Equals(
+			ECFieldElement other)
+		{
+			return ToBigInteger().Equals(other.ToBigInteger());
+		}
+
+		public override int GetHashCode()
+		{
+			return ToBigInteger().GetHashCode();
+		}
+
+		public override string ToString()
+		{
+			return this.ToBigInteger().ToString(2);
+		}
+	}
+
+	public class FpFieldElement
+		: ECFieldElement
+	{
+		private readonly BigInteger q, x;
+
+		public FpFieldElement(
+			BigInteger	q,
+			BigInteger	x)
+		{
+			if (x.CompareTo(q) >= 0)
+				throw new ArgumentException("x value too large in field element");
+
+			this.q = q;
+			this.x = x;
+		}
+
+		public override BigInteger ToBigInteger()
+		{
+			return x;
+		}
+
+		/**
+		 * return the field name for this field.
+		 *
+		 * @return the string "Fp".
+		 */
+		public override string FieldName
+		{
+			get { return "Fp"; }
+		}
+
+		public override int FieldSize
+		{
+			get { return q.BitLength; }
+		}
+
+		public BigInteger Q
+		{
+			get { return q; }
+		}
+
+		public override ECFieldElement Add(
+			ECFieldElement b)
+		{
+			return new FpFieldElement(q, x.Add(b.ToBigInteger()).Mod(q));
+		}
+
+		public override ECFieldElement Subtract(
+			ECFieldElement b)
+		{
+			return new FpFieldElement(q, x.Subtract(b.ToBigInteger()).Mod(q));
+		}
+
+		public override ECFieldElement Multiply(
+			ECFieldElement b)
+		{
+			return new FpFieldElement(q, x.Multiply(b.ToBigInteger()).Mod(q));
+		}
+
+		public override ECFieldElement Divide(
+			ECFieldElement b)
+		{
+			return new FpFieldElement(q, x.Multiply(b.ToBigInteger().ModInverse(q)).Mod(q));
+		}
+
+		public override ECFieldElement Negate()
+		{
+			return new FpFieldElement(q, x.Negate().Mod(q));
+		}
+
+		public override ECFieldElement Square()
+		{
+			return new FpFieldElement(q, x.Multiply(x).Mod(q));
+		}
+
+		public override ECFieldElement Invert()
+		{
+			return new FpFieldElement(q, x.ModInverse(q));
+		}
+
+		// D.1.4 91
+		/**
+		 * return a sqrt root - the routine verifies that the calculation
+		 * returns the right value - if none exists it returns null.
+		 */
+		public override ECFieldElement Sqrt()
+		{
+			if (!q.TestBit(0))
+				throw Platform.CreateNotImplementedException("even value of q");
+
+			// p mod 4 == 3
+			if (q.TestBit(1))
+			{
+				// TODO Can this be optimised (inline the Square?)
+				// z = g^(u+1) + p, p = 4u + 3
+				ECFieldElement z = new FpFieldElement(q, x.ModPow(q.ShiftRight(2).Add(BigInteger.One), q));
+
+                return this.Equals(z.Square()) ? z : null;
+			}
+
+			// p mod 4 == 1
+			BigInteger qMinusOne = q.Subtract(BigInteger.One);
+
+			BigInteger legendreExponent = qMinusOne.ShiftRight(1);
+			if (!(x.ModPow(legendreExponent, q).Equals(BigInteger.One)))
+				return null;
+
+			BigInteger u = qMinusOne.ShiftRight(2);
+			BigInteger k = u.ShiftLeft(1).Add(BigInteger.One);
+
+			BigInteger Q = this.x;
+			BigInteger fourQ = Q.ShiftLeft(2).Mod(q);
+
+			BigInteger U, V;
+			do
+			{
+				Random rand = new Random();
+				BigInteger P;
+				do
+				{
+					P = new BigInteger(q.BitLength, rand);
+				}
+				while (P.CompareTo(q) >= 0
+					|| !(P.Multiply(P).Subtract(fourQ).ModPow(legendreExponent, q).Equals(qMinusOne)));
+
+				BigInteger[] result = fastLucasSequence(q, P, Q, k);
+				U = result[0];
+				V = result[1];
+
+				if (V.Multiply(V).Mod(q).Equals(fourQ))
+				{
+					// Integer division by 2, mod q
+					if (V.TestBit(0))
+					{
+						V = V.Add(q);
+					}
+
+					V = V.ShiftRight(1);
+
+					Debug.Assert(V.Multiply(V).Mod(q).Equals(x));
+
+					return new FpFieldElement(q, V);
+				}
+			}
+			while (U.Equals(BigInteger.One) || U.Equals(qMinusOne));
+
+			return null;
+
+
+//			BigInteger qMinusOne = q.Subtract(BigInteger.One);
+//
+//			BigInteger legendreExponent = qMinusOne.ShiftRight(1);
+//			if (!(x.ModPow(legendreExponent, q).Equals(BigInteger.One)))
+//				return null;
+//
+//			Random rand = new Random();
+//			BigInteger fourX = x.ShiftLeft(2);
+//
+//			BigInteger r;
+//			do
+//			{
+//				r = new BigInteger(q.BitLength, rand);
+//			}
+//			while (r.CompareTo(q) >= 0
+//				|| !(r.Multiply(r).Subtract(fourX).ModPow(legendreExponent, q).Equals(qMinusOne)));
+//
+//			BigInteger n1 = qMinusOne.ShiftRight(2);
+//			BigInteger n2 = n1.Add(BigInteger.One);
+//
+//			BigInteger wOne = WOne(r, x, q);
+//			BigInteger wSum = W(n1, wOne, q).Add(W(n2, wOne, q)).Mod(q);
+//			BigInteger twoR = r.ShiftLeft(1);
+//
+//			BigInteger root = twoR.ModPow(q.Subtract(BigInteger.Two), q)
+//				.Multiply(x).Mod(q)
+//				.Multiply(wSum).Mod(q);
+//
+//			return new FpFieldElement(q, root);
+		}
+
+//		private static BigInteger W(BigInteger n, BigInteger wOne, BigInteger p)
+//		{
+//			if (n.Equals(BigInteger.One))
+//				return wOne;
+//
+//			bool isEven = !n.TestBit(0);
+//			n = n.ShiftRight(1);
+//			if (isEven)
+//			{
+//				BigInteger w = W(n, wOne, p);
+//				return w.Multiply(w).Subtract(BigInteger.Two).Mod(p);
+//			}
+//			BigInteger w1 = W(n.Add(BigInteger.One), wOne, p);
+//			BigInteger w2 = W(n, wOne, p);
+//			return w1.Multiply(w2).Subtract(wOne).Mod(p);
+//		}
+//
+//		private BigInteger WOne(BigInteger r, BigInteger x, BigInteger p)
+//		{
+//			return r.Multiply(r).Multiply(x.ModPow(q.Subtract(BigInteger.Two), q)).Subtract(BigInteger.Two).Mod(p);
+//		}
+
+		private static BigInteger[] fastLucasSequence(
+			BigInteger	p,
+			BigInteger	P,
+			BigInteger	Q,
+			BigInteger	k)
+		{
+			// TODO Research and apply "common-multiplicand multiplication here"
+
+			int n = k.BitLength;
+			int s = k.GetLowestSetBit();
+
+			Debug.Assert(k.TestBit(s));
+
+			BigInteger Uh = BigInteger.One;
+			BigInteger Vl = BigInteger.Two;
+			BigInteger Vh = P;
+			BigInteger Ql = BigInteger.One;
+			BigInteger Qh = BigInteger.One;
+
+			for (int j = n - 1; j >= s + 1; --j)
+			{
+				Ql = Ql.Multiply(Qh).Mod(p);
+
+				if (k.TestBit(j))
+				{
+					Qh = Ql.Multiply(Q).Mod(p);
+					Uh = Uh.Multiply(Vh).Mod(p);
+					Vl = Vh.Multiply(Vl).Subtract(P.Multiply(Ql)).Mod(p);
+					Vh = Vh.Multiply(Vh).Subtract(Qh.ShiftLeft(1)).Mod(p);
+				}
+				else
+				{
+					Qh = Ql;
+					Uh = Uh.Multiply(Vl).Subtract(Ql).Mod(p);
+					Vh = Vh.Multiply(Vl).Subtract(P.Multiply(Ql)).Mod(p);
+					Vl = Vl.Multiply(Vl).Subtract(Ql.ShiftLeft(1)).Mod(p);
+				}
+			}
+
+			Ql = Ql.Multiply(Qh).Mod(p);
+			Qh = Ql.Multiply(Q).Mod(p);
+			Uh = Uh.Multiply(Vl).Subtract(Ql).Mod(p);
+			Vl = Vh.Multiply(Vl).Subtract(P.Multiply(Ql)).Mod(p);
+			Ql = Ql.Multiply(Qh).Mod(p);
+
+			for (int j = 1; j <= s; ++j)
+			{
+				Uh = Uh.Multiply(Vl).Mod(p);
+				Vl = Vl.Multiply(Vl).Subtract(Ql.ShiftLeft(1)).Mod(p);
+				Ql = Ql.Multiply(Ql).Mod(p);
+			}
+
+			return new BigInteger[]{ Uh, Vl };
+		}
+
+//		private static BigInteger[] verifyLucasSequence(
+//			BigInteger	p,
+//			BigInteger	P,
+//			BigInteger	Q,
+//			BigInteger	k)
+//		{
+//			BigInteger[] actual = fastLucasSequence(p, P, Q, k);
+//			BigInteger[] plus1 = fastLucasSequence(p, P, Q, k.Add(BigInteger.One));
+//			BigInteger[] plus2 = fastLucasSequence(p, P, Q, k.Add(BigInteger.Two));
+//
+//			BigInteger[] check = stepLucasSequence(p, P, Q, actual, plus1);
+//
+//			Debug.Assert(check[0].Equals(plus2[0]));
+//			Debug.Assert(check[1].Equals(plus2[1]));
+//
+//			return actual;
+//		}
+//
+//		private static BigInteger[] stepLucasSequence(
+//			BigInteger		p,
+//			BigInteger		P,
+//			BigInteger		Q,
+//			BigInteger[]	backTwo,
+//			BigInteger[]	backOne)
+//		{
+//			return new BigInteger[]
+//			{
+//				P.Multiply(backOne[0]).Subtract(Q.Multiply(backTwo[0])).Mod(p),
+//				P.Multiply(backOne[1]).Subtract(Q.Multiply(backTwo[1])).Mod(p)
+//			};
+//		}
+
+		public override bool Equals(
+			object obj)
+		{
+			if (obj == this)
+				return true;
+
+			FpFieldElement other = obj as FpFieldElement;
+
+			if (other == null)
+				return false;
+
+			return Equals(other);
+		}
+
+		protected bool Equals(
+			FpFieldElement other)
+		{
+			return q.Equals(other.q) && base.Equals(other);
+		}
+
+		public override int GetHashCode()
+		{
+			return q.GetHashCode() ^ base.GetHashCode();
+		}
+	}
+
+//	/**
+//	 * Class representing the Elements of the finite field
+//	 * <code>F<sub>2<sup>m</sup></sub></code> in polynomial basis (PB)
+//	 * representation. Both trinomial (Tpb) and pentanomial (Ppb) polynomial
+//	 * basis representations are supported. Gaussian normal basis (GNB)
+//	 * representation is not supported.
+//	 */
+//	public class F2mFieldElement
+//		: ECFieldElement
+//	{
+//		/**
+//		 * Indicates gaussian normal basis representation (GNB). Number chosen
+//		 * according to X9.62. GNB is not implemented at present.
+//		 */
+//		public const int Gnb = 1;
+//
+//		/**
+//		 * Indicates trinomial basis representation (Tpb). Number chosen
+//		 * according to X9.62.
+//		 */
+//		public const int Tpb = 2;
+//
+//		/**
+//		 * Indicates pentanomial basis representation (Ppb). Number chosen
+//		 * according to X9.62.
+//		 */
+//		public const int Ppb = 3;
+//
+//		/**
+//		 * Tpb or Ppb.
+//		 */
+//		private int representation;
+//
+//		/**
+//		 * The exponent <code>m</code> of <code>F<sub>2<sup>m</sup></sub></code>.
+//		 */
+//		private int m;
+//
+//		/**
+//		 * Tpb: The integer <code>k</code> where <code>x<sup>m</sup> +
+//		 * x<sup>k</sup> + 1</code> represents the reduction polynomial
+//		 * <code>f(z)</code>.<br/>
+//		 * Ppb: The integer <code>k1</code> where <code>x<sup>m</sup> +
+//		 * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+//		 * represents the reduction polynomial <code>f(z)</code>.<br/>
+//		 */
+//		private int k1;
+//
+//		/**
+//		 * Tpb: Always set to <code>0</code><br/>
+//		 * Ppb: The integer <code>k2</code> where <code>x<sup>m</sup> +
+//		 * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+//		 * represents the reduction polynomial <code>f(z)</code>.<br/>
+//		 */
+//		private int k2;
+//
+//		/**
+//			* Tpb: Always set to <code>0</code><br/>
+//			* Ppb: The integer <code>k3</code> where <code>x<sup>m</sup> +
+//			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+//			* represents the reduction polynomial <code>f(z)</code>.<br/>
+//			*/
+//		private int k3;
+//
+//		/**
+//			* Constructor for Ppb.
+//			* @param m  The exponent <code>m</code> of
+//			* <code>F<sub>2<sup>m</sup></sub></code>.
+//			* @param k1 The integer <code>k1</code> where <code>x<sup>m</sup> +
+//			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+//			* represents the reduction polynomial <code>f(z)</code>.
+//			* @param k2 The integer <code>k2</code> where <code>x<sup>m</sup> +
+//			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+//			* represents the reduction polynomial <code>f(z)</code>.
+//			* @param k3 The integer <code>k3</code> where <code>x<sup>m</sup> +
+//			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+//			* represents the reduction polynomial <code>f(z)</code>.
+//			* @param x The BigInteger representing the value of the field element.
+//			*/
+//		public F2mFieldElement(
+//			int			m,
+//			int			k1,
+//			int			k2,
+//			int			k3,
+//			BigInteger	x)
+//			: base(x)
+//		{
+//			if ((k2 == 0) && (k3 == 0))
+//			{
+//				this.representation = Tpb;
+//			}
+//			else
+//			{
+//				if (k2 >= k3)
+//					throw new ArgumentException("k2 must be smaller than k3");
+//				if (k2 <= 0)
+//					throw new ArgumentException("k2 must be larger than 0");
+//
+//				this.representation = Ppb;
+//			}
+//
+//			if (x.SignValue < 0)
+//				throw new ArgumentException("x value cannot be negative");
+//
+//			this.m = m;
+//			this.k1 = k1;
+//			this.k2 = k2;
+//			this.k3 = k3;
+//		}
+//
+//		/**
+//			* Constructor for Tpb.
+//			* @param m  The exponent <code>m</code> of
+//			* <code>F<sub>2<sup>m</sup></sub></code>.
+//			* @param k The integer <code>k</code> where <code>x<sup>m</sup> +
+//			* x<sup>k</sup> + 1</code> represents the reduction
+//			* polynomial <code>f(z)</code>.
+//			* @param x The BigInteger representing the value of the field element.
+//			*/
+//		public F2mFieldElement(
+//			int			m,
+//			int			k,
+//			BigInteger	x)
+//			: this(m, k, 0, 0, x)
+//		{
+//			// Set k1 to k, and set k2 and k3 to 0
+//		}
+//
+//		public override string FieldName
+//		{
+//			get { return "F2m"; }
+//		}
+//
+//		/**
+//		* Checks, if the ECFieldElements <code>a</code> and <code>b</code>
+//		* are elements of the same field <code>F<sub>2<sup>m</sup></sub></code>
+//		* (having the same representation).
+//		* @param a field element.
+//		* @param b field element to be compared.
+//		* @throws ArgumentException if <code>a</code> and <code>b</code>
+//		* are not elements of the same field
+//		* <code>F<sub>2<sup>m</sup></sub></code> (having the same
+//		* representation).
+//		*/
+//		public static void CheckFieldElements(
+//			ECFieldElement	a,
+//			ECFieldElement	b)
+//		{
+//			if (!(a is F2mFieldElement) || !(b is F2mFieldElement))
+//			{
+//				throw new ArgumentException("Field elements are not "
+//					+ "both instances of F2mFieldElement");
+//			}
+//
+//			if ((a.x.SignValue < 0) || (b.x.SignValue < 0))
+//			{
+//				throw new ArgumentException(
+//					"x value may not be negative");
+//			}
+//
+//			F2mFieldElement aF2m = (F2mFieldElement)a;
+//			F2mFieldElement bF2m = (F2mFieldElement)b;
+//
+//			if ((aF2m.m != bF2m.m) || (aF2m.k1 != bF2m.k1)
+//				|| (aF2m.k2 != bF2m.k2) || (aF2m.k3 != bF2m.k3))
+//			{
+//				throw new ArgumentException("Field elements are not "
+//					+ "elements of the same field F2m");
+//			}
+//
+//			if (aF2m.representation != bF2m.representation)
+//			{
+//				// Should never occur
+//				throw new ArgumentException(
+//					"One of the field "
+//					+ "elements are not elements has incorrect representation");
+//			}
+//		}
+//
+//		/**
+//			* Computes <code>z * a(z) mod f(z)</code>, where <code>f(z)</code> is
+//			* the reduction polynomial of <code>this</code>.
+//			* @param a The polynomial <code>a(z)</code> to be multiplied by
+//			* <code>z mod f(z)</code>.
+//			* @return <code>z * a(z) mod f(z)</code>
+//			*/
+//		private BigInteger multZModF(
+//			BigInteger a)
+//		{
+//			// Left-shift of a(z)
+//			BigInteger az = a.ShiftLeft(1);
+//			if (az.TestBit(this.m))
+//			{
+//				// If the coefficient of z^m in a(z) Equals 1, reduction
+//				// modulo f(z) is performed: Add f(z) to to a(z):
+//				// Step 1: Unset mth coeffient of a(z)
+//				az = az.ClearBit(this.m);
+//
+//				// Step 2: Add r(z) to a(z), where r(z) is defined as
+//				// f(z) = z^m + r(z), and k1, k2, k3 are the positions of
+//				// the non-zero coefficients in r(z)
+//				az = az.FlipBit(0);
+//				az = az.FlipBit(this.k1);
+//				if (this.representation == Ppb)
+//				{
+//					az = az.FlipBit(this.k2);
+//					az = az.FlipBit(this.k3);
+//				}
+//			}
+//			return az;
+//		}
+//
+//		public override ECFieldElement Add(
+//			ECFieldElement b)
+//		{
+//			// No check performed here for performance reasons. Instead the
+//			// elements involved are checked in ECPoint.F2m
+//			// checkFieldElements(this, b);
+//			if (b.x.SignValue == 0)
+//				return this;
+//
+//			return new F2mFieldElement(this.m, this.k1, this.k2, this.k3, this.x.Xor(b.x));
+//		}
+//
+//		public override ECFieldElement Subtract(
+//			ECFieldElement b)
+//		{
+//			// Addition and subtraction are the same in F2m
+//			return Add(b);
+//		}
+//
+//		public override ECFieldElement Multiply(
+//			ECFieldElement b)
+//		{
+//			// Left-to-right shift-and-add field multiplication in F2m
+//			// Input: Binary polynomials a(z) and b(z) of degree at most m-1
+//			// Output: c(z) = a(z) * b(z) mod f(z)
+//
+//			// No check performed here for performance reasons. Instead the
+//			// elements involved are checked in ECPoint.F2m
+//			// checkFieldElements(this, b);
+//			BigInteger az = this.x;
+//			BigInteger bz = b.x;
+//			BigInteger cz;
+//
+//			// Compute c(z) = a(z) * b(z) mod f(z)
+//			if (az.TestBit(0))
+//			{
+//				cz = bz;
+//			}
+//			else
+//			{
+//				cz = BigInteger.Zero;
+//			}
+//
+//			for (int i = 1; i < this.m; i++)
+//			{
+//				// b(z) := z * b(z) mod f(z)
+//				bz = multZModF(bz);
+//
+//				if (az.TestBit(i))
+//				{
+//					// If the coefficient of x^i in a(z) Equals 1, b(z) is added
+//					// to c(z)
+//					cz = cz.Xor(bz);
+//				}
+//			}
+//			return new F2mFieldElement(m, this.k1, this.k2, this.k3, cz);
+//		}
+//
+//
+//		public override ECFieldElement Divide(
+//			ECFieldElement b)
+//		{
+//			// There may be more efficient implementations
+//			ECFieldElement bInv = b.Invert();
+//			return Multiply(bInv);
+//		}
+//
+//		public override ECFieldElement Negate()
+//		{
+//			// -x == x holds for all x in F2m
+//			return this;
+//		}
+//
+//		public override ECFieldElement Square()
+//		{
+//			// Naive implementation, can probably be speeded up using modular
+//			// reduction
+//			return Multiply(this);
+//		}
+//
+//		public override ECFieldElement Invert()
+//		{
+//			// Inversion in F2m using the extended Euclidean algorithm
+//			// Input: A nonzero polynomial a(z) of degree at most m-1
+//			// Output: a(z)^(-1) mod f(z)
+//
+//			// u(z) := a(z)
+//			BigInteger uz = this.x;
+//			if (uz.SignValue <= 0)
+//			{
+//				throw new ArithmeticException("x is zero or negative, " +
+//					"inversion is impossible");
+//			}
+//
+//			// v(z) := f(z)
+//			BigInteger vz = BigInteger.One.ShiftLeft(m);
+//			vz = vz.SetBit(0);
+//			vz = vz.SetBit(this.k1);
+//			if (this.representation == Ppb)
+//			{
+//				vz = vz.SetBit(this.k2);
+//				vz = vz.SetBit(this.k3);
+//			}
+//
+//			// g1(z) := 1, g2(z) := 0
+//			BigInteger g1z = BigInteger.One;
+//			BigInteger g2z = BigInteger.Zero;
+//
+//			// while u != 1
+//			while (uz.SignValue != 0)
+//			{
+//				// j := deg(u(z)) - deg(v(z))
+//				int j = uz.BitLength - vz.BitLength;
+//
+//				// If j < 0 then: u(z) <-> v(z), g1(z) <-> g2(z), j := -j
+//				if (j < 0)
+//				{
+//					BigInteger uzCopy = uz;
+//					uz = vz;
+//					vz = uzCopy;
+//
+//					BigInteger g1zCopy = g1z;
+//					g1z = g2z;
+//					g2z = g1zCopy;
+//
+//					j = -j;
+//				}
+//
+//				// u(z) := u(z) + z^j * v(z)
+//				// Note, that no reduction modulo f(z) is required, because
+//				// deg(u(z) + z^j * v(z)) <= max(deg(u(z)), j + deg(v(z)))
+//				// = max(deg(u(z)), deg(u(z)) - deg(v(z)) + deg(v(z))
+//				// = deg(u(z))
+//				uz = uz.Xor(vz.ShiftLeft(j));
+//
+//				// g1(z) := g1(z) + z^j * g2(z)
+//				g1z = g1z.Xor(g2z.ShiftLeft(j));
+//				//                if (g1z.BitLength() > this.m) {
+//				//                    throw new ArithmeticException(
+//				//                            "deg(g1z) >= m, g1z = " + g1z.ToString(2));
+//				//                }
+//			}
+//			return new F2mFieldElement(this.m, this.k1, this.k2, this.k3, g2z);
+//		}
+//
+//		public override ECFieldElement Sqrt()
+//		{
+//			throw new ArithmeticException("Not implemented");
+//		}
+//
+//		/**
+//			* @return the representation of the field
+//			* <code>F<sub>2<sup>m</sup></sub></code>, either of
+//			* {@link F2mFieldElement.Tpb} (trinomial
+//			* basis representation) or
+//			* {@link F2mFieldElement.Ppb} (pentanomial
+//			* basis representation).
+//			*/
+//		public int Representation
+//		{
+//			get { return this.representation; }
+//		}
+//
+//		/**
+//			* @return the degree <code>m</code> of the reduction polynomial
+//			* <code>f(z)</code>.
+//			*/
+//		public int M
+//		{
+//			get { return this.m; }
+//		}
+//
+//		/**
+//			* @return Tpb: The integer <code>k</code> where <code>x<sup>m</sup> +
+//			* x<sup>k</sup> + 1</code> represents the reduction polynomial
+//			* <code>f(z)</code>.<br/>
+//			* Ppb: The integer <code>k1</code> where <code>x<sup>m</sup> +
+//			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+//			* represents the reduction polynomial <code>f(z)</code>.<br/>
+//			*/
+//		public int K1
+//		{
+//			get { return this.k1; }
+//		}
+//
+//		/**
+//			* @return Tpb: Always returns <code>0</code><br/>
+//			* Ppb: The integer <code>k2</code> where <code>x<sup>m</sup> +
+//			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+//			* represents the reduction polynomial <code>f(z)</code>.<br/>
+//			*/
+//		public int K2
+//		{
+//			get { return this.k2; }
+//		}
+//
+//		/**
+//			* @return Tpb: Always set to <code>0</code><br/>
+//			* Ppb: The integer <code>k3</code> where <code>x<sup>m</sup> +
+//			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+//			* represents the reduction polynomial <code>f(z)</code>.<br/>
+//			*/
+//		public int K3
+//		{
+//			get { return this.k3; }
+//		}
+//
+//		public override bool Equals(
+//			object obj)
+//		{
+//			if (obj == this)
+//				return true;
+//
+//			F2mFieldElement other = obj as F2mFieldElement;
+//
+//			if (other == null)
+//				return false;
+//
+//			return Equals(other);
+//		}
+//
+//		protected bool Equals(
+//			F2mFieldElement other)
+//		{
+//			return m == other.m
+//				&& k1 == other.k1
+//				&& k2 == other.k2
+//				&& k3 == other.k3
+//				&& representation == other.representation
+//				&& base.Equals(other);
+//		}
+//
+//		public override int GetHashCode()
+//		{
+//			return base.GetHashCode() ^ m ^ k1 ^ k2 ^ k3;
+//		}
+//	}
+
+	/**
+	 * Class representing the Elements of the finite field
+	 * <code>F<sub>2<sup>m</sup></sub></code> in polynomial basis (PB)
+	 * representation. Both trinomial (Tpb) and pentanomial (Ppb) polynomial
+	 * basis representations are supported. Gaussian normal basis (GNB)
+	 * representation is not supported.
+	 */
+	public class F2mFieldElement
+		: ECFieldElement
+	{
+		/**
+		 * Indicates gaussian normal basis representation (GNB). Number chosen
+		 * according to X9.62. GNB is not implemented at present.
+		 */
+		public const int Gnb = 1;
+
+		/**
+		 * Indicates trinomial basis representation (Tpb). Number chosen
+		 * according to X9.62.
+		 */
+		public const int Tpb = 2;
+
+		/**
+		 * Indicates pentanomial basis representation (Ppb). Number chosen
+		 * according to X9.62.
+		 */
+		public const int Ppb = 3;
+
+		/**
+		 * Tpb or Ppb.
+		 */
+		private int representation;
+
+		/**
+		 * The exponent <code>m</code> of <code>F<sub>2<sup>m</sup></sub></code>.
+		 */
+		private int m;
+
+		/**
+		 * Tpb: The integer <code>k</code> where <code>x<sup>m</sup> +
+		 * x<sup>k</sup> + 1</code> represents the reduction polynomial
+		 * <code>f(z)</code>.<br/>
+		 * Ppb: The integer <code>k1</code> where <code>x<sup>m</sup> +
+		 * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+		 * represents the reduction polynomial <code>f(z)</code>.<br/>
+		 */
+		private int k1;
+
+		/**
+		 * Tpb: Always set to <code>0</code><br/>
+		 * Ppb: The integer <code>k2</code> where <code>x<sup>m</sup> +
+		 * x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+		 * represents the reduction polynomial <code>f(z)</code>.<br/>
+		 */
+		private int k2;
+
+		/**
+			* Tpb: Always set to <code>0</code><br/>
+			* Ppb: The integer <code>k3</code> where <code>x<sup>m</sup> +
+			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+			* represents the reduction polynomial <code>f(z)</code>.<br/>
+			*/
+		private int k3;
+
+		/**
+		 * The <code>IntArray</code> holding the bits.
+		 */
+		private IntArray x;
+
+		/**
+		 * The number of <code>int</code>s required to hold <code>m</code> bits.
+		 */
+		private readonly int t;
+
+		/**
+			* Constructor for Ppb.
+			* @param m  The exponent <code>m</code> of
+			* <code>F<sub>2<sup>m</sup></sub></code>.
+			* @param k1 The integer <code>k1</code> where <code>x<sup>m</sup> +
+			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+			* represents the reduction polynomial <code>f(z)</code>.
+			* @param k2 The integer <code>k2</code> where <code>x<sup>m</sup> +
+			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+			* represents the reduction polynomial <code>f(z)</code>.
+			* @param k3 The integer <code>k3</code> where <code>x<sup>m</sup> +
+			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+			* represents the reduction polynomial <code>f(z)</code>.
+			* @param x The BigInteger representing the value of the field element.
+			*/
+		public F2mFieldElement(
+			int			m,
+			int			k1,
+			int			k2,
+			int			k3,
+			BigInteger	x)
+		{
+			// t = m / 32 rounded up to the next integer
+			this.t = (m + 31) >> 5;
+			this.x = new IntArray(x, t);
+
+			if ((k2 == 0) && (k3 == 0))
+			{
+				this.representation = Tpb;
+			}
+			else
+			{
+				if (k2 >= k3)
+					throw new ArgumentException("k2 must be smaller than k3");
+				if (k2 <= 0)
+					throw new ArgumentException("k2 must be larger than 0");
+
+				this.representation = Ppb;
+			}
+
+			if (x.SignValue < 0)
+				throw new ArgumentException("x value cannot be negative");
+
+			this.m = m;
+			this.k1 = k1;
+			this.k2 = k2;
+			this.k3 = k3;
+		}
+
+		/**
+			* Constructor for Tpb.
+			* @param m  The exponent <code>m</code> of
+			* <code>F<sub>2<sup>m</sup></sub></code>.
+			* @param k The integer <code>k</code> where <code>x<sup>m</sup> +
+			* x<sup>k</sup> + 1</code> represents the reduction
+			* polynomial <code>f(z)</code>.
+			* @param x The BigInteger representing the value of the field element.
+			*/
+		public F2mFieldElement(
+			int			m,
+			int			k,
+			BigInteger	x)
+			: this(m, k, 0, 0, x)
+		{
+			// Set k1 to k, and set k2 and k3 to 0
+		}
+
+		private F2mFieldElement(int m, int k1, int k2, int k3, IntArray x)
+		{
+			t = (m + 31) >> 5;
+			this.x = x;
+			this.m = m;
+			this.k1 = k1;
+			this.k2 = k2;
+			this.k3 = k3;
+
+			if ((k2 == 0) && (k3 == 0))
+			{
+				this.representation = Tpb;
+			}
+			else
+			{
+				this.representation = Ppb;
+			}
+		}
+
+		public override BigInteger ToBigInteger()
+		{
+			return x.ToBigInteger();
+		}
+
+		public override string FieldName
+		{
+			get { return "F2m"; }
+		}
+
+		public override int FieldSize
+		{
+			get { return m; }
+		}
+
+		/**
+		* Checks, if the ECFieldElements <code>a</code> and <code>b</code>
+		* are elements of the same field <code>F<sub>2<sup>m</sup></sub></code>
+		* (having the same representation).
+		* @param a field element.
+		* @param b field element to be compared.
+		* @throws ArgumentException if <code>a</code> and <code>b</code>
+		* are not elements of the same field
+		* <code>F<sub>2<sup>m</sup></sub></code> (having the same
+		* representation).
+		*/
+		public static void CheckFieldElements(
+			ECFieldElement	a,
+			ECFieldElement	b)
+		{
+			if (!(a is F2mFieldElement) || !(b is F2mFieldElement))
+			{
+				throw new ArgumentException("Field elements are not "
+					+ "both instances of F2mFieldElement");
+			}
+
+			F2mFieldElement aF2m = (F2mFieldElement)a;
+			F2mFieldElement bF2m = (F2mFieldElement)b;
+
+			if ((aF2m.m != bF2m.m) || (aF2m.k1 != bF2m.k1)
+				|| (aF2m.k2 != bF2m.k2) || (aF2m.k3 != bF2m.k3))
+			{
+				throw new ArgumentException("Field elements are not "
+					+ "elements of the same field F2m");
+			}
+
+			if (aF2m.representation != bF2m.representation)
+			{
+				// Should never occur
+				throw new ArgumentException(
+					"One of the field "
+					+ "elements are not elements has incorrect representation");
+			}
+		}
+
+		public override ECFieldElement Add(
+			ECFieldElement b)
+		{
+			// No check performed here for performance reasons. Instead the
+			// elements involved are checked in ECPoint.F2m
+			// checkFieldElements(this, b);
+			IntArray iarrClone = (IntArray) this.x.Copy();
+			F2mFieldElement bF2m = (F2mFieldElement) b;
+			iarrClone.AddShifted(bF2m.x, 0);
+			return new F2mFieldElement(m, k1, k2, k3, iarrClone);
+		}
+
+		public override ECFieldElement Subtract(
+			ECFieldElement b)
+		{
+			// Addition and subtraction are the same in F2m
+			return Add(b);
+		}
+
+		public override ECFieldElement Multiply(
+			ECFieldElement b)
+		{
+			// Right-to-left comb multiplication in the IntArray
+			// Input: Binary polynomials a(z) and b(z) of degree at most m-1
+			// Output: c(z) = a(z) * b(z) mod f(z)
+
+			// No check performed here for performance reasons. Instead the
+			// elements involved are checked in ECPoint.F2m
+			// checkFieldElements(this, b);
+			F2mFieldElement bF2m = (F2mFieldElement) b;
+			IntArray mult = x.Multiply(bF2m.x, m);
+			mult.Reduce(m, new int[]{k1, k2, k3});
+			return new F2mFieldElement(m, k1, k2, k3, mult);
+		}
+
+		public override ECFieldElement Divide(
+			ECFieldElement b)
+		{
+			// There may be more efficient implementations
+			ECFieldElement bInv = b.Invert();
+			return Multiply(bInv);
+		}
+
+		public override ECFieldElement Negate()
+		{
+			// -x == x holds for all x in F2m
+			return this;
+		}
+
+		public override ECFieldElement Square()
+		{
+			IntArray squared = x.Square(m);
+			squared.Reduce(m, new int[]{k1, k2, k3});
+			return new F2mFieldElement(m, k1, k2, k3, squared);
+		}
+
+		public override ECFieldElement Invert()
+		{
+			// Inversion in F2m using the extended Euclidean algorithm
+			// Input: A nonzero polynomial a(z) of degree at most m-1
+			// Output: a(z)^(-1) mod f(z)
+
+			// u(z) := a(z)
+            IntArray uz = (IntArray)this.x.Copy();
+
+			// v(z) := f(z)
+			IntArray vz = new IntArray(t);
+			vz.SetBit(m);
+			vz.SetBit(0);
+			vz.SetBit(this.k1);
+			if (this.representation == Ppb)
+			{
+				vz.SetBit(this.k2);
+				vz.SetBit(this.k3);
+			}
+
+			// g1(z) := 1, g2(z) := 0
+			IntArray g1z = new IntArray(t);
+			g1z.SetBit(0);
+			IntArray g2z = new IntArray(t);
+
+			// while u != 0
+			while (uz.GetUsedLength() > 0)
+//            while (uz.bitLength() > 1)
+			{
+				// j := deg(u(z)) - deg(v(z))
+				int j = uz.BitLength - vz.BitLength;
+
+				// If j < 0 then: u(z) <-> v(z), g1(z) <-> g2(z), j := -j
+				if (j < 0)
+				{
+                    IntArray uzCopy = uz;
+					uz = vz;
+					vz = uzCopy;
+
+                    IntArray g1zCopy = g1z;
+					g1z = g2z;
+					g2z = g1zCopy;
+
+					j = -j;
+				}
+
+				// u(z) := u(z) + z^j * v(z)
+				// Note, that no reduction modulo f(z) is required, because
+				// deg(u(z) + z^j * v(z)) <= max(deg(u(z)), j + deg(v(z)))
+				// = max(deg(u(z)), deg(u(z)) - deg(v(z)) + deg(v(z))
+				// = deg(u(z))
+				// uz = uz.xor(vz.ShiftLeft(j));
+				// jInt = n / 32
+				int jInt = j >> 5;
+				// jInt = n % 32
+				int jBit = j & 0x1F;
+				IntArray vzShift = vz.ShiftLeft(jBit);
+				uz.AddShifted(vzShift, jInt);
+
+				// g1(z) := g1(z) + z^j * g2(z)
+//                g1z = g1z.xor(g2z.ShiftLeft(j));
+				IntArray g2zShift = g2z.ShiftLeft(jBit);
+				g1z.AddShifted(g2zShift, jInt);
+			}
+			return new F2mFieldElement(this.m, this.k1, this.k2, this.k3, g2z);
+		}
+
+		public override ECFieldElement Sqrt()
+		{
+			throw new ArithmeticException("Not implemented");
+		}
+
+		/**
+			* @return the representation of the field
+			* <code>F<sub>2<sup>m</sup></sub></code>, either of
+			* {@link F2mFieldElement.Tpb} (trinomial
+			* basis representation) or
+			* {@link F2mFieldElement.Ppb} (pentanomial
+			* basis representation).
+			*/
+		public int Representation
+		{
+			get { return this.representation; }
+		}
+
+		/**
+			* @return the degree <code>m</code> of the reduction polynomial
+			* <code>f(z)</code>.
+			*/
+		public int M
+		{
+			get { return this.m; }
+		}
+
+		/**
+			* @return Tpb: The integer <code>k</code> where <code>x<sup>m</sup> +
+			* x<sup>k</sup> + 1</code> represents the reduction polynomial
+			* <code>f(z)</code>.<br/>
+			* Ppb: The integer <code>k1</code> where <code>x<sup>m</sup> +
+			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+			* represents the reduction polynomial <code>f(z)</code>.<br/>
+			*/
+		public int K1
+		{
+			get { return this.k1; }
+		}
+
+		/**
+			* @return Tpb: Always returns <code>0</code><br/>
+			* Ppb: The integer <code>k2</code> where <code>x<sup>m</sup> +
+			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+			* represents the reduction polynomial <code>f(z)</code>.<br/>
+			*/
+		public int K2
+		{
+			get { return this.k2; }
+		}
+
+		/**
+			* @return Tpb: Always set to <code>0</code><br/>
+			* Ppb: The integer <code>k3</code> where <code>x<sup>m</sup> +
+			* x<sup>k3</sup> + x<sup>k2</sup> + x<sup>k1</sup> + 1</code>
+			* represents the reduction polynomial <code>f(z)</code>.<br/>
+			*/
+		public int K3
+		{
+			get { return this.k3; }
+		}
+
+		public override bool Equals(
+			object obj)
+		{
+			if (obj == this)
+				return true;
+
+			F2mFieldElement other = obj as F2mFieldElement;
+
+			if (other == null)
+				return false;
+
+			return Equals(other);
+		}
+
+		protected bool Equals(
+			F2mFieldElement other)
+		{
+			return m == other.m
+				&& k1 == other.k1
+				&& k2 == other.k2
+				&& k3 == other.k3
+				&& representation == other.representation
+				&& base.Equals(other);
+		}
+
+		public override int GetHashCode()
+		{
+			return m.GetHashCode()
+				^	k1.GetHashCode()
+				^	k2.GetHashCode()
+				^	k3.GetHashCode()
+				^	representation.GetHashCode()
+				^	base.GetHashCode();
+		}
+	}
+}
diff --git a/crypto/src/math/ec/ECPoint.cs b/crypto/src/math/ec/ECPoint.cs
new file mode 100644
index 000000000..6a06be4d1
--- /dev/null
+++ b/crypto/src/math/ec/ECPoint.cs
@@ -0,0 +1,572 @@
+using System;
+using System.Collections;
+using System.Diagnostics;
+
+using Org.BouncyCastle.Asn1.X9;
+
+using Org.BouncyCastle.Math.EC.Multiplier;
+
+namespace Org.BouncyCastle.Math.EC
+{
+    /**
+     * base class for points on elliptic curves.
+     */
+    public abstract class ECPoint
+    {
+        internal readonly ECCurve			curve;
+        internal readonly ECFieldElement	x, y;
+        internal readonly bool				withCompression;
+        internal ECMultiplier				multiplier = null;
+        internal PreCompInfo				preCompInfo = null;
+
+        protected internal ECPoint(
+            ECCurve			curve,
+            ECFieldElement	x,
+            ECFieldElement	y,
+            bool			withCompression)
+        {
+            if (curve == null)
+                throw new ArgumentNullException("curve");
+
+            this.curve = curve;
+            this.x = x;
+            this.y = y;
+            this.withCompression = withCompression;
+        }
+
+        public ECCurve Curve
+        {
+            get { return curve; }
+        }
+
+        public ECFieldElement X
+        {
+            get { return x; }
+        }
+
+        public ECFieldElement Y
+        {
+            get { return y; }
+        }
+
+        public bool IsInfinity
+        {
+            get { return x == null && y == null; }
+        }
+
+        public bool IsCompressed
+        {
+            get { return withCompression; }
+        }
+
+        public override bool Equals(
+            object obj)
+        {
+            if (obj == this)
+                return true;
+
+            ECPoint o = obj as ECPoint;
+
+            if (o == null)
+                return false;
+
+            if (this.IsInfinity)
+                return o.IsInfinity;
+
+            return x.Equals(o.x) && y.Equals(o.y);
+        }
+
+        public override int GetHashCode()
+        {
+            if (this.IsInfinity)
+                return 0;
+
+            return x.GetHashCode() ^ y.GetHashCode();
+        }
+
+//		/**
+//		 * Mainly for testing. Explicitly set the <code>ECMultiplier</code>.
+//		 * @param multiplier The <code>ECMultiplier</code> to be used to multiply
+//		 * this <code>ECPoint</code>.
+//		 */
+//		internal void SetECMultiplier(
+//			ECMultiplier multiplier)
+//		{
+//			this.multiplier = multiplier;
+//		}
+
+        /**
+         * Sets the <code>PreCompInfo</code>. Used by <code>ECMultiplier</code>s
+         * to save the precomputation for this <code>ECPoint</code> to store the
+         * precomputation result for use by subsequent multiplication.
+         * @param preCompInfo The values precomputed by the
+         * <code>ECMultiplier</code>.
+         */
+        internal void SetPreCompInfo(
+            PreCompInfo preCompInfo)
+        {
+            this.preCompInfo = preCompInfo;
+        }
+
+        public virtual byte[] GetEncoded()
+        {
+            return GetEncoded(withCompression);
+        }
+
+        public abstract byte[] GetEncoded(bool compressed);
+
+        public abstract ECPoint Add(ECPoint b);
+        public abstract ECPoint Subtract(ECPoint b);
+        public abstract ECPoint Negate();
+        public abstract ECPoint Twice();
+        public abstract ECPoint Multiply(BigInteger b);
+
+        /**
+        * Sets the appropriate <code>ECMultiplier</code>, unless already set. 
+        */
+        internal virtual void AssertECMultiplier()
+        {
+            if (this.multiplier == null)
+            {
+                lock (this)
+                {
+                    if (this.multiplier == null)
+                    {
+                        this.multiplier = new FpNafMultiplier();
+                    }
+                }
+            }
+        }
+    }
+
+    public abstract class ECPointBase
+        : ECPoint
+    {
+        protected internal ECPointBase(
+            ECCurve			curve,
+            ECFieldElement	x,
+            ECFieldElement	y,
+            bool			withCompression)
+            : base(curve, x, y, withCompression)
+        {
+        }
+
+        protected internal abstract bool YTilde { get; }
+
+        /**
+         * return the field element encoded with point compression. (S 4.3.6)
+         */
+        public override byte[] GetEncoded(bool compressed)
+        {
+            if (this.IsInfinity)
+                return new byte[1];
+
+            // Note: some of the tests rely on calculating byte length from the field element
+            // (since the test cases use mismatching fields for curve/elements)
+            int byteLength = X9IntegerConverter.GetByteLength(x);
+            byte[] X = X9IntegerConverter.IntegerToBytes(this.X.ToBigInteger(), byteLength);
+            byte[] PO;
+
+            if (compressed)
+            {
+                PO = new byte[1 + X.Length];
+
+                PO[0] = (byte)(YTilde ? 0x03 : 0x02);
+            }
+            else
+            {
+                byte[] Y = X9IntegerConverter.IntegerToBytes(this.Y.ToBigInteger(), byteLength);
+                PO = new byte[1 + X.Length + Y.Length];
+
+                PO[0] = 0x04;
+
+                Y.CopyTo(PO, 1 + X.Length);
+            }
+
+            X.CopyTo(PO, 1);
+
+            return PO;
+        }
+
+        /**
+         * Multiplies this <code>ECPoint</code> by the given number.
+         * @param k The multiplicator.
+         * @return <code>k * this</code>.
+         */
+        public override ECPoint Multiply(
+            BigInteger k)
+        {
+            if (k.SignValue < 0)
+                throw new ArgumentException("The multiplicator cannot be negative", "k");
+
+            if (this.IsInfinity)
+                return this;
+
+            if (k.SignValue == 0)
+                return this.curve.Infinity;
+
+            AssertECMultiplier();
+            return this.multiplier.Multiply(this, k, preCompInfo);
+        }
+    }
+
+    /**
+     * Elliptic curve points over Fp
+     */
+    public class FpPoint
+        : ECPointBase
+    {
+        /**
+         * Create a point which encodes with point compression.
+         *
+         * @param curve the curve to use
+         * @param x affine x co-ordinate
+         * @param y affine y co-ordinate
+         */
+        public FpPoint(
+            ECCurve			curve,
+            ECFieldElement	x,
+            ECFieldElement	y)
+            : this(curve, x, y, false)
+        {
+        }
+
+        /**
+         * Create a point that encodes with or without point compresion.
+         *
+         * @param curve the curve to use
+         * @param x affine x co-ordinate
+         * @param y affine y co-ordinate
+         * @param withCompression if true encode with point compression
+         */
+        public FpPoint(
+            ECCurve			curve,
+            ECFieldElement	x,
+            ECFieldElement	y,
+            bool			withCompression)
+            : base(curve, x, y, withCompression)
+        {
+            if ((x == null) != (y == null))
+                throw new ArgumentException("Exactly one of the field elements is null");
+        }
+
+        protected internal override bool YTilde
+        {
+            get
+            {
+                return this.Y.ToBigInteger().TestBit(0);
+            }
+        }
+
+        // B.3 pg 62
+        public override ECPoint Add(
+            ECPoint b)
+        {
+            if (this.IsInfinity)
+                return b;
+
+            if (b.IsInfinity)
+                return this;
+
+            // Check if b = this or b = -this
+            if (this.x.Equals(b.x))
+            {
+                if (this.y.Equals(b.y))
+                {
+                    // this = b, i.e. this must be doubled
+                    return this.Twice();
+                }
+
+                Debug.Assert(this.y.Equals(b.y.Negate()));
+
+                // this = -b, i.e. the result is the point at infinity
+                return this.curve.Infinity;
+            }
+
+            ECFieldElement gamma = b.y.Subtract(this.y).Divide(b.x.Subtract(this.x));
+
+            ECFieldElement x3 = gamma.Square().Subtract(this.x).Subtract(b.x);
+            ECFieldElement y3 = gamma.Multiply(this.x.Subtract(x3)).Subtract(this.y);
+
+            return new FpPoint(curve, x3, y3, withCompression);
+        }
+
+        // B.3 pg 62
+        public override ECPoint Twice()
+        {
+            // Twice identity element (point at infinity) is identity
+            if (this.IsInfinity)
+                return this;
+
+            // if y1 == 0, then (x1, y1) == (x1, -y1)
+            // and hence this = -this and thus 2(x1, y1) == infinity
+            if (this.y.ToBigInteger().SignValue == 0)
+                return this.curve.Infinity;
+
+            ECFieldElement TWO = this.curve.FromBigInteger(BigInteger.Two);
+            ECFieldElement THREE = this.curve.FromBigInteger(BigInteger.Three);
+            ECFieldElement gamma = this.x.Square().Multiply(THREE).Add(curve.a).Divide(y.Multiply(TWO));
+
+            ECFieldElement x3 = gamma.Square().Subtract(this.x.Multiply(TWO));
+            ECFieldElement y3 = gamma.Multiply(this.x.Subtract(x3)).Subtract(this.y);
+
+            return new FpPoint(curve, x3, y3, this.withCompression);
+        }
+
+        // D.3.2 pg 102 (see Note:)
+        public override ECPoint Subtract(
+            ECPoint b)
+        {
+            if (b.IsInfinity)
+                return this;
+
+            // Add -b
+            return Add(b.Negate());
+        }
+
+        public override ECPoint Negate()
+        {
+            return new FpPoint(this.curve, this.x, this.y.Negate(), this.withCompression);
+        }
+
+        /**
+         * Sets the default <code>ECMultiplier</code>, unless already set. 
+         */
+        internal override void AssertECMultiplier()
+        {
+            if (this.multiplier == null)
+            {
+                lock (this)
+                {
+                    if (this.multiplier == null)
+                    {
+                        this.multiplier = new WNafMultiplier();
+                    }
+                }
+            }
+        }
+    }
+
+    /**
+     * Elliptic curve points over F2m
+     */
+    public class F2mPoint
+        : ECPointBase
+    {
+        /**
+         * @param curve base curve
+         * @param x x point
+         * @param y y point
+         */
+        public F2mPoint(
+            ECCurve			curve,
+            ECFieldElement	x,
+            ECFieldElement	y)
+            :  this(curve, x, y, false)
+        {
+        }
+
+        /**
+         * @param curve base curve
+         * @param x x point
+         * @param y y point
+         * @param withCompression true if encode with point compression.
+         */
+        public F2mPoint(
+            ECCurve			curve,
+            ECFieldElement	x,
+            ECFieldElement	y,
+            bool			withCompression)
+            : base(curve, x, y, withCompression)
+        {
+            if ((x != null && y == null) || (x == null && y != null))
+            {
+                throw new ArgumentException("Exactly one of the field elements is null");
+            }
+
+            if (x != null)
+            {
+                // Check if x and y are elements of the same field
+                F2mFieldElement.CheckFieldElements(this.x, this.y);
+
+                // Check if x and a are elements of the same field
+                F2mFieldElement.CheckFieldElements(this.x, this.curve.A);
+            }
+        }
+
+        /**
+         * Constructor for point at infinity
+         */
+        [Obsolete("Use ECCurve.Infinity property")]
+        public F2mPoint(
+            ECCurve curve)
+            : this(curve, null, null)
+        {
+        }
+
+        protected internal override bool YTilde
+        {
+            get
+            {
+                // X9.62 4.2.2 and 4.3.6:
+                // if x = 0 then ypTilde := 0, else ypTilde is the rightmost
+                // bit of y * x^(-1)
+                return this.X.ToBigInteger().SignValue != 0
+                    && this.Y.Multiply(this.X.Invert()).ToBigInteger().TestBit(0);
+            }
+        }
+
+        /**
+         * Check, if two <code>ECPoint</code>s can be added or subtracted.
+         * @param a The first <code>ECPoint</code> to check.
+         * @param b The second <code>ECPoint</code> to check.
+         * @throws IllegalArgumentException if <code>a</code> and <code>b</code>
+         * cannot be added.
+         */
+        private static void CheckPoints(
+            ECPoint	a,
+            ECPoint	b)
+        {
+            // Check, if points are on the same curve
+            if (!a.curve.Equals(b.curve))
+                throw new ArgumentException("Only points on the same curve can be added or subtracted");
+
+//			F2mFieldElement.CheckFieldElements(a.x, b.x);
+        }
+
+        /* (non-Javadoc)
+         * @see org.bouncycastle.math.ec.ECPoint#add(org.bouncycastle.math.ec.ECPoint)
+         */
+        public override ECPoint Add(ECPoint b)
+        {
+            CheckPoints(this, b);
+            return AddSimple((F2mPoint) b);
+        }
+
+        /**
+         * Adds another <code>ECPoints.F2m</code> to <code>this</code> without
+         * checking if both points are on the same curve. Used by multiplication
+         * algorithms, because there all points are a multiple of the same point
+         * and hence the checks can be omitted.
+         * @param b The other <code>ECPoints.F2m</code> to add to
+         * <code>this</code>.
+         * @return <code>this + b</code>
+         */
+        internal F2mPoint AddSimple(F2mPoint b)
+        {
+            if (this.IsInfinity)
+                return b;
+
+            if (b.IsInfinity)
+                return this;
+
+            F2mFieldElement x2 = (F2mFieldElement) b.X;
+            F2mFieldElement y2 = (F2mFieldElement) b.Y;
+
+            // Check if b == this or b == -this
+            if (this.x.Equals(x2))
+            {
+                // this == b, i.e. this must be doubled
+                if (this.y.Equals(y2))
+                    return (F2mPoint) this.Twice();
+
+                // this = -other, i.e. the result is the point at infinity
+                return (F2mPoint) this.curve.Infinity;
+            }
+
+            ECFieldElement xSum = this.x.Add(x2);
+
+            F2mFieldElement lambda
+                = (F2mFieldElement)(this.y.Add(y2)).Divide(xSum);
+
+            F2mFieldElement x3
+                = (F2mFieldElement)lambda.Square().Add(lambda).Add(xSum).Add(this.curve.A);
+
+            F2mFieldElement y3
+                = (F2mFieldElement)lambda.Multiply(this.x.Add(x3)).Add(x3).Add(this.y);
+
+            return new F2mPoint(curve, x3, y3, withCompression);
+        }
+
+        /* (non-Javadoc)
+         * @see org.bouncycastle.math.ec.ECPoint#subtract(org.bouncycastle.math.ec.ECPoint)
+         */
+        public override ECPoint Subtract(
+            ECPoint b)
+        {
+            CheckPoints(this, b);
+            return SubtractSimple((F2mPoint) b);
+        }
+
+        /**
+         * Subtracts another <code>ECPoints.F2m</code> from <code>this</code>
+         * without checking if both points are on the same curve. Used by
+         * multiplication algorithms, because there all points are a multiple
+         * of the same point and hence the checks can be omitted.
+         * @param b The other <code>ECPoints.F2m</code> to subtract from
+         * <code>this</code>.
+         * @return <code>this - b</code>
+         */
+        internal F2mPoint SubtractSimple(
+            F2mPoint b)
+        {
+            if (b.IsInfinity)
+                return this;
+
+            // Add -b
+            return AddSimple((F2mPoint) b.Negate());
+        }
+
+        /* (non-Javadoc)
+         * @see Org.BouncyCastle.Math.EC.ECPoint#twice()
+         */
+        public override ECPoint Twice()
+        {
+            // Twice identity element (point at infinity) is identity
+            if (this.IsInfinity)
+                return this;
+
+            // if x1 == 0, then (x1, y1) == (x1, x1 + y1)
+            // and hence this = -this and thus 2(x1, y1) == infinity
+            if (this.x.ToBigInteger().SignValue == 0)
+                return this.curve.Infinity;
+
+            F2mFieldElement lambda = (F2mFieldElement) this.x.Add(this.y.Divide(this.x));
+            F2mFieldElement x2 = (F2mFieldElement)lambda.Square().Add(lambda).Add(this.curve.A);
+            ECFieldElement ONE = this.curve.FromBigInteger(BigInteger.One);
+            F2mFieldElement y2 = (F2mFieldElement)this.x.Square().Add(
+                x2.Multiply(lambda.Add(ONE)));
+
+            return new F2mPoint(this.curve, x2, y2, withCompression);
+        }
+
+        public override ECPoint Negate()
+        {
+            return new F2mPoint(curve, this.x, this.x.Add(this.y), withCompression);
+        }
+
+        /**
+         * Sets the appropriate <code>ECMultiplier</code>, unless already set. 
+         */
+        internal override void AssertECMultiplier()
+        {
+            if (this.multiplier == null)
+            {
+                lock (this)
+                {
+                    if (this.multiplier == null)
+                    {
+                        if (((F2mCurve) this.curve).IsKoblitz)
+                        {
+                            this.multiplier = new WTauNafMultiplier();
+                        }
+                        else
+                        {
+                            this.multiplier = new WNafMultiplier();
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
diff --git a/crypto/src/math/ec/IntArray.cs b/crypto/src/math/ec/IntArray.cs
new file mode 100644
index 000000000..1089966a8
--- /dev/null
+++ b/crypto/src/math/ec/IntArray.cs
@@ -0,0 +1,485 @@
+using System;
+using System.Text;
+
+namespace Org.BouncyCastle.Math.EC
+{
+	internal class IntArray
+    {
+        // TODO make m fixed for the IntArray, and hence compute T once and for all
+
+		// TODO Use uint's internally?
+		private int[] m_ints;
+
+		public IntArray(int intLen)
+		{
+			m_ints = new int[intLen];
+		}
+
+		private IntArray(int[] ints)
+		{
+			m_ints = ints;
+		}
+
+		public IntArray(BigInteger bigInt)
+			: this(bigInt, 0)
+		{
+		}
+
+		public IntArray(BigInteger bigInt, int minIntLen)
+		{
+			if (bigInt.SignValue == -1)
+				throw new ArgumentException("Only positive Integers allowed", "bigint");
+
+			if (bigInt.SignValue == 0)
+			{
+				m_ints = new int[] { 0 };
+				return;
+			}
+
+			byte[] barr = bigInt.ToByteArrayUnsigned();
+			int barrLen = barr.Length;
+
+			int intLen = (barrLen + 3) / 4;
+			m_ints = new int[System.Math.Max(intLen, minIntLen)];
+
+			int rem = barrLen % 4;
+			int barrI = 0;
+
+			if (0 < rem)
+			{
+				int temp = (int) barr[barrI++];
+				while (barrI < rem)
+				{
+					temp = temp << 8 | (int) barr[barrI++];
+				}
+				m_ints[--intLen] = temp;
+			}
+
+			while (intLen > 0)
+			{
+				int temp = (int) barr[barrI++];
+				for (int i = 1; i < 4; i++)
+				{
+					temp = temp << 8 | (int) barr[barrI++];
+				}
+				m_ints[--intLen] = temp;
+			}
+		}
+
+		public int GetUsedLength()
+		{
+			int highestIntPos = m_ints.Length;
+
+			if (highestIntPos < 1)
+				return 0;
+
+			// Check if first element will act as sentinel
+			if (m_ints[0] != 0)
+			{
+				while (m_ints[--highestIntPos] == 0)
+				{
+				}
+				return highestIntPos + 1;
+			}
+
+			do
+			{
+				if (m_ints[--highestIntPos] != 0)
+				{
+					return highestIntPos + 1;
+				}
+			}
+			while (highestIntPos > 0);
+
+			return 0;
+		}
+
+		public int BitLength
+		{
+			get
+			{
+				// JDK 1.5: see Integer.numberOfLeadingZeros()
+				int intLen = GetUsedLength();
+				if (intLen == 0)
+					return 0;
+
+				int last = intLen - 1;
+				uint highest = (uint) m_ints[last];
+				int bits = (last << 5) + 1;
+
+				// A couple of binary search steps
+				if (highest > 0x0000ffff)
+				{
+					if (highest > 0x00ffffff)
+					{
+						bits += 24;
+						highest >>= 24;
+					}
+					else
+					{
+						bits += 16;
+						highest >>= 16;
+					}
+				}
+				else if (highest > 0x000000ff)
+				{
+					bits += 8;
+					highest >>= 8;
+				}
+
+				while (highest > 1)
+				{
+					++bits;
+					highest >>= 1;
+				}
+
+				return bits;
+			}
+		}
+
+		private int[] resizedInts(int newLen)
+		{
+			int[] newInts = new int[newLen];
+			int oldLen = m_ints.Length;
+			int copyLen = oldLen < newLen ? oldLen : newLen;
+			Array.Copy(m_ints, 0, newInts, 0, copyLen);
+			return newInts;
+		}
+
+		public BigInteger ToBigInteger()
+		{
+			int usedLen = GetUsedLength();
+			if (usedLen == 0)
+			{
+				return BigInteger.Zero;
+			}
+
+			int highestInt = m_ints[usedLen - 1];
+			byte[] temp = new byte[4];
+			int barrI = 0;
+			bool trailingZeroBytesDone = false;
+			for (int j = 3; j >= 0; j--)
+			{
+				byte thisByte = (byte)((int)((uint) highestInt >> (8 * j)));
+				if (trailingZeroBytesDone || (thisByte != 0))
+				{
+					trailingZeroBytesDone = true;
+					temp[barrI++] = thisByte;
+				}
+			}
+
+			int barrLen = 4 * (usedLen - 1) + barrI;
+			byte[] barr = new byte[barrLen];
+			for (int j = 0; j < barrI; j++)
+			{
+				barr[j] = temp[j];
+			}
+			// Highest value int is done now
+
+			for (int iarrJ = usedLen - 2; iarrJ >= 0; iarrJ--)
+			{
+				for (int j = 3; j >= 0; j--)
+				{
+					barr[barrI++] = (byte)((int)((uint)m_ints[iarrJ] >> (8 * j)));
+				}
+			}
+			return new BigInteger(1, barr);
+		}
+
+		public void ShiftLeft()
+		{
+			int usedLen = GetUsedLength();
+			if (usedLen == 0)
+			{
+				return;
+			}
+			if (m_ints[usedLen - 1] < 0)
+			{
+				// highest bit of highest used byte is set, so shifting left will
+				// make the IntArray one byte longer
+				usedLen++;
+				if (usedLen > m_ints.Length)
+				{
+					// make the m_ints one byte longer, because we need one more
+					// byte which is not available in m_ints
+					m_ints = resizedInts(m_ints.Length + 1);
+				}
+			}
+
+			bool carry = false;
+			for (int i = 0; i < usedLen; i++)
+			{
+				// nextCarry is true if highest bit is set
+				bool nextCarry = m_ints[i] < 0;
+				m_ints[i] <<= 1;
+				if (carry)
+				{
+					// set lowest bit
+					m_ints[i] |= 1;
+				}
+				carry = nextCarry;
+			}
+		}
+
+		public IntArray ShiftLeft(int n)
+		{
+			int usedLen = GetUsedLength();
+			if (usedLen == 0)
+			{
+				return this;
+			}
+
+			if (n == 0)
+			{
+				return this;
+			}
+
+			if (n > 31)
+			{
+				throw new ArgumentException("shiftLeft() for max 31 bits "
+					+ ", " + n + "bit shift is not possible", "n");
+			}
+
+			int[] newInts = new int[usedLen + 1];
+
+			int nm32 = 32 - n;
+			newInts[0] = m_ints[0] << n;
+			for (int i = 1; i < usedLen; i++)
+			{
+				newInts[i] = (m_ints[i] << n) | (int)((uint)m_ints[i - 1] >> nm32);
+			}
+			newInts[usedLen] = (int)((uint)m_ints[usedLen - 1] >> nm32);
+
+			return new IntArray(newInts);
+		}
+
+		public void AddShifted(IntArray other, int shift)
+		{
+			int usedLenOther = other.GetUsedLength();
+			int newMinUsedLen = usedLenOther + shift;
+			if (newMinUsedLen > m_ints.Length)
+			{
+				m_ints = resizedInts(newMinUsedLen);
+				//Console.WriteLine("Resize required");
+			}
+
+			for (int i = 0; i < usedLenOther; i++)
+			{
+				m_ints[i + shift] ^= other.m_ints[i];
+			}
+		}
+
+		public int Length
+		{
+			get { return m_ints.Length; }
+		}
+
+		public bool TestBit(int n)
+		{
+			// theInt = n / 32
+			int theInt = n >> 5;
+			// theBit = n % 32
+			int theBit = n & 0x1F;
+			int tester = 1 << theBit;
+			return ((m_ints[theInt] & tester) != 0);
+		}
+
+		public void FlipBit(int n)
+		{
+			// theInt = n / 32
+			int theInt = n >> 5;
+			// theBit = n % 32
+			int theBit = n & 0x1F;
+			int flipper = 1 << theBit;
+			m_ints[theInt] ^= flipper;
+		}
+
+		public void SetBit(int n)
+		{
+			// theInt = n / 32
+			int theInt = n >> 5;
+			// theBit = n % 32
+			int theBit = n & 0x1F;
+			int setter = 1 << theBit;
+			m_ints[theInt] |= setter;
+		}
+
+		public IntArray Multiply(IntArray other, int m)
+		{
+			// Lenght of c is 2m bits rounded up to the next int (32 bit)
+			int t = (m + 31) >> 5;
+			if (m_ints.Length < t)
+			{
+				m_ints = resizedInts(t);
+			}
+
+			IntArray b = new IntArray(other.resizedInts(other.Length + 1));
+			IntArray c = new IntArray((m + m + 31) >> 5);
+			// IntArray c = new IntArray(t + t);
+			int testBit = 1;
+			for (int k = 0; k < 32; k++)
+			{
+				for (int j = 0; j < t; j++)
+				{
+					if ((m_ints[j] & testBit) != 0)
+					{
+						// The kth bit of m_ints[j] is set
+						c.AddShifted(b, j);
+					}
+				}
+				testBit <<= 1;
+				b.ShiftLeft();
+			}
+			return c;
+		}
+
+		// public IntArray multiplyLeftToRight(IntArray other, int m) {
+		// // Lenght of c is 2m bits rounded up to the next int (32 bit)
+		// int t = (m + 31) / 32;
+		// if (m_ints.Length < t) {
+		// m_ints = resizedInts(t);
+		// }
+		//
+		// IntArray b = new IntArray(other.resizedInts(other.getLength() + 1));
+		// IntArray c = new IntArray((m + m + 31) / 32);
+		// // IntArray c = new IntArray(t + t);
+		// int testBit = 1 << 31;
+		// for (int k = 31; k >= 0; k--) {
+		// for (int j = 0; j < t; j++) {
+		// if ((m_ints[j] & testBit) != 0) {
+		// // The kth bit of m_ints[j] is set
+		// c.addShifted(b, j);
+		// }
+		// }
+		// testBit >>>= 1;
+		// if (k > 0) {
+		// c.shiftLeft();
+		// }
+		// }
+		// return c;
+		// }
+
+		// TODO note, redPol.Length must be 3 for TPB and 5 for PPB
+		public void Reduce(int m, int[] redPol)
+		{
+			for (int i = m + m - 2; i >= m; i--)
+			{
+				if (TestBit(i))
+				{
+					int bit = i - m;
+					FlipBit(bit);
+					FlipBit(i);
+					int l = redPol.Length;
+					while (--l >= 0)
+					{
+						FlipBit(redPol[l] + bit);
+					}
+				}
+			}
+			m_ints = resizedInts((m + 31) >> 5);
+		}
+
+		public IntArray Square(int m)
+		{
+			// TODO make the table static readonly
+			int[] table = { 0x0, 0x1, 0x4, 0x5, 0x10, 0x11, 0x14, 0x15, 0x40,
+									0x41, 0x44, 0x45, 0x50, 0x51, 0x54, 0x55 };
+
+			int t = (m + 31) >> 5;
+			if (m_ints.Length < t)
+			{
+				m_ints = resizedInts(t);
+			}
+
+			IntArray c = new IntArray(t + t);
+
+			// TODO twice the same code, put in separate private method
+			for (int i = 0; i < t; i++)
+			{
+				int v0 = 0;
+				for (int j = 0; j < 4; j++)
+				{
+					v0 = (int)((uint) v0 >> 8);
+					int u = (int)((uint)m_ints[i] >> (j * 4)) & 0xF;
+					int w = table[u] << 24;
+					v0 |= w;
+				}
+				c.m_ints[i + i] = v0;
+
+				v0 = 0;
+				int upper = (int)((uint) m_ints[i] >> 16);
+				for (int j = 0; j < 4; j++)
+				{
+					v0 = (int)((uint) v0 >> 8);
+					int u = (int)((uint)upper >> (j * 4)) & 0xF;
+					int w = table[u] << 24;
+					v0 |= w;
+				}
+				c.m_ints[i + i + 1] = v0;
+			}
+			return c;
+		}
+
+		public override bool Equals(object o)
+		{
+			if (!(o is IntArray))
+			{
+				return false;
+			}
+			IntArray other = (IntArray) o;
+			int usedLen = GetUsedLength();
+			if (other.GetUsedLength() != usedLen)
+			{
+				return false;
+			}
+			for (int i = 0; i < usedLen; i++)
+			{
+				if (m_ints[i] != other.m_ints[i])
+				{
+					return false;
+				}
+			}
+			return true;
+		}
+
+		public override int GetHashCode()
+		{
+			int i = GetUsedLength();
+			int hc = i;
+			while (--i >= 0)
+			{
+				hc *= 17;
+				hc ^= m_ints[i];
+			}
+			return hc;
+		}
+
+		internal IntArray Copy()
+		{
+			return new IntArray((int[]) m_ints.Clone());
+		}
+
+		public override string ToString()
+		{
+			int usedLen = GetUsedLength();
+			if (usedLen == 0)
+			{
+				return "0";
+			}
+
+			StringBuilder sb = new StringBuilder(Convert.ToString(m_ints[usedLen - 1], 2));
+			for (int iarrJ = usedLen - 2; iarrJ >= 0; iarrJ--)
+			{
+				string hexString = Convert.ToString(m_ints[iarrJ], 2);
+
+				// Add leading zeroes, except for highest significant int
+				for (int i = hexString.Length; i < 8; i++)
+				{
+					hexString = "0" + hexString;
+				}
+				sb.Append(hexString);
+			}
+			return sb.ToString();
+		}
+	}
+}
diff --git a/crypto/src/math/ec/abc/SimpleBigDecimal.cs b/crypto/src/math/ec/abc/SimpleBigDecimal.cs
new file mode 100644
index 000000000..d5664dbfd
--- /dev/null
+++ b/crypto/src/math/ec/abc/SimpleBigDecimal.cs
@@ -0,0 +1,241 @@
+using System;
+using System.Text;
+
+namespace Org.BouncyCastle.Math.EC.Abc
+{
+	/**
+	* Class representing a simple version of a big decimal. A
+	* <code>SimpleBigDecimal</code> is basically a
+	* {@link java.math.BigInteger BigInteger} with a few digits on the right of
+	* the decimal point. The number of (binary) digits on the right of the decimal
+	* point is called the <code>scale</code> of the <code>SimpleBigDecimal</code>.
+	* Unlike in {@link java.math.BigDecimal BigDecimal}, the scale is not adjusted
+	* automatically, but must be set manually. All <code>SimpleBigDecimal</code>s
+	* taking part in the same arithmetic operation must have equal scale. The
+	* result of a multiplication of two <code>SimpleBigDecimal</code>s returns a
+	* <code>SimpleBigDecimal</code> with double scale.
+	*/
+	internal class SimpleBigDecimal
+		//	: Number
+	{
+		//	private static final long serialVersionUID = 1L;
+
+		private readonly BigInteger	bigInt;
+		private readonly int		scale;
+
+		/**
+		* Returns a <code>SimpleBigDecimal</code> representing the same numerical
+		* value as <code>value</code>.
+		* @param value The value of the <code>SimpleBigDecimal</code> to be
+		* created. 
+		* @param scale The scale of the <code>SimpleBigDecimal</code> to be
+		* created. 
+		* @return The such created <code>SimpleBigDecimal</code>.
+		*/
+		public static SimpleBigDecimal GetInstance(BigInteger val, int scale)
+		{
+			return new SimpleBigDecimal(val.ShiftLeft(scale), scale);
+		}
+
+		/**
+		* Constructor for <code>SimpleBigDecimal</code>. The value of the
+		* constructed <code>SimpleBigDecimal</code> Equals <code>bigInt / 
+		* 2<sup>scale</sup></code>.
+		* @param bigInt The <code>bigInt</code> value parameter.
+		* @param scale The scale of the constructed <code>SimpleBigDecimal</code>.
+		*/
+		public SimpleBigDecimal(BigInteger bigInt, int scale)
+		{
+			if (scale < 0)
+				throw new ArgumentException("scale may not be negative");
+
+			this.bigInt = bigInt;
+			this.scale = scale;
+		}
+
+		private SimpleBigDecimal(SimpleBigDecimal limBigDec)
+		{
+			bigInt = limBigDec.bigInt;
+			scale = limBigDec.scale;
+		}
+
+		private void CheckScale(SimpleBigDecimal b)
+		{
+			if (scale != b.scale)
+				throw new ArgumentException("Only SimpleBigDecimal of same scale allowed in arithmetic operations");
+		}
+
+		public SimpleBigDecimal AdjustScale(int newScale)
+		{
+			if (newScale < 0)
+				throw new ArgumentException("scale may not be negative");
+
+			if (newScale == scale)
+				return this;
+
+			return new SimpleBigDecimal(bigInt.ShiftLeft(newScale - scale), newScale);
+		}
+
+		public SimpleBigDecimal Add(SimpleBigDecimal b)
+		{
+			CheckScale(b);
+			return new SimpleBigDecimal(bigInt.Add(b.bigInt), scale);
+		}
+
+		public SimpleBigDecimal Add(BigInteger b)
+		{
+			return new SimpleBigDecimal(bigInt.Add(b.ShiftLeft(scale)), scale);
+		}
+
+		public SimpleBigDecimal Negate()
+		{
+			return new SimpleBigDecimal(bigInt.Negate(), scale);
+		}
+
+		public SimpleBigDecimal Subtract(SimpleBigDecimal b)
+		{
+			return Add(b.Negate());
+		}
+
+		public SimpleBigDecimal Subtract(BigInteger b)
+		{
+			return new SimpleBigDecimal(bigInt.Subtract(b.ShiftLeft(scale)), scale);
+		}
+
+		public SimpleBigDecimal Multiply(SimpleBigDecimal b)
+		{
+			CheckScale(b);
+			return new SimpleBigDecimal(bigInt.Multiply(b.bigInt), scale + scale);
+		}
+
+		public SimpleBigDecimal Multiply(BigInteger b)
+		{
+			return new SimpleBigDecimal(bigInt.Multiply(b), scale);
+		}
+
+		public SimpleBigDecimal Divide(SimpleBigDecimal b)
+		{
+			CheckScale(b);
+			BigInteger dividend = bigInt.ShiftLeft(scale);
+			return new SimpleBigDecimal(dividend.Divide(b.bigInt), scale);
+		}
+
+		public SimpleBigDecimal Divide(BigInteger b)
+		{
+			return new SimpleBigDecimal(bigInt.Divide(b), scale);
+		}
+
+		public SimpleBigDecimal ShiftLeft(int n)
+		{
+			return new SimpleBigDecimal(bigInt.ShiftLeft(n), scale);
+		}
+
+		public int CompareTo(SimpleBigDecimal val)
+		{
+			CheckScale(val);
+			return bigInt.CompareTo(val.bigInt);
+		}
+
+		public int CompareTo(BigInteger val)
+		{
+			return bigInt.CompareTo(val.ShiftLeft(scale));
+		}
+
+		public BigInteger Floor()
+		{
+			return bigInt.ShiftRight(scale);
+		}
+
+		public BigInteger Round()
+		{
+			SimpleBigDecimal oneHalf = new SimpleBigDecimal(BigInteger.One, 1);
+			return Add(oneHalf.AdjustScale(scale)).Floor();
+		}
+
+		public int IntValue
+		{
+			get { return Floor().IntValue; }
+		}
+
+		public long LongValue
+		{
+			get { return Floor().LongValue; }
+		}
+
+//		public double doubleValue()
+//		{
+//			return new Double(ToString()).doubleValue();
+//		}
+//
+//		public float floatValue()
+//		{
+//			return new Float(ToString()).floatValue();
+//		}
+
+		public int Scale
+		{
+			get { return scale; }
+		}
+
+		public override string ToString()
+		{
+			if (scale == 0)
+				return bigInt.ToString();
+
+			BigInteger floorBigInt = Floor();
+	        
+			BigInteger fract = bigInt.Subtract(floorBigInt.ShiftLeft(scale));
+			if (bigInt.SignValue < 0)
+			{
+				fract = BigInteger.One.ShiftLeft(scale).Subtract(fract);
+			}
+
+			if ((floorBigInt.SignValue == -1) && (!(fract.Equals(BigInteger.Zero))))
+			{
+				floorBigInt = floorBigInt.Add(BigInteger.One);
+			}
+			string leftOfPoint = floorBigInt.ToString();
+
+			char[] fractCharArr = new char[scale];
+				string fractStr = fract.ToString(2);
+			int fractLen = fractStr.Length;
+			int zeroes = scale - fractLen;
+			for (int i = 0; i < zeroes; i++)
+			{
+				fractCharArr[i] = '0';
+			}
+			for (int j = 0; j < fractLen; j++)
+			{
+				fractCharArr[zeroes + j] = fractStr[j];
+			}
+			string rightOfPoint = new string(fractCharArr);
+
+			StringBuilder sb = new StringBuilder(leftOfPoint);
+			sb.Append(".");
+			sb.Append(rightOfPoint);
+
+			return sb.ToString();
+		}
+
+		public override bool Equals(
+			object obj)
+		{
+			if (this == obj)
+				return true;
+
+			SimpleBigDecimal other = obj as SimpleBigDecimal;
+
+			if (other == null)
+				return false;
+
+			return bigInt.Equals(other.bigInt)
+				&& scale == other.scale;
+		}
+
+		public override int GetHashCode()
+		{
+			return bigInt.GetHashCode() ^ scale;
+		}
+
+	}
+}
diff --git a/crypto/src/math/ec/abc/Tnaf.cs b/crypto/src/math/ec/abc/Tnaf.cs
new file mode 100644
index 000000000..225fc3075
--- /dev/null
+++ b/crypto/src/math/ec/abc/Tnaf.cs
@@ -0,0 +1,834 @@
+using System;
+
+namespace Org.BouncyCastle.Math.EC.Abc
+{
+	/**
+	* Class holding methods for point multiplication based on the window
+	* &#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;
+		}
+	}
+}