ECC: Binary field perf. opt.

13 files changed, 543 insertions, 287 deletions

diff --git a/crypto/src/math/ec/custom/sec/SecT113Field.cs b/crypto/src/math/ec/custom/sec/SecT113Field.cs
index 3c9e0938d..56738a219 100644
--- a/crypto/src/math/ec/custom/sec/SecT113Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT113Field.cs
@@ -87,14 +87,14 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Multiply(ulong[] x, ulong[] y, ulong[] z)
         {
-            ulong[] tt = Nat128.CreateExt64();
+            ulong[] tt = new ulong[8];
             ImplMultiply(x, y, tt);
             Reduce(tt, z);
         }
 
         public static void MultiplyAddToExt(ulong[] x, ulong[] y, ulong[] zz)
         {
-            ulong[] tt = Nat128.CreateExt64();
+            ulong[] tt = new ulong[8];
             ImplMultiply(x, y, tt);
             AddExt(zz, tt, zz);
         }
@@ -180,11 +180,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             g1  = ((g0 >> 57) ^ (g1 << 7)) & M57;
             g0 &= M57;
 
+            ulong[] u = zz;
             ulong[] H = new ulong[6];
 
-            ImplMulw(f0, g0, H, 0);               // H(0)       57/56 bits                                
-            ImplMulw(f1, g1, H, 2);               // H(INF)     57/54 bits                                
-            ImplMulw(f0 ^ f1, g0 ^ g1, H, 4);     // H(1)       57/56 bits
+            ImplMulw(u, f0, g0, H, 0);              // H(0)       57/56 bits                                
+            ImplMulw(u, f1, g1, H, 2);              // H(INF)     57/54 bits                                
+            ImplMulw(u, f0 ^ f1, g0 ^ g1, H, 4);    // H(1)       57/56 bits
 
             ulong r  = H[1] ^ H[2];
             ulong z0 = H[0],
@@ -198,12 +199,11 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             zz[3] = (z3 >> 21);
         }
 
-        protected static void ImplMulw(ulong x, ulong y, ulong[] z, int zOff)
+        protected static void ImplMulw(ulong[] u, ulong x, ulong y, ulong[] z, int zOff)
         {
             Debug.Assert(x >> 57 == 0);
             Debug.Assert(y >> 57 == 0);
 
-            ulong[] u = new ulong[8];
             //u[0] = 0;
             u[1] = y;
             u[2] = u[1] << 1;
@@ -237,8 +237,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         protected static void ImplSquare(ulong[] x, ulong[] zz)
         {
-            Interleave.Expand64To128(x[0], zz, 0);
-            Interleave.Expand64To128(x[1], zz, 2);
+            Interleave.Expand64To128(x, 0, 2, zz, 0);
         }
     }
 }
diff --git a/crypto/src/math/ec/custom/sec/SecT131Field.cs b/crypto/src/math/ec/custom/sec/SecT131Field.cs
index db703d9e0..adf4f0448 100644
--- a/crypto/src/math/ec/custom/sec/SecT131Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT131Field.cs
@@ -93,14 +93,14 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Multiply(ulong[] x, ulong[] y, ulong[] z)
         {
-            ulong[] tt = Nat192.CreateExt64();
+            ulong[] tt = new ulong[8];
             ImplMultiply(x, y, tt);
             Reduce(tt, z);
         }
 
         public static void MultiplyAddToExt(ulong[] x, ulong[] y, ulong[] zz)
         {
-            ulong[] tt = Nat192.CreateExt64();
+            ulong[] tt = new ulong[8];
             ImplMultiply(x, y, tt);
             AddExt(zz, tt, zz);
         }
@@ -214,21 +214,22 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             g1  = ((g0 >> 44) ^ (g1 << 20)) & M44;
             g0 &= M44;
 
+            ulong[] u = zz;
             ulong[] H = new ulong[10];
 
-            ImplMulw(f0, g0, H, 0);               // H(0)       44/43 bits
-            ImplMulw(f2, g2, H, 2);               // H(INF)     44/41 bits
+            ImplMulw(u, f0, g0, H, 0);              // H(0)       44/43 bits
+            ImplMulw(u, f2, g2, H, 2);              // H(INF)     44/41 bits
 
             ulong t0 = f0 ^ f1 ^ f2;
             ulong t1 = g0 ^ g1 ^ g2;
 
-            ImplMulw(t0, t1, H, 4);               // H(1)       44/43 bits
+            ImplMulw(u, t0, t1, H, 4);              // H(1)       44/43 bits
         
             ulong t2 = (f1 << 1) ^ (f2 << 2);
             ulong t3 = (g1 << 1) ^ (g2 << 2);
 
-            ImplMulw(f0 ^ t2, g0 ^ t3, H, 6);     // H(t)       44/45 bits
-            ImplMulw(t0 ^ t2, t1 ^ t3, H, 8);     // H(t + 1)   44/45 bits
+            ImplMulw(u, f0 ^ t2, g0 ^ t3, H, 6);    // H(t)       44/45 bits
+            ImplMulw(u, t0 ^ t2, t1 ^ t3, H, 8);    // H(t + 1)   44/45 bits
 
             ulong t4 = H[6] ^ H[8];
             ulong t5 = H[7] ^ H[9];
@@ -301,12 +302,11 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplCompactExt(zz);
         }
 
-        protected static void ImplMulw(ulong x, ulong y, ulong[] z, int zOff)
+        protected static void ImplMulw(ulong[] u, ulong x, ulong y, ulong[] z, int zOff)
         {
             Debug.Assert(x >> 45 == 0);
             Debug.Assert(y >> 45 == 0);
 
-            ulong[] u = new ulong[8];
             //u[0] = 0;
             u[1] = y;
             u[2] = u[1] << 1;
@@ -318,20 +318,23 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
             uint j = (uint)x;
             ulong g, h = 0, l = u[j & 7]
-                              ^ u[(j >> 3) & 7] << 3
-                              ^ u[(j >> 6) & 7] << 6;
-            int k = 33;
+                              ^ u[(j >>  3) & 7] <<  3
+                              ^ u[(j >>  6) & 7] <<  6
+                              ^ u[(j >>  9) & 7] <<  9
+                              ^ u[(j >> 12) & 7] << 12;
+            int k = 30;
             do
             {
                 j  = (uint)(x >> k);
                 g  = u[j & 7]
-                   ^ u[(j >> 3) & 7] << 3
-                   ^ u[(j >> 6) & 7] << 6
-                   ^ u[(j >> 9) & 7] << 9;
-                l ^= (g <<  k);
+                   ^ u[(j >>  3) & 7] <<  3
+                   ^ u[(j >>  6) & 7] <<  6
+                   ^ u[(j >>  9) & 7] <<  9
+                   ^ u[(j >> 12) & 7] << 12;
+                l ^= (g << k);
                 h ^= (g >> -k);
             }
-            while ((k -= 12) > 0);
+            while ((k -= 15) > 0);
 
             Debug.Assert(h >> 25 == 0);
 
@@ -341,8 +344,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         protected static void ImplSquare(ulong[] x, ulong[] zz)
         {
-            Interleave.Expand64To128(x[0], zz, 0);
-            Interleave.Expand64To128(x[1], zz, 2);
+            Interleave.Expand64To128(x, 0, 2, zz, 0);
             zz[4] = Interleave.Expand8to16((uint)x[2]);
         }
     }
diff --git a/crypto/src/math/ec/custom/sec/SecT163Field.cs b/crypto/src/math/ec/custom/sec/SecT163Field.cs
index b7f60d860..79079ac0b 100644
--- a/crypto/src/math/ec/custom/sec/SecT163Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT163Field.cs
@@ -106,14 +106,14 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Multiply(ulong[] x, ulong[] y, ulong[] z)
         {
-            ulong[] tt = Nat192.CreateExt64();
+            ulong[] tt = new ulong[8];
             ImplMultiply(x, y, tt);
             Reduce(tt, z);
         }
 
         public static void MultiplyAddToExt(ulong[] x, ulong[] y, ulong[] zz)
         {
-            ulong[] tt = Nat192.CreateExt64();
+            ulong[] tt = new ulong[8];
             ImplMultiply(x, y, tt);
             AddExt(zz, tt, zz);
         }
@@ -225,21 +225,22 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             g1  = ((g0 >> 55) ^ (g1 <<  9)) & M55;
             g0 &= M55;
 
+            ulong[] u = zz;
             ulong[] H = new ulong[10];
 
-            ImplMulw(f0, g0, H, 0);               // H(0)       55/54 bits
-            ImplMulw(f2, g2, H, 2);               // H(INF)     55/50 bits
+            ImplMulw(u, f0, g0, H, 0);              // H(0)       55/54 bits
+            ImplMulw(u, f2, g2, H, 2);              // H(INF)     55/50 bits
 
             ulong t0 = f0 ^ f1 ^ f2;
             ulong t1 = g0 ^ g1 ^ g2;
 
-            ImplMulw(t0, t1, H, 4);               // H(1)       55/54 bits
+            ImplMulw(u, t0, t1, H, 4);              // H(1)       55/54 bits
         
             ulong t2 = (f1 << 1) ^ (f2 << 2);
             ulong t3 = (g1 << 1) ^ (g2 << 2);
 
-            ImplMulw(f0 ^ t2, g0 ^ t3, H, 6);     // H(t)       55/56 bits
-            ImplMulw(t0 ^ t2, t1 ^ t3, H, 8);     // H(t + 1)   55/56 bits
+            ImplMulw(u, f0 ^ t2, g0 ^ t3, H, 6);    // H(t)       55/56 bits
+            ImplMulw(u, t0 ^ t2, t1 ^ t3, H, 8);    // H(t + 1)   55/56 bits
 
             ulong t4 = H[6] ^ H[8];
             ulong t5 = H[7] ^ H[9];
@@ -312,12 +313,11 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplCompactExt(zz);
         }
 
-        protected static void ImplMulw(ulong x, ulong y, ulong[] z, int zOff)
+        protected static void ImplMulw(ulong[] u, ulong x, ulong y, ulong[] z, int zOff)
         {
             Debug.Assert(x >> 56 == 0);
             Debug.Assert(y >> 56 == 0);
 
-            ulong[] u = new ulong[8];
             //u[0] = 0;
             u[1] = y;
             u[2] = u[1] << 1;
@@ -349,9 +349,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         protected static void ImplSquare(ulong[] x, ulong[] zz)
         {
-            Interleave.Expand64To128(x[0], zz, 0);
-            Interleave.Expand64To128(x[1], zz, 2);
-            Interleave.Expand64To128(x[2], zz, 4);
+            Interleave.Expand64To128(x, 0, 3, zz, 0);
         }
     }
 }
diff --git a/crypto/src/math/ec/custom/sec/SecT193Field.cs b/crypto/src/math/ec/custom/sec/SecT193Field.cs
index 3ad9b0af2..1a4739b69 100644
--- a/crypto/src/math/ec/custom/sec/SecT193Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT193Field.cs
@@ -239,10 +239,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplExpand(x, f);
             ImplExpand(y, g);
 
-            ImplMulwAcc(f[0], g[0], zz, 0);
-            ImplMulwAcc(f[1], g[1], zz, 1);
-            ImplMulwAcc(f[2], g[2], zz, 2);
-            ImplMulwAcc(f[3], g[3], zz, 3);
+            ulong[] u = new ulong[8];
+
+            ImplMulwAcc(u, f[0], g[0], zz, 0);
+            ImplMulwAcc(u, f[1], g[1], zz, 1);
+            ImplMulwAcc(u, f[2], g[2], zz, 2);
+            ImplMulwAcc(u, f[3], g[3], zz, 3);
 
             // U *= (1 - t^n)
             for (int i = 5; i > 0; --i)
@@ -250,8 +252,8 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
                 zz[i] ^= zz[i - 1];
             }
 
-            ImplMulwAcc(f[0] ^ f[1], g[0] ^ g[1], zz, 1);
-            ImplMulwAcc(f[2] ^ f[3], g[2] ^ g[3], zz, 3);
+            ImplMulwAcc(u, f[0] ^ f[1], g[0] ^ g[1], zz, 1);
+            ImplMulwAcc(u, f[2] ^ f[3], g[2] ^ g[3], zz, 3);
 
             // V *= (1 - t^2n)
             for (int i = 7; i > 1; --i)
@@ -263,10 +265,10 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             {
                 ulong c0 = f[0] ^ f[2], c1 = f[1] ^ f[3];
                 ulong d0 = g[0] ^ g[2], d1 = g[1] ^ g[3];
-                ImplMulwAcc(c0 ^ c1, d0 ^ d1, zz, 3);
+                ImplMulwAcc(u, c0 ^ c1, d0 ^ d1, zz, 3);
                 ulong[] t = new ulong[3];
-                ImplMulwAcc(c0, d0, t, 0);
-                ImplMulwAcc(c1, d1, t, 1);
+                ImplMulwAcc(u, c0, d0, t, 0);
+                ImplMulwAcc(u, c1, d1, t, 1);
                 ulong t0 = t[0], t1 = t[1], t2 = t[2];
                 zz[2] ^= t0;
                 zz[3] ^= t0 ^ t1;
@@ -277,12 +279,11 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplCompactExt(zz);
         }
 
-        protected static void ImplMulwAcc(ulong x, ulong y, ulong[] z, int zOff)
+        protected static void ImplMulwAcc(ulong[] u, ulong x, ulong y, ulong[] z, int zOff)
         {
             Debug.Assert(x >> 49 == 0);
             Debug.Assert(y >> 49 == 0);
 
-            ulong[] u = new ulong[8];
             //u[0] = 0;
             u[1] = y;
             u[2] = u[1] << 1;
@@ -317,9 +318,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         protected static void ImplSquare(ulong[] x, ulong[] zz)
         {
-            Interleave.Expand64To128(x[0], zz, 0);
-            Interleave.Expand64To128(x[1], zz, 2);
-            Interleave.Expand64To128(x[2], zz, 4);
+            Interleave.Expand64To128(x, 0, 3, zz, 0);
             zz[6] = (x[3] & M01);
         }
     }
diff --git a/crypto/src/math/ec/custom/sec/SecT233Field.cs b/crypto/src/math/ec/custom/sec/SecT233Field.cs
index d7916c57d..1ebac2eac 100644
--- a/crypto/src/math/ec/custom/sec/SecT233Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT233Field.cs
@@ -251,10 +251,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplExpand(x, f);
             ImplExpand(y, g);
 
-            ImplMulwAcc(f[0], g[0], zz, 0);
-            ImplMulwAcc(f[1], g[1], zz, 1);
-            ImplMulwAcc(f[2], g[2], zz, 2);
-            ImplMulwAcc(f[3], g[3], zz, 3);
+            ulong[] u = new ulong[8];
+
+            ImplMulwAcc(u, f[0], g[0], zz, 0);
+            ImplMulwAcc(u, f[1], g[1], zz, 1);
+            ImplMulwAcc(u, f[2], g[2], zz, 2);
+            ImplMulwAcc(u, f[3], g[3], zz, 3);
 
             // U *= (1 - t^n)
             for (int i = 5; i > 0; --i)
@@ -262,8 +264,8 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
                 zz[i] ^= zz[i - 1];
             }
 
-            ImplMulwAcc(f[0] ^ f[1], g[0] ^ g[1], zz, 1);
-            ImplMulwAcc(f[2] ^ f[3], g[2] ^ g[3], zz, 3);
+            ImplMulwAcc(u, f[0] ^ f[1], g[0] ^ g[1], zz, 1);
+            ImplMulwAcc(u, f[2] ^ f[3], g[2] ^ g[3], zz, 3);
 
             // V *= (1 - t^2n)
             for (int i = 7; i > 1; --i)
@@ -275,10 +277,10 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             {
                 ulong c0 = f[0] ^ f[2], c1 = f[1] ^ f[3];
                 ulong d0 = g[0] ^ g[2], d1 = g[1] ^ g[3];
-                ImplMulwAcc(c0 ^ c1, d0 ^ d1, zz, 3);
+                ImplMulwAcc(u, c0 ^ c1, d0 ^ d1, zz, 3);
                 ulong[] t = new ulong[3];
-                ImplMulwAcc(c0, d0, t, 0);
-                ImplMulwAcc(c1, d1, t, 1);
+                ImplMulwAcc(u, c0, d0, t, 0);
+                ImplMulwAcc(u, c1, d1, t, 1);
                 ulong t0 = t[0], t1 = t[1], t2 = t[2];
                 zz[2] ^= t0;
                 zz[3] ^= t0 ^ t1;
@@ -289,12 +291,11 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplCompactExt(zz);
         }
 
-        protected static void ImplMulwAcc(ulong x, ulong y, ulong[] z, int zOff)
+        protected static void ImplMulwAcc(ulong[] u, ulong x, ulong y, ulong[] z, int zOff)
         {
             Debug.Assert(x >> 59 == 0);
             Debug.Assert(y >> 59 == 0);
 
-            ulong[] u = new ulong[8];
             //u[0] = 0;
             u[1] = y;
             u[2] = u[1] << 1;
@@ -326,10 +327,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         protected static void ImplSquare(ulong[] x, ulong[] zz)
         {
-            Interleave.Expand64To128(x[0], zz, 0);
-            Interleave.Expand64To128(x[1], zz, 2);
-            Interleave.Expand64To128(x[2], zz, 4);
-            Interleave.Expand64To128(x[3], zz, 6);
+            Interleave.Expand64To128(x, 0, 4, zz, 0);
         }
     }
 }
diff --git a/crypto/src/math/ec/custom/sec/SecT239Field.cs b/crypto/src/math/ec/custom/sec/SecT239Field.cs
index eab929359..ce2e3ba84 100644
--- a/crypto/src/math/ec/custom/sec/SecT239Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT239Field.cs
@@ -260,10 +260,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplExpand(x, f);
             ImplExpand(y, g);
 
-            ImplMulwAcc(f[0], g[0], zz, 0);
-            ImplMulwAcc(f[1], g[1], zz, 1);
-            ImplMulwAcc(f[2], g[2], zz, 2);
-            ImplMulwAcc(f[3], g[3], zz, 3);
+            ulong[] u = new ulong[8];
+
+            ImplMulwAcc(u, f[0], g[0], zz, 0);
+            ImplMulwAcc(u, f[1], g[1], zz, 1);
+            ImplMulwAcc(u, f[2], g[2], zz, 2);
+            ImplMulwAcc(u, f[3], g[3], zz, 3);
 
             // U *= (1 - t^n)
             for (int i = 5; i > 0; --i)
@@ -271,8 +273,8 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
                 zz[i] ^= zz[i - 1];
             }
 
-            ImplMulwAcc(f[0] ^ f[1], g[0] ^ g[1], zz, 1);
-            ImplMulwAcc(f[2] ^ f[3], g[2] ^ g[3], zz, 3);
+            ImplMulwAcc(u, f[0] ^ f[1], g[0] ^ g[1], zz, 1);
+            ImplMulwAcc(u, f[2] ^ f[3], g[2] ^ g[3], zz, 3);
 
             // V *= (1 - t^2n)
             for (int i = 7; i > 1; --i)
@@ -284,10 +286,10 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             {
                 ulong c0 = f[0] ^ f[2], c1 = f[1] ^ f[3];
                 ulong d0 = g[0] ^ g[2], d1 = g[1] ^ g[3];
-                ImplMulwAcc(c0 ^ c1, d0 ^ d1, zz, 3);
+                ImplMulwAcc(u, c0 ^ c1, d0 ^ d1, zz, 3);
                 ulong[] t = new ulong[3];
-                ImplMulwAcc(c0, d0, t, 0);
-                ImplMulwAcc(c1, d1, t, 1);
+                ImplMulwAcc(u, c0, d0, t, 0);
+                ImplMulwAcc(u, c1, d1, t, 1);
                 ulong t0 = t[0], t1 = t[1], t2 = t[2];
                 zz[2] ^= t0;
                 zz[3] ^= t0 ^ t1;
@@ -298,12 +300,11 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplCompactExt(zz);
         }
 
-        protected static void ImplMulwAcc(ulong x, ulong y, ulong[] z, int zOff)
+        protected static void ImplMulwAcc(ulong[] u, ulong x, ulong y, ulong[] z, int zOff)
         {
             Debug.Assert(x >> 60 == 0);
             Debug.Assert(y >> 60 == 0);
 
-            ulong[] u = new ulong[8];
             //u[0] = 0;
             u[1] = y;
             u[2] = u[1] << 1;
@@ -337,10 +338,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         protected static void ImplSquare(ulong[] x, ulong[] zz)
         {
-            Interleave.Expand64To128(x[0], zz, 0);
-            Interleave.Expand64To128(x[1], zz, 2);
-            Interleave.Expand64To128(x[2], zz, 4);
-            Interleave.Expand64To128(x[3], zz, 6);
+            Interleave.Expand64To128(x, 0, 4, zz, 0);
         }
     }
 }
diff --git a/crypto/src/math/ec/custom/sec/SecT283Field.cs b/crypto/src/math/ec/custom/sec/SecT283Field.cs
index 4e2cee0f8..61a1c9afd 100644
--- a/crypto/src/math/ec/custom/sec/SecT283Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT283Field.cs
@@ -10,7 +10,8 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
         private const ulong M27 = ulong.MaxValue >> 37;
         private const ulong M57 = ulong.MaxValue >> 7;
 
-        private static readonly ulong[] ROOT_Z = new ulong[]{ 0x0C30C30C30C30808UL, 0x30C30C30C30C30C3UL, 0x820820820820830CUL, 0x0820820820820820UL, 0x2082082UL };
+        private static readonly ulong[] ROOT_Z = new ulong[]{ 0x0C30C30C30C30808UL, 0x30C30C30C30C30C3UL,
+            0x820820820820830CUL, 0x0820820820820820UL, 0x2082082UL };
 
         public static void Add(ulong[] x, ulong[] y, ulong[] z)
         {
@@ -263,32 +264,33 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplExpand(x, a);
             ImplExpand(y, b);
 
+            ulong[] u = zz;
             ulong[] p = new ulong[26];
 
-            ImplMulw(a[0], b[0], p, 0);                                 // m1
-            ImplMulw(a[1], b[1], p, 2);                                 // m2
-            ImplMulw(a[2], b[2], p, 4);                                 // m3
-            ImplMulw(a[3], b[3], p, 6);                                 // m4
-            ImplMulw(a[4], b[4], p, 8);                                 // m5
+            ImplMulw(u, a[0], b[0], p, 0);                  // m1
+            ImplMulw(u, a[1], b[1], p, 2);                  // m2
+            ImplMulw(u, a[2], b[2], p, 4);                  // m3
+            ImplMulw(u, a[3], b[3], p, 6);                  // m4
+            ImplMulw(u, a[4], b[4], p, 8);                  // m5
 
             ulong u0 = a[0] ^ a[1], v0 = b[0] ^ b[1];
             ulong u1 = a[0] ^ a[2], v1 = b[0] ^ b[2];
             ulong u2 = a[2] ^ a[4], v2 = b[2] ^ b[4];
             ulong u3 = a[3] ^ a[4], v3 = b[3] ^ b[4];
 
-            ImplMulw(u1 ^ a[3], v1 ^ b[3], p, 18);                      // m10
-            ImplMulw(u2 ^ a[1], v2 ^ b[1], p, 20);                      // m11
+            ImplMulw(u, u1 ^ a[3], v1 ^ b[3], p, 18);       // m10
+            ImplMulw(u, u2 ^ a[1], v2 ^ b[1], p, 20);       // m11
 
             ulong A4 = u0 ^ u3  , B4 = v0 ^ v3;
             ulong A5 = A4 ^ a[2], B5 = B4 ^ b[2];
 
-            ImplMulw(A4, B4, p, 22);                                    // m12
-            ImplMulw(A5, B5, p, 24);                                    // m13
+            ImplMulw(u, A4, B4, p, 22);                     // m12
+            ImplMulw(u, A5, B5, p, 24);                     // m13
 
-            ImplMulw(u0, v0, p, 10);                                    // m6
-            ImplMulw(u1, v1, p, 12);                                    // m7
-            ImplMulw(u2, v2, p, 14);                                    // m8
-            ImplMulw(u3, v3, p, 16);                                    // m9
+            ImplMulw(u, u0, v0, p, 10);                     // m6
+            ImplMulw(u, u1, v1, p, 12);                     // m7
+            ImplMulw(u, u2, v2, p, 14);                     // m8
+            ImplMulw(u, u3, v3, p, 16);                     // m9
 
 
             // Original method, corresponding to formula (16)
@@ -375,12 +377,11 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplCompactExt(zz);
         }
 
-        protected static void ImplMulw(ulong x, ulong y, ulong[] z, int zOff)
+        protected static void ImplMulw(ulong[] u, ulong x, ulong y, ulong[] z, int zOff)
         {
             Debug.Assert(x >> 57 == 0);
             Debug.Assert(y >> 57 == 0);
 
-            ulong[] u = new ulong[8];
             //u[0] = 0;
             u[1] = y;
             u[2] = u[1] << 1;
@@ -414,10 +415,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         protected static void ImplSquare(ulong[] x, ulong[] zz)
         {
-            Interleave.Expand64To128(x[0], zz, 0);
-            Interleave.Expand64To128(x[1], zz, 2);
-            Interleave.Expand64To128(x[2], zz, 4);
-            Interleave.Expand64To128(x[3], zz, 6);
+            Interleave.Expand64To128(x, 0, 4, zz, 0);
             zz[8] = Interleave.Expand32to64((uint)x[4]);
         }
     }
diff --git a/crypto/src/math/ec/custom/sec/SecT409Field.cs b/crypto/src/math/ec/custom/sec/SecT409Field.cs
index 2e5609542..c35d3cef0 100644
--- a/crypto/src/math/ec/custom/sec/SecT409Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT409Field.cs
@@ -271,9 +271,8 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             zz[10] = (z10 >> 50) ^ (z11 <<  9);
             zz[11] = (z11 >> 55) ^ (z12 <<  4)
                                  ^ (z13 << 63);
-            zz[12] = (z12 >> 60)
-                   ^ (z13 >> 1);
-            zz[13] = 0;
+            zz[12] = (z13 >>  1);
+            //zz[13] = 0;
         }
 
         protected static void ImplExpand(ulong[] x, ulong[] z)
@@ -294,19 +293,69 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             ImplExpand(x, a);
             ImplExpand(y, b);
 
+            ulong[] u = new ulong[8];
             for (int i = 0; i < 7; ++i)
             {
-                ImplMulwAcc(a, b[i], zz, i);
+                ImplMulwAcc(u, a[i], b[i], zz, i << 1);
             }
 
+            ulong v0 = zz[0], v1 = zz[1];
+            v0 ^= zz[ 2]; zz[1] = v0 ^ v1; v1 ^= zz[ 3];
+            v0 ^= zz[ 4]; zz[2] = v0 ^ v1; v1 ^= zz[ 5];
+            v0 ^= zz[ 6]; zz[3] = v0 ^ v1; v1 ^= zz[ 7];
+            v0 ^= zz[ 8]; zz[4] = v0 ^ v1; v1 ^= zz[ 9];
+            v0 ^= zz[10]; zz[5] = v0 ^ v1; v1 ^= zz[11];
+            v0 ^= zz[12]; zz[6] = v0 ^ v1; v1 ^= zz[13];
+
+            ulong w = v0 ^ v1;
+            zz[ 7] = zz[0] ^ w;
+            zz[ 8] = zz[1] ^ w;
+            zz[ 9] = zz[2] ^ w;
+            zz[10] = zz[3] ^ w;
+            zz[11] = zz[4] ^ w;
+            zz[12] = zz[5] ^ w;
+            zz[13] = zz[6] ^ w;
+
+            ImplMulwAcc(u, a[0] ^ a[1], b[0] ^ b[1], zz,  1);
+
+            ImplMulwAcc(u, a[0] ^ a[2], b[0] ^ b[2], zz,  2);
+
+            ImplMulwAcc(u, a[0] ^ a[3], b[0] ^ b[3], zz,  3);
+            ImplMulwAcc(u, a[1] ^ a[2], b[1] ^ b[2], zz,  3);
+
+            ImplMulwAcc(u, a[0] ^ a[4], b[0] ^ b[4], zz,  4);
+            ImplMulwAcc(u, a[1] ^ a[3], b[1] ^ b[3], zz,  4);
+
+            ImplMulwAcc(u, a[0] ^ a[5], b[0] ^ b[5], zz,  5);
+            ImplMulwAcc(u, a[1] ^ a[4], b[1] ^ b[4], zz,  5);
+            ImplMulwAcc(u, a[2] ^ a[3], b[2] ^ b[3], zz,  5);
+
+            ImplMulwAcc(u, a[0] ^ a[6], b[0] ^ b[6], zz,  6);
+            ImplMulwAcc(u, a[1] ^ a[5], b[1] ^ b[5], zz,  6);
+            ImplMulwAcc(u, a[2] ^ a[4], b[2] ^ b[4], zz,  6);
+
+            ImplMulwAcc(u, a[1] ^ a[6], b[1] ^ b[6], zz,  7);
+            ImplMulwAcc(u, a[2] ^ a[5], b[2] ^ b[5], zz,  7);
+            ImplMulwAcc(u, a[3] ^ a[4], b[3] ^ b[4], zz,  7);
+
+            ImplMulwAcc(u, a[2] ^ a[6], b[2] ^ b[6], zz,  8);
+            ImplMulwAcc(u, a[3] ^ a[5], b[3] ^ b[5], zz,  8);
+
+            ImplMulwAcc(u, a[3] ^ a[6], b[3] ^ b[6], zz,  9);
+            ImplMulwAcc(u, a[4] ^ a[5], b[4] ^ b[5], zz,  9);
+
+            ImplMulwAcc(u, a[4] ^ a[6], b[4] ^ b[6], zz, 10);
+
+            ImplMulwAcc(u, a[5] ^ a[6], b[5] ^ b[6], zz, 11);
+
             ImplCompactExt(zz);
         }
 
-        protected static void ImplMulwAcc(ulong[] xs, ulong y, ulong[] z, int zOff)
+        protected static void ImplMulwAcc(ulong[] u, ulong x, ulong y, ulong[] z, int zOff)
         {
+            Debug.Assert(x >> 59 == 0);
             Debug.Assert(y >> 59 == 0);
 
-            ulong[] u = new ulong[8];
             //u[0] = 0;
             u[1] = y;
             u[2] = u[1] << 1;
@@ -316,41 +365,29 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             u[6] = u[3] << 1;
             u[7] = u[6] ^  y;
 
-            for (int i = 0; i < 7; ++i)
+            uint j = (uint)x;
+            ulong g, h = 0, l = u[j & 7]
+                              ^ (u[(j >> 3) & 7] << 3);
+            int k = 54;
+            do
             {
-                ulong x = xs[i];
-
-                Debug.Assert(x >> 59 == 0);
-
-                uint j = (uint)x;
-                ulong g, h = 0, l = u[j & 7]
-                                  ^ (u[(j >> 3) & 7] << 3);
-                int k = 54;
-                do
-                {
-                    j  = (uint)(x >> k);
-                    g  = u[j & 7]
-                       ^ u[(j >> 3) & 7] << 3;
-                    l ^= (g <<  k);
-                    h ^= (g >> -k);
-                }
-                while ((k -= 6) > 0);
-
-                Debug.Assert(h >> 53 == 0);
-
-                z[zOff + i    ] ^= l & M59;
-                z[zOff + i + 1] ^= (l >> 59) ^ (h << 5);
+                j  = (uint)(x >> k);
+                g  = u[j & 7]
+                   ^ u[(j >> 3) & 7] << 3;
+                l ^= (g << k);
+                h ^= (g >> -k);
             }
+            while ((k -= 6) > 0);
+
+            Debug.Assert(h >> 53 == 0);
+
+            z[zOff    ] ^= l & M59;
+            z[zOff + 1] ^= (l >> 59) ^ (h << 5);
         }
 
         protected static void ImplSquare(ulong[] x, ulong[] zz)
         {
-            Interleave.Expand64To128(x[0], zz, 0);
-            Interleave.Expand64To128(x[1], zz, 2);
-            Interleave.Expand64To128(x[2], zz, 4);
-            Interleave.Expand64To128(x[3], zz, 6);
-            Interleave.Expand64To128(x[4], zz, 8);
-            Interleave.Expand64To128(x[5], zz, 10);
+            Interleave.Expand64To128(x, 0, 6, zz, 0);
             zz[12] = Interleave.Expand32to64((uint)x[6]);
         }
     }
diff --git a/crypto/src/math/ec/custom/sec/SecT571Field.cs b/crypto/src/math/ec/custom/sec/SecT571Field.cs
index 0d9b337fc..1b8bb763e 100644
--- a/crypto/src/math/ec/custom/sec/SecT571Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT571Field.cs
@@ -9,10 +9,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
     {
         private const ulong M59 = ulong.MaxValue >> 5;
 
-        private const ulong RM = 0xEF7BDEF7BDEF7BDEUL;
-
-        private static readonly ulong[] ROOT_Z = new ulong[]{ 0x2BE1195F08CAFB99UL, 0x95F08CAF84657C23UL, 0xCAF84657C232BE11UL, 0x657C232BE1195F08UL,
-            0xF84657C2308CAF84UL, 0x7C232BE1195F08CAUL, 0xBE1195F08CAF8465UL, 0x5F08CAF84657C232UL, 0x784657C232BE119UL };
+        private static readonly ulong[] ROOT_Z = new ulong[]{ 0x2BE1195F08CAFB99UL, 0x95F08CAF84657C23UL,
+            0xCAF84657C232BE11UL, 0x657C232BE1195F08UL, 0xF84657C2308CAF84UL, 0x7C232BE1195F08CAUL,
+            0xBE1195F08CAF8465UL, 0x5F08CAF84657C232UL, 0x784657C232BE119UL };
 
         public static void Add(ulong[] x, ulong[] y, ulong[] z)
         {
@@ -30,6 +29,14 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             }
         }
 
+        public static void AddBothTo(ulong[] x, ulong[] y, ulong[] z)
+        {
+            for (int i = 0; i < 9; ++i)
+            {
+                z[i] ^= x[i] ^ y[i];
+            }
+        }
+
         private static void AddBothTo(ulong[] x, int xOff, ulong[] y, int yOff, ulong[] z, int zOff)
         {
             for (int i = 0; i < 9; ++i)
@@ -148,6 +155,46 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             AddExt(zz, tt, zz);
         }
 
+        public static void MultiplyPrecomp(ulong[] x, ulong[] precomp, ulong[] z)
+        {
+            ulong[] tt = Nat576.CreateExt64();
+            ImplMultiplyPrecomp(x, precomp, tt);
+            Reduce(tt, z);
+        }
+
+        public static void MultiplyPrecompAddToExt(ulong[] x, ulong[] precomp, ulong[] zz)
+        {
+            ulong[] tt = Nat576.CreateExt64();
+            ImplMultiplyPrecomp(x, precomp, tt);
+            AddExt(zz, tt, zz);
+        }
+
+        public static ulong[] PrecompMultiplicand(ulong[] x)
+        {
+            /*
+             * Precompute table of all 4-bit products of x (first section)
+             */
+            int len = 9 << 4;
+            ulong[] t = new ulong[len << 1];
+            Array.Copy(x, 0, t, 9, 9);
+            //Reduce5(t, 9);
+            int tOff = 0;
+            for (int i = 7; i > 0; --i)
+            {
+                tOff += 18;
+                Nat.ShiftUpBit64(9, t, tOff >> 1, 0UL, t, tOff);
+                Reduce5(t, tOff);
+                Add(t, 9, t, tOff, t, tOff + 9);
+            }
+
+            /*
+             * Second section with all 4-bit products of x shifted 4 bits
+             */
+            Nat.ShiftUpBits64(len, t, 0, 4, 0UL, t, len);
+
+            return t;
+        }
+
         public static void Reduce(ulong[] xx, ulong[] z)
         {
             ulong xx09 = xx[9];
@@ -239,32 +286,91 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         protected static void ImplMultiply(ulong[] x, ulong[] y, ulong[] zz)
         {
-            //for (int i = 0; i < 9; ++i)
-            //{
-            //    ImplMulwAcc(x, y[i], zz, i);
-            //}
+            //ulong[] precomp = PrecompMultiplicand(y);
 
-            /*
-             * Precompute table of all 4-bit products of y
-             */
-            ulong[] T0 = new ulong[9 << 4];
-            Array.Copy(y, 0, T0, 9, 9);
-    //        Reduce5(T0, 9);
-            int tOff = 0;
-            for (int i = 7; i > 0; --i)
+            //ImplMultiplyPrecomp(x, precomp, zz);
+
+            ulong[] u = new ulong[16];
+            for (int i = 0; i < 9; ++i)
             {
-                tOff += 18;
-                Nat.ShiftUpBit64(9, T0, tOff >> 1, 0UL, T0, tOff);
-                Reduce5(T0, tOff);
-                Add(T0, 9, T0, tOff, T0, tOff + 9);
+                ImplMulwAcc(u, x[i], y[i], zz, i << 1);
             }
 
-            /*
-             * Second table with all 4-bit products of B shifted 4 bits
-             */
-            ulong[] T1 = new ulong[T0.Length];
-            Nat.ShiftUpBits64(T0.Length, T0, 0, 4, 0L, T1, 0);
+            ulong v0 = zz[0], v1 = zz[1];
+            v0 ^= zz[ 2]; zz[1] = v0 ^ v1; v1 ^= zz[ 3];
+            v0 ^= zz[ 4]; zz[2] = v0 ^ v1; v1 ^= zz[ 5];
+            v0 ^= zz[ 6]; zz[3] = v0 ^ v1; v1 ^= zz[ 7];
+            v0 ^= zz[ 8]; zz[4] = v0 ^ v1; v1 ^= zz[ 9];
+            v0 ^= zz[10]; zz[5] = v0 ^ v1; v1 ^= zz[11];
+            v0 ^= zz[12]; zz[6] = v0 ^ v1; v1 ^= zz[13];
+            v0 ^= zz[14]; zz[7] = v0 ^ v1; v1 ^= zz[15];
+            v0 ^= zz[16]; zz[8] = v0 ^ v1; v1 ^= zz[17];
+
+            ulong w = v0 ^ v1;
+            zz[ 9] = zz[0] ^ w;
+            zz[10] = zz[1] ^ w;
+            zz[11] = zz[2] ^ w;
+            zz[12] = zz[3] ^ w;
+            zz[13] = zz[4] ^ w;
+            zz[14] = zz[5] ^ w;
+            zz[15] = zz[6] ^ w;
+            zz[16] = zz[7] ^ w;
+            zz[17] = zz[8] ^ w;
+
+            ImplMulwAcc(u, x[0] ^ x[1], y[0] ^ y[1], zz, 1);
+
+            ImplMulwAcc(u, x[0] ^ x[2], y[0] ^ y[2], zz, 2);
+
+            ImplMulwAcc(u, x[0] ^ x[3], y[0] ^ y[3], zz, 3);
+            ImplMulwAcc(u, x[1] ^ x[2], y[1] ^ y[2], zz, 3);
+
+            ImplMulwAcc(u, x[0] ^ x[4], y[0] ^ y[4], zz, 4);
+            ImplMulwAcc(u, x[1] ^ x[3], y[1] ^ y[3], zz, 4);
+
+            ImplMulwAcc(u, x[0] ^ x[5], y[0] ^ y[5], zz, 5);
+            ImplMulwAcc(u, x[1] ^ x[4], y[1] ^ y[4], zz, 5);
+            ImplMulwAcc(u, x[2] ^ x[3], y[2] ^ y[3], zz, 5);
+
+            ImplMulwAcc(u, x[0] ^ x[6], y[0] ^ y[6], zz, 6);
+            ImplMulwAcc(u, x[1] ^ x[5], y[1] ^ y[5], zz, 6);
+            ImplMulwAcc(u, x[2] ^ x[4], y[2] ^ y[4], zz, 6);
+
+            ImplMulwAcc(u, x[0] ^ x[7], y[0] ^ y[7], zz, 7);
+            ImplMulwAcc(u, x[1] ^ x[6], y[1] ^ y[6], zz, 7);
+            ImplMulwAcc(u, x[2] ^ x[5], y[2] ^ y[5], zz, 7);
+            ImplMulwAcc(u, x[3] ^ x[4], y[3] ^ y[4], zz, 7);
+
+            ImplMulwAcc(u, x[0] ^ x[8], y[0] ^ y[8], zz, 8);
+            ImplMulwAcc(u, x[1] ^ x[7], y[1] ^ y[7], zz, 8);
+            ImplMulwAcc(u, x[2] ^ x[6], y[2] ^ y[6], zz, 8);
+            ImplMulwAcc(u, x[3] ^ x[5], y[3] ^ y[5], zz, 8);
+
+            ImplMulwAcc(u, x[1] ^ x[8], y[1] ^ y[8], zz, 9);
+            ImplMulwAcc(u, x[2] ^ x[7], y[2] ^ y[7], zz, 9);
+            ImplMulwAcc(u, x[3] ^ x[6], y[3] ^ y[6], zz, 9);
+            ImplMulwAcc(u, x[4] ^ x[5], y[4] ^ y[5], zz, 9);
+
+            ImplMulwAcc(u, x[2] ^ x[8], y[2] ^ y[8], zz, 10);
+            ImplMulwAcc(u, x[3] ^ x[7], y[3] ^ y[7], zz, 10);
+            ImplMulwAcc(u, x[4] ^ x[6], y[4] ^ y[6], zz, 10);
+
+            ImplMulwAcc(u, x[3] ^ x[8], y[3] ^ y[8], zz, 11);
+            ImplMulwAcc(u, x[4] ^ x[7], y[4] ^ y[7], zz, 11);
+            ImplMulwAcc(u, x[5] ^ x[6], y[5] ^ y[6], zz, 11);
+
+            ImplMulwAcc(u, x[4] ^ x[8], y[4] ^ y[8], zz, 12);
+            ImplMulwAcc(u, x[5] ^ x[7], y[5] ^ y[7], zz, 12);
+
+            ImplMulwAcc(u, x[5] ^ x[8], y[5] ^ y[8], zz, 13);
+            ImplMulwAcc(u, x[6] ^ x[7], y[6] ^ y[7], zz, 13);
+
+            ImplMulwAcc(u, x[6] ^ x[8], y[6] ^ y[8], zz, 14);
+
+            ImplMulwAcc(u, x[7] ^ x[8], y[7] ^ y[8], zz, 15);
+        }
 
+        protected static void ImplMultiplyPrecomp(ulong[] x, ulong[] precomp, ulong[] zz)
+        {
             uint MASK = 0xF;
 
             /*
@@ -278,9 +384,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
                     uint aVal = (uint)(x[j] >> k);
                     uint u = aVal & MASK;
                     uint v = (aVal >> 4) & MASK;
-                    AddBothTo(T0, (int)(9 * u), T1, (int)(9 * v), zz, j - 1);
+                    AddBothTo(precomp, (int)(9 * u), precomp, (int)(9 * (v + 16)), zz, j - 1);
                 }
-                Nat.ShiftUpBits64(16, zz, 0, 8, 0L);
+                Nat.ShiftUpBits64(16, zz, 0, 8, 0UL);
             }
 
             for (int k = 56; k >= 0; k -= 8)
@@ -290,70 +396,54 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
                     uint aVal = (uint)(x[j] >> k);
                     uint u = aVal & MASK;
                     uint v = (aVal >> 4) & MASK;
-                    AddBothTo(T0, (int)(9 * u), T1, (int)(9 * v), zz, j);
+                    AddBothTo(precomp, (int)(9 * u), precomp, (int)(9 * (v + 16)), zz, j);
                 }
                 if (k > 0)
                 {
-                    Nat.ShiftUpBits64(18, zz, 0, 8, 0L);
+                    Nat.ShiftUpBits64(18, zz, 0, 8, 0UL);
                 }
             }
         }
 
-        protected static void ImplMulwAcc(ulong[] xs, ulong y, ulong[] z, int zOff)
+        protected static void ImplMulwAcc(ulong[] u, ulong x, ulong y, ulong[] z, int zOff)
         {
-            ulong[] u = new ulong[32];
             //u[0] = 0;
             u[1] = y;
-            for (int i = 2; i < 32; i += 2)
+            for (int i = 2; i < 16; i += 2)
             {
                 u[i    ] = u[i >> 1] << 1;
                 u[i + 1] = u[i     ] ^  y;
             }
 
-            ulong l = 0;
-            for (int i = 0; i < 9; ++i)
+            uint j = (uint)x;
+            ulong g, h = 0, l = u[j & 15]
+                              ^ u[(j >> 4) & 15] << 4;
+            int k = 56;
+            do
             {
-                ulong x = xs[i];
-
-                uint j = (uint)x;
-
-                l ^= u[j & 31];
-
-                ulong g, h = 0;
-                int k = 60;
-                do
-                {
-                    j  = (uint)(x >> k);
-                    g  = u[j & 31];
-                    l ^= (g <<  k);
-                    h ^= (g >> -k);
-                }
-                while ((k -= 5) > 0);
+                j  = (uint)(x >> k);
+                g  = u[j & 15]
+                   ^ u[(j >> 4) & 15] << 4;
+                l ^= (g << k);
+                h ^= (g >> -k);
+            }
+            while ((k -= 8) > 0);
 
-                for (int p = 0; p < 4; ++p)
-                {
-                    x = (x & RM) >> 1;
-                    h ^= x & (ulong)(((long)y << p) >> 63);
-                }
+            for (int p = 0; p < 7; ++p)
+            {
+                x = (x & 0xFEFEFEFEFEFEFEFEUL) >> 1;
+                h ^= x & (ulong)((long)(y << p) >> 63);
+            }
 
-                z[zOff + i] ^= l;
+            Debug.Assert(h >> 63 == 0);
 
-                l = h;
-            }
-            z[zOff + 9] ^= l;
+            z[zOff    ] ^= l;
+            z[zOff + 1] ^= h;
         }
 
         protected static void ImplSquare(ulong[] x, ulong[] zz)
         {
-            Interleave.Expand64To128(x[0], zz, 0);
-            Interleave.Expand64To128(x[1], zz, 2);
-            Interleave.Expand64To128(x[2], zz, 4);
-            Interleave.Expand64To128(x[3], zz, 6);
-            Interleave.Expand64To128(x[4], zz, 8);
-            Interleave.Expand64To128(x[5], zz, 10);
-            Interleave.Expand64To128(x[6], zz, 12);
-            Interleave.Expand64To128(x[7], zz, 14);
-            Interleave.Expand64To128(x[8], zz, 16);
+            Interleave.Expand64To128(x, 0, 9, zz, 0);
         }
     }
 }
diff --git a/crypto/src/math/ec/custom/sec/SecT571K1Point.cs b/crypto/src/math/ec/custom/sec/SecT571K1Point.cs
index a6a97ee5e..0f132c84c 100644
--- a/crypto/src/math/ec/custom/sec/SecT571K1Point.cs
+++ b/crypto/src/math/ec/custom/sec/SecT571K1Point.cs
@@ -1,5 +1,7 @@
 using System;
 
+using Org.BouncyCastle.Math.Raw;
+
 namespace Org.BouncyCastle.Math.EC.Custom.Sec
 {
     internal class SecT571K1Point
@@ -79,8 +81,8 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
             ECCurve curve = this.Curve;
 
-            ECFieldElement X1 = this.RawXCoord;
-            ECFieldElement X2 = b.RawXCoord;
+            SecT571FieldElement X1 = (SecT571FieldElement)this.RawXCoord;
+            SecT571FieldElement X2 = (SecT571FieldElement)b.RawXCoord;
 
             if (X1.IsZero)
             {
@@ -90,82 +92,118 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
                 return b.Add(this);
             }
 
-            ECFieldElement L1 = this.RawYCoord, Z1 = this.RawZCoords[0];
-            ECFieldElement L2 = b.RawYCoord, Z2 = b.RawZCoords[0];
+            SecT571FieldElement L1 = (SecT571FieldElement)this.RawYCoord, Z1 = (SecT571FieldElement)this.RawZCoords[0];
+            SecT571FieldElement L2 = (SecT571FieldElement)b.RawYCoord, Z2 = (SecT571FieldElement)b.RawZCoords[0];
 
-            bool Z1IsOne = Z1.IsOne;
-            ECFieldElement U2 = X2, S2 = L2;
-            if (!Z1IsOne)
+            ulong[] t1 = Nat576.Create64();
+            ulong[] t2 = Nat576.Create64();
+            ulong[] t3 = Nat576.Create64();
+            ulong[] t4 = Nat576.Create64();
+
+            ulong[] Z1Precomp = Z1.IsOne ? null : SecT571Field.PrecompMultiplicand(Z1.x);
+            ulong[] U2, S2;
+            if (Z1Precomp == null)
             {
-                U2 = U2.Multiply(Z1);
-                S2 = S2.Multiply(Z1);
+                U2 = X2.x;
+                S2 = L2.x;
+            }
+            else
+            {
+                SecT571Field.MultiplyPrecomp(X2.x, Z1Precomp, U2 = t2);
+                SecT571Field.MultiplyPrecomp(L2.x, Z1Precomp, S2 = t4);
             }
 
-            bool Z2IsOne = Z2.IsOne;
-            ECFieldElement U1 = X1, S1 = L1;
-            if (!Z2IsOne)
+            ulong[] Z2Precomp = Z2.IsOne ? null : SecT571Field.PrecompMultiplicand(Z2.x);
+            ulong[] U1, S1;
+            if (Z2Precomp == null)
+            {
+                U1 = X1.x;
+                S1 = L1.x;
+            }
+            else
             {
-                U1 = U1.Multiply(Z2);
-                S1 = S1.Multiply(Z2);
+                SecT571Field.MultiplyPrecomp(X1.x, Z2Precomp, U1 = t1);
+                SecT571Field.MultiplyPrecomp(L1.x, Z2Precomp, S1 = t3);
             }
 
-            ECFieldElement A = S1.Add(S2);
-            ECFieldElement B = U1.Add(U2);
+            ulong[] A = t3;
+            SecT571Field.Add(S1, S2, A);
 
-            if (B.IsZero)
+            ulong[] B = t4;
+            SecT571Field.Add(U1, U2, B);
+
+            if (Nat576.IsZero64(B))
             {
-                if (A.IsZero)
+                if (Nat576.IsZero64(A))
                     return Twice();
 
                 return curve.Infinity;
             }
 
-            ECFieldElement X3, L3, Z3;
+            SecT571FieldElement X3, L3, Z3;
             if (X2.IsZero)
             {
                 // TODO This can probably be optimized quite a bit
                 ECPoint p = this.Normalize();
-                X1 = p.XCoord;
+                X1 = (SecT571FieldElement)p.XCoord;
                 ECFieldElement Y1 = p.YCoord;
 
                 ECFieldElement Y2 = L2;
                 ECFieldElement L = Y1.Add(Y2).Divide(X1);
 
-                X3 = L.Square().Add(L).Add(X1);
+                X3 = (SecT571FieldElement)L.Square().Add(L).Add(X1);
                 if (X3.IsZero)
                 {
                     return new SecT571K1Point(curve, X3, curve.B, IsCompressed);
                 }
 
                 ECFieldElement Y3 = L.Multiply(X1.Add(X3)).Add(X3).Add(Y1);
-                L3 = Y3.Divide(X3).Add(X3);
-                Z3 = curve.FromBigInteger(BigInteger.One);
+                L3 = (SecT571FieldElement)Y3.Divide(X3).Add(X3);
+                Z3 = (SecT571FieldElement)curve.FromBigInteger(BigInteger.One);
             }
             else
             {
-                B = B.Square();
+                SecT571Field.Square(B, B);
+
+                ulong[] APrecomp = SecT571Field.PrecompMultiplicand(A);
+
+                ulong[] AU1 = t1;
+                ulong[] AU2 = t2;
 
-                ECFieldElement AU1 = A.Multiply(U1);
-                ECFieldElement AU2 = A.Multiply(U2);
+                SecT571Field.MultiplyPrecomp(U1, APrecomp, AU1);
+                SecT571Field.MultiplyPrecomp(U2, APrecomp, AU2);
+
+                X3 = new SecT571FieldElement(t1);
+                SecT571Field.Multiply(AU1, AU2, X3.x);
 
-                X3 = AU1.Multiply(AU2);
                 if (X3.IsZero)
                 {
                     return new SecT571K1Point(curve, X3, curve.B, IsCompressed);
                 }
 
-                ECFieldElement ABZ2 = A.Multiply(B);
-                if (!Z2IsOne)
+                Z3 = new SecT571FieldElement(t3);
+                SecT571Field.MultiplyPrecomp(B, APrecomp, Z3.x);
+
+                if (Z2Precomp != null)
                 {
-                    ABZ2 = ABZ2.Multiply(Z2);
+                    SecT571Field.MultiplyPrecomp(Z3.x, Z2Precomp, Z3.x);
                 }
 
-                L3 = AU2.Add(B).SquarePlusProduct(ABZ2, L1.Add(Z1));
+                //L3 = AU2.Add(B).SquarePlusProduct(ABZ2, L1.Add(Z1));
+                ulong[] tt = Nat576.CreateExt64();
+
+                SecT571Field.Add(AU2, B, t4);
+                SecT571Field.SquareAddToExt(t4, tt);
+
+                SecT571Field.Add(L1.x, Z1.x, t4);
+                SecT571Field.MultiplyAddToExt(t4, Z3.x, tt);
+
+                L3 = new SecT571FieldElement(t4);
+                SecT571Field.Reduce(tt, L3.x);
 
-                Z3 = ABZ2;
-                if (!Z1IsOne)
+                if (Z1Precomp != null)
                 {
-                    Z3 = Z3.Multiply(Z1);
+                    SecT571Field.MultiplyPrecomp(Z3.x, Z1Precomp, Z3.x);
                 }
             }
 
diff --git a/crypto/src/math/ec/custom/sec/SecT571R1Point.cs b/crypto/src/math/ec/custom/sec/SecT571R1Point.cs
index 8e4d5d52f..6a82a5ef5 100644
--- a/crypto/src/math/ec/custom/sec/SecT571R1Point.cs
+++ b/crypto/src/math/ec/custom/sec/SecT571R1Point.cs
@@ -1,5 +1,7 @@
 using System;
 
+using Org.BouncyCastle.Math.Raw;
+
 namespace Org.BouncyCastle.Math.EC.Custom.Sec
 {
     internal class SecT571R1Point
@@ -79,8 +81,8 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
             ECCurve curve = this.Curve;
 
-            ECFieldElement X1 = this.RawXCoord;
-            ECFieldElement X2 = b.RawXCoord;
+            SecT571FieldElement X1 = (SecT571FieldElement)this.RawXCoord;
+            SecT571FieldElement X2 = (SecT571FieldElement)b.RawXCoord;
 
             if (X1.IsZero)
             {
@@ -90,82 +92,117 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
                 return b.Add(this);
             }
 
-            ECFieldElement L1 = this.RawYCoord, Z1 = this.RawZCoords[0];
-            ECFieldElement L2 = b.RawYCoord, Z2 = b.RawZCoords[0];
+            SecT571FieldElement L1 = (SecT571FieldElement)this.RawYCoord, Z1 = (SecT571FieldElement)this.RawZCoords[0];
+            SecT571FieldElement L2 = (SecT571FieldElement)b.RawYCoord, Z2 = (SecT571FieldElement)b.RawZCoords[0];
 
-            bool Z1IsOne = Z1.IsOne;
-            ECFieldElement U2 = X2, S2 = L2;
-            if (!Z1IsOne)
+            ulong[] t1 = Nat576.Create64();
+            ulong[] t2 = Nat576.Create64();
+            ulong[] t3 = Nat576.Create64();
+            ulong[] t4 = Nat576.Create64();
+
+            ulong[] Z1Precomp = Z1.IsOne ? null : SecT571Field.PrecompMultiplicand(Z1.x);
+            ulong[] U2, S2;
+            if (Z1Precomp == null)
+            {
+                U2 = X2.x;
+                S2 = L2.x;
+            }
+            else
             {
-                U2 = U2.Multiply(Z1);
-                S2 = S2.Multiply(Z1);
+                SecT571Field.MultiplyPrecomp(X2.x, Z1Precomp, U2 = t2);
+                SecT571Field.MultiplyPrecomp(L2.x, Z1Precomp, S2 = t4);
             }
 
-            bool Z2IsOne = Z2.IsOne;
-            ECFieldElement U1 = X1, S1 = L1;
-            if (!Z2IsOne)
+            ulong[] Z2Precomp = Z2.IsOne ? null : SecT571Field.PrecompMultiplicand(Z2.x);
+            ulong[] U1, S1;
+            if (Z2Precomp == null)
+            {
+                U1 = X1.x;
+                S1 = L1.x;
+            }
+            else
             {
-                U1 = U1.Multiply(Z2);
-                S1 = S1.Multiply(Z2);
+                SecT571Field.MultiplyPrecomp(X1.x, Z2Precomp, U1 = t1);
+                SecT571Field.MultiplyPrecomp(L1.x, Z2Precomp, S1 = t3);
             }
 
-            ECFieldElement A = S1.Add(S2);
-            ECFieldElement B = U1.Add(U2);
+            ulong[] A = t3;
+            SecT571Field.Add(S1, S2, A);
 
-            if (B.IsZero)
+            ulong[] B = t4;
+            SecT571Field.Add(U1, U2, B);
+
+            if (Nat576.IsZero64(B))
             {
-                if (A.IsZero)
+                if (Nat576.IsZero64(A))
                     return Twice();
 
                 return curve.Infinity;
             }
 
-            ECFieldElement X3, L3, Z3;
+            SecT571FieldElement X3, L3, Z3;
             if (X2.IsZero)
             {
                 // TODO This can probably be optimized quite a bit
                 ECPoint p = this.Normalize();
-                X1 = p.XCoord;
+                X1 = (SecT571FieldElement)p.XCoord;
                 ECFieldElement Y1 = p.YCoord;
 
                 ECFieldElement Y2 = L2;
                 ECFieldElement L = Y1.Add(Y2).Divide(X1);
 
-                X3 = L.Square().Add(L).Add(X1).AddOne();
+                X3 = (SecT571FieldElement)L.Square().Add(L).Add(X1).AddOne();
                 if (X3.IsZero)
                 {
                     return new SecT571R1Point(curve, X3, SecT571R1Curve.SecT571R1_B_SQRT, IsCompressed);
                 }
 
                 ECFieldElement Y3 = L.Multiply(X1.Add(X3)).Add(X3).Add(Y1);
-                L3 = Y3.Divide(X3).Add(X3);
-                Z3 = curve.FromBigInteger(BigInteger.One);
+                L3 = (SecT571FieldElement)Y3.Divide(X3).Add(X3);
+                Z3 = (SecT571FieldElement)curve.FromBigInteger(BigInteger.One);
             }
             else
             {
-                B = B.Square();
+                SecT571Field.Square(B, B);
+
+                ulong[] APrecomp = SecT571Field.PrecompMultiplicand(A);
+
+                ulong[] AU1 = t1;
+                ulong[] AU2 = t2;
 
-                ECFieldElement AU1 = A.Multiply(U1);
-                ECFieldElement AU2 = A.Multiply(U2);
+                SecT571Field.MultiplyPrecomp(U1, APrecomp, AU1);
+                SecT571Field.MultiplyPrecomp(U2, APrecomp, AU2);
+
+                X3 = new SecT571FieldElement(t1);
+                SecT571Field.Multiply(AU1, AU2, X3.x);
 
-                X3 = AU1.Multiply(AU2);
                 if (X3.IsZero)
                 {
                     return new SecT571R1Point(curve, X3, SecT571R1Curve.SecT571R1_B_SQRT, IsCompressed);
                 }
 
-                ECFieldElement ABZ2 = A.Multiply(B);
-                if (!Z2IsOne)
+                Z3 = new SecT571FieldElement(t3);
+                SecT571Field.MultiplyPrecomp(B, APrecomp, Z3.x);
+
+                if (Z2Precomp != null)
                 {
-                    ABZ2 = ABZ2.Multiply(Z2);
+                    SecT571Field.MultiplyPrecomp(Z3.x, Z2Precomp, Z3.x);
                 }
 
-                L3 = AU2.Add(B).SquarePlusProduct(ABZ2, L1.Add(Z1));
+                ulong[] tt = Nat576.CreateExt64();
+
+                SecT571Field.Add(AU2, B, t4);
+                SecT571Field.SquareAddToExt(t4, tt);
+
+                SecT571Field.Add(L1.x, Z1.x, t4);
+                SecT571Field.MultiplyAddToExt(t4, Z3.x, tt);
+
+                L3 = new SecT571FieldElement(t4);
+                SecT571Field.Reduce(tt, L3.x);
 
-                Z3 = ABZ2;
-                if (!Z1IsOne)
+                if (Z1Precomp != null)
                 {
-                    Z3 = Z3.Multiply(Z1);
+                    SecT571Field.MultiplyPrecomp(Z3.x, Z1Precomp, Z3.x);
                 }
             }
 
@@ -179,29 +216,66 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
             ECCurve curve = this.Curve;
 
-            ECFieldElement X1 = this.RawXCoord;
+            SecT571FieldElement X1 = (SecT571FieldElement)this.RawXCoord;
             if (X1.IsZero)
             {
                 // A point with X == 0 is its own additive inverse
                 return curve.Infinity;
             }
 
-            ECFieldElement L1 = this.RawYCoord, Z1 = this.RawZCoords[0];
+            SecT571FieldElement L1 = (SecT571FieldElement)this.RawYCoord, Z1 = (SecT571FieldElement)this.RawZCoords[0];
+
+            ulong[] t1 = Nat576.Create64();
+            ulong[] t2 = Nat576.Create64();
+
+            ulong[] Z1Precomp = Z1.IsOne ? null : SecT571Field.PrecompMultiplicand(Z1.x);
+            ulong[] L1Z1, Z1Sq;
+            if (Z1Precomp == null)
+            {
+                L1Z1 = L1.x;
+                Z1Sq = Z1.x;
+            }
+            else
+            {
+                SecT571Field.MultiplyPrecomp(L1.x, Z1Precomp, L1Z1 = t1);
+                SecT571Field.Square(Z1.x, Z1Sq = t2);
+            }
+
+            ulong[] T = Nat576.Create64();
+            SecT571Field.Square(L1.x, T);
+            SecT571Field.AddBothTo(L1Z1, Z1Sq, T);
 
-            bool Z1IsOne = Z1.IsOne;
-            ECFieldElement L1Z1 = Z1IsOne ? L1 : L1.Multiply(Z1);
-            ECFieldElement Z1Sq = Z1IsOne ? Z1 : Z1.Square();
-            ECFieldElement T = L1.Square().Add(L1Z1).Add(Z1Sq);
-            if (T.IsZero)
+            if (Nat576.IsZero64(T))
             {
-                return new SecT571R1Point(curve, T, SecT571R1Curve.SecT571R1_B_SQRT, IsCompressed);
+                return new SecT571R1Point(curve, new SecT571FieldElement(T), SecT571R1Curve.SecT571R1_B_SQRT, IsCompressed);
             }
 
-            ECFieldElement X3 = T.Square();
-            ECFieldElement Z3 = Z1IsOne ? T : T.Multiply(Z1Sq);
+            ulong[] tt = Nat576.CreateExt64();
+            SecT571Field.MultiplyAddToExt(T, L1Z1, tt);
+
+            SecT571FieldElement X3 = new SecT571FieldElement(t1);
+            SecT571Field.Square(T, X3.x);
+
+            SecT571FieldElement Z3 = new SecT571FieldElement(T);
+            if (Z1Precomp != null)
+            {
+                SecT571Field.Multiply(Z3.x, Z1Sq, Z3.x);
+            }
+
+            ulong[] X1Z1;
+            if (Z1Precomp == null)
+            {
+                X1Z1 = X1.x;
+            }
+            else
+            {
+                SecT571Field.MultiplyPrecomp(X1.x, Z1Precomp, X1Z1 = t2);
+            }
 
-            ECFieldElement X1Z1 = Z1IsOne ? X1 : X1.Multiply(Z1);
-            ECFieldElement L3 = X1Z1.SquarePlusProduct(T, L1Z1).Add(X3).Add(Z3);
+            SecT571Field.SquareAddToExt(X1Z1, tt);
+            SecT571Field.Reduce(tt, t2);
+            SecT571Field.AddBothTo(X3.x, Z3.x, t2);
+            SecT571FieldElement L3 = new SecT571FieldElement(t2);
 
             return new SecT571R1Point(curve, X3, L3, new ECFieldElement[] { Z3 }, IsCompressed);
         }
diff --git a/crypto/src/math/raw/Interleave.cs b/crypto/src/math/raw/Interleave.cs
index 591ba3f15..49d3768d7 100644
--- a/crypto/src/math/raw/Interleave.cs
+++ b/crypto/src/math/raw/Interleave.cs
@@ -93,6 +93,23 @@ namespace Org.BouncyCastle.Math.Raw
             z[zOff + 1] = (x >> 1) & M64;
         }
 
+        internal static void Expand64To128(ulong[] xs, int xsOff, int xsLen, ulong[] zs, int zsOff)
+        {
+            for (int i = 0; i < xsLen; ++i)
+            {
+                // "shuffle" low half to even bits and high half to odd bits
+                ulong x = xs[xsOff + i], t;
+                t = (x ^ (x >> 16)) & 0x00000000FFFF0000UL; x ^= (t ^ (t << 16));
+                t = (x ^ (x >>  8)) & 0x0000FF000000FF00UL; x ^= (t ^ (t <<  8));
+                t = (x ^ (x >>  4)) & 0x00F000F000F000F0UL; x ^= (t ^ (t <<  4));
+                t = (x ^ (x >>  2)) & 0x0C0C0C0C0C0C0C0CUL; x ^= (t ^ (t <<  2));
+                t = (x ^ (x >>  1)) & 0x2222222222222222UL; x ^= (t ^ (t <<  1));
+
+                zs[zsOff++] = (x     ) & M64;
+                zs[zsOff++] = (x >> 1) & M64;
+            }
+        }
+
         internal static void Expand64To128Rev(ulong x, ulong[] z, int zOff)
         {
             // "shuffle" low half to even bits and high half to odd bits
diff --git a/crypto/src/math/raw/Nat.cs b/crypto/src/math/raw/Nat.cs
index d67de0a5c..effe46454 100644
--- a/crypto/src/math/raw/Nat.cs
+++ b/crypto/src/math/raw/Nat.cs
@@ -1406,5 +1406,13 @@ namespace Org.BouncyCastle.Math.Raw
                 z[i] = 0;
             }
         }
+
+        public static void Zero64(int len, ulong[] z)
+        {
+            for (int i = 0; i < len; ++i)
+            {
+                z[i] = 0UL;
+            }
+        }
     }
 }