Add support for delayed modular reduction

3 files changed, 314 insertions, 35 deletions

diff --git a/crypto/src/math/ec/ECFieldElement.cs b/crypto/src/math/ec/ECFieldElement.cs
index 66fadddec..13ae32e5d 100644
--- a/crypto/src/math/ec/ECFieldElement.cs
+++ b/crypto/src/math/ec/ECFieldElement.cs
@@ -35,6 +35,26 @@ namespace Org.BouncyCastle.Math.EC
             get { return 0 == ToBigInteger().SignValue; }
         }
 
+        public virtual ECFieldElement MultiplyMinusProduct(ECFieldElement b, ECFieldElement x, ECFieldElement y)
+        {
+            return Multiply(b).Subtract(x.Multiply(y));
+        }
+
+        public virtual ECFieldElement MultiplyPlusProduct(ECFieldElement b, ECFieldElement x, ECFieldElement y)
+        {
+            return Multiply(b).Add(x.Multiply(y));
+        }
+
+        public virtual ECFieldElement SquareMinusProduct(ECFieldElement x, ECFieldElement y)
+        {
+            return Square().Subtract(x.Multiply(y));
+        }
+
+        public virtual ECFieldElement SquarePlusProduct(ECFieldElement x, ECFieldElement y)
+        {
+            return Square().Add(x.Multiply(y));
+        }
+
         public virtual bool TestBitZero()
         {
             return ToBigInteger().TestBit(0);
@@ -162,6 +182,27 @@ namespace Org.BouncyCastle.Math.EC
             return new FpFieldElement(q, r, ModMult(x, b.ToBigInteger()));
         }
 
+        public override ECFieldElement MultiplyMinusProduct(ECFieldElement b, ECFieldElement x, ECFieldElement y)
+        {
+            BigInteger ax = this.x, bx = b.ToBigInteger(), xx = x.ToBigInteger(), yx = y.ToBigInteger();
+            BigInteger ab = ax.Multiply(bx);
+            BigInteger xy = xx.Multiply(yx);
+            return new FpFieldElement(q, r, ModReduce(ab.Subtract(xy)));
+        }
+
+        public override ECFieldElement MultiplyPlusProduct(ECFieldElement b, ECFieldElement x, ECFieldElement y)
+        {
+            BigInteger ax = this.x, bx = b.ToBigInteger(), xx = x.ToBigInteger(), yx = y.ToBigInteger();
+            BigInteger ab = ax.Multiply(bx);
+            BigInteger xy = xx.Multiply(yx);
+            BigInteger sum = ab.Add(xy);
+            if (r != null && r.SignValue < 0 && sum.BitLength > (q.BitLength << 1))
+            {
+                sum = sum.Subtract(q.ShiftLeft(q.BitLength));
+            }
+            return new FpFieldElement(q, r, ModReduce(sum));
+        }
+
         public override ECFieldElement Divide(
             ECFieldElement b)
         {
@@ -178,13 +219,33 @@ namespace Org.BouncyCastle.Math.EC
             return new FpFieldElement(q, r, ModMult(x, x));
         }
 
+        public override ECFieldElement SquareMinusProduct(ECFieldElement x, ECFieldElement y)
+        {
+            BigInteger ax = this.x, xx = x.ToBigInteger(), yx = y.ToBigInteger();
+            BigInteger aa = ax.Multiply(ax);
+            BigInteger xy = xx.Multiply(yx);
+            return new FpFieldElement(q, r, ModReduce(aa.Subtract(xy)));
+        }
+
+        public override ECFieldElement SquarePlusProduct(ECFieldElement x, ECFieldElement y)
+        {
+            BigInteger ax = this.x, xx = x.ToBigInteger(), yx = y.ToBigInteger();
+            BigInteger aa = ax.Multiply(ax);
+            BigInteger xy = xx.Multiply(yx);
+            BigInteger sum = aa.Add(xy);
+            if (r != null && r.SignValue < 0 && sum.BitLength > (q.BitLength << 1))
+            {
+                sum = sum.Subtract(q.ShiftLeft(q.BitLength));
+            }
+            return new FpFieldElement(q, r, ModReduce(sum));
+        }
+
         public override ECFieldElement Invert()
         {
             // TODO Modular inversion can be faster for a (Generalized) Mersenne Prime.
             return new FpFieldElement(q, r, ModInverse(x));
         }
 
-        // D.1.4 91
         /**
          * return a sqrt root - the routine verifies that the calculation
          * returns the right value - if none exists it returns null.
@@ -1134,6 +1195,29 @@ namespace Org.BouncyCastle.Math.EC
             return new F2mFieldElement(m, ks, x.ModMultiply(((F2mFieldElement)b).x, m, ks));
         }
 
+        public override ECFieldElement MultiplyMinusProduct(ECFieldElement b, ECFieldElement x, ECFieldElement y)
+        {
+            return MultiplyPlusProduct(b, x, y);
+        }
+
+        public override ECFieldElement MultiplyPlusProduct(ECFieldElement b, ECFieldElement x, ECFieldElement y)
+        {
+            LongArray ax = this.x, bx = ((F2mFieldElement)b).x, xx = ((F2mFieldElement)x).x, yx = ((F2mFieldElement)y).x;
+
+            LongArray ab = ax.Multiply(bx, m, ks);
+            LongArray xy = xx.Multiply(yx, m, ks);
+
+            if (ab == ax || ab == bx)
+            {
+                ab = (LongArray)ab.Copy();
+            }
+
+            ab.AddShiftedByWords(xy, 0);
+            ab.Reduce(m, ks);
+
+            return new F2mFieldElement(m, ks, ab);
+        }
+
         public override ECFieldElement Divide(
             ECFieldElement b)
         {
@@ -1153,6 +1237,29 @@ namespace Org.BouncyCastle.Math.EC
             return new F2mFieldElement(m, ks, x.ModSquare(m, ks));
         }
 
+        public override ECFieldElement SquareMinusProduct(ECFieldElement x, ECFieldElement y)
+        {
+            return SquarePlusProduct(x, y);
+        }
+
+        public override ECFieldElement SquarePlusProduct(ECFieldElement x, ECFieldElement y)
+        {
+            LongArray ax = this.x, xx = ((F2mFieldElement)x).x, yx = ((F2mFieldElement)y).x;
+
+            LongArray aa = ax.Square(m, ks);
+            LongArray xy = xx.Multiply(yx, m, ks);
+
+            if (aa == ax)
+            {
+                aa = (LongArray)aa.Copy();
+            }
+
+            aa.AddShiftedByWords(xy, 0);
+            aa.Reduce(m, ks);
+
+            return new F2mFieldElement(m, ks, aa);
+        }
+
         public override ECFieldElement Invert()
         {
             return new F2mFieldElement(this.m, this.ks, this.x.ModInverse(m, ks));
diff --git a/crypto/src/math/ec/ECPoint.cs b/crypto/src/math/ec/ECPoint.cs
index f2b0cdc27..257e0fd5d 100644
--- a/crypto/src/math/ec/ECPoint.cs
+++ b/crypto/src/math/ec/ECPoint.cs
@@ -617,7 +617,7 @@ namespace Org.BouncyCastle.Math.EC
                     ECFieldElement A = u.Square().Multiply(w).Subtract(vCubed).Subtract(Two(vSquaredV2));
 
                     ECFieldElement X3 = v.Multiply(A);
-                    ECFieldElement Y3 = vSquaredV2.Subtract(A).Multiply(u).Subtract(vCubed.Multiply(u2));
+                    ECFieldElement Y3 = vSquaredV2.Subtract(A).MultiplyMinusProduct(u, u2, vCubed);
                     ECFieldElement Z3 = vCubed.Multiply(w);
 
                     return new FpPoint(curve, X3, Y3, new ECFieldElement[] { Z3 }, IsCompressed);
@@ -714,7 +714,7 @@ namespace Org.BouncyCastle.Math.EC
                         ECFieldElement V = HSquared.Multiply(U1);
 
                         X3 = R.Square().Add(G).Subtract(Two(V));
-                        Y3 = V.Subtract(X3).Multiply(R).Subtract(S1.Multiply(G));
+                        Y3 = V.Subtract(X3).MultiplyMinusProduct(R, G, S1);
 
                         Z3 = H;
                         if (!Z1IsOne)
@@ -1038,7 +1038,6 @@ namespace Org.BouncyCastle.Math.EC
             return a.Add(b).Square().Subtract(aSquared).Subtract(bSquared);
         }
 
-        // D.3.2 pg 102 (see Note:)
         public override ECPoint Subtract(
             ECPoint b)
         {
@@ -1068,7 +1067,7 @@ namespace Org.BouncyCastle.Math.EC
         protected virtual ECFieldElement CalculateJacobianModifiedW(ECFieldElement Z, ECFieldElement ZSquared)
         {
             ECFieldElement a4 = this.Curve.A;
-            if (a4.IsZero)
+            if (a4.IsZero || Z.IsOne)
                 return a4;
 
             if (ZSquared == null)
@@ -1334,13 +1333,23 @@ namespace Org.BouncyCastle.Math.EC
                     ECFieldElement Y1 = this.RawYCoord, Z1 = this.RawZCoords[0];
                     ECFieldElement Y2 = b.RawYCoord, Z2 = b.RawZCoords[0];
 
+                    bool Z1IsOne = Z1.IsOne;
+                    ECFieldElement U1 = Y2, V1 = X2;
+                    if (!Z1IsOne)
+                    {
+                        U1 = U1.Multiply(Z1);
+                        V1 = V1.Multiply(Z1);
+                    }
+
                     bool Z2IsOne = Z2.IsOne;
+                    ECFieldElement U2 = Y1, V2 = X1;
+                    if (!Z2IsOne)
+                    {
+                        U2 = U2.Multiply(Z2);
+                        V2 = V2.Multiply(Z2);
+                    }
 
-                    ECFieldElement U1 = Z1.Multiply(Y2);
-                    ECFieldElement U2 = Z2IsOne ? Y1 : Y1.Multiply(Z2);
                     ECFieldElement U = U1.Add(U2);
-                    ECFieldElement V1 = Z1.Multiply(X2);
-                    ECFieldElement V2 = Z2IsOne ? X1 : X1.Multiply(Z2);
                     ECFieldElement V = V1.Add(V2);
 
                     if (V.IsZero)
@@ -1355,15 +1364,13 @@ namespace Org.BouncyCastle.Math.EC
 
                     ECFieldElement VSq = V.Square();
                     ECFieldElement VCu = VSq.Multiply(V);
-                    ECFieldElement W = Z2IsOne ? Z1 : Z1.Multiply(Z2);
+                    ECFieldElement W = Z1IsOne ? Z2 : Z2IsOne ? Z1 : Z1.Multiply(Z2);
                     ECFieldElement uv = U.Add(V);
-                    // TODO Delayed modular reduction for sum of products
-                    ECFieldElement A = uv.Multiply(U).Add(VSq.Multiply(curve.A)).Multiply(W).Add(VCu);
+                    ECFieldElement A = uv.MultiplyPlusProduct(U, VSq, curve.A).Multiply(W).Add(VCu);
 
                     ECFieldElement X3 = V.Multiply(A);
                     ECFieldElement VSqZ2 = Z2IsOne ? VSq : VSq.Multiply(Z2);
-                    // TODO Delayed modular reduction for sum of products
-                    ECFieldElement Y3 = U.Multiply(X1).Add(Y1.Multiply(V)).Multiply(VSqZ2).Add(A.Multiply(uv));
+                    ECFieldElement Y3 = U.MultiplyPlusProduct(X1, V, Y1).MultiplyPlusProduct(VSqZ2, uv, A);
                     ECFieldElement Z3 = VCu.Multiply(W);
 
                     return new F2mPoint(curve, X3, Y3, new ECFieldElement[] { Z3 }, IsCompressed);
@@ -1450,8 +1457,7 @@ namespace Org.BouncyCastle.Math.EC
                             ABZ2 = ABZ2.Multiply(Z2);
                         }
 
-                        // TODO Delayed modular reduction for sum of products
-                        L3 = AU2.Add(B).Square().Add(ABZ2.Multiply(L1.Add(Z1)));
+                        L3 = AU2.Add(B).SquarePlusProduct(ABZ2, L1.Add(Z1));
 
                         Z3 = ABZ2;
                         if (!Z1IsOne)
@@ -1559,8 +1565,7 @@ namespace Org.BouncyCastle.Math.EC
                     ECFieldElement L1 = Y1.Divide(X1).Add(X1);
 
                     ECFieldElement X3 = L1.Square().Add(L1).Add(curve.A);
-                    // TODO Delayed modular reduction for sum of products
-                    ECFieldElement Y3 = X1.Square().Add(X3.Multiply(L1.AddOne()));
+                    ECFieldElement Y3 = X1.SquarePlusProduct(X3, L1.AddOne());
 
                     return new F2mPoint(curve, X3, Y3, IsCompressed);
                 }
@@ -1577,12 +1582,10 @@ namespace Org.BouncyCastle.Math.EC
                     ECFieldElement V = X1Z1;
                     ECFieldElement vSquared = V.Square();
                     ECFieldElement sv = S.Add(V);
-                    // TODO Delayed modular reduction for sum of products
-                    ECFieldElement h = sv.Multiply(S).Add(curve.A.Multiply(vSquared));
+                    ECFieldElement h = sv.MultiplyPlusProduct(S, vSquared, curve.A);
 
                     ECFieldElement X3 = V.Multiply(h);
-                    // TODO Delayed modular reduction for sum of products
-                    ECFieldElement Y3 = h.Multiply(sv).Add(X1Sq.Square().Multiply(V));
+                    ECFieldElement Y3 = X1Sq.Square().MultiplyPlusProduct(V, h, sv);
                     ECFieldElement Z3 = V.Multiply(vSquared);
 
                     return new F2mPoint(curve, X3, Y3, new ECFieldElement[] { Z3 }, IsCompressed);
@@ -1610,19 +1613,17 @@ namespace Org.BouncyCastle.Math.EC
                     if (b.BitLength < (curve.FieldSize >> 1))
                     {
                         ECFieldElement t1 = L1.Add(X1).Square();
-                        ECFieldElement t4;
+                        ECFieldElement t2;
                         if (b.IsOne)
                         {
-                            t4 = aZ1Sq.Add(Z1Sq).Square();
+                            t2 = aZ1Sq.Add(Z1Sq).Square();
                         }
                         else
                         {
-                            // TODO t2/t3 can be calculated with one square if we pre-compute sqrt(b)
-                            ECFieldElement t2 = aZ1Sq.Square();
-                            ECFieldElement t3 = b.Multiply(Z1Sq.Square());
-                            t4 = t2.Add(t3);
+                            // TODO Can be calculated with one square if we pre-compute sqrt(b)
+                            t2 = aZ1Sq.SquarePlusProduct(b, Z1Sq.Square());
                         }
-                        L3 = t1.Add(T).Add(Z1Sq).Multiply(t1).Add(t4).Add(X3);
+                        L3 = t1.Add(T).Add(Z1Sq).Multiply(t1).Add(t2).Add(X3);
                         if (a.IsZero)
                         {
                             L3 = L3.Add(Z3);
@@ -1635,8 +1636,7 @@ namespace Org.BouncyCastle.Math.EC
                     else
                     {
                         ECFieldElement X1Z1 = Z1IsOne ? X1 : X1.Multiply(Z1);
-                        // TODO Delayed modular reduction for sum of products
-                        L3 = X1Z1.Square().Add(T.Multiply(L1Z1)).Add(X3).Add(Z3);
+                        L3 = X1Z1.SquarePlusProduct(T, L1Z1).Add(X3).Add(Z3);
                     }
 
                     return new F2mPoint(curve, X3, L3, new ECFieldElement[] { Z3 }, IsCompressed);
@@ -1687,8 +1687,7 @@ namespace Org.BouncyCastle.Math.EC
 
                     ECFieldElement T = curve.A.Multiply(Z1Sq).Add(L1Sq).Add(L1Z1);
                     ECFieldElement L2plus1 = L2.AddOne();
-                    // TODO Delayed modular reduction for sum of products
-                    ECFieldElement A = curve.A.Add(L2plus1).Multiply(Z1Sq).Add(L1Sq).Multiply(T).Add(X1Sq.Multiply(Z1Sq));
+                    ECFieldElement A = curve.A.Add(L2plus1).Multiply(Z1Sq).Add(L1Sq).MultiplyPlusProduct(T, X1Sq, Z1Sq);
                     ECFieldElement X2Z1Sq = X2.Multiply(Z1Sq);
                     ECFieldElement B = X2Z1Sq.Add(T).Square();
 
@@ -1709,8 +1708,7 @@ namespace Org.BouncyCastle.Math.EC
 
                     ECFieldElement X3 = A.Square().Multiply(X2Z1Sq);
                     ECFieldElement Z3 = A.Multiply(B).Multiply(Z1Sq);
-                    // TODO Delayed modular reduction for sum of products
-                    ECFieldElement L3 = A.Add(B).Square().Multiply(T).Add(L2plus1.Multiply(Z3));
+                    ECFieldElement L3 = A.Add(B).Square().MultiplyPlusProduct(T, L2plus1, Z3);
 
                     return new F2mPoint(curve, X3, L3, new ECFieldElement[] { Z3 }, IsCompressed);
                 }
diff --git a/crypto/src/math/ec/LongArray.cs b/crypto/src/math/ec/LongArray.cs
index 29f971c46..c4e3dacbc 100644
--- a/crypto/src/math/ec/LongArray.cs
+++ b/crypto/src/math/ec/LongArray.cs
@@ -1335,6 +1335,158 @@ namespace Org.BouncyCastle.Math.EC
             return ReduceResult(c, ci[1], cLen, m, ks);
         }
 
+        public LongArray ModReduce(int m, int[] ks)
+        {
+            long[] buf = Arrays.Clone(m_ints);
+            int rLen = ReduceInPlace(buf, 0, buf.Length, m, ks);
+            return new LongArray(buf, 0, rLen);
+        }
+
+        public LongArray Multiply(LongArray other, int m, int[] ks)
+        {
+            /*
+             * Find out the degree of each argument and handle the zero cases
+             */
+            int aDeg = Degree();
+            if (aDeg == 0)
+            {
+                return this;
+            }
+            int bDeg = other.Degree();
+            if (bDeg == 0)
+            {
+                return other;
+            }
+
+            /*
+             * Swap if necessary so that A is the smaller argument
+             */
+            LongArray A = this, B = other;
+            if (aDeg > bDeg)
+            {
+                A = other; B = this;
+                int tmp = aDeg; aDeg = bDeg; bDeg = tmp;
+            }
+
+            /*
+             * Establish the word lengths of the arguments and result
+             */
+            int aLen = (int)((uint)(aDeg + 63) >> 6);
+            int bLen = (int)((uint)(bDeg + 63) >> 6);
+            int cLen = (int)((uint)(aDeg + bDeg + 62) >> 6);
+
+            if (aLen == 1)
+            {
+                long a0 = A.m_ints[0];
+                if (a0 == 1L)
+                {
+                    return B;
+                }
+
+                /*
+                 * Fast path for small A, with performance dependent only on the number of set bits
+                 */
+                long[] c0 = new long[cLen];
+                MultiplyWord(a0, B.m_ints, bLen, c0, 0);
+
+                /*
+                 * Reduce the raw answer against the reduction coefficients
+                 */
+                //return ReduceResult(c0, 0, cLen, m, ks);
+                return new LongArray(c0, 0, cLen);
+            }
+
+            /*
+             * Determine if B will get bigger during shifting
+             */
+            int bMax = (int)((uint)(bDeg + 7 + 63) >> 6);
+
+            /*
+             * Lookup table for the offset of each B in the tables
+             */
+            int[] ti = new int[16];
+
+            /*
+             * Precompute table of all 4-bit products of B
+             */
+            long[] T0 = new long[bMax << 4];
+            int tOff = bMax;
+            ti[1] = tOff;
+            Array.Copy(B.m_ints, 0, T0, tOff, bLen);
+            for (int i = 2; i < 16; ++i)
+            {
+                ti[i] = (tOff += bMax);
+                if ((i & 1) == 0)
+                {
+                    ShiftUp(T0, (int)((uint)tOff >> 1), T0, tOff, bMax, 1);
+                }
+                else
+                {
+                    Add(T0, bMax, T0, tOff - bMax, T0, tOff, bMax);
+                }
+            }
+
+            /*
+             * Second table with all 4-bit products of B shifted 4 bits
+             */
+            long[] T1 = new long[T0.Length];
+            ShiftUp(T0, 0, T1, 0, T0.Length, 4);
+            //        ShiftUp(T0, bMax, T1, bMax, tOff, 4);
+
+            long[] a = A.m_ints;
+            long[] c = new long[cLen << 3];
+
+            int MASK = 0xF;
+
+            /*
+             * Lopez-Dahab (Modified) algorithm
+             */
+
+            for (int aPos = 0; aPos < aLen; ++aPos)
+            {
+                long aVal = a[aPos];
+                int cOff = aPos;
+                for (; ; )
+                {
+                    int u = (int)aVal & MASK;
+                    aVal = (long)((ulong)aVal >> 4);
+                    int v = (int)aVal & MASK;
+                    AddBoth(c, cOff, T0, ti[u], T1, ti[v], bMax);
+                    aVal = (long)((ulong)aVal >> 4);
+                    if (aVal == 0L)
+                    {
+                        break;
+                    }
+                    cOff += cLen;
+                }
+            }
+
+            {
+                int cOff = c.Length;
+                while ((cOff -= cLen) != 0)
+                {
+                    AddShiftedUp(c, cOff - cLen, c, cOff, cLen, 8);
+                }
+            }
+
+            /*
+             * Finally the raw answer is collected, reduce it against the reduction coefficients
+             */
+            //return ReduceResult(c, 0, cLen, m, ks);
+            return new LongArray(c, 0, cLen);
+        }
+
+        public void Reduce(int m, int[] ks)
+        {
+            long[] buf = m_ints;
+            int rLen = ReduceInPlace(buf, 0, buf.Length, m, ks);
+            if (rLen < buf.Length)
+            {
+                m_ints = new long[rLen];
+                Array.Copy(buf, 0, m_ints, 0, rLen);
+            }
+        }
+
         private static LongArray ReduceResult(long[] buf, int off, int len, int m, int[] ks)
         {
             int rLen = ReduceInPlace(buf, off, len, m, ks);
@@ -1546,6 +1698,28 @@ namespace Org.BouncyCastle.Math.EC
             return new LongArray(r, 0, len);
         }
 
+        public LongArray Square(int m, int[] ks)
+        {
+            int len = GetUsedLength();
+            if (len == 0)
+            {
+                return this;
+            }
+
+            int _2len = len << 1;
+            long[] r = new long[_2len];
+
+            int pos = 0;
+            while (pos < _2len)
+            {
+                long mi = m_ints[(uint)pos >> 1];
+                r[pos++] = Interleave2_32to64((int)mi);
+                r[pos++] = Interleave2_32to64((int)((ulong)mi >> 32));
+            }
+
+            return new LongArray(r, 0, r.Length);
+        }
+
         private static void SquareInPlace(long[] x, int xLen, int m, int[] ks)
         {
             int pos = xLen << 1;