Optimize Curve25519 point operations

2 files changed, 214 insertions, 111 deletions

diff --git a/crypto/src/math/ec/custom/djb/Curve25519Field.cs b/crypto/src/math/ec/custom/djb/Curve25519Field.cs
index 084ca96af..809e51b80 100644
--- a/crypto/src/math/ec/custom/djb/Curve25519Field.cs
+++ b/crypto/src/math/ec/custom/djb/Curve25519Field.cs
@@ -10,7 +10,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
         // 2^255 - 2^4 - 2^1 - 1
         internal static readonly uint[] P = new uint[]{ 0xFFFFFFED, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF,
             0xFFFFFFFF, 0x7FFFFFFF };
-        private const int P7 = 0x7FFFFFFF;
+        private const uint P7 = 0x7FFFFFFF;
         private static readonly uint[] PExt = new uint[]{ 0x00000169, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
             0x00000000, 0x00000000, 0x00000000, 0xFFFFFFED, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF,
             0xFFFFFFFF, 0x3FFFFFFF };
@@ -21,7 +21,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
             Nat256.Add(x, y, z);
             if (Nat256.Gte(z, P))
             {
-                AddPInvTo(z);
+                SubPFrom(z);
             }
         }
 
@@ -39,7 +39,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
             Nat.Inc(8, x, z);
             if (Nat256.Gte(z, P))
             {
-                AddPInvTo(z);
+                SubPFrom(z);
             }
         }
 
@@ -73,6 +73,15 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
             Reduce(tt, z);
         }
 
+        public static void MultiplyAddToExt(uint[] x, uint[] y, uint[] zz)
+        {
+            Nat256.MulAddTo(x, y, zz);
+            if (Nat.Gte(16, zz, PExt))
+            {
+                SubPExtFrom(zz);
+            }
+        }
+
         public static void Negate(uint[] x, uint[] z)
         {
             if (Nat256.IsZero(x))
@@ -92,13 +101,29 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
             uint xx07 = xx[7];
             Nat.ShiftUpBit(8, xx, 8, xx07, z, 0);
             uint c = Nat256.MulByWordAddTo(PInv, xx, z) << 1;
-            uint z07 = z[7];
-            z[7] = z07 & P7;
-            c += (z07 >> 31) - (xx07 >> 31);
-            Nat.AddWordTo(8, c * PInv, z);
-            if (Nat256.Gte(z, P))
+            uint z7 = z[7];
+            c += (z7 >> 31) - (xx07 >> 31);
+            z7 &= P7;
+            z7 += Nat.AddWordTo(7, c * PInv, z);
+            z[7] = z7;
+            if (z7 >= P7 && Nat256.Gte(z, P))
+            {
+                SubPFrom(z);
+            }
+        }
+
+        public static void Reduce27(uint x, uint[] z)
+        {
+            Debug.Assert(x >> 26 == 0);
+
+            uint z7 = z[7];
+            uint c = (x << 1 | z7 >> 31);
+            z7 &= P7;
+            z7 += Nat.AddWordTo(7, c * PInv, z);
+            z[7] = z7;
+            if (z7 >= P7 && Nat256.Gte(z, P))
             {
-                AddPInvTo(z);
+                SubPFrom(z);
             }
         }
 
@@ -111,7 +136,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
 
         public static void SquareN(uint[] x, int n, uint[] z)
         {
-    //        assert n > 0;
+            Debug.Assert(n > 0);
 
             uint[] tt = Nat256.CreateExt();
             Nat256.Square(x, tt);
@@ -129,7 +154,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
             int c = Nat256.Sub(x, y, z);
             if (c != 0)
             {
-                SubPInvFrom(z);
+                AddPTo(z);
             }
         }
 
@@ -147,66 +172,82 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
             Nat.ShiftUpBit(8, x, 0, z);
             if (Nat256.Gte(z, P))
             {
-                AddPInvTo(z);
+                SubPFrom(z);
             }
         }
 
-        private static void AddPExtTo(uint[] zz)
+        private static uint AddPTo(uint[] z)
         {
-            ulong c = (ulong)zz[0] + PExt[0];
-            zz[0] = (uint)c;
+            long c = (long)z[0] - PInv;
+            z[0] = (uint)c;
             c >>= 32;
-
-            int i = 1 - (int)c;
-            i = (i << 3) - i;
-
-            while (++i < 16)
+            if (c != 0)
             {
-                c += (ulong)zz[i] + PExt[i];
-                zz[i] = (uint)c;
-                c >>= 32;
+                c = Nat.DecAt(7, z, 1);
             }
+            c += (long)z[7] + (P7 + 1);
+            z[7] = (uint)c;
+            c >>= 32;
+            return (uint)c;
         }
 
-        private static void SubPExtFrom(uint[] zz)
+        private static uint AddPExtTo(uint[] zz)
         {
-            long c = (long)zz[0] - PExt[0];
+            long c = (long)zz[0] + PExt[0];
             zz[0] = (uint)c;
             c >>= 32;
-
-            int i = 1 + (int)c;
-            i = (i << 3) - i;
-
-            while (++i < 16)
+            if (c != 0)
             {
-                c += (long)zz[i] - PExt[i];
-                zz[i] = (uint)c;
-                c >>= 32;
+                c = Nat.IncAt(8, zz, 1);
             }
+            c += (long)zz[8] - PInv;
+            zz[8] = (uint)c;
+            c >>= 32;
+            if (c != 0)
+            {
+                c = Nat.DecAt(15, zz, 9);
+            }
+            c += (long)zz[15] + (PExt[15] + 1);
+            zz[15] = (uint)c;
+            c >>= 32;
+            return (uint)c;
         }
 
-        private static void AddPInvTo(uint[] z)
+        private static int SubPFrom(uint[] z)
         {
-            ulong c = (ulong)z[0] + PInv;
+            long c = (long)z[0] + PInv;
             z[0] = (uint)c;
             c >>= 32;
             if (c != 0)
             {
-                Nat.IncAt(8, z, 1);
+                c = Nat.IncAt(7, z, 1);
             }
-            z[7] &= P7;
+            c += (long)z[7] - (P7 + 1);
+            z[7] = (uint)c;
+            c >>= 32;
+            return (int)c;
         }
 
-        private static void SubPInvFrom(uint[] z)
+        private static int SubPExtFrom(uint[] zz)
         {
-            long c = (long)z[0] - PInv;
-            z[0] = (uint)c;
+            long c = (long)zz[0] - PExt[0];
+            zz[0] = (uint)c;
             c >>= 32;
             if (c != 0)
             {
-                Nat.DecAt(8, z, 1);
+                c = Nat.DecAt(8, zz, 1);
             }
-            z[7] &= P7;
+            c += (long)zz[8] + PInv;
+            zz[8] = (uint)c;
+            c >>= 32;
+            if (c != 0)
+            {
+                c = Nat.IncAt(15, zz, 9);
+            }
+            c += (long)zz[15] - (PExt[15] + 1);
+            zz[15] = (uint)c;
+            c >>= 32;
+            return (int)c;
         }
     }
 }
diff --git a/crypto/src/math/ec/custom/djb/Curve25519Point.cs b/crypto/src/math/ec/custom/djb/Curve25519Point.cs
index 65b6792eb..f3da59d16 100644
--- a/crypto/src/math/ec/custom/djb/Curve25519Point.cs
+++ b/crypto/src/math/ec/custom/djb/Curve25519Point.cs
@@ -5,7 +5,7 @@ using Org.BouncyCastle.Math.EC.Custom.Sec;
 namespace Org.BouncyCastle.Math.EC.Custom.Djb
 {
     internal class Curve25519Point
-        :   ECPointBase
+        : ECPointBase
     {
         /**
          * Create a point which encodes with point compression.
@@ -14,10 +14,10 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
          * @param x affine x co-ordinate
          * @param y affine y co-ordinate
          * 
-         * @deprecated Use ECCurve.createPoint to construct points
+         * @deprecated Use ECCurve.CreatePoint to construct points
          */
         public Curve25519Point(ECCurve curve, ECFieldElement x, ECFieldElement y)
-            :   this(curve, x, y, false)
+            : this(curve, x, y, false)
         {
         }
 
@@ -32,7 +32,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
          * @deprecated per-point compression property will be removed, refer {@link #getEncoded(bool)}
          */
         public Curve25519Point(ECCurve curve, ECFieldElement x, ECFieldElement y, bool withCompression)
-            :   base(curve, x, y, withCompression)
+            : base(curve, x, y, withCompression)
         {
             if ((x == null) != (y == null))
                 throw new ArgumentException("Exactly one of the field elements is null");
@@ -74,48 +74,65 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
 
             ECCurve curve = this.Curve;
 
-            ECFieldElement X1 = this.RawXCoord, Y1 = this.RawYCoord;
-            ECFieldElement X2 = b.RawXCoord, Y2 = b.RawYCoord;
+            Curve25519FieldElement X1 = (Curve25519FieldElement)this.RawXCoord, Y1 = (Curve25519FieldElement)this.RawYCoord,
+                Z1 = (Curve25519FieldElement)this.RawZCoords[0];
+            Curve25519FieldElement X2 = (Curve25519FieldElement)b.RawXCoord, Y2 = (Curve25519FieldElement)b.RawYCoord,
+                Z2 = (Curve25519FieldElement)b.RawZCoords[0];
 
-            ECFieldElement Z1 = this.RawZCoords[0];
-            ECFieldElement Z2 = b.RawZCoords[0];
+            uint c;
+            uint[] tt1 = Nat256.CreateExt();
+            uint[] t2 = Nat256.Create();
+            uint[] t3 = Nat256.Create();
+            uint[] t4 = Nat256.Create();
 
             bool Z1IsOne = Z1.IsOne;
-
-            ECFieldElement Z1Squared, U2, S2;
+            uint[] U2, S2;
             if (Z1IsOne)
             {
-                Z1Squared = Z1; U2 = X2; S2 = Y2;
+                U2 = X2.x;
+                S2 = Y2.x;
             }
             else
             {
-                Z1Squared = Z1.Square();
-                U2 = Z1Squared.Multiply(X2);
-                ECFieldElement Z1Cubed = Z1Squared.Multiply(Z1);
-                S2 = Z1Cubed.Multiply(Y2);
+                S2 = t3;
+                Curve25519Field.Square(Z1.x, S2);
+
+                U2 = t2;
+                Curve25519Field.Multiply(S2, X2.x, U2);
+
+                Curve25519Field.Multiply(S2, Z1.x, S2);
+                Curve25519Field.Multiply(S2, Y2.x, S2);
             }
 
             bool Z2IsOne = Z2.IsOne;
-            ECFieldElement Z2Squared, U1, S1;
+            uint[] U1, S1;
             if (Z2IsOne)
             {
-                Z2Squared = Z2; U1 = X1; S1 = Y1;
+                U1 = X1.x;
+                S1 = Y1.x;
             }
             else
             {
-                Z2Squared = Z2.Square();
-                U1 = Z2Squared.Multiply(X1); 
-                ECFieldElement Z2Cubed = Z2Squared.Multiply(Z2);
-                S1 = Z2Cubed.Multiply(Y1);
+                S1 = t4;
+                Curve25519Field.Square(Z2.x, S1);
+
+                U1 = tt1;
+                Curve25519Field.Multiply(S1, X1.x, U1);
+
+                Curve25519Field.Multiply(S1, Z2.x, S1);
+                Curve25519Field.Multiply(S1, Y1.x, S1);
             }
 
-            ECFieldElement H = U1.Subtract(U2);
-            ECFieldElement R = S1.Subtract(S2);
+            uint[] H = Nat256.Create();
+            Curve25519Field.Subtract(U1, U2, H);
+
+            uint[] R = t2;
+            Curve25519Field.Subtract(S1, S2, R);
 
             // Check if b == this or b == -this
-            if (H.IsZero)
+            if (Nat256.IsZero(H))
             {
-                if (R.IsZero)
+                if (Nat256.IsZero(R))
                 {
                     // this == b, i.e. this must be doubled
                     return this.Twice();
@@ -125,29 +142,46 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
                 return curve.Infinity;
             }
 
-            ECFieldElement HSquared = H.Square();
-            ECFieldElement G = HSquared.Multiply(H);
-            ECFieldElement V = HSquared.Multiply(U1);
+            uint[] HSquared = Nat256.Create();
+            Curve25519Field.Square(H, HSquared);
+
+            uint[] G = Nat256.Create();
+            Curve25519Field.Multiply(HSquared, H, G);
+
+            uint[] V = t3;
+            Curve25519Field.Multiply(HSquared, U1, V);
+
+            Curve25519Field.Negate(G, G);
+            Nat256.Mul(S1, G, tt1);
 
-            ECFieldElement X3 = R.Square().Add(G).Subtract(Two(V));
-            ECFieldElement Y3 = V.Subtract(X3).MultiplyMinusProduct(R, G, S1);
+            c = Nat256.AddBothTo(V, V, G);
+            Curve25519Field.Reduce27(c, G);
 
-            ECFieldElement Z3 = H;
+            Curve25519FieldElement X3 = new Curve25519FieldElement(t4);
+            Curve25519Field.Square(R, X3.x);
+            Curve25519Field.Subtract(X3.x, G, X3.x);
+
+            Curve25519FieldElement Y3 = new Curve25519FieldElement(G);
+            Curve25519Field.Subtract(V, X3.x, Y3.x);
+            Curve25519Field.MultiplyAddToExt(Y3.x, R, tt1);
+            Curve25519Field.Reduce(tt1, Y3.x);
+
+            Curve25519FieldElement Z3 = new Curve25519FieldElement(H);
             if (!Z1IsOne)
             {
-                Z3 = Z3.Multiply(Z1);
+                Curve25519Field.Multiply(Z3.x, Z1.x, Z3.x);
             }
             if (!Z2IsOne)
             {
-                Z3 = Z3.Multiply(Z2);
+                Curve25519Field.Multiply(Z3.x, Z2.x, Z3.x);
             }
 
-            ECFieldElement Z3Squared = (Z3 == H) ? HSquared : null;
+            uint[] Z3Squared = (Z1IsOne && Z2IsOne) ? HSquared : null;
 
             // TODO If the result will only be used in a subsequent addition, we don't need W3
-            ECFieldElement W3 = CalculateJacobianModifiedW(Z3, Z3Squared);
+            Curve25519FieldElement W3 = CalculateJacobianModifiedW((Curve25519FieldElement)Z3, Z3Squared);
 
-            ECFieldElement[] zs = new ECFieldElement[]{ Z3, W3 };
+            ECFieldElement[] zs = new ECFieldElement[] { Z3, W3 };
 
             return new Curve25519Point(curve, X3, Y3, zs, IsCompressed);
         }
@@ -160,7 +194,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
             ECCurve curve = this.Curve;
 
             ECFieldElement Y1 = this.RawYCoord;
-            if (Y1.IsZero) 
+            if (Y1.IsZero)
                 return curve.Infinity;
 
             return TwiceJacobianModified(true);
@@ -176,7 +210,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
                 return Twice();
 
             ECFieldElement Y1 = this.RawYCoord;
-            if (Y1.IsZero) 
+            if (Y1.IsZero)
                 return b;
 
             return TwiceJacobianModified(false).Add(b);
@@ -190,16 +224,6 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
             return TwiceJacobianModified(false).Add(this);
         }
 
-        protected virtual ECFieldElement Two(ECFieldElement x)
-        {
-            return x.Add(x);
-        }
-
-        protected virtual ECFieldElement Three(ECFieldElement x)
-        {
-            return Two(x).Add(x);
-        }
-
         public override ECPoint Subtract(ECPoint b)
         {
             if (b.IsInfinity)
@@ -216,47 +240,85 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
             return new Curve25519Point(Curve, RawXCoord, RawYCoord.Negate(), RawZCoords, IsCompressed);
         }
 
-        protected virtual ECFieldElement CalculateJacobianModifiedW(ECFieldElement Z, ECFieldElement ZSquared)
+        protected virtual Curve25519FieldElement CalculateJacobianModifiedW(Curve25519FieldElement Z, uint[] ZSquared)
         {
-            ECFieldElement a4 = this.Curve.A;
+            Curve25519FieldElement a4 = (Curve25519FieldElement)this.Curve.A;
             if (Z.IsOne)
                 return a4;
 
+            Curve25519FieldElement W = new Curve25519FieldElement();
             if (ZSquared == null)
             {
-                ZSquared = Z.Square();
+                ZSquared = W.x;
+                Curve25519Field.Square(Z.x, ZSquared);
             }
-
-            return ZSquared.Square().Multiply(a4);
+            Curve25519Field.Square(ZSquared, W.x);
+            Curve25519Field.Multiply(W.x, a4.x, W.x);
+            return W;
         }
 
-        protected virtual ECFieldElement GetJacobianModifiedW()
+        protected virtual Curve25519FieldElement GetJacobianModifiedW()
         {
             ECFieldElement[] ZZ = this.RawZCoords;
-            ECFieldElement W = ZZ[1];
+            Curve25519FieldElement W = (Curve25519FieldElement)ZZ[1];
             if (W == null)
             {
                 // NOTE: Rarely, TwicePlus will result in the need for a lazy W1 calculation here
-                ZZ[1] = W = CalculateJacobianModifiedW(ZZ[0], null);
+                ZZ[1] = W = CalculateJacobianModifiedW((Curve25519FieldElement)ZZ[0], null);
             }
             return W;
         }
 
         protected virtual Curve25519Point TwiceJacobianModified(bool calculateW)
         {
-            ECFieldElement X1 = this.RawXCoord, Y1 = this.RawYCoord, Z1 = this.RawZCoords[0], W1 = GetJacobianModifiedW();
-
-            ECFieldElement X1Squared = X1.Square();
-            ECFieldElement M = Three(X1Squared).Add(W1);
-            ECFieldElement _2Y1 = Two(Y1);
-            ECFieldElement _2Y1Squared = _2Y1.Multiply(Y1);
-            ECFieldElement S = Two(X1.Multiply(_2Y1Squared));
-            ECFieldElement X3 = M.Square().Subtract(Two(S));
-            ECFieldElement _4T = _2Y1Squared.Square();
-            ECFieldElement _8T = Two(_4T);
-            ECFieldElement Y3 = M.Multiply(S.Subtract(X3)).Subtract(_8T);
-            ECFieldElement W3 = calculateW ? Two(_8T.Multiply(W1)) : null;
-            ECFieldElement Z3 = Z1.IsOne ? _2Y1 : _2Y1.Multiply(Z1);
+            Curve25519FieldElement X1 = (Curve25519FieldElement)this.RawXCoord, Y1 = (Curve25519FieldElement)this.RawYCoord,
+                Z1 = (Curve25519FieldElement)this.RawZCoords[0], W1 = GetJacobianModifiedW();
+
+            uint c;
+
+            uint[] M = Nat256.Create();
+            Curve25519Field.Square(X1.x, M);
+            c = Nat256.AddBothTo(M, M, M);
+            c += Nat256.AddTo(W1.x, M);
+            Curve25519Field.Reduce27(c, M);
+
+            uint[] _2Y1 = Nat256.Create();
+            Curve25519Field.Twice(Y1.x, _2Y1);
+
+            uint[] _2Y1Squared = Nat256.Create();
+            Curve25519Field.Multiply(_2Y1, Y1.x, _2Y1Squared);
+
+            uint[] S = Nat256.Create();
+            Curve25519Field.Multiply(_2Y1Squared, X1.x, S);
+            Curve25519Field.Twice(S, S);
+
+            uint[] _8T = Nat256.Create();
+            Curve25519Field.Square(_2Y1Squared, _8T);
+            Curve25519Field.Twice(_8T, _8T);
+
+            Curve25519FieldElement X3 = new Curve25519FieldElement(_2Y1Squared);
+            Curve25519Field.Square(M, X3.x);
+            Curve25519Field.Subtract(X3.x, S, X3.x);
+            Curve25519Field.Subtract(X3.x, S, X3.x);
+
+            Curve25519FieldElement Y3 = new Curve25519FieldElement(S);
+            Curve25519Field.Subtract(S, X3.x, Y3.x);
+            Curve25519Field.Multiply(Y3.x, M, Y3.x);
+            Curve25519Field.Subtract(Y3.x, _8T, Y3.x);
+
+            Curve25519FieldElement Z3 = new Curve25519FieldElement(_2Y1);
+            if (!Nat256.IsOne(Z1.x))
+            {
+                Curve25519Field.Multiply(Z3.x, Z1.x, Z3.x);
+            }
+
+            Curve25519FieldElement W3 = null;
+            if (calculateW)
+            {
+                W3 = new Curve25519FieldElement(_8T);
+                Curve25519Field.Multiply(W3.x, W1.x, W3.x);
+                Curve25519Field.Twice(W3.x, W3.x);
+            }
 
             return new Curve25519Point(this.Curve, X3, Y3, new ECFieldElement[] { Z3, W3 }, IsCompressed);
         }