Optimized Sqrt and Trace for custom binary curves

20 files changed, 327 insertions, 9 deletions

diff --git a/crypto/src/math/ec/custom/sec/SecT113Field.cs b/crypto/src/math/ec/custom/sec/SecT113Field.cs
index 640c6e787..49773b66d 100644
--- a/crypto/src/math/ec/custom/sec/SecT113Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT113Field.cs
@@ -102,6 +102,16 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             z[zOff + 1]  = z1 & M49;
         }
 
+        public static void Sqrt(ulong[] x, ulong[] z)
+        {
+            ulong u0 = Interleave.Unshuffle(x[0]), u1 = Interleave.Unshuffle(x[1]);
+            ulong e0 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            ulong c0  = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            z[0] = e0 ^ (c0 << 57) ^ (c0 <<  5);
+            z[1] =      (c0 >>  7) ^ (c0 >> 59); 
+        }
+
         public static void Square(ulong[] x, ulong[] z)
         {
             ulong[] tt = Nat128.CreateExt64();
@@ -131,6 +141,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             }
         }
 
+        public static uint Trace(ulong[] x)
+        {
+            // Non-zero-trace bits: 0
+            return (uint)(x[0]) & 1U;
+        }
+
         protected static void ImplMultiply(ulong[] x, ulong[] y, ulong[] zz)
         {
             /*
diff --git a/crypto/src/math/ec/custom/sec/SecT113FieldElement.cs b/crypto/src/math/ec/custom/sec/SecT113FieldElement.cs
index f217e28cb..9ba25d987 100644
--- a/crypto/src/math/ec/custom/sec/SecT113FieldElement.cs
+++ b/crypto/src/math/ec/custom/sec/SecT113FieldElement.cs
@@ -159,7 +159,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public override ECFieldElement Sqrt()
         {
-            return SquarePow(M - 1);
+            ulong[] z = Nat128.Create64();
+            SecT113Field.Sqrt(x, z);
+            return new SecT113FieldElement(z);
         }
 
         public virtual int Representation
diff --git a/crypto/src/math/ec/custom/sec/SecT131Field.cs b/crypto/src/math/ec/custom/sec/SecT131Field.cs
index 47f97078c..1b6697afe 100644
--- a/crypto/src/math/ec/custom/sec/SecT131Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT131Field.cs
@@ -10,6 +10,8 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
         private const ulong M03 = ulong.MaxValue >> 61;
         private const ulong M44 = ulong.MaxValue >> 20;
 
+        private static readonly ulong[] ROOT_Z = new ulong[]{ 0x26BC4D789AF13523UL, 0x26BC4D789AF135E2UL, 0x6UL };
+
         public static void Add(ulong[] x, ulong[] y, ulong[] z)
         {
             z[0] = x[0] ^ y[0];
@@ -109,6 +111,25 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             z[zOff + 2]  = z2 & M03;
         }
 
+        public static void Sqrt(ulong[] x, ulong[] z)
+        {
+            ulong[] odd = Nat192.Create64();
+
+            ulong u0, u1;
+            u0 = Interleave.Unshuffle(x[0]); u1 = Interleave.Unshuffle(x[1]);
+            ulong e0 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            odd[0]   = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            u0 = Interleave.Unshuffle(x[2]);
+            ulong e1 = (u0 & 0x00000000FFFFFFFFUL);
+            odd[1]   = (u0 >> 32);
+
+            Multiply(odd, ROOT_Z, z);
+
+            z[0] ^= e0;
+            z[1] ^= e1;
+        }
+
         public static void Square(ulong[] x, ulong[] z)
         {
             ulong[] tt = Nat.Create64(5);
@@ -138,6 +159,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             }
         }
 
+        public static uint Trace(ulong[] x)
+        {
+            // Non-zero-trace bits: 0, 123, 129
+            return (uint)(x[0] ^ (x[1] >> 59) ^ (x[2] >> 1)) & 1U;
+        }
+
         protected static void ImplCompactExt(ulong[] zz)
         {
             ulong z0 = zz[0], z1 = zz[1], z2 = zz[2], z3 = zz[3], z4 = zz[4], z5 = zz[5];
diff --git a/crypto/src/math/ec/custom/sec/SecT131FieldElement.cs b/crypto/src/math/ec/custom/sec/SecT131FieldElement.cs
index 0ea00ea07..e0ecc100f 100644
--- a/crypto/src/math/ec/custom/sec/SecT131FieldElement.cs
+++ b/crypto/src/math/ec/custom/sec/SecT131FieldElement.cs
@@ -159,7 +159,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public override ECFieldElement Sqrt()
         {
-            return SquarePow(M - 1);
+            ulong[] z = Nat192.Create64();
+            SecT131Field.Sqrt(x, z);
+            return new SecT131FieldElement(z);
         }
 
         public virtual int Representation
diff --git a/crypto/src/math/ec/custom/sec/SecT163Field.cs b/crypto/src/math/ec/custom/sec/SecT163Field.cs
index f921a5bc7..b1e9aa725 100644
--- a/crypto/src/math/ec/custom/sec/SecT163Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT163Field.cs
@@ -10,6 +10,8 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
         private const ulong M35 = ulong.MaxValue >> 29;
         private const ulong M55 = ulong.MaxValue >> 9;
 
+        private static readonly ulong[] ROOT_Z = new ulong[]{ 0xB6DB6DB6DB6DB6B0UL, 0x492492492492DB6DUL, 0x492492492UL };
+
         public static void Add(ulong[] x, ulong[] y, ulong[] z)
         {
             z[0] = x[0] ^ y[0];
@@ -122,6 +124,25 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             z[zOff + 2]  = z2 & M35;
         }
 
+        public static void Sqrt(ulong[] x, ulong[] z)
+        {
+            ulong[] odd = Nat192.Create64();
+
+            ulong u0, u1;
+            u0 = Interleave.Unshuffle(x[0]); u1 = Interleave.Unshuffle(x[1]);
+            ulong e0 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            odd[0]   = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            u0 = Interleave.Unshuffle(x[2]);
+            ulong e1 = (u0 & 0x00000000FFFFFFFFUL);
+            odd[1]   = (u0 >> 32);
+
+            Multiply(odd, ROOT_Z, z);
+
+            z[0] ^= e0;
+            z[1] ^= e1;
+        }
+
         public static void Square(ulong[] x, ulong[] z)
         {
             ulong[] tt = Nat192.CreateExt64();
@@ -151,6 +172,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             }
         }
 
+        public static uint Trace(ulong[] x)
+        {
+            // Non-zero-trace bits: 0, 157
+            return (uint)(x[0] ^ (x[2] >> 29)) & 1U;
+        }
+
         protected static void ImplCompactExt(ulong[] zz)
         {
             ulong z0 = zz[0], z1 = zz[1], z2 = zz[2], z3 = zz[3], z4 = zz[4], z5 = zz[5];
diff --git a/crypto/src/math/ec/custom/sec/SecT163FieldElement.cs b/crypto/src/math/ec/custom/sec/SecT163FieldElement.cs
index c7a0b5639..8953fb529 100644
--- a/crypto/src/math/ec/custom/sec/SecT163FieldElement.cs
+++ b/crypto/src/math/ec/custom/sec/SecT163FieldElement.cs
@@ -159,7 +159,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public override ECFieldElement Sqrt()
         {
-            return SquarePow(M - 1);
+            ulong[] z = Nat192.Create64();
+            SecT163Field.Sqrt(x, z);
+            return new SecT163FieldElement(z);
         }
 
         public virtual int Representation
diff --git a/crypto/src/math/ec/custom/sec/SecT193Field.cs b/crypto/src/math/ec/custom/sec/SecT193Field.cs
index 5154f1e0a..41acb4f94 100644
--- a/crypto/src/math/ec/custom/sec/SecT193Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT193Field.cs
@@ -132,6 +132,23 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             z[zOff + 3]  = z3 & M01;
         }
 
+        public static void Sqrt(ulong[] x, ulong[] z)
+        {
+            ulong u0, u1;
+            u0 = Interleave.Unshuffle(x[0]); u1 = Interleave.Unshuffle(x[1]);
+            ulong e0 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            ulong c0 = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            u0 = Interleave.Unshuffle(x[2]);
+            ulong e1 = (u0 & 0x00000000FFFFFFFFUL) ^ (x[3] << 32);
+            ulong c1 = (u0 >> 32);
+
+            z[0] = e0 ^ (c0 << 8);
+            z[1] = e1 ^ (c1 << 8) ^ (c0 >> 56) ^ (c0 << 33);
+            z[2] =                  (c1 >> 56) ^ (c1 << 33) ^ (c0 >> 31);
+            z[3] =                                            (c1 >> 31);
+        }
+
         public static void Square(ulong[] x, ulong[] z)
         {
             ulong[] tt = Nat256.CreateExt64();
@@ -161,6 +178,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             }
         }
 
+        public static uint Trace(ulong[] x)
+        {
+            // Non-zero-trace bits: 0
+            return (uint)(x[0]) & 1U;
+        }
+
         protected static void ImplCompactExt(ulong[] zz)
         {
             ulong z0 = zz[0], z1 = zz[1], z2 = zz[2], z3 = zz[3], z4 = zz[4], z5 = zz[5], z6 = zz[6], z7 = zz[7];
diff --git a/crypto/src/math/ec/custom/sec/SecT193FieldElement.cs b/crypto/src/math/ec/custom/sec/SecT193FieldElement.cs
index eba4d10e6..a1150b3f9 100644
--- a/crypto/src/math/ec/custom/sec/SecT193FieldElement.cs
+++ b/crypto/src/math/ec/custom/sec/SecT193FieldElement.cs
@@ -159,7 +159,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public override ECFieldElement Sqrt()
         {
-            return SquarePow(M - 1);
+            ulong[] z = Nat256.Create64();
+            SecT193Field.Sqrt(x, z);
+            return new SecT193FieldElement(z);
         }
 
         public virtual int Representation
diff --git a/crypto/src/math/ec/custom/sec/SecT233Field.cs b/crypto/src/math/ec/custom/sec/SecT233Field.cs
index a2f73fd5d..870dade50 100644
--- a/crypto/src/math/ec/custom/sec/SecT233Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT233Field.cs
@@ -128,6 +128,41 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             z[zOff + 3]  = z3 & M41;
         }
 
+        public static void Sqrt(ulong[] x, ulong[] z)
+        {
+            ulong u0, u1;
+            u0 = Interleave.Unshuffle(x[0]); u1 = Interleave.Unshuffle(x[1]);
+            ulong e0 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            ulong c0 = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            u0 = Interleave.Unshuffle(x[2]); u1 = Interleave.Unshuffle(x[3]);
+            ulong e1 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            ulong c1 = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            ulong c2;
+            c2  = (c1 >> 27);
+            c1 ^= (c0 >> 27) | (c1 << 37);
+            c0 ^=              (c0 << 37);
+
+            ulong[] tt = Nat256.CreateExt64();
+
+            int[] shifts = { 32, 117, 191 };
+            for (int i = 0; i < shifts.Length; ++i)
+            {
+                int w = shifts[i] >> 6, s = shifts[i] & 63;
+                Debug.Assert(s != 0);
+                tt[w    ] ^= (c0 << s);
+                tt[w + 1] ^= (c1 << s) | (c0 >> -s);
+                tt[w + 2] ^= (c2 << s) | (c1 >> -s);
+                tt[w + 3] ^=             (c2 >> -s);
+            }
+
+            Reduce(tt, z);
+
+            z[0] ^= e0;
+            z[1] ^= e1;
+        }
+
         public static void Square(ulong[] x, ulong[] z)
         {
             ulong[] tt = Nat256.CreateExt64();
@@ -157,6 +192,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             }
         }
 
+        public static uint Trace(ulong[] x)
+        {
+            // Non-zero-trace bits: 0, 159
+            return (uint)(x[0] ^ (x[2] >> 31)) & 1U;
+        }
+
         protected static void ImplCompactExt(ulong[] zz)
         {
             ulong z0 = zz[0], z1 = zz[1], z2 = zz[2], z3 = zz[3], z4 = zz[4], z5 = zz[5], z6 = zz[6], z7 = zz[7];
diff --git a/crypto/src/math/ec/custom/sec/SecT233FieldElement.cs b/crypto/src/math/ec/custom/sec/SecT233FieldElement.cs
index a9041efde..91b8e2f1c 100644
--- a/crypto/src/math/ec/custom/sec/SecT233FieldElement.cs
+++ b/crypto/src/math/ec/custom/sec/SecT233FieldElement.cs
@@ -159,7 +159,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public override ECFieldElement Sqrt()
         {
-            return SquarePow(M - 1);
+            ulong[] z = Nat256.Create64();
+            SecT233Field.Sqrt(x, z);
+            return new SecT233FieldElement(z);
         }
 
         public virtual int Representation
diff --git a/crypto/src/math/ec/custom/sec/SecT239Field.cs b/crypto/src/math/ec/custom/sec/SecT239Field.cs
index 6b8ad696f..2e6ed2ad6 100644
--- a/crypto/src/math/ec/custom/sec/SecT239Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT239Field.cs
@@ -136,6 +136,42 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             z[zOff + 3]  = z3 & M47;
         }
 
+        public static void Sqrt(ulong[] x, ulong[] z)
+        {
+            ulong u0, u1;
+            u0 = Interleave.Unshuffle(x[0]); u1 = Interleave.Unshuffle(x[1]);
+            ulong e0 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            ulong c0 = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            u0 = Interleave.Unshuffle(x[2]); u1 = Interleave.Unshuffle(x[3]);
+            ulong e1 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            ulong c1 = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            ulong c2, c3;
+            c3  = (c1 >> 49);
+            c2  = (c0 >> 49) | (c1 << 15);
+            c1 ^=              (c0 << 15);
+
+            ulong[] tt = Nat256.CreateExt64();
+
+            int[] shifts = { 39, 120 };
+            for (int i = 0; i < shifts.Length; ++i)
+            {
+                int w = shifts[i] >> 6, s = shifts[i] & 63;
+                Debug.Assert(s != 0);
+                tt[w    ] ^= (c0 << s);
+                tt[w + 1] ^= (c1 << s) | (c0 >> -s);
+                tt[w + 2] ^= (c2 << s) | (c1 >> -s);
+                tt[w + 3] ^= (c3 << s) | (c2 >> -s);
+                tt[w + 4] ^=             (c3 >> -s);
+            }
+
+            Reduce(tt, z);
+
+            z[0] ^= e0;
+            z[1] ^= e1;
+        }
+
         public static void Square(ulong[] x, ulong[] z)
         {
             ulong[] tt = Nat256.CreateExt64();
@@ -165,6 +201,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             }
         }
 
+        public static uint Trace(ulong[] x)
+        {
+            // Non-zero-trace bits: 0, 81, 162
+            return (uint)(x[0] ^ (x[1] >> 17) ^ (x[2] >> 34)) & 1U;
+        }
+
         protected static void ImplCompactExt(ulong[] zz)
         {
             ulong z0 = zz[0], z1 = zz[1], z2 = zz[2], z3 = zz[3], z4 = zz[4], z5 = zz[5], z6 = zz[6], z7 = zz[7];
diff --git a/crypto/src/math/ec/custom/sec/SecT239FieldElement.cs b/crypto/src/math/ec/custom/sec/SecT239FieldElement.cs
index de074c55f..a32ffc5d2 100644
--- a/crypto/src/math/ec/custom/sec/SecT239FieldElement.cs
+++ b/crypto/src/math/ec/custom/sec/SecT239FieldElement.cs
@@ -159,7 +159,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public override ECFieldElement Sqrt()
         {
-            return SquarePow(M - 1);
+            ulong[] z = Nat256.Create64();
+            SecT239Field.Sqrt(x, z);
+            return new SecT239FieldElement(z);
         }
 
         public virtual int Representation
diff --git a/crypto/src/math/ec/custom/sec/SecT283Field.cs b/crypto/src/math/ec/custom/sec/SecT283Field.cs
index 903ea02ff..22b7eaaab 100644
--- a/crypto/src/math/ec/custom/sec/SecT283Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT283Field.cs
@@ -10,6 +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 };
+
         public static void Add(ulong[] x, ulong[] y, ulong[] z)
         {
             z[0] = x[0] ^ y[0];
@@ -129,6 +131,30 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             z[zOff + 4]  = z4 & M27;
         }
 
+        public static void Sqrt(ulong[] x, ulong[] z)
+        {
+            ulong[] odd = Nat320.Create64();
+
+            ulong u0, u1;
+            u0 = Interleave.Unshuffle(x[0]); u1 = Interleave.Unshuffle(x[1]);
+            ulong e0 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            odd[0]   = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            u0 = Interleave.Unshuffle(x[2]); u1 = Interleave.Unshuffle(x[3]);
+            ulong e1 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            odd[1]   = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            u0 = Interleave.Unshuffle(x[4]);
+            ulong e2 = (u0 & 0x00000000FFFFFFFFUL);
+            odd[2]   = (u0 >> 32);
+
+            Multiply(odd, ROOT_Z, z);
+
+            z[0] ^= e0;
+            z[1] ^= e1;
+            z[2] ^= e2;
+        }
+
         public static void Square(ulong[] x, ulong[] z)
         {
             ulong[] tt = Nat.Create64(9);
@@ -158,6 +184,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             }
         }
 
+        public static uint Trace(ulong[] x)
+        {
+            // Non-zero-trace bits: 0, 271
+            return (uint)(x[0] ^ (x[4] >> 15)) & 1U;
+        }
+
         protected static void ImplCompactExt(ulong[] zz)
         {
             ulong z0 = zz[0], z1 = zz[1], z2 = zz[2], z3 = zz[3], z4 = zz[4];
diff --git a/crypto/src/math/ec/custom/sec/SecT283FieldElement.cs b/crypto/src/math/ec/custom/sec/SecT283FieldElement.cs
index e02108f73..adfd4e0ed 100644
--- a/crypto/src/math/ec/custom/sec/SecT283FieldElement.cs
+++ b/crypto/src/math/ec/custom/sec/SecT283FieldElement.cs
@@ -159,7 +159,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public override ECFieldElement Sqrt()
         {
-            return SquarePow(M - 1);
+            ulong[] z = Nat320.Create64();
+            SecT283Field.Sqrt(x, z);
+            return new SecT283FieldElement(z);
         }
 
         public virtual int Representation
diff --git a/crypto/src/math/ec/custom/sec/SecT409Field.cs b/crypto/src/math/ec/custom/sec/SecT409Field.cs
index 84eada96e..861b77aa1 100644
--- a/crypto/src/math/ec/custom/sec/SecT409Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT409Field.cs
@@ -165,6 +165,36 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             z[zOff + 6]  = z6 & M25;
         }
 
+        public static void Sqrt(ulong[] x, ulong[] z)
+        {
+            ulong u0, u1;
+            u0 = Interleave.Unshuffle(x[0]); u1 = Interleave.Unshuffle(x[1]);
+            ulong e0 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            ulong c0 = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            u0 = Interleave.Unshuffle(x[2]); u1 = Interleave.Unshuffle(x[3]);
+            ulong e1 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            ulong c1 = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            u0 = Interleave.Unshuffle(x[4]); u1 = Interleave.Unshuffle(x[5]);
+            ulong e2 = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+            ulong c2 = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+
+            u0 = Interleave.Unshuffle(x[6]);
+            ulong e3 = (u0 & 0x00000000FFFFFFFFUL);
+            ulong c3 = (u0 >> 32);
+
+            z[0] = e0 ^ (c0 << 44);
+            z[1] = e1 ^ (c1 << 44) ^ (c0 >> 20);
+            z[2] = e2 ^ (c2 << 44) ^ (c1 >> 20);
+            z[3] = e3 ^ (c3 << 44) ^ (c2 >> 20) ^ (c0 << 13);
+            z[4] =                   (c3 >> 20) ^ (c1 << 13) ^ (c0 >> 51);
+            z[5] =                                (c2 << 13) ^ (c1 >> 51);
+            z[6] =                                (c3 << 13) ^ (c2 >> 51);
+
+            Debug.Assert((c3 >> 51) == 0);
+        }
+
         public static void Square(ulong[] x, ulong[] z)
         {
             ulong[] tt = Nat.Create64(13);
@@ -194,6 +224,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             }
         }
 
+        public static uint Trace(ulong[] x)
+        {
+            // Non-zero-trace bits: 0
+            return (uint)(x[0]) & 1U;
+        }
+
         protected static void ImplCompactExt(ulong[] zz)
         {
             ulong z00 = zz[ 0], z01 = zz[ 1], z02 = zz[ 2], z03 = zz[ 3], z04 = zz[ 4], z05 = zz[ 5], z06 = zz[ 6];
diff --git a/crypto/src/math/ec/custom/sec/SecT409FieldElement.cs b/crypto/src/math/ec/custom/sec/SecT409FieldElement.cs
index 581ea73df..f954f46e7 100644
--- a/crypto/src/math/ec/custom/sec/SecT409FieldElement.cs
+++ b/crypto/src/math/ec/custom/sec/SecT409FieldElement.cs
@@ -159,7 +159,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public override ECFieldElement Sqrt()
         {
-            return SquarePow(M - 1);
+            ulong[] z = Nat448.Create64();
+            SecT409Field.Sqrt(x, z);
+            return new SecT409FieldElement(z);
         }
 
         public virtual int Representation
diff --git a/crypto/src/math/ec/custom/sec/SecT571Field.cs b/crypto/src/math/ec/custom/sec/SecT571Field.cs
index fc84e336b..98f4f7fc2 100644
--- a/crypto/src/math/ec/custom/sec/SecT571Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecT571Field.cs
@@ -11,6 +11,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         private const ulong RM = 0xEF7BDEF7BDEF7BDEUL;
 
+        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)
         {
             for (int i = 0; i < 9; ++i)
@@ -156,6 +159,28 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             z[zOff + 8]  = z8 & M59;
         }
 
+        public static void Sqrt(ulong[] x, ulong[] z)
+        {
+            ulong[] evn = Nat576.Create64(), odd = Nat576.Create64();
+
+            int pos = 0;
+            for (int i = 0; i < 4; ++i)
+            {
+                ulong u0 = Interleave.Unshuffle(x[pos++]);
+                ulong u1 = Interleave.Unshuffle(x[pos++]);
+                evn[i] = (u0 & 0x00000000FFFFFFFFUL) | (u1 << 32);
+                odd[i] = (u0 >> 32) | (u1 & 0xFFFFFFFF00000000UL);
+            }
+            {
+                ulong u0 = Interleave.Unshuffle(x[pos]);
+                evn[4] = (u0 & 0x00000000FFFFFFFFUL);
+                odd[4] = (u0 >> 32);
+            }
+
+            Multiply(odd, ROOT_Z, z);
+            Add(z, evn, z);
+        }
+
         public static void Square(ulong[] x, ulong[] z)
         {
             ulong[] tt = Nat576.CreateExt64();
@@ -185,6 +210,12 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
             }
         }
 
+        public static uint Trace(ulong[] x)
+        {
+            // Non-zero-trace bits: 0, 561, 569
+            return (uint)(x[0] ^ (x[8] >> 49) ^ (x[8] >> 57)) & 1U;
+        }
+
         protected static void ImplMultiply(ulong[] x, ulong[] y, ulong[] zz)
         {
             //for (int i = 0; i < 9; ++i)
diff --git a/crypto/src/math/ec/custom/sec/SecT571FieldElement.cs b/crypto/src/math/ec/custom/sec/SecT571FieldElement.cs
index 5d5458412..c43b8dc3a 100644
--- a/crypto/src/math/ec/custom/sec/SecT571FieldElement.cs
+++ b/crypto/src/math/ec/custom/sec/SecT571FieldElement.cs
@@ -159,7 +159,9 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public override ECFieldElement Sqrt()
         {
-            return SquarePow(M - 1);
+            ulong[] z = Nat576.Create64();
+            SecT571Field.Sqrt(x, z);
+            return new SecT571FieldElement(z);
         }
 
         public virtual int Representation
diff --git a/crypto/src/math/raw/Interleave.cs b/crypto/src/math/raw/Interleave.cs
index a45ee1e08..d21840644 100644
--- a/crypto/src/math/raw/Interleave.cs
+++ b/crypto/src/math/raw/Interleave.cs
@@ -91,5 +91,17 @@ namespace Org.BouncyCastle.Math.Raw
             z[zOff    ] = (x     ) & M64;
             z[zOff + 1] = (x >> 1) & M64;
         }
+
+        internal static ulong Unshuffle(ulong x)
+        {
+            // "unshuffle" even bits to low half and odd bits to high half
+            ulong t;
+            t = (x ^ (x >>  1)) & 0x2222222222222222UL; x ^= (t ^ (t <<  1));
+            t = (x ^ (x >>  2)) & 0x0C0C0C0C0C0C0C0CUL; x ^= (t ^ (t <<  2));
+            t = (x ^ (x >>  4)) & 0x00F000F000F000F0UL; x ^= (t ^ (t <<  4));
+            t = (x ^ (x >>  8)) & 0x0000FF000000FF00UL; x ^= (t ^ (t <<  8));
+            t = (x ^ (x >> 16)) & 0x00000000FFFF0000UL; x ^= (t ^ (t << 16));
+            return x;
+        }
     }
 }
diff --git a/crypto/test/src/math/ec/test/ECPointTest.cs b/crypto/test/src/math/ec/test/ECPointTest.cs
index 3c10170f7..57f838f9f 100644
--- a/crypto/test/src/math/ec/test/ECPointTest.cs
+++ b/crypto/test/src/math/ec/test/ECPointTest.cs
@@ -437,6 +437,19 @@ namespace Org.BouncyCastle.Math.EC.Tests
                     }
                 }
             }
+            else if (ECAlgorithms.IsF2mCurve(c))
+            {
+                int m = c.FieldSize;
+                BigInteger x = new BigInteger(m, secRand);
+                ECFieldElement fe = c.FromBigInteger(x);
+                for (int i = 0; i < 100; ++i)
+                {
+                    ECFieldElement sq = fe.Square();
+                    ECFieldElement check = sq.Sqrt();
+                    Assert.AreEqual(fe, check);
+                    fe = sq;
+                }
+            }
         }
 
         private void ImplAddSubtractMultiplyTwiceEncodingTestAllCoords(X9ECParameters x9ECParameters)