Port of WNaf precomp optimization from Java

3 files changed, 140 insertions, 36 deletions

diff --git a/crypto/src/math/ec/ECAlgorithms.cs b/crypto/src/math/ec/ECAlgorithms.cs
index 3c911b173..a1349a9e0 100644
--- a/crypto/src/math/ec/ECAlgorithms.cs
+++ b/crypto/src/math/ec/ECAlgorithms.cs
@@ -117,6 +117,11 @@ namespace Org.BouncyCastle.Math.EC
 
         public static void MontgomeryTrick(ECFieldElement[] zs, int off, int len)
         {
+            MontgomeryTrick(zs, off, len, null);
+        }
+
+        public static void MontgomeryTrick(ECFieldElement[] zs, int off, int len, ECFieldElement scale)
+        {
             /*
              * Uses the "Montgomery Trick" to invert many field elements, with only a single actual
              * field inversion. See e.g. the paper:
@@ -133,7 +138,14 @@ namespace Org.BouncyCastle.Math.EC
                 c[i] = c[i - 1].Multiply(zs[off + i]);
             }
 
-            ECFieldElement u = c[--i].Invert();
+            --i;
+
+            if (scale != null)
+            {
+                c[i] = c[i].Multiply(scale);
+            }
+
+            ECFieldElement u = c[i].Invert();
 
             while (i > 0)
             {
diff --git a/crypto/src/math/ec/ECCurve.cs b/crypto/src/math/ec/ECCurve.cs
index eaa3e0c3d..339d37f7c 100644
--- a/crypto/src/math/ec/ECCurve.cs
+++ b/crypto/src/math/ec/ECCurve.cs
@@ -221,26 +221,56 @@ namespace Org.BouncyCastle.Math.EC
          */
         public virtual void NormalizeAll(ECPoint[] points)
         {
-            CheckPoints(points);
+            NormalizeAll(points, 0, points.Length, null);
+        }
+
+        /**
+         * Normalization ensures that any projective coordinate is 1, and therefore that the x, y
+         * coordinates reflect those of the equivalent point in an affine coordinate system. Where more
+         * than one point is to be normalized, this method will generally be more efficient than
+         * normalizing each point separately. An (optional) z-scaling factor can be applied; effectively
+         * each z coordinate is scaled by this value prior to normalization (but only one
+         * actual multiplication is needed).
+         * 
+         * @param points
+         *            An array of points that will be updated in place with their normalized versions,
+         *            where necessary
+         * @param off
+         *            The start of the range of points to normalize
+         * @param len
+         *            The length of the range of points to normalize
+         * @param iso
+         *            The (optional) z-scaling factor - can be null
+         */
+        public virtual void NormalizeAll(ECPoint[] points, int off, int len, ECFieldElement iso)
+        {
+            CheckPoints(points, off, len);
 
-            if (this.CoordinateSystem == ECCurve.COORD_AFFINE)
+            switch (this.CoordinateSystem)
             {
-                return;
+                case ECCurve.COORD_AFFINE:
+                case ECCurve.COORD_LAMBDA_AFFINE:
+                {
+                    if (iso != null)
+                        throw new ArgumentException("not valid for affine coordinates", "iso");
+
+                    return;
+                }
             }
 
             /*
              * Figure out which of the points actually need to be normalized
              */
-            ECFieldElement[] zs = new ECFieldElement[points.Length];
-            int[] indices = new int[points.Length];
+            ECFieldElement[] zs = new ECFieldElement[len];
+            int[] indices = new int[len];
             int count = 0;
-            for (int i = 0; i < points.Length; ++i)
+            for (int i = 0; i < len; ++i)
             {
-                ECPoint p = points[i];
-                if (null != p && !p.IsNormalized())
+                ECPoint p = points[off + i];
+                if (null != p && (iso != null || !p.IsNormalized()))
                 {
                     zs[count] = p.GetZCoord(0);
-                    indices[count++] = i;
+                    indices[count++] = off + i;
                 }
             }
 
@@ -249,7 +279,7 @@ namespace Org.BouncyCastle.Math.EC
                 return;
             }
 
-            ECAlgorithms.MontgomeryTrick(zs, 0, count);
+            ECAlgorithms.MontgomeryTrick(zs, 0, count, iso);
 
             for (int j = 0; j < count; ++j)
             {
@@ -298,12 +328,19 @@ namespace Org.BouncyCastle.Math.EC
 
         protected virtual void CheckPoints(ECPoint[] points)
         {
+            CheckPoints(points, 0, points.Length);
+        }
+
+        protected virtual void CheckPoints(ECPoint[] points, int off, int len)
+        {
             if (points == null)
                 throw new ArgumentNullException("points");
+            if (off < 0 || len < 0 || (off > (points.Length - len)))
+                throw new ArgumentException("invalid range specified", "points");
 
-            for (int i = 0; i < points.Length; ++i)
+            for (int i = 0; i < len; ++i)
             {
-                ECPoint point = points[i];
+                ECPoint point = points[off + i];
                 if (null != point && this != point.Curve)
                     throw new ArgumentException("entries must be null or on this curve", "points");
             }
diff --git a/crypto/src/math/ec/multiplier/WNafUtilities.cs b/crypto/src/math/ec/multiplier/WNafUtilities.cs
index 865b9073e..5491297d7 100644
--- a/crypto/src/math/ec/multiplier/WNafUtilities.cs
+++ b/crypto/src/math/ec/multiplier/WNafUtilities.cs
@@ -10,6 +10,7 @@ namespace Org.BouncyCastle.Math.EC.Multiplier
 
         private static readonly byte[] EMPTY_BYTES = new byte[0];
         private static readonly int[] EMPTY_INTS = new int[0];
+        private static readonly ECPoint[] EMPTY_POINTS = new ECPoint[0];
 
         public static int[] GenerateCompactNaf(BigInteger k)
         {
@@ -368,46 +369,100 @@ namespace Org.BouncyCastle.Math.EC.Multiplier
         {
             ECCurve c = p.Curve;
             WNafPreCompInfo wnafPreCompInfo = GetWNafPreCompInfo(c.GetPreCompInfo(p, PRECOMP_NAME));
-            
+
+            int iniPreCompLen = 0, reqPreCompLen = 1 << System.Math.Max(0, width - 2);
+
             ECPoint[] preComp = wnafPreCompInfo.PreComp;
             if (preComp == null)
             {
-                preComp = new ECPoint[]{ p };
+                preComp = EMPTY_POINTS;
+            }
+            else
+            {
+                iniPreCompLen = preComp.Length;
             }
 
-            int preCompLen = preComp.Length;
-            int reqPreCompLen = 1 << System.Math.Max(0, width - 2);
-
-            if (preCompLen < reqPreCompLen)
+            if (iniPreCompLen < reqPreCompLen)
             {
                 preComp = ResizeTable(preComp, reqPreCompLen);
-                if (reqPreCompLen == 2)
+
+                if (reqPreCompLen == 1)
                 {
-                    preComp[1] = preComp[0].ThreeTimes();
+                    preComp[0] = p.Normalize();
                 }
                 else
                 {
-                    ECPoint twiceP = wnafPreCompInfo.Twice;
-                    if (twiceP == null)
+                    int curPreCompLen = iniPreCompLen;
+                    if (curPreCompLen == 0)
                     {
-                        twiceP = preComp[0].Twice();
-                        wnafPreCompInfo.Twice = twiceP;
+                        preComp[0] = p;
+                        curPreCompLen = 1;
                     }
 
-                    for (int i = preCompLen; i < reqPreCompLen; i++)
+                    ECFieldElement iso = null;
+
+                    if (reqPreCompLen == 2)
                     {
-                        /*
-                         * 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]);
+                        preComp[1] = p.ThreeTimes();
+                    }
+                    else
+                    {
+                        ECPoint twiceP = wnafPreCompInfo.Twice, last = preComp[curPreCompLen - 1];
+                        if (twiceP == null)
+                        {
+                            twiceP = preComp[0].Twice();
+                            wnafPreCompInfo.Twice = twiceP;
+
+                            /*
+                             * For Fp curves with Jacobian projective coordinates, use a (quasi-)isomorphism
+                             * where 'twiceP' is "affine", so that the subsequent additions are cheaper. This
+                             * also requires scaling the initial point's X, Y coordinates, and reversing the
+                             * isomorphism as part of the subsequent normalization.
+                             * 
+                             *  NOTE: The correctness of this optimization depends on:
+                             *      1) additions do not use the curve's A, B coefficients.
+                             *      2) no special cases (i.e. Q +/- Q) when calculating 1P, 3P, 5P, ...
+                             */
+                            if (ECAlgorithms.IsFpCurve(c) && c.FieldSize >= 64)
+                            {
+                                switch (c.CoordinateSystem)
+                                {
+                                    case ECCurve.COORD_JACOBIAN:
+                                    case ECCurve.COORD_JACOBIAN_CHUDNOVSKY:
+                                    case ECCurve.COORD_JACOBIAN_MODIFIED:
+                                    {
+                                        iso = twiceP.GetZCoord(0);
+                                        twiceP = c.CreatePoint(twiceP.XCoord.ToBigInteger(),
+                                            twiceP.YCoord.ToBigInteger());
+
+                                        ECFieldElement iso2 = iso.Square(), iso3 = iso2.Multiply(iso);
+                                        last = last.ScaleX(iso2).ScaleY(iso3);
+
+                                        if (iniPreCompLen == 0)
+                                        {
+                                            preComp[0] = last;
+                                        }
+                                        break;
+                                    }
+                                }
+                            }
+                        }
+
+                        while (curPreCompLen < reqPreCompLen)
+                        {
+                            /*
+                             * Compute the new ECPoints for the precomputation array. The values 1, 3,
+                             * 5, ..., 2^(width-1)-1 times p are computed
+                             */
+                            preComp[curPreCompLen++] = last = last.Add(twiceP);
+                        }
                     }
-                }
 
-                /*
-                 * Having oft-used operands in affine form makes operations faster.
-                 */
-                c.NormalizeAll(preComp);
+                    /*
+                     * Having oft-used operands in affine form makes operations faster.
+                     */
+                    c.NormalizeAll(preComp, iniPreCompLen, reqPreCompLen - iniPreCompLen, iso);
+                }
             }
 
             wnafPreCompInfo.PreComp = preComp;