Implementation of homogeneous coordinates for Fp

Various changes to point methods to deal with non-affine points
Changes in client code and tests to apply point normalization

3 files changed, 1324 insertions, 966 deletions

diff --git a/crypto/src/math/ec/ECCurve.cs b/crypto/src/math/ec/ECCurve.cs
index 38daa719c..d369cb2b7 100644
--- a/crypto/src/math/ec/ECCurve.cs
+++ b/crypto/src/math/ec/ECCurve.cs
@@ -99,6 +99,8 @@ namespace Org.BouncyCastle.Math.EC
 
         protected abstract ECCurve CloneCurve();
 
+        protected internal abstract ECPoint CreateRawPoint(ECFieldElement x, ECFieldElement y, bool withCompression);
+
         protected virtual ECMultiplier CreateDefaultMultiplier()
         {
             return new WNafMultiplier();
@@ -145,7 +147,7 @@ namespace Org.BouncyCastle.Math.EC
             // TODO Default behaviour could be improved if the two curves have the same coordinate system by copying any Z coordinates.
             p = p.Normalize();
 
-            return CreatePoint(p.X.ToBigInteger(), p.Y.ToBigInteger(), p.IsCompressed);
+            return CreatePoint(p.XCoord.ToBigInteger(), p.YCoord.ToBigInteger(), p.IsCompressed);
         }
 
         /**
@@ -375,6 +377,20 @@ namespace Org.BouncyCastle.Math.EC
             return new FpCurve(m_q, m_r, m_a, m_b);
         }
 
+        public override bool SupportsCoordinateSystem(int coord)
+        {
+            switch (coord)
+            {
+                case COORD_AFFINE:
+                case COORD_HOMOGENEOUS:
+                //case COORD_JACOBIAN:
+                //case COORD_JACOBIAN_MODIFIED:
+                    return true;
+                default:
+                    return false;
+            }
+        }
+
         public virtual BigInteger Q
         {
             get { return m_q; }
@@ -395,6 +411,11 @@ namespace Org.BouncyCastle.Math.EC
             return new FpFieldElement(this.m_q, this.m_r, x);
         }
 
+        protected internal override ECPoint CreateRawPoint(ECFieldElement x, ECFieldElement y, bool withCompression)
+        {
+            return new FpPoint(this, x, y, withCompression);
+        }
+
         public override ECPoint CreatePoint(
             BigInteger	X1,
             BigInteger	Y1,
@@ -710,6 +731,11 @@ namespace Org.BouncyCastle.Math.EC
             return base.CreateDefaultMultiplier();
         }
 
+        protected internal override ECPoint CreateRawPoint(ECFieldElement x, ECFieldElement y, bool withCompression)
+        {
+            return new F2mPoint(this, x, y, withCompression);
+        }
+
         public override ECPoint Infinity
         {
             get { return m_infinity; }
diff --git a/crypto/src/math/ec/ECPoint.cs b/crypto/src/math/ec/ECPoint.cs
index 7a4450ac1..86134f80c 100644
--- a/crypto/src/math/ec/ECPoint.cs
+++ b/crypto/src/math/ec/ECPoint.cs
@@ -1,6 +1,7 @@
 using System;
 using System.Collections;
 using System.Diagnostics;
+using System.Text;
 
 using Org.BouncyCastle.Math.EC.Multiplier;
 
@@ -51,24 +52,18 @@ namespace Org.BouncyCastle.Math.EC
 
         protected internal PreCompInfo m_preCompInfo = null;
 
-        protected ECPoint(
-            ECCurve			curve,
-            ECFieldElement	x,
-            ECFieldElement	y,
-            bool			withCompression)
+        protected ECPoint(ECCurve curve, ECFieldElement	x, ECFieldElement y, bool withCompression)
+            : this(curve, x, y, GetInitialZCoords(curve), withCompression)
         {
-            this.m_curve = curve;
-            this.m_x = x;
-            this.m_y = y;
-            this.m_withCompression = withCompression;
         }
 
-        protected ECPoint(ECCurve curve, ECFieldElement x, ECFieldElement y, ECFieldElement[] zs)
+        internal ECPoint(ECCurve curve, ECFieldElement x, ECFieldElement y, ECFieldElement[] zs, bool withCompression)
         {
             this.m_curve = curve;
             this.m_x = x;
             this.m_y = y;
             this.m_zs = zs;
+            this.m_withCompression = withCompression;
         }
 
         public virtual ECCurve Curve
@@ -85,13 +80,85 @@ namespace Org.BouncyCastle.Math.EC
             }
         }
 
+        /**
+         * Normalizes this point, and then returns the affine x-coordinate.
+         * 
+         * Note: normalization can be expensive, this method is deprecated in favour
+         * of caller-controlled normalization.
+         */
+        [Obsolete("Use AffineXCoord, or Normalize() and XCoord, instead")]
         public virtual ECFieldElement X
         {
-            get { return m_x; }
+            get { return Normalize().XCoord; }
         }
 
+        /**
+         * Normalizes this point, and then returns the affine y-coordinate.
+         * 
+         * Note: normalization can be expensive, this method is deprecated in favour
+         * of caller-controlled normalization.
+         */
+        [Obsolete("Use AffineYCoord, or Normalize() and YCoord, instead")]
         public virtual ECFieldElement Y
         {
+            get { return Normalize().YCoord; }
+        }
+
+        /**
+         * Returns the affine x-coordinate after checking that this point is normalized.
+         * 
+         * @return The affine x-coordinate of this point
+         * @throws IllegalStateException if the point is not normalized
+         */
+        public virtual ECFieldElement AffineXCoord
+        {
+            get
+            {
+                CheckNormalized();
+                return XCoord;
+            }
+        }
+
+        /**
+         * Returns the affine y-coordinate after checking that this point is normalized
+         * 
+         * @return The affine y-coordinate of this point
+         * @throws IllegalStateException if the point is not normalized
+         */
+        public virtual ECFieldElement AffineYCoord
+        {
+            get
+            {
+                CheckNormalized();
+                return YCoord;
+            }
+        }
+
+        /**
+         * Returns the x-coordinate.
+         * 
+         * Caution: depending on the curve's coordinate system, this may not be the same value as in an
+         * affine coordinate system; use Normalize() to get a point where the coordinates have their
+         * affine values, or use AffineXCoord if you expect the point to already have been normalized.
+         * 
+         * @return the x-coordinate of this point
+         */
+        public virtual ECFieldElement XCoord
+        {
+            get { return m_x; }
+        }
+
+        /**
+         * Returns the y-coordinate.
+         * 
+         * Caution: depending on the curve's coordinate system, this may not be the same value as in an
+         * affine coordinate system; use Normalize() to get a point where the coordinates have their
+         * affine values, or use AffineYCoord if you expect the point to already have been normalized.
+         * 
+         * @return the y-coordinate of this point
+         */
+        public virtual ECFieldElement YCoord
+        {
             get { return m_y; }
         }
 
@@ -112,6 +179,16 @@ namespace Org.BouncyCastle.Math.EC
             return copy;
         }
 
+        protected virtual ECFieldElement RawXCoord
+        {
+            get { return m_x; }
+        }
+
+        protected virtual ECFieldElement RawYCoord
+        {
+            get { return m_y; }
+        }
+
         protected virtual void CheckNormalized()
         {
             if (!IsNormalized())
@@ -124,8 +201,8 @@ namespace Org.BouncyCastle.Math.EC
 
             return coord == ECCurve.COORD_AFFINE
                 || coord == ECCurve.COORD_LAMBDA_AFFINE
-                || IsInfinity;
-                //|| zs[0].isOne();
+                || IsInfinity
+                || GetZCoord(0).IsOne;
         }
 
         /**
@@ -163,27 +240,30 @@ namespace Org.BouncyCastle.Math.EC
 
         internal virtual ECPoint Normalize(ECFieldElement zInv)
         {
-            throw new InvalidOperationException("not a projective coordinate system");
-
-            //switch (this.CurveCoordinateSystem)
-            //{
-            //    case ECCurve.COORD_HOMOGENEOUS:
-            //    case ECCurve.COORD_LAMBDA_PROJECTIVE:
-            //    {
-            //        return CreateScaledPoint(zInv, zInv);
-            //    }
-            //    case ECCurve.COORD_JACOBIAN:
-            //    case ECCurve.COORD_JACOBIAN_CHUDNOVSKY:
-            //    case ECCurve.COORD_JACOBIAN_MODIFIED:
-            //    {
-            //        ECFieldElement zInv2 = zInv.Square(), zInv3 = zInv2.Multiply(zInv);
-            //        return CreateScaledPoint(zInv2, zInv3);
-            //    }
-            //    default:
-            //    {
-            //        throw new InvalidOperationException("not a projective coordinate system");
-            //    }
-            //}
+            switch (this.CurveCoordinateSystem)
+            {
+                case ECCurve.COORD_HOMOGENEOUS:
+                case ECCurve.COORD_LAMBDA_PROJECTIVE:
+                {
+                    return CreateScaledPoint(zInv, zInv);
+                }
+                case ECCurve.COORD_JACOBIAN:
+                case ECCurve.COORD_JACOBIAN_CHUDNOVSKY:
+                case ECCurve.COORD_JACOBIAN_MODIFIED:
+                {
+                    ECFieldElement zInv2 = zInv.Square(), zInv3 = zInv2.Multiply(zInv);
+                    return CreateScaledPoint(zInv2, zInv3);
+                }
+                default:
+                {
+                    throw new InvalidOperationException("not a projective coordinate system");
+                }
+            }
+        }
+
+        protected virtual ECPoint CreateScaledPoint(ECFieldElement sx, ECFieldElement sy)
+        {
+            return Curve.CreateRawPoint(RawXCoord.Multiply(sx), RawYCoord.Multiply(sy), IsCompressed);
         }
 
         public bool IsInfinity
@@ -208,26 +288,87 @@ namespace Org.BouncyCastle.Math.EC
             if (null == other)
                 return false;
 
+            ECCurve c1 = this.Curve, c2 = other.Curve;
+            bool n1 = (null == c1), n2 = (null == c2);
             bool i1 = IsInfinity, i2 = other.IsInfinity;
+
             if (i1 || i2)
             {
-                return i1 && i2;
+                return (i1 && i2) && (n1 || n2 || c1.Equals(c2));
+            }
+
+            ECPoint p1 = this, p2 = other;
+            if (n1 && n2)
+            {
+                // Points with null curve are in affine form, so already normalized
+            }
+            else if (n1)
+            {
+                p2 = p2.Normalize();
+            }
+            else if (n2)
+            {
+                p1 = p1.Normalize();
+            }
+            else if (!c1.Equals(c2))
+            {
+                return false;
+            }
+            else
+            {
+                // TODO Consider just requiring already normalized, to avoid silent performance degradation
+
+                ECPoint[] points = new ECPoint[] { this, c1.ImportPoint(p2) };
+
+                // TODO This is a little strong, really only requires coZNormalizeAll to get Zs equal
+                c1.NormalizeAll(points);
+
+                p1 = points[0];
+                p2 = points[1];
             }
 
-            return X.Equals(other.X) && Y.Equals(other.Y);
+            return p1.XCoord.Equals(p2.XCoord) && p1.YCoord.Equals(p2.YCoord);
         }
 
         public override int GetHashCode()
         {
-            int hc = 0;
-            if (!IsInfinity)
+            ECCurve c = this.Curve;
+            int hc = (null == c) ? 0 : ~c.GetHashCode();
+
+            if (!this.IsInfinity)
             {
-                hc ^= X.GetHashCode() * 17;
-                hc ^= Y.GetHashCode() * 257;
+                // TODO Consider just requiring already normalized, to avoid silent performance degradation
+
+                ECPoint p = Normalize();
+
+                hc ^= p.XCoord.GetHashCode() * 17;
+                hc ^= p.YCoord.GetHashCode() * 257;
             }
+
             return hc;
         }
 
+        public override string ToString()
+        {
+            if (this.IsInfinity)
+            {
+                return "INF";
+            }
+
+            StringBuilder sb = new StringBuilder();
+            sb.Append('(');
+            sb.Append(RawXCoord);
+            sb.Append(',');
+            sb.Append(RawYCoord);
+            for (int i = 0; i < m_zs.Length; ++i)
+            {
+                sb.Append(',');
+                sb.Append(m_zs[i]);
+            }
+            sb.Append(')');
+            return sb.ToString();
+        }
+
         public virtual byte[] GetEncoded()
         {
             return GetEncoded(m_withCompression);
@@ -280,6 +421,11 @@ namespace Org.BouncyCastle.Math.EC
         {
         }
 
+        protected internal ECPointBase(ECCurve curve, ECFieldElement x, ECFieldElement y, ECFieldElement[] zs, bool withCompression)
+            : base(curve, x, y, zs, withCompression)
+        {
+        }
+
         /**
          * return the field element encoded with point compression. (S 4.3.6)
          */
@@ -292,7 +438,7 @@ namespace Org.BouncyCastle.Math.EC
 
             ECPoint normed = Normalize();
 
-            byte[] X = normed.X.GetEncoded();
+            byte[] X = normed.XCoord.GetEncoded();
 
             if (compressed)
             {
@@ -302,7 +448,7 @@ namespace Org.BouncyCastle.Math.EC
                 return PO;
             }
 
-            byte[] Y = normed.Y.GetEncoded();
+            byte[] Y = normed.YCoord.GetEncoded();
 
             {
                 byte[] PO = new byte[X.Length + Y.Length + 1];
@@ -374,9 +520,14 @@ namespace Org.BouncyCastle.Math.EC
                 throw new ArgumentException("Exactly one of the field elements is null");
         }
 
+        internal FpPoint(ECCurve curve, ECFieldElement x, ECFieldElement y, ECFieldElement[] zs, bool withCompression)
+            : base(curve, x, y, zs, withCompression)
+        {
+        }
+
         protected internal override bool CompressionYTilde
         {
-            get { return this.Y.TestBitZero(); }
+            get { return this.AffineYCoord.TestBitZero(); }
         }
 
         // B.3 pg 62
@@ -396,28 +547,84 @@ namespace Org.BouncyCastle.Math.EC
                 return Twice();
             }
 
-            ECFieldElement X1 = this.X, Y1 = this.Y;
-            ECFieldElement X2 = b.X, Y2 = b.Y;
+            ECCurve curve = this.Curve;
+            int coord = curve.CoordinateSystem;
 
-            ECFieldElement dx = X2.Subtract(X1), dy = Y2.Subtract(Y1);
+            ECFieldElement X1 = this.XCoord, Y1 = this.YCoord;
+            ECFieldElement X2 = b.XCoord, Y2 = b.YCoord;
 
-            if (dx.IsZero)
+            switch (coord)
             {
-                if (dy.IsZero)
+                case ECCurve.COORD_AFFINE:
                 {
-                    // this == b, i.e. this must be doubled
-                    return Twice();
+                    ECFieldElement dx = X2.Subtract(X1), dy = Y2.Subtract(Y1);
+
+                    if (dx.IsZero)
+                    {
+                        if (dy.IsZero)
+                        {
+                            // this == b, i.e. this must be doubled
+                            return Twice();
+                        }
+
+                        // this == -b, i.e. the result is the point at infinity
+                        return Curve.Infinity;
+                    }
+
+                    ECFieldElement gamma = dy.Divide(dx);
+                    ECFieldElement X3 = gamma.Square().Subtract(X1).Subtract(X2);
+                    ECFieldElement Y3 = gamma.Multiply(X1.Subtract(X3)).Subtract(Y1);
+
+                    return new FpPoint(Curve, X3, Y3, IsCompressed);
                 }
 
-                // this == -b, i.e. the result is the point at infinity
-                return Curve.Infinity;
-            }
+                case ECCurve.COORD_HOMOGENEOUS:
+                {
+                    ECFieldElement Z1 = this.GetZCoord(0);
+                    ECFieldElement Z2 = b.GetZCoord(0);
+
+                    bool Z1IsOne = Z1.IsOne;
+                    bool Z2IsOne = Z2.IsOne;
+
+                    ECFieldElement u1 = Z1IsOne ? Y2 : Y2.Multiply(Z1);
+                    ECFieldElement u2 = Z2IsOne ? Y1 : Y1.Multiply(Z2);
+                    ECFieldElement u = u1.Subtract(u2);
+                    ECFieldElement v1 = Z1IsOne ? X2 : X2.Multiply(Z1);
+                    ECFieldElement v2 = Z2IsOne ? X1 : X1.Multiply(Z2);
+                    ECFieldElement v = v1.Subtract(v2);
+
+                    // Check if b == this or b == -this
+                    if (v.IsZero)
+                    {
+                        if (u.IsZero)
+                        {
+                            // this == b, i.e. this must be doubled
+                            return this.Twice();
+                        }
+
+                        // this == -b, i.e. the result is the point at infinity
+                        return curve.Infinity;
+                    }
+
+                    // TODO Optimize for when w == 1
+                    ECFieldElement w = Z1IsOne ? Z2 : Z2IsOne ? Z1 : Z1.Multiply(Z2);
+                    ECFieldElement vSquared = v.Square();
+                    ECFieldElement vCubed = vSquared.Multiply(v);
+                    ECFieldElement vSquaredV2 = vSquared.Multiply(v2);
+                    ECFieldElement A = u.Square().Multiply(w).Subtract(vCubed).Subtract(Two(vSquaredV2));
 
-            ECFieldElement gamma = dy.Divide(dx);
-            ECFieldElement X3 = gamma.Square().Subtract(X1).Subtract(X2);
-            ECFieldElement Y3 = gamma.Multiply(X1.Subtract(X3)).Subtract(Y1);
+                    ECFieldElement X3 = v.Multiply(A);
+                    ECFieldElement Y3 = vSquaredV2.Subtract(A).Multiply(u).Subtract(vCubed.Multiply(u2));
+                    ECFieldElement Z3 = vCubed.Multiply(w);
 
-            return new FpPoint(Curve, X3, Y3, IsCompressed);
+                    return new FpPoint(curve, X3, Y3, new ECFieldElement[] { Z3 }, IsCompressed);
+                }
+
+                default:
+                {
+                    throw new InvalidOperationException("unsupported coordinate system");
+                }
+            }
         }
 
         // B.3 pg 62
@@ -428,20 +635,65 @@ namespace Org.BouncyCastle.Math.EC
                 return this;
             }
 
-            ECFieldElement Y1 = this.Y;
+            ECCurve curve = this.Curve;
+
+            ECFieldElement Y1 = this.YCoord;
             if (Y1.IsZero) 
             {
-                return Curve.Infinity;
+                return curve.Infinity;
             }
 
-            ECFieldElement X1 = this.X;
+            int coord = curve.CoordinateSystem;
+
+            ECFieldElement X1 = this.XCoord;
+
+            switch (coord)
+            {
+                case ECCurve.COORD_AFFINE:
+                {
+                    ECFieldElement X1Squared = X1.Square();
+                    ECFieldElement gamma = Three(X1Squared).Add(this.Curve.A).Divide(Two(Y1));
+                    ECFieldElement X3 = gamma.Square().Subtract(Two(X1));
+                    ECFieldElement Y3 = gamma.Multiply(X1.Subtract(X3)).Subtract(Y1);
 
-            ECFieldElement X1Squared = X1.Square();
-            ECFieldElement gamma = Three(X1Squared).Add(this.Curve.A).Divide(Two(Y1));
-            ECFieldElement X3 = gamma.Square().Subtract(Two(X1));
-            ECFieldElement Y3 = gamma.Multiply(X1.Subtract(X3)).Subtract(Y1);
+                    return new FpPoint(Curve, X3, Y3, IsCompressed);
+                }
+
+                case ECCurve.COORD_HOMOGENEOUS:
+                {
+                    ECFieldElement Z1 = this.GetZCoord(0);
+
+                    bool Z1IsOne = Z1.IsOne;
+
+                    // TODO Optimize for small negative a4 and -3
+                    ECFieldElement w = curve.A;
+                    if (!w.IsZero && !Z1IsOne)
+                    {
+                        w = w.Multiply(Z1.Square());
+                    }
+                    w = w.Add(Three(X1.Square()));
+
+                    ECFieldElement s = Z1IsOne ? Y1 : Y1.Multiply(Z1);
+                    ECFieldElement t = Z1IsOne ? Y1.Square() : s.Multiply(Y1);
+                    ECFieldElement B = X1.Multiply(t);
+                    ECFieldElement _4B = Four(B);
+                    ECFieldElement h = w.Square().Subtract(Two(_4B));
+
+                    ECFieldElement _2s = Two(s);
+                    ECFieldElement X3 = h.Multiply(_2s);
+                    ECFieldElement _2t = Two(t);
+                    ECFieldElement Y3 = _4B.Subtract(h).Multiply(w).Subtract(Two(_2t.Square()));
+                    ECFieldElement _4sSquared = Z1IsOne ? Two(_2t) : _2s.Square();
+                    ECFieldElement Z3 = Two(_4sSquared).Multiply(s);
+
+                    return new FpPoint(curve, X3, Y3, new ECFieldElement[] { Z3 }, IsCompressed);
+                }
 
-            return new FpPoint(Curve, X3, Y3, IsCompressed);
+                default:
+                {
+                    throw new InvalidOperationException("unsupported coordinate system");
+                }
+            }
         }
 
         public override ECPoint TwicePlus(ECPoint b)
@@ -459,79 +711,106 @@ namespace Org.BouncyCastle.Math.EC
                 return Twice();
             }
 
-            ECFieldElement Y1 = this.Y;
+            ECFieldElement Y1 = this.YCoord;
             if (Y1.IsZero)
             {
                 return b;
             }
 
-            ECFieldElement X1 = this.X;
-            ECFieldElement X2 = b.X, Y2 = b.Y;
-
-            ECFieldElement dx = X2.Subtract(X1), dy = Y2.Subtract(Y1);
+            ECCurve curve = this.Curve;
+            int coord = curve.CoordinateSystem;
 
-            if (dx.IsZero)
+            switch (coord)
             {
-                if (dy.IsZero)
+                case ECCurve.COORD_AFFINE:
                 {
-                    // this == b i.e. the result is 3P
-                    return ThreeTimes();
-                }
+                    ECFieldElement X1 = this.XCoord;
+                    ECFieldElement X2 = b.XCoord, Y2 = b.YCoord;
 
-                // this == -b, i.e. the result is P
-                return this;
-            }
+                    ECFieldElement dx = X2.Subtract(X1), dy = Y2.Subtract(Y1);
 
-            /*
-             * Optimized calculation of 2P + Q, as described in "Trading Inversions for
-             * Multiplications in Elliptic Curve Cryptography", by Ciet, Joye, Lauter, Montgomery.
-             */
+                    if (dx.IsZero)
+                    {
+                        if (dy.IsZero)
+                        {
+                            // this == b i.e. the result is 3P
+                            return ThreeTimes();
+                        }
 
-            ECFieldElement X = dx.Square(), Y = dy.Square();
-            ECFieldElement d = X.Multiply(Two(X1).Add(X2)).Subtract(Y);
-            if (d.IsZero)
-            {
-                return Curve.Infinity;
-            }
+                        // this == -b, i.e. the result is P
+                        return this;
+                    }
 
-            ECFieldElement D = d.Multiply(dx);
-            ECFieldElement I = D.Invert();
-            ECFieldElement L1 = d.Multiply(I).Multiply(dy);
-            ECFieldElement L2 = Two(Y1).Multiply(X).Multiply(dx).Multiply(I).Subtract(L1);
-            ECFieldElement X4 = (L2.Subtract(L1)).Multiply(L1.Add(L2)).Add(X2);
-            ECFieldElement Y4 = (X1.Subtract(X4)).Multiply(L2).Subtract(Y1);
+                    /*
+                     * Optimized calculation of 2P + Q, as described in "Trading Inversions for
+                     * Multiplications in Elliptic Curve Cryptography", by Ciet, Joye, Lauter, Montgomery.
+                     */
 
-            return new FpPoint(Curve, X4, Y4, IsCompressed);
+                    ECFieldElement X = dx.Square(), Y = dy.Square();
+                    ECFieldElement d = X.Multiply(Two(X1).Add(X2)).Subtract(Y);
+                    if (d.IsZero)
+                    {
+                        return Curve.Infinity;
+                    }
+
+                    ECFieldElement D = d.Multiply(dx);
+                    ECFieldElement I = D.Invert();
+                    ECFieldElement L1 = d.Multiply(I).Multiply(dy);
+                    ECFieldElement L2 = Two(Y1).Multiply(X).Multiply(dx).Multiply(I).Subtract(L1);
+                    ECFieldElement X4 = (L2.Subtract(L1)).Multiply(L1.Add(L2)).Add(X2);
+                    ECFieldElement Y4 = (X1.Subtract(X4)).Multiply(L2).Subtract(Y1);
+
+                    return new FpPoint(Curve, X4, Y4, IsCompressed);
+                }
+                default:
+                {
+                    return Twice().Add(b);
+                }
+            }
         }
 
         public override ECPoint ThreeTimes()
         {
-            if (IsInfinity || this.Y.IsZero)
+            if (IsInfinity || this.YCoord.IsZero)
             {
                 return this;
             }
 
-            ECFieldElement X1 = this.X, Y1 = this.Y;
-
-            ECFieldElement _2Y1 = Two(Y1);
-            ECFieldElement X = _2Y1.Square();
-            ECFieldElement Z = Three(X1.Square()).Add(Curve.A);
-            ECFieldElement Y = Z.Square();
+            ECCurve curve = this.Curve;
+            int coord = curve.CoordinateSystem;
 
-            ECFieldElement d = Three(X1).Multiply(X).Subtract(Y);
-            if (d.IsZero)
+            switch (coord)
             {
-                return Curve.Infinity;
-            }
+                case ECCurve.COORD_AFFINE:
+                {
+                    ECFieldElement X1 = this.XCoord, Y1 = this.YCoord;
 
-            ECFieldElement D = d.Multiply(_2Y1);
-            ECFieldElement I = D.Invert();
-            ECFieldElement L1 = d.Multiply(I).Multiply(Z);
-            ECFieldElement L2 = X.Square().Multiply(I).Subtract(L1);
+                    ECFieldElement _2Y1 = Two(Y1);
+                    ECFieldElement X = _2Y1.Square();
+                    ECFieldElement Z = Three(X1.Square()).Add(Curve.A);
+                    ECFieldElement Y = Z.Square();
 
-            ECFieldElement X4 = (L2.Subtract(L1)).Multiply(L1.Add(L2)).Add(X1);
-            ECFieldElement Y4 = (X1.Subtract(X4)).Multiply(L2).Subtract(Y1);
-            return new FpPoint(Curve, X4, Y4, IsCompressed);
+                    ECFieldElement d = Three(X1).Multiply(X).Subtract(Y);
+                    if (d.IsZero)
+                    {
+                        return Curve.Infinity;
+                    }
+
+                    ECFieldElement D = d.Multiply(_2Y1);
+                    ECFieldElement I = D.Invert();
+                    ECFieldElement L1 = d.Multiply(I).Multiply(Z);
+                    ECFieldElement L2 = X.Square().Multiply(I).Subtract(L1);
+
+                    ECFieldElement X4 = (L2.Subtract(L1)).Multiply(L1.Add(L2)).Add(X1);
+                    ECFieldElement Y4 = (X1.Subtract(X4)).Multiply(L2).Subtract(Y1);
+                    return new FpPoint(Curve, X4, Y4, IsCompressed);
+                }
+                default:
+                {
+                    // NOTE: Be careful about recursions between twicePlus and threeTimes
+                    return Twice().Add(this);
+                }
+            }
         }
 
         protected virtual ECFieldElement Two(ECFieldElement x)
@@ -583,14 +862,14 @@ namespace Org.BouncyCastle.Math.EC
             }
 
             ECCurve curve = this.Curve;
-            //int coord = curve.CoordinateSystem;
+            int coord = curve.CoordinateSystem;
 
-            //if (ECCurve.COORD_AFFINE != coord)
-            //{
-            //    return new FpPoint(curve, X, Y.Negate(), this.m_zs, IsCompressed);
-            //}
+            if (ECCurve.COORD_AFFINE != coord)
+            {
+                return new FpPoint(curve, XCoord, YCoord.Negate(), this.m_zs, IsCompressed);
+            }
 
-            return new FpPoint(curve, X, Y.Negate(), IsCompressed);
+            return new FpPoint(curve, XCoord, YCoord.Negate(), IsCompressed);
         }
     }
 
@@ -637,10 +916,18 @@ namespace Org.BouncyCastle.Math.EC
                 F2mFieldElement.CheckFieldElements(x, y);
 
                 // Check if x and a are elements of the same field
-                F2mFieldElement.CheckFieldElements(x, curve.A);
+                if (curve != null)
+                {
+                    F2mFieldElement.CheckFieldElements(x, curve.A);
+                }
             }
         }
 
+        internal F2mPoint(ECCurve curve, ECFieldElement x, ECFieldElement y, ECFieldElement[] zs, bool withCompression)
+            : base(curve, x, y, zs, withCompression)
+        {
+        }
+
         /**
          * Constructor for point at infinity
          */
@@ -655,10 +942,28 @@ namespace Org.BouncyCastle.Math.EC
         {
             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.IsZero && this.Y.Divide(this.X).TestBitZero();
+                ECFieldElement X = this.RawXCoord;
+                if (X.IsZero)
+                {
+                    return false;
+                }
+
+                ECFieldElement Y = this.RawYCoord;
+
+                switch (this.CurveCoordinateSystem)
+                {
+                    case ECCurve.COORD_LAMBDA_AFFINE:
+                    case ECCurve.COORD_LAMBDA_PROJECTIVE:
+                    {
+                        // Y is actually Lambda (X + Y/X) here
+                        return Y.Subtract(X).TestBitZero();
+                    }
+                    default:
+                    {
+                        return Y.Divide(X).TestBitZero();
+                    }
+                }
+
             }
         }
 
@@ -706,30 +1011,30 @@ namespace Org.BouncyCastle.Math.EC
             if (b.IsInfinity)
                 return this;
 
-            F2mFieldElement x2 = (F2mFieldElement) b.X;
-            F2mFieldElement y2 = (F2mFieldElement) b.Y;
+            F2mFieldElement x2 = (F2mFieldElement) b.XCoord;
+            F2mFieldElement y2 = (F2mFieldElement) b.YCoord;
 
             // Check if b == this or b == -this
-            if (this.X.Equals(x2))
+            if (this.XCoord.Equals(x2))
             {
                 // this == b, i.e. this must be doubled
-                if (this.Y.Equals(y2))
+                if (this.YCoord.Equals(y2))
                     return (F2mPoint) this.Twice();
 
                 // this = -other, i.e. the result is the point at infinity
                 return (F2mPoint) Curve.Infinity;
             }
 
-            ECFieldElement xSum = this.X.Add(x2);
+            ECFieldElement xSum = this.XCoord.Add(x2);
 
             F2mFieldElement lambda
-                = (F2mFieldElement)(this.Y.Add(y2)).Divide(xSum);
+                = (F2mFieldElement)(this.YCoord.Add(y2)).Divide(xSum);
 
             F2mFieldElement x3
                 = (F2mFieldElement)lambda.Square().Add(lambda).Add(xSum).Add(Curve.A);
 
             F2mFieldElement y3
-                = (F2mFieldElement)lambda.Multiply(this.X.Add(x3)).Add(x3).Add(this.Y);
+                = (F2mFieldElement)lambda.Multiply(this.XCoord.Add(x3)).Add(x3).Add(this.YCoord);
 
             return new F2mPoint(Curve, x3, y3, IsCompressed);
         }
@@ -763,6 +1068,39 @@ namespace Org.BouncyCastle.Math.EC
             return AddSimple((F2mPoint) b.Negate());
         }
 
+        public virtual F2mPoint Tau()
+        {
+            if (this.IsInfinity)
+            {
+                return this;
+            }
+
+            ECCurve curve = this.Curve;
+            int coord = curve.CoordinateSystem;
+
+            ECFieldElement X1 = this.XCoord;
+
+            switch (coord)
+            {
+                case ECCurve.COORD_AFFINE:
+                case ECCurve.COORD_LAMBDA_AFFINE:
+                {
+                    ECFieldElement Y1 = this.YCoord;
+                    return new F2mPoint(curve, X1.Square(), Y1.Square(), IsCompressed);
+                }
+                case ECCurve.COORD_HOMOGENEOUS:
+                case ECCurve.COORD_LAMBDA_PROJECTIVE:
+                {
+                    ECFieldElement Y1 = this.YCoord, Z1 = this.GetZCoord(0);
+                    return new F2mPoint(curve, X1.Square(), Y1.Square(), new ECFieldElement[] { Z1.Square() }, IsCompressed);
+                }
+                default:
+                {
+                    throw new InvalidOperationException("unsupported coordinate system");
+                }
+            }
+        }
+
         /* (non-Javadoc)
          * @see Org.BouncyCastle.Math.EC.ECPoint#twice()
          */
@@ -774,15 +1112,15 @@ namespace Org.BouncyCastle.Math.EC
 
             // if x1 == 0, then (x1, y1) == (x1, x1 + y1)
             // and hence this = -this and thus 2(x1, y1) == infinity
-            if (this.X.IsZero)
+            if (this.XCoord.IsZero)
             {
                 return Curve.Infinity;
             }
 
-            F2mFieldElement lambda = (F2mFieldElement) this.X.Add(this.Y.Divide(this.X));
+            F2mFieldElement lambda = (F2mFieldElement) this.XCoord.Add(this.YCoord.Divide(this.XCoord));
             F2mFieldElement x2 = (F2mFieldElement)lambda.Square().Add(lambda).Add(Curve.A);
             ECFieldElement ONE = Curve.FromBigInteger(BigInteger.One);
-            F2mFieldElement y2 = (F2mFieldElement)this.X.Square().Add(
+            F2mFieldElement y2 = (F2mFieldElement)this.XCoord.Square().Add(
                 x2.Multiply(lambda.Add(ONE)));
 
             return new F2mPoint(Curve, x2, y2, IsCompressed);
@@ -795,13 +1133,13 @@ namespace Org.BouncyCastle.Math.EC
                 return this;
             }
 
-            ECFieldElement X1 = this.X;
+            ECFieldElement X1 = this.XCoord;
             if (X1.IsZero)
             {
                 return this;
             }
 
-            return new F2mPoint(Curve, X1, X1.Add(this.Y), IsCompressed);
+            return new F2mPoint(Curve, X1, X1.Add(this.YCoord), IsCompressed);
         }
     }
 }
diff --git a/crypto/src/math/ec/abc/Tnaf.cs b/crypto/src/math/ec/abc/Tnaf.cs
index 225fc3075..0ba414e68 100644
--- a/crypto/src/math/ec/abc/Tnaf.cs
+++ b/crypto/src/math/ec/abc/Tnaf.cs
@@ -2,833 +2,827 @@ 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;
-		}
-	}
+    /**
+    * 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)
+        {
+            return p.Tau();
+        }
+
+        /**
+        * 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;
+        }
+    }
 }