'safegcd' modular inversion

31 files changed, 770 insertions, 846 deletions

diff --git a/crypto/src/bcpg/RsaSecretBcpgKey.cs b/crypto/src/bcpg/RsaSecretBcpgKey.cs
index 5c04d9f85..783f083ce 100644
--- a/crypto/src/bcpg/RsaSecretBcpgKey.cs
+++ b/crypto/src/bcpg/RsaSecretBcpgKey.cs
@@ -1,6 +1,7 @@
 using System;
 
 using Org.BouncyCastle.Math;
+using Org.BouncyCastle.Utilities;
 
 namespace Org.BouncyCastle.Bcpg
 {
@@ -21,7 +22,7 @@ namespace Org.BouncyCastle.Bcpg
 
 			this.expP = d.Value.Remainder(p.Value.Subtract(BigInteger.One));
 			this.expQ = d.Value.Remainder(q.Value.Subtract(BigInteger.One));
-			this.crt = q.Value.ModInverse(p.Value);
+			this.crt = BigIntegers.ModOddInverse(p.Value, q.Value);
 		}
 
 		public RsaSecretBcpgKey(
@@ -44,11 +45,11 @@ namespace Org.BouncyCastle.Bcpg
 			this.d = new MPInteger(d);
 			this.p = new MPInteger(p);
 			this.q = new MPInteger(q);
-			this.u = new MPInteger(p.ModInverse(q));
+			this.u = new MPInteger(BigIntegers.ModOddInverse(q, p));
 
 			this.expP = d.Remainder(p.Subtract(BigInteger.One));
 			this.expQ = d.Remainder(q.Subtract(BigInteger.One));
-			this.crt = q.ModInverse(p);
+			this.crt = BigIntegers.ModOddInverse(p, q);
 		}
 
 		public BigInteger Modulus
diff --git a/crypto/src/crypto/engines/RSABlindedEngine.cs b/crypto/src/crypto/engines/RSABlindedEngine.cs
index 98108f9b4..637bf3cc0 100644
--- a/crypto/src/crypto/engines/RSABlindedEngine.cs
+++ b/crypto/src/crypto/engines/RSABlindedEngine.cs
@@ -132,7 +132,7 @@ namespace Org.BouncyCastle.Crypto.Engines
                     BigInteger blindedInput = r.ModPow(e, m).Multiply(input).Mod(m);
                     BigInteger blindedResult = core.ProcessBlock(blindedInput);
 
-                    BigInteger rInv = r.ModInverse(m);
+                    BigInteger rInv = BigIntegers.ModOddInverse(m, r);
                     result = blindedResult.Multiply(rInv).Mod(m);
 
                     // defence against Arjen Lenstra’s CRT attack
diff --git a/crypto/src/crypto/engines/RSABlindingEngine.cs b/crypto/src/crypto/engines/RSABlindingEngine.cs
index 1289456a6..11bb8d9d9 100644
--- a/crypto/src/crypto/engines/RSABlindingEngine.cs
+++ b/crypto/src/crypto/engines/RSABlindingEngine.cs
@@ -2,6 +2,7 @@ using System;
 
 using Org.BouncyCastle.Crypto.Parameters;
 using Org.BouncyCastle.Math;
+using Org.BouncyCastle.Utilities;
 
 namespace Org.BouncyCastle.Crypto.Engines
 {
@@ -139,7 +140,7 @@ namespace Org.BouncyCastle.Crypto.Engines
 		{
 			BigInteger m = key.Modulus;
 			BigInteger msg = blindedMsg;
-			BigInteger blindFactorInverse = blindingFactor.ModInverse(m);
+			BigInteger blindFactorInverse = BigIntegers.ModOddInverse(m, blindingFactor);
 			msg = msg.Multiply(blindFactorInverse);
 			msg = msg.Mod(m);
 
diff --git a/crypto/src/crypto/generators/RsaKeyPairGenerator.cs b/crypto/src/crypto/generators/RsaKeyPairGenerator.cs
index 449976550..78c9eb18f 100644
--- a/crypto/src/crypto/generators/RsaKeyPairGenerator.cs
+++ b/crypto/src/crypto/generators/RsaKeyPairGenerator.cs
@@ -127,7 +127,7 @@ namespace Org.BouncyCastle.Crypto.Generators
                 //
                 BigInteger dP = d.Remainder(pSub1);
                 BigInteger dQ = d.Remainder(qSub1);
-                BigInteger qInv = q.ModInverse(p);
+                BigInteger qInv = BigIntegers.ModOddInverse(p, q);
 
                 return new AsymmetricCipherKeyPair(
                     new RsaKeyParameters(false, n, e),
diff --git a/crypto/src/crypto/parameters/ECDomainParameters.cs b/crypto/src/crypto/parameters/ECDomainParameters.cs
index ce17ab6df..b5ca183de 100644
--- a/crypto/src/crypto/parameters/ECDomainParameters.cs
+++ b/crypto/src/crypto/parameters/ECDomainParameters.cs
@@ -89,7 +89,7 @@ namespace Org.BouncyCastle.Crypto.Parameters
                 {
                     if (hInv == null)
                     {
-                        hInv = h.ModInverse(n);
+                        hInv = BigIntegers.ModOddInverseVar(n, h);
                     }
                     return hInv;
                 }
diff --git a/crypto/src/crypto/signers/DsaSigner.cs b/crypto/src/crypto/signers/DsaSigner.cs
index 32a33d084..1f5d9b937 100644
--- a/crypto/src/crypto/signers/DsaSigner.cs
+++ b/crypto/src/crypto/signers/DsaSigner.cs
@@ -1,9 +1,9 @@
 using System;
 
-using Org.BouncyCastle.Crypto.Digests;
 using Org.BouncyCastle.Crypto.Parameters;
 using Org.BouncyCastle.Math;
 using Org.BouncyCastle.Security;
+using Org.BouncyCastle.Utilities;
 
 namespace Org.BouncyCastle.Crypto.Signers
 {
@@ -104,7 +104,7 @@ namespace Org.BouncyCastle.Crypto.Signers
 
             BigInteger r = parameters.G.ModPow(k, parameters.P).Mod(q);
 
-            k = k.ModInverse(q).Multiply(m.Add(x.Multiply(r)));
+            k = BigIntegers.ModOddInverse(q, k).Multiply(m.Add(x.Multiply(r)));
 
             BigInteger s = k.Mod(q);
 
@@ -132,7 +132,7 @@ namespace Org.BouncyCastle.Crypto.Signers
                 return false;
             }
 
-            BigInteger w = s.ModInverse(q);
+            BigInteger w = BigIntegers.ModOddInverseVar(q, s);
 
             BigInteger u1 = m.Multiply(w).Mod(q);
             BigInteger u2 = r.Multiply(w).Mod(q);
diff --git a/crypto/src/crypto/signers/ECDsaSigner.cs b/crypto/src/crypto/signers/ECDsaSigner.cs
index 7b04d1076..0a265d96e 100644
--- a/crypto/src/crypto/signers/ECDsaSigner.cs
+++ b/crypto/src/crypto/signers/ECDsaSigner.cs
@@ -6,6 +6,7 @@ using Org.BouncyCastle.Math;
 using Org.BouncyCastle.Math.EC;
 using Org.BouncyCastle.Math.EC.Multiplier;
 using Org.BouncyCastle.Security;
+using Org.BouncyCastle.Utilities;
 
 namespace Org.BouncyCastle.Crypto.Signers
 {
@@ -123,7 +124,7 @@ namespace Org.BouncyCastle.Crypto.Signers
                 }
                 while (r.SignValue == 0);
 
-                s = k.ModInverse(n).Multiply(e.Add(d.Multiply(r))).Mod(n);
+                s = BigIntegers.ModOddInverse(n, k).Multiply(e.Add(d.Multiply(r))).Mod(n);
             }
             while (s.SignValue == 0);
 
@@ -148,7 +149,7 @@ namespace Org.BouncyCastle.Crypto.Signers
             }
 
             BigInteger e = CalculateE(n, message);
-            BigInteger c = s.ModInverse(n);
+            BigInteger c = BigIntegers.ModOddInverseVar(n, s);
 
             BigInteger u1 = e.Multiply(c).Mod(n);
             BigInteger u2 = r.Multiply(c).Mod(n);
diff --git a/crypto/src/crypto/signers/ECGOST3410Signer.cs b/crypto/src/crypto/signers/ECGOST3410Signer.cs
index ac6b080a5..7df43f0a0 100644
--- a/crypto/src/crypto/signers/ECGOST3410Signer.cs
+++ b/crypto/src/crypto/signers/ECGOST3410Signer.cs
@@ -145,7 +145,7 @@ namespace Org.BouncyCastle.Crypto.Signers
                 return false;
             }
 
-            BigInteger v = e.ModInverse(n);
+            BigInteger v = BigIntegers.ModOddInverseVar(n, e);
 
             BigInteger z1 = s.Multiply(v).Mod(n);
             BigInteger z2 = (n.Subtract(r)).Multiply(v).Mod(n);
diff --git a/crypto/src/crypto/signers/SM2Signer.cs b/crypto/src/crypto/signers/SM2Signer.cs
index 2597cbf3d..c344a220a 100644
--- a/crypto/src/crypto/signers/SM2Signer.cs
+++ b/crypto/src/crypto/signers/SM2Signer.cs
@@ -6,6 +6,7 @@ using Org.BouncyCastle.Math;
 using Org.BouncyCastle.Math.EC;
 using Org.BouncyCastle.Math.EC.Multiplier;
 using Org.BouncyCastle.Security;
+using Org.BouncyCastle.Utilities;
 using Org.BouncyCastle.Utilities.Encoders;
 
 namespace Org.BouncyCastle.Crypto.Signers
@@ -164,7 +165,7 @@ namespace Org.BouncyCastle.Crypto.Signers
                 while (r.SignValue == 0 || r.Add(k).Equals(n));
 
                 // A6
-                BigInteger dPlus1ModN = d.Add(BigInteger.One).ModInverse(n);
+                BigInteger dPlus1ModN = BigIntegers.ModOddInverse(n, d.Add(BigIntegers.One));
 
                 s = k.Subtract(r.Multiply(d)).Mod(n);
                 s = dPlus1ModN.Multiply(s).Mod(n);
diff --git a/crypto/src/math/ec/ECFieldElement.cs b/crypto/src/math/ec/ECFieldElement.cs
index ef10dbf90..ed530b6b7 100644
--- a/crypto/src/math/ec/ECFieldElement.cs
+++ b/crypto/src/math/ec/ECFieldElement.cs
@@ -1,7 +1,6 @@
 using System;
 using System.Diagnostics;
 
-using Org.BouncyCastle.Math.Raw;
 using Org.BouncyCastle.Utilities;
 
 namespace Org.BouncyCastle.Math.EC
@@ -434,13 +433,7 @@ namespace Org.BouncyCastle.Math.EC
 
         protected virtual BigInteger ModInverse(BigInteger x)
         {
-            int bits = FieldSize;
-            int len = (bits + 31) >> 5;
-            uint[] p = Nat.FromBigInteger(bits, q);
-            uint[] n = Nat.FromBigInteger(bits, x);
-            uint[] z = Nat.Create(len);
-            Mod.Invert(p, n, z);
-            return Nat.ToBigInteger(len, z);
+            return BigIntegers.ModOddInverse(q, x);
         }
 
         protected virtual BigInteger ModMult(BigInteger x1, BigInteger x2)
diff --git a/crypto/src/math/ec/custom/djb/Curve25519Field.cs b/crypto/src/math/ec/custom/djb/Curve25519Field.cs
index 4e4cfbaa5..0006acd94 100644
--- a/crypto/src/math/ec/custom/djb/Curve25519Field.cs
+++ b/crypto/src/math/ec/custom/djb/Curve25519Field.cs
@@ -70,56 +70,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Djb
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^255 - 21
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 250 1s } { 1 0s } { 1 1s } { 1 0s } { 2 1s }
-             *
-             * Therefore we need an addition chain containing 1, 2, 250 (the lengths of the repunits)
-             * We use: [1], [2], 3, 5, 10, 15, 25, 50, 75, 125, [250]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat256.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x3 = Nat256.Create();
-            Square(x2, x3);
-            Multiply(x3, x1, x3);
-            uint[] x5 = x3;
-            SquareN(x3, 2, x5);
-            Multiply(x5, x2, x5);
-            uint[] x10 = Nat256.Create();
-            SquareN(x5, 5, x10);
-            Multiply(x10, x5, x10);
-            uint[] x15 = Nat256.Create();
-            SquareN(x10, 5, x15);
-            Multiply(x15, x5, x15);
-            uint[] x25 = x5;
-            SquareN(x15, 10, x25);
-            Multiply(x25, x10, x25);
-            uint[] x50 = x10;
-            SquareN(x25, 25, x50);
-            Multiply(x50, x25, x50);
-            uint[] x75 = x15;
-            SquareN(x50, 25, x75);
-            Multiply(x75, x25, x75);
-            uint[] x125 = x25;
-            SquareN(x75, 50, x125);
-            Multiply(x125, x50, x125);
-            uint[] x250 = x50;
-            SquareN(x125, 125, x250);
-            Multiply(x250, x125, x250);
-
-            uint[] t = x250;
-            SquareN(t, 2, t);
-            Multiply(t, x1, t);
-            SquareN(t, 3, t);
-            Multiply(t, x2, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/gm/SM2P256V1Field.cs b/crypto/src/math/ec/custom/gm/SM2P256V1Field.cs
index 55596b844..38743189a 100644
--- a/crypto/src/math/ec/custom/gm/SM2P256V1Field.cs
+++ b/crypto/src/math/ec/custom/gm/SM2P256V1Field.cs
@@ -57,60 +57,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.GM
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^256 - 2^224 - 2^96 + 2^64 - 3
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 31 1s } { 1 0s } { 128 1s } { 32 0s } { 62 1s } { 1 0s } { 1 1s }
-             *
-             * We use an addition chain for the beginning: [1], 2, [4], 6, 12, 24, 30, [31] 
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat256.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x4 = Nat256.Create();
-            SquareN(x2, 2, x4);
-            Multiply(x4, x2, x4);
-            uint[] x6 = Nat256.Create();
-            SquareN(x4, 2, x6);
-            Multiply(x6, x2, x6);
-            uint[] x12 = x2;
-            SquareN(x6, 6, x12);
-            Multiply(x12, x6, x12);
-            uint[] x24 = Nat256.Create();
-            SquareN(x12, 12, x24);
-            Multiply(x24, x12, x24);
-            uint[] x30 = x12;
-            SquareN(x24, 6, x30);
-            Multiply(x30, x6, x30);
-            uint[] x31 = x6;
-            Square(x30, x31);
-            Multiply(x31, x1, x31);
-
-            uint[] t = x24;
-            SquareN(x31, 32, t);
-            Multiply(t, x31, t);
-            SquareN(t, 31, t);
-            Multiply(t, x31, t);
-            SquareN(t, 31, t);
-            Multiply(t, x31, t);
-            SquareN(t, 31, t);
-            Multiply(t, x31, t);
-            SquareN(t, 4, t);
-            Multiply(t, x4, t);
-            SquareN(t, 63, t);
-            Multiply(t, x31, t);
-            SquareN(t, 31, t);
-            Multiply(t, x31, t);
-            SquareN(t, 2, t);
-
-            // NOTE that x1 and z could be the same array
-            Multiply(x1, t, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static void Half(uint[] x, uint[] z)
diff --git a/crypto/src/math/ec/custom/sec/SecP128R1Field.cs b/crypto/src/math/ec/custom/sec/SecP128R1Field.cs
index 23ea361a0..03a07f79b 100644
--- a/crypto/src/math/ec/custom/sec/SecP128R1Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP128R1Field.cs
@@ -70,51 +70,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^128 - 2^97 - 3
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 30 1s } { 1 0s } { 95 1s } { 1 0s } { 1 1s }
-             *
-             * We use an addition chain for the beginning: [1], 2, 3, [5], 10, 20, [30]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat128.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x3 = Nat128.Create();
-            Square(x2, x3);
-            Multiply(x3, x1, x3);
-            uint[] x5 = x3;
-            SquareN(x3, 2, x5);
-            Multiply(x5, x2, x5);
-            uint[] x10 = x2;
-            SquareN(x5, 5, x10);
-            Multiply(x10, x5, x10);
-            uint[] x20 = Nat128.Create();
-            SquareN(x10, 10, x20);
-            Multiply(x20, x10, x20);
-            uint[] x30 = x20;
-            SquareN(x20, 10, x30);
-            Multiply(x30, x10, x30);
-
-            uint[] t = x10;
-            SquareN(x30, 31, t);
-            Multiply(t, x30, t);
-            SquareN(t, 30, t);
-            Multiply(t, x30, t);
-            SquareN(t, 30, t);
-            Multiply(t, x30, t);
-            SquareN(t, 5, t);
-            Multiply(t, x5, t);
-            SquareN(t, 2, t);
-
-            // NOTE that x1 and z could be the same array
-            Multiply(x1, t, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/sec/SecP160R1Field.cs b/crypto/src/math/ec/custom/sec/SecP160R1Field.cs
index 139cd80d6..31c957301 100644
--- a/crypto/src/math/ec/custom/sec/SecP160R1Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP160R1Field.cs
@@ -74,58 +74,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^160 - 2^31 - 3
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 128 1s } { 1 0s } { 29 1s } { 1 0s } { 1 1s }
-             *
-             * Therefore we need an addition chain containing 1, 29, 128 (the lengths of the repunits)
-             * We use: [1], 2, 3, 6, 12, 24, 27, [29], 32, 64, [128]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat160.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x3 = Nat160.Create();
-            Square(x2, x3);
-            Multiply(x3, x1, x3);
-            uint[] x6 = Nat160.Create();
-            SquareN(x3, 3, x6);
-            Multiply(x6, x3, x6);
-            uint[] x12 = Nat160.Create();
-            SquareN(x6, 6, x12);
-            Multiply(x12, x6, x12);
-            uint[] x24 = x6;
-            SquareN(x12, 12, x24);
-            Multiply(x24, x12, x24);
-            uint[] x27 = x12;
-            SquareN(x24, 3, x27);
-            Multiply(x27, x3, x27);
-            uint[] x29 = x24;
-            SquareN(x27, 2, x29);
-            Multiply(x29, x2, x29);
-            uint[] x32 = x2;
-            SquareN(x29, 3, x32);
-            Multiply(x32, x3, x32);
-            uint[] x64 = x3;
-            SquareN(x32, 32, x64);
-            Multiply(x64, x32, x64);
-            uint[] x128 = x27;
-            SquareN(x64, 64, x128);
-            Multiply(x128, x64, x128);
-
-            uint[] t = x128;
-            SquareN(t, 30, t);
-            Multiply(t, x29, t);
-            SquareN(t, 2, t);
-
-            // NOTE that x1 and z could be the same array
-            Multiply(x1, t, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/sec/SecP160R2Field.cs b/crypto/src/math/ec/custom/sec/SecP160R2Field.cs
index bc36d9de1..55f02e438 100644
--- a/crypto/src/math/ec/custom/sec/SecP160R2Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP160R2Field.cs
@@ -74,70 +74,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^160 - 2^32 - 2^14 - 2^12 - 2^9 - 2^8 - 2^7 - 2^3 - 2^2 - 3
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 127 1s } { 1 0s } { 17 1s } "010110001110001"
-             *
-             * Therefore we need an addition chain containing 1, 2, 3, 17, 127 (the lengths of the repunits)
-             * We use: 1, 2, 3, 6, 12, 15, [17], 34, 68, 102, 119, 125, [127]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat160.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x3 = Nat160.Create();
-            Square(x2, x3);
-            Multiply(x3, x1, x3);
-            uint[] x6 = Nat160.Create();
-            SquareN(x3, 3, x6);
-            Multiply(x6, x3, x6);
-            uint[] x12 = Nat160.Create();
-            SquareN(x6, 6, x12);
-            Multiply(x12, x6, x12);
-            uint[] x15 = x12;
-            SquareN(x12, 3, x15);
-            Multiply(x15, x3, x15);
-            uint[] x17 = x15;
-            SquareN(x15, 2, x17);
-            Multiply(x17, x2, x17);
-            uint[] x34 = Nat160.Create();
-            SquareN(x17, 17, x34);
-            Multiply(x34, x17, x34);
-            uint[] x68 = Nat160.Create();
-            SquareN(x34, 34, x68);
-            Multiply(x68, x34, x68);
-            uint[] x102 = x68;
-            SquareN(x68, 34, x102);
-            Multiply(x102, x34, x102);
-            uint[] x119 = x34;
-            SquareN(x102, 17, x119);
-            Multiply(x119, x17, x119);
-            uint[] x125 = x102;
-            SquareN(x119, 6, x125);
-            Multiply(x125, x6, x125);
-            uint[] x127 = x6;
-            SquareN(x125, 2, x127);
-            Multiply(x127, x2, x127);
-
-            uint[] t = x127;
-            SquareN(t, 18, t);
-            Multiply(t, x17, t);
-            SquareN(t, 2, t);
-            Multiply(t, x1, t);
-            SquareN(t, 3, t);
-            Multiply(t, x2, t);
-            SquareN(t, 6, t);
-            Multiply(t, x3, t);
-            SquareN(t, 4, t);
-
-            // NOTE that x1 and z could be the same array
-            Multiply(x1, t, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/sec/SecP192K1Field.cs b/crypto/src/math/ec/custom/sec/SecP192K1Field.cs
index 30d53f7dc..23bd732bd 100644
--- a/crypto/src/math/ec/custom/sec/SecP192K1Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP192K1Field.cs
@@ -75,67 +75,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^192 - 2^32 - 2^12 - 2^8 - 2^7 - 2^6 - 2^3 - 3
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 159 1s } { 1 0s } { 19 1s } { 1 0s } { 3 1s } "000110101"
-             *
-             * Therefore we need an addition chain containing 1, 2, 3, 19, 159 (the lengths of the repunits)
-             * We use: [1], [2], [3], 6, 12, 18, [19], 38, 76, 152, 158, [159]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat192.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x3 = Nat192.Create();
-            Square(x2, x3);
-            Multiply(x3, x1, x3);
-            uint[] x6 = Nat192.Create();
-            SquareN(x3, 3, x6);
-            Multiply(x6, x3, x6);
-            uint[] x12 = Nat192.Create();
-            SquareN(x6, 6, x12);
-            Multiply(x12, x6, x12);
-            uint[] x18 = x12;
-            SquareN(x12, 6, x18);
-            Multiply(x18, x6, x18);
-            uint[] x19 = x18;
-            Square(x18, x19);
-            Multiply(x19, x1, x19);
-            uint[] x38 = Nat192.Create();
-            SquareN(x19, 19, x38);
-            Multiply(x38, x19, x38);
-            uint[] x76 = Nat192.Create();
-            SquareN(x38, 38, x76);
-            Multiply(x76, x38, x76);
-            uint[] x152 = x38;
-            SquareN(x76, 76, x152);
-            Multiply(x152, x76, x152);
-            uint[] x158 = x76;
-            SquareN(x152, 6, x158);
-            Multiply(x158, x6, x158);
-            uint[] x159 = x6;
-            Square(x158, x159);
-            Multiply(x159, x1, x159);
-
-            uint[] t = x159;
-            SquareN(t, 20, t);
-            Multiply(t, x19, t);
-            SquareN(t, 4, t);
-            Multiply(t, x3, t);
-            SquareN(t, 5, t);
-            Multiply(t, x2, t);
-            SquareN(t, 2, t);
-            Multiply(t, x1, t);
-            SquareN(t, 2, t);
-
-            // NOTE that x1 and z could be the same array
-            Multiply(x1, t, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/sec/SecP192R1Field.cs b/crypto/src/math/ec/custom/sec/SecP192R1Field.cs
index 2061d1359..a4fb4bb76 100644
--- a/crypto/src/math/ec/custom/sec/SecP192R1Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP192R1Field.cs
@@ -74,58 +74,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^192 - 2^64 - 1
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 127 1s } { 1 0s } { 62 1s } { 1 0s } { 1 1s }
-             *
-             * Therefore we need an addition chain containing 1, 62, 127 (the lengths of the repunits)
-             * We use: [1], 2, 3, 6, 12, 24, 30, 32, [62], 65, [127]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat192.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x3 = Nat192.Create();
-            Square(x2, x3);
-            Multiply(x3, x1, x3);
-            uint[] x6 = Nat192.Create();
-            SquareN(x3, 3, x6);
-            Multiply(x6, x3, x6);
-            uint[] x12 = Nat192.Create();
-            SquareN(x6, 6, x12);
-            Multiply(x12, x6, x12);
-            uint[] x24 = Nat192.Create();
-            SquareN(x12, 12, x24);
-            Multiply(x24, x12, x24);
-            uint[] x30 = x12;
-            SquareN(x24, 6, x30);
-            Multiply(x30, x6, x30);
-            uint[] x32 = x6;
-            SquareN(x30, 2, x32);
-            Multiply(x32, x2, x32);
-            uint[] x62 = x2;
-            SquareN(x32, 30, x62);
-            Multiply(x62, x30, x62);
-            uint[] x65 = x24;
-            SquareN(x62, 3, x65);
-            Multiply(x65, x3, x65);
-            uint[] x127 = x3;
-            SquareN(x65, 62, x127);
-            Multiply(x127, x62, x127);
-
-            uint[] t = x127;
-            SquareN(t, 63, t);
-            Multiply(t, x62, t);
-            SquareN(t, 2, t);
-
-            // NOTE that x1 and z could be the same array
-            Multiply(x1, t, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/sec/SecP224K1Field.cs b/crypto/src/math/ec/custom/sec/SecP224K1Field.cs
index d20ac63f3..5d4237708 100644
--- a/crypto/src/math/ec/custom/sec/SecP224K1Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP224K1Field.cs
@@ -76,67 +76,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^224 - 2^32 - 2^12 - 2^11 - 2^9 - 2^7 - 2^4 - 5
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 191 1s } { 1 0s } { 19 1s } "0010101101011"
-             *
-             * Therefore we need an addition chain containing 1, 2, 19, 191 (the lengths of the repunits)
-             * We use: [1], [2], 4, 5, 9, 10, [19], 38, 76, 152, 190 [191]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat224.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x4 = Nat224.Create();
-            SquareN(x2, 2, x4);
-            Multiply(x4, x2, x4);
-            uint[] x5 = Nat224.Create();
-            Square(x4, x5);
-            Multiply(x5, x1, x5);
-            uint[] x9 = x5;
-            SquareN(x5, 4, x9);
-            Multiply(x9, x4, x9);
-            uint[] x10 = x4;
-            Square(x9, x10);
-            Multiply(x10, x1, x10);
-            uint[] x19 = x10;
-            SquareN(x10, 9, x19);
-            Multiply(x19, x9, x19);
-            uint[] x38 = x9;
-            SquareN(x19, 19, x38);
-            Multiply(x38, x19, x38);
-            uint[] x76 = Nat224.Create();
-            SquareN(x38, 38, x76);
-            Multiply(x76, x38, x76);
-            uint[] x152 = Nat224.Create();
-            SquareN(x76, 76, x152);
-            Multiply(x152, x76, x152);
-            uint[] x190 = x76;
-            SquareN(x152, 38, x190);
-            Multiply(x190, x38, x190);
-            uint[] x191 = x38;
-            Square(x190, x191);
-            Multiply(x191, x1, x191);
-
-            uint[] t = x191;
-            SquareN(t, 20, t);
-            Multiply(t, x19, t);
-            SquareN(t, 3, t);
-            Multiply(t, x1, t);
-            SquareN(t, 2, t);
-            Multiply(t, x1, t);
-            SquareN(t, 3, t);
-            Multiply(t, x2, t);
-            SquareN(t, 2, t);
-            Multiply(t, x1, t);
-            SquareN(t, 3, t);
-            Multiply(t, x2, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/sec/SecP224R1Field.cs b/crypto/src/math/ec/custom/sec/SecP224R1Field.cs
index 06d451c2b..dde291d5e 100644
--- a/crypto/src/math/ec/custom/sec/SecP224R1Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP224R1Field.cs
@@ -75,54 +75,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^224 - 2^96 - 1
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 127 1s } { 1 0s } { 96 1s }
-             *
-             * Therefore we need an addition chain containing 96, 127 (the lengths of the repunits)
-             * We use: 1, 2, 3, 6, 12, 24, 48, [96], 120, 126, [127]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat224.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x3 = x2;
-            Square(x2, x3);
-            Multiply(x3, x1, x3);
-            uint[] x6 = Nat224.Create();
-            SquareN(x3, 3, x6);
-            Multiply(x6, x3, x6);
-            uint[] x12 = x3;
-            SquareN(x6, 6, x12);
-            Multiply(x12, x6, x12);
-            uint[] x24 = Nat224.Create();
-            SquareN(x12, 12, x24);
-            Multiply(x24, x12, x24);
-            uint[] x48 = x12;
-            SquareN(x24, 24, x48);
-            Multiply(x48, x24, x48);
-            uint[] x96 = Nat224.Create();
-            SquareN(x48, 48, x96);
-            Multiply(x96, x48, x96);
-            uint[] x120 = x48;
-            SquareN(x96, 24, x120);
-            Multiply(x120, x24, x120);
-            uint[] x126 = x24;
-            SquareN(x120, 6, x126);
-            Multiply(x126, x6, x126);
-            uint[] x127 = x6;
-            Square(x126, x127);
-            Multiply(x127, x1, x127);
-
-            uint[] t = x127;
-            SquareN(t, 97, t);
-            Multiply(t, x96, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/sec/SecP256K1Field.cs b/crypto/src/math/ec/custom/sec/SecP256K1Field.cs
index 2193c94e6..acdb1f362 100644
--- a/crypto/src/math/ec/custom/sec/SecP256K1Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP256K1Field.cs
@@ -76,65 +76,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^256 - 2^32 - 2^9 - 2^8 - 2^7 - 2^6 - 2^4 - 3
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 223 1s } { 1 0s } { 22 1s } { 4 0s } { 1 1s } { 1 0s } { 2 1s } { 1 0s } { 1 1s }
-             *
-             * Therefore we need an addition chain containing 1, 2, 22, 223 (the lengths of the repunits)
-             * We use: [1], [2], 3, 6, 9, 11, [22], 44, 88, 176, 220, [223]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat256.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x3 = Nat256.Create();
-            Square(x2, x3);
-            Multiply(x3, x1, x3);
-            uint[] x6 = Nat256.Create();
-            SquareN(x3, 3, x6);
-            Multiply(x6, x3, x6);
-            uint[] x9 = x6;
-            SquareN(x6, 3, x9);
-            Multiply(x9, x3, x9);
-            uint[] x11 = x9;
-            SquareN(x9, 2, x11);
-            Multiply(x11, x2, x11);
-            uint[] x22 = Nat256.Create();
-            SquareN(x11, 11, x22);
-            Multiply(x22, x11, x22);
-            uint[] x44 = x11;
-            SquareN(x22, 22, x44);
-            Multiply(x44, x22, x44);
-            uint[] x88 = Nat256.Create();
-            SquareN(x44, 44, x88);
-            Multiply(x88, x44, x88);
-            uint[] x176 = Nat256.Create();
-            SquareN(x88, 88, x176);
-            Multiply(x176, x88, x176);
-            uint[] x220 = x88;
-            SquareN(x176, 44, x220);
-            Multiply(x220, x44, x220);
-            uint[] x223 = x44;
-            SquareN(x220, 3, x223);
-            Multiply(x223, x3, x223);
-
-            uint[] t = x223;
-            SquareN(t, 23, t);
-            Multiply(t, x22, t);
-            SquareN(t, 5, t);
-            Multiply(t, x1, t);
-            SquareN(t, 3, t);
-            Multiply(t, x2, t);
-            SquareN(t, 2, t);
-
-            // NOTE that x1 and z could be the same array
-            Multiply(x1, t, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/sec/SecP256R1Field.cs b/crypto/src/math/ec/custom/sec/SecP256R1Field.cs
index eadc7ee58..668efc895 100644
--- a/crypto/src/math/ec/custom/sec/SecP256R1Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP256R1Field.cs
@@ -70,60 +70,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^256 - 2^224 + 2^192 + 2^96 - 3
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 32 1s } { 31 0s } { 1 1s } { 96 0s } { 94 1s } { 1 0s } { 1 1s }
-             *
-             * Therefore we need an addition chain containing 1, 32, 94 (the lengths of the repunits)
-             * We use: [1], 2, 4, 8, 16, [32], 64, 80, 88, 92, [94]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat256.Create();
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x4 = Nat256.Create();
-            SquareN(x2, 2, x4);
-            Multiply(x4, x2, x4);
-            uint[] x8 = Nat256.Create();
-            SquareN(x4, 4, x8);
-            Multiply(x8, x4, x8);
-            uint[] x16 = Nat256.Create();
-            SquareN(x8, 8, x16);
-            Multiply(x16, x8, x16);
-            uint[] x32 = Nat256.Create();
-            SquareN(x16, 16, x32);
-            Multiply(x32, x16, x32);
-            uint[] x64 = Nat256.Create();
-            SquareN(x32, 32, x64);
-            Multiply(x64, x32, x64);
-            uint[] x80 = x64;
-            SquareN(x64, 16, x80);
-            Multiply(x80, x16, x80);
-            uint[] x88 = x16;
-            SquareN(x80, 8, x88);
-            Multiply(x88, x8, x88);
-            uint[] x92 = x8;
-            SquareN(x88, 4, x92);
-            Multiply(x92, x4, x92);
-            uint[] x94 = x4;
-            SquareN(x92, 2, x94);
-            Multiply(x94, x2, x94);
-
-            uint[] t = x32;
-            SquareN(t, 32, t);
-            Multiply(t, x1, t);
-            SquareN(t, 190, t);
-            Multiply(t, x94, t);
-            SquareN(t, 2, t);
-
-            // NOTE that x1 and z could be the same array
-            Multiply(x1, t, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/sec/SecP384R1Field.cs b/crypto/src/math/ec/custom/sec/SecP384R1Field.cs
index 9b20db1b0..cddb46895 100644
--- a/crypto/src/math/ec/custom/sec/SecP384R1Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP384R1Field.cs
@@ -77,66 +77,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^384 - 2^128 - 2^96 + 2^32 - 3
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 255 1s } { 1 0s } { 32 1s } { 64 0s } { 30 1s } { 1 0s } { 1 1s }
-             *
-             * Therefore we need an addition chain containing 1, 30, 32, 255 (the lengths of the repunits)
-             * We use: [1], 2, 3, 6, 12, 24, [30], [32], 62, 124, 248, 254, [255]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat.Create(12);
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x3 = Nat.Create(12);
-            Square(x2, x3);
-            Multiply(x3, x1, x3);
-            uint[] x6 = Nat.Create(12);
-            SquareN(x3, 3, x6);
-            Multiply(x6, x3, x6);
-            uint[] x12 = x3;
-            SquareN(x6, 6, x12);
-            Multiply(x12, x6, x12);
-            uint[] x24 = Nat.Create(12);
-            SquareN(x12, 12, x24);
-            Multiply(x24, x12, x24);
-            uint[] x30 = x12;
-            SquareN(x24, 6, x30);
-            Multiply(x30, x6, x30);
-            uint[] x32 = x24;
-            SquareN(x30, 2, x32);
-            Multiply(x32, x2, x32);
-            uint[] x62 = x2;
-            SquareN(x32, 30, x62);
-            Multiply(x62, x30, x62);
-            uint[] x124 = Nat.Create(12);
-            SquareN(x62, 62, x124);
-            Multiply(x124, x62, x124);
-            uint[] x248 = x62;
-            SquareN(x124, 124, x248);
-            Multiply(x248, x124, x248);
-            uint[] x254 = x124;
-            SquareN(x248, 6, x254);
-            Multiply(x254, x6, x254);
-            uint[] x255 = x6;
-            Square(x254, x255);
-            Multiply(x255, x1, x255);
-
-            uint[] t = x255;
-            SquareN(t, 33, t);
-            Multiply(t, x32, t);
-            SquareN(t, 94, t);
-            Multiply(t, x30, t);
-            SquareN(t, 2, t);
-
-            // NOTE that x1 and z could be the same array
-            Multiply(x1, t, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/custom/sec/SecP521R1Field.cs b/crypto/src/math/ec/custom/sec/SecP521R1Field.cs
index 10b98fc21..0f1922f36 100644
--- a/crypto/src/math/ec/custom/sec/SecP521R1Field.cs
+++ b/crypto/src/math/ec/custom/sec/SecP521R1Field.cs
@@ -56,62 +56,7 @@ namespace Org.BouncyCastle.Math.EC.Custom.Sec
 
         public static void Inv(uint[] x, uint[] z)
         {
-            /*
-             * Raise this element to the exponent 2^521 - 3
-             *
-             * Breaking up the exponent's binary representation into "repunits", we get:
-             * { 519 1s } { 1 0s } { 1 1s }
-             *
-             * Therefore we need an addition chain containing 1, 519 (the lengths of the repunits)
-             * We use: [1], 2, 4, 8, 16, 32, 64, 128, 256, 512, 516, 518, [519]
-             */
-
-            if (0 != IsZero(x))
-                throw new ArgumentException("cannot be 0", "x");
-
-            uint[] x1 = x;
-            uint[] x2 = Nat.Create(17);
-            Square(x1, x2);
-            Multiply(x2, x1, x2);
-            uint[] x4 = Nat.Create(17);
-            SquareN(x2, 2, x4);
-            Multiply(x4, x2, x4);
-            uint[] x8 = Nat.Create(17);
-            SquareN(x4, 4, x8);
-            Multiply(x8, x4, x8);
-            uint[] x16 = Nat.Create(17);
-            SquareN(x8, 8, x16);
-            Multiply(x16, x8, x16);
-            uint[] x32 = x8;
-            SquareN(x16, 16, x32);
-            Multiply(x32, x16, x32);
-            uint[] x64 = x16;
-            SquareN(x32, 32, x64);
-            Multiply(x64, x32, x64);
-            uint[] x128 = x32;
-            SquareN(x64, 64, x128);
-            Multiply(x128, x64, x128);
-            uint[] x256 = x64;
-            SquareN(x128, 128, x256);
-            Multiply(x256, x128, x256);
-            uint[] x512 = x128;
-            SquareN(x256, 256, x512);
-            Multiply(x512, x256, x512);
-            uint[] x516 = x256;
-            SquareN(x512, 4, x516);
-            Multiply(x516, x4, x516);
-            uint[] x518 = x4;
-            SquareN(x516, 2, x518);
-            Multiply(x518, x2, x518);
-            uint[] x519 = x2;
-            Square(x518, x519);
-            Multiply(x519, x1, x519);
-
-            uint[] t = x519;
-            SquareN(t, 2, t);
-
-            // NOTE that x1 and z could be the same array
-            Multiply(x1, t, z);
+            Mod.CheckedModOddInverse(P, x, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/rfc7748/X25519Field.cs b/crypto/src/math/ec/rfc7748/X25519Field.cs
index 6843e274a..ffede563b 100644
--- a/crypto/src/math/ec/rfc7748/X25519Field.cs
+++ b/crypto/src/math/ec/rfc7748/X25519Field.cs
@@ -1,6 +1,8 @@
 using System;
 using System.Diagnostics;
 
+using Org.BouncyCastle.Math.Raw;
+
 namespace Org.BouncyCastle.Math.EC.Rfc7748
 {
     public abstract class X25519Field
@@ -11,6 +13,8 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
         private const int M25 = 0x01FFFFFF;
         private const int M26 = 0x03FFFFFF;
 
+        private static readonly uint[] P32 = new uint[]{ 0xFFFFFFEDU, 0xFFFFFFFFU, 0xFFFFFFFFU, 0xFFFFFFFFU,
+            0xFFFFFFFFU, 0xFFFFFFFFU, 0xFFFFFFFFU, 0x7FFFFFFFU };
         private static readonly int[] RootNegOne = { 0x020EA0B0, 0x0386C9D2, 0x00478C4E, 0x0035697F, 0x005E8630,
             0x01FBD7A7, 0x0340264F, 0x01F0B2B4, 0x00027E0E, 0x00570649 };
 
@@ -128,6 +132,14 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
             }
         }
 
+        [CLSCompliantAttribute(false)]
+        public static void Decode(uint[] x, int xOff, int[] z)
+        {
+            Decode128(x, xOff, z, 0);
+            Decode128(x, xOff + 4, z, 5);
+            z[9] &= M24;
+        }
+
         public static void Decode(byte[] x, int xOff, int[] z)
         {
             Decode128(x, xOff, z, 0);
@@ -135,6 +147,17 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
             z[9] &= M24;
         }
 
+        private static void Decode128(uint[] x, int xOff, int[] z, int zOff)
+        {
+            uint t0 = x[xOff + 0], t1 = x[xOff + 1], t2 = x[xOff + 2], t3 = x[xOff + 3];
+
+            z[zOff + 0] = (int)t0 & M26;
+            z[zOff + 1] = (int)((t1 <<  6) | (t0 >> 26)) & M26;
+            z[zOff + 2] = (int)((t2 << 12) | (t1 >> 20)) & M25;
+            z[zOff + 3] = (int)((t3 << 19) | (t2 >> 13)) & M26;
+            z[zOff + 4] = (int)(t3 >> 7);
+        }
+
         private static void Decode128(byte[] bs, int off, int[] z, int zOff)
         {
             uint t0 = Decode32(bs, off + 0);
@@ -158,12 +181,30 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
             return n;
         }
 
+        [CLSCompliantAttribute(false)]
+        public static void Encode(int[] x, uint[] z, int zOff)
+        {
+            Encode128(x, 0, z, zOff);
+            Encode128(x, 5, z, zOff + 4);
+        }
+
         public static void Encode(int[] x, byte[] z, int zOff)
         {
             Encode128(x, 0, z, zOff);
             Encode128(x, 5, z, zOff + 16);
         }
 
+        private static void Encode128(int[] x, int xOff, uint[] z, int zOff)
+        {
+            uint x0 = (uint)x[xOff + 0], x1 = (uint)x[xOff + 1], x2 = (uint)x[xOff + 2], x3 = (uint)x[xOff + 3],
+                x4 = (uint)x[xOff + 4];
+
+            z[zOff + 0] =  x0        | (x1 << 26);
+            z[zOff + 1] = (x1 >>  6) | (x2 << 20);
+            z[zOff + 2] = (x2 >> 12) | (x3 << 13);
+            z[zOff + 3] = (x3 >> 19) | (x4 <<  7);
+        }
+
         private static void Encode128(int[] x, int xOff, byte[] bs, int off)
         {
             uint x0 = (uint)x[xOff + 0], x1 = (uint)x[xOff + 1], x2 = (uint)x[xOff + 2];
@@ -185,15 +226,36 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
 
         public static void Inv(int[] x, int[] z)
         {
-            // z = x^(p-2) = x^7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEB
-            // (250 1s) (1 0s) (1 1s) (1 0s) (2 1s)
-            // Addition chain: [1] [2] 3 5 10 15 25 50 75 125 [250]
+            //int[] x2 = Create();
+            //int[] t = Create();
+            //PowPm5d8(x, x2, t);
+            //Sqr(t, 3, t);
+            //Mul(t, x2, z);
 
-            int[] x2 = Create();
             int[] t = Create();
-            PowPm5d8(x, x2, t);
-            Sqr(t, 3, t);
-            Mul(t, x2, z);
+            uint[] u = new uint[8];
+
+            Copy(x, 0, t, 0);
+            Normalize(t);
+            Encode(t, u, 0);
+
+            Mod.ModOddInverse(P32, u, u);
+
+            Decode(u, 0, z);
+        }
+
+        public static void InvVar(int[] x, int[] z)
+        {
+            int[] t = Create();
+            uint[] u = new uint[8];
+
+            Copy(x, 0, t, 0);
+            Normalize(t);
+            Encode(t, u, 0);
+
+            Mod.ModOddInverseVar(P32, u, u);
+
+            Decode(u, 0, z);
         }
 
         public static int IsZero(int[] x)
diff --git a/crypto/src/math/ec/rfc7748/X448Field.cs b/crypto/src/math/ec/rfc7748/X448Field.cs
index 240518cde..4d3be5cda 100644
--- a/crypto/src/math/ec/rfc7748/X448Field.cs
+++ b/crypto/src/math/ec/rfc7748/X448Field.cs
@@ -1,6 +1,8 @@
 using System;
 using System.Diagnostics;
 
+using Org.BouncyCastle.Math.Raw;
+
 namespace Org.BouncyCastle.Math.EC.Rfc7748
 {
     [CLSCompliantAttribute(false)]
@@ -10,6 +12,10 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
 
         private const uint M28 = 0x0FFFFFFFU;
 
+        private static readonly uint[] P32 = new uint[]{ 0xFFFFFFFFU, 0xFFFFFFFFU, 0xFFFFFFFFU, 0xFFFFFFFFU,
+            0xFFFFFFFFU, 0xFFFFFFFFU, 0xFFFFFFFFU, 0xFFFFFFFEU, 0xFFFFFFFFU, 0xFFFFFFFFU, 0xFFFFFFFFU, 0xFFFFFFFFU,
+            0xFFFFFFFFU, 0xFFFFFFFFU };
+
         protected X448Field() {}
 
         public static void Add(uint[] x, uint[] y, uint[] z)
@@ -129,6 +135,12 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
             }
         }
 
+        public static void Decode(uint[] x, int xOff, uint[] z)
+        {
+            Decode224(x, xOff, z, 0);
+            Decode224(x, xOff + 7, z, 8);
+        }
+
         public static void Decode(byte[] x, int xOff, uint[] z)
         {
             Decode56(x, xOff, z, 0);
@@ -141,6 +153,21 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
             Decode56(x, xOff + 49, z, 14);
         }
 
+        private static void Decode224(uint[] x, int xOff, uint[] z, int zOff)
+        {
+            uint x0 = x[xOff + 0], x1 = x[xOff + 1], x2 = x[xOff + 2], x3 = x[xOff + 3];
+            uint x4 = x[xOff + 4], x5 = x[xOff + 5], x6 = x[xOff + 6];
+
+            z[zOff + 0] = x0 & M28;
+            z[zOff + 1] = (x0 >> 28 | x1 <<  4) & M28;
+            z[zOff + 2] = (x1 >> 24 | x2 <<  8) & M28;
+            z[zOff + 3] = (x2 >> 20 | x3 << 12) & M28;
+            z[zOff + 4] = (x3 >> 16 | x4 << 16) & M28;
+            z[zOff + 5] = (x4 >> 12 | x5 << 20) & M28;
+            z[zOff + 6] = (x5 >>  8 | x6 << 24) & M28;
+            z[zOff + 7] = x6 >> 4;
+        }
+
         private static uint Decode24(byte[] bs, int off)
         {
             uint n = bs[off];
@@ -166,6 +193,12 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
             z[zOff + 1] = (lo >> 28) | (hi << 4);
         }
 
+        public static void Encode(uint[] x, uint[] z, int zOff)
+        {
+            Encode224(x, 0, z, zOff);
+            Encode224(x, 8, z, zOff + 7);
+        }
+
         public static void Encode(uint[] x, byte[] z, int zOff)
         {
             Encode56(x, 0, z, zOff);
@@ -178,6 +211,20 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
             Encode56(x, 14, z, zOff + 49);
         }
 
+        private static void Encode224(uint[] x, int xOff, uint[] z, int zOff)
+        {
+            uint x0 = x[xOff + 0], x1 = x[xOff + 1], x2 = x[xOff + 2], x3 = x[xOff + 3];
+            uint x4 = x[xOff + 4], x5 = x[xOff + 5], x6 = x[xOff + 6], x7 = x[xOff + 7];
+
+            z[zOff + 0] =  x0        | (x1 << 28);
+            z[zOff + 1] = (x1 >>  4) | (x2 << 24);
+            z[zOff + 2] = (x2 >>  8) | (x3 << 20);
+            z[zOff + 3] = (x3 >> 12) | (x4 << 16);
+            z[zOff + 4] = (x4 >> 16) | (x5 << 12);
+            z[zOff + 5] = (x5 >> 20) | (x6 <<  8);
+            z[zOff + 6] = (x6 >> 24) | (x7 <<  4);
+        }
+
         private static void Encode24(uint n, byte[] bs, int off)
         {
             bs[  off] = (byte)(n      );
@@ -202,14 +249,35 @@ namespace Org.BouncyCastle.Math.EC.Rfc7748
 
         public static void Inv(uint[] x, uint[] z)
         {
-            // z = x^(p-2) = x^(2^448 - 2^224 - 3)
-            // (223 1s) (1 0s) (222 1s) (1 0s) (1 1s)
-            // Addition chain: [1] 2 3 6 9 18 19 37 74 111 [222] [223]
+            //uint[] t = Create();
+            //PowPm3d4(x, t);
+            //Sqr(t, 2, t);
+            //Mul(t, x, z);
 
             uint[] t = Create();
-            PowPm3d4(x, t);
-            Sqr(t, 2, t);
-            Mul(t, x, z);
+            uint[] u = new uint[14];
+
+            Copy(x, 0, t, 0);
+            Normalize(t);
+            Encode(t, u, 0);
+
+            Mod.ModOddInverse(P32, u, u);
+
+            Decode(u, 0, z);
+        }
+
+        public static void InvVar(uint[] x, uint[] z)
+        {
+            uint[] t = Create();
+            uint[] u = new uint[14];
+
+            Copy(x, 0, t, 0);
+            Normalize(t);
+            Encode(t, u, 0);
+
+            Mod.ModOddInverseVar(P32, u, u);
+
+            Decode(u, 0, z);
         }
 
         public static int IsZero(uint[] x)
diff --git a/crypto/src/math/ec/rfc8032/Ed25519.cs b/crypto/src/math/ec/rfc8032/Ed25519.cs
index 95ba43472..3a39ae53d 100644
--- a/crypto/src/math/ec/rfc8032/Ed25519.cs
+++ b/crypto/src/math/ec/rfc8032/Ed25519.cs
@@ -754,7 +754,7 @@ namespace Org.BouncyCastle.Math.EC.Rfc8032
 
                         X25519Field.Add(q.z, q.z, x);
                         // TODO[ed25519] Batch inversion
-                        X25519Field.Inv(x, y);
+                        X25519Field.InvVar(x, y);
                         X25519Field.Mul(q.x, y, x);
                         X25519Field.Mul(q.y, y, y);
 
diff --git a/crypto/src/math/ec/rfc8032/Ed448.cs b/crypto/src/math/ec/rfc8032/Ed448.cs
index 12f24c1ff..c1202dc02 100644
--- a/crypto/src/math/ec/rfc8032/Ed448.cs
+++ b/crypto/src/math/ec/rfc8032/Ed448.cs
@@ -702,7 +702,7 @@ namespace Org.BouncyCastle.Math.EC.Rfc8032
                     {
                         PointExt q = points[i];
                         // TODO[ed448] Batch inversion
-                        X448Field.Inv(q.z, q.z);
+                        X448Field.InvVar(q.z, q.z);
                         X448Field.Mul(q.x, q.z, q.x);
                         X448Field.Mul(q.y, q.z, q.y);
 
diff --git a/crypto/src/math/raw/Mod.cs b/crypto/src/math/raw/Mod.cs
index 197b5c82b..1a1524c11 100644
--- a/crypto/src/math/raw/Mod.cs
+++ b/crypto/src/math/raw/Mod.cs
@@ -7,74 +7,190 @@ using Org.BouncyCastle.Utilities;
 
 namespace Org.BouncyCastle.Math.Raw
 {
+    /*
+     * Modular inversion as implemented in this class is based on the paper "Fast constant-time gcd
+     * computation and modular inversion" by Daniel J. Bernstein and Bo-Yin Yang.
+     */
+
     internal abstract class Mod
     {
         private static readonly SecureRandom RandomSource = new SecureRandom();
 
-        public static void Invert(uint[] p, uint[] x, uint[] z)
+        private const int M30 = 0x3FFFFFFF;
+        private const ulong M32UL = 0xFFFFFFFFUL;
+
+        [Obsolete("Will be removed")]
+        public static void Add(uint[] p, uint[] x, uint[] y, uint[] z)
         {
             int len = p.Length;
-            if (Nat.IsZero(len, x))
-                throw new ArgumentException("cannot be 0", "x");
-            if (Nat.IsOne(len, x))
+            uint c = Nat.Add(len, x, y, z);
+            if (c != 0)
             {
-                Array.Copy(x, 0, z, 0, len);
-                return;
+                Nat.SubFrom(len, p, z);
             }
+        }
 
-            uint[] u = Nat.Copy(len, x);
-            uint[] a = Nat.Create(len);
-            a[0] = 1;
-            int ac = 0;
+        public static void CheckedModOddInverse(uint[] m, uint[] x, uint[] z)
+        {
+            if (0 == ModOddInverse(m, x, z))
+                throw new ArithmeticException("Inverse does not exist.");
+        }
 
-            if ((u[0] & 1) == 0)
-            {
-                InversionStep(p, u, len, a, ref ac);
-            }
-            if (Nat.IsOne(len, u))
+        public static void CheckedModOddInverseVar(uint[] m, uint[] x, uint[] z)
+        {
+            if (!ModOddInverseVar(m, x, z))
+                throw new ArithmeticException("Inverse does not exist.");
+        }
+
+        public static uint Inverse32(uint d)
+        {
+            Debug.Assert((d & 1) == 1);
+
+            //int x = d + (((d + 1) & 4) << 1);   // d.x == 1 mod 2**4
+            uint x = d;                         // d.x == 1 mod 2**3
+            x *= 2 - d * x;                     // d.x == 1 mod 2**6
+            x *= 2 - d * x;                     // d.x == 1 mod 2**12
+            x *= 2 - d * x;                     // d.x == 1 mod 2**24
+            x *= 2 - d * x;                     // d.x == 1 mod 2**48
+            Debug.Assert(d * x == 1);
+            return x;
+        }
+
+        [Obsolete("Use 'CheckedModOddInverseVar' instead")]
+        public static void Invert(uint[] m, uint[] x, uint[] z)
+        {
+            CheckedModOddInverseVar(m, x, z);
+        }
+
+        public static uint ModOddInverse(uint[] m, uint[] x, uint[] z)
+        {
+            int len32 = m.Length;
+            Debug.Assert(len32 > 0);
+            Debug.Assert((m[0] & 1) != 0);
+            Debug.Assert(m[len32 - 1] != 0);
+
+            int bits = (len32 << 5) - Integers.NumberOfLeadingZeros((int)m[len32 - 1]);
+            int len30 = (bits + 29) / 30;
+            int m0Inv30x4 = -(int)Inverse32(m[0]) << 2;
+
+            int[] t = new int[4];
+            int[] D = new int[len30];
+            int[] E = new int[len30];
+            int[] F = new int[len30];
+            int[] G = new int[len30];
+            int[] M = new int[len30];
+
+            E[0] = 1;
+            Encode30(bits, x, 0, G, 0);
+            Encode30(bits, m, 0, M, 0);
+            Array.Copy(M, 0, F, 0, len30);
+
+            int eta = -1;
+            int maxDivsteps = GetMaximumDivsteps(bits);
+
+            for (int divSteps = 0; divSteps < maxDivsteps; divSteps += 30)
             {
-                InversionResult(p, ac, a, z);
-                return;
+                eta = Divsteps30(eta, F[0], G[0], t);
+                UpdateDE30(len30, D, E, t, m0Inv30x4, M);
+                UpdateFG30(len30, F, G, t);
             }
 
-            uint[] v = Nat.Copy(len, p);
-            uint[] b = Nat.Create(len);
-            int bc = 0;
+            int signF = F[len30 - 1] >> 31;
+            Debug.Assert(-1 == signF | 0 == signF);
 
-            int uvLen = len;
+            CNegate30(len30, signF, F);
+            CNegate30(len30, signF, D);
 
-            for (;;)
+            Decode30(bits, D, 0, z, 0);
+
+            int signD = D[len30 - 1] >> 31;
+            Debug.Assert(-1 == signD | 0 == signD);
+
+            signD += (int)Nat.CAdd(len32, signD, z, m, z);
+            Debug.Assert(0 == signD & 0 != Nat.LessThan(len32, z, m));
+
+            return (uint)(EqualTo(len30, F, 1) & EqualToZero(len30, G));
+        }
+
+        public static bool ModOddInverseVar(uint[] m, uint[] x, uint[] z)
+        {
+            int len32 = m.Length;
+            Debug.Assert(len32 > 0);
+            Debug.Assert((m[0] & 1) != 0);
+            Debug.Assert(m[len32 - 1] != 0);
+
+            int bits = (len32 << 5) - Integers.NumberOfLeadingZeros((int)m[len32 - 1]);
+            int len30 = (bits + 29) / 30;
+            int m0Inv30x4 = -(int)Inverse32(m[0]) << 2;
+
+            int[] t = new int[4];
+            int[] D = new int[len30];
+            int[] E = new int[len30];
+            int[] F = new int[len30];
+            int[] G = new int[len30];
+            int[] M = new int[len30];
+
+            E[0] = 1;
+            Encode30(bits, x, 0, G, 0);
+            Encode30(bits, m, 0, M, 0);
+            Array.Copy(M, 0, F, 0, len30);
+
+            int clzG = Integers.NumberOfLeadingZeros(G[len30 - 1] | 1) - (len30 * 30 + 2 - bits);
+            int eta = -1 - clzG;
+            int lenDE = len30, lenFG = len30;
+            int maxDivsteps = GetMaximumDivsteps(bits);
+
+            int divsteps = 0;
+            while (!IsZero(lenFG, G))
             {
-                while (u[uvLen - 1] == 0 && v[uvLen - 1] == 0)
-                {
-                    --uvLen;
-                }
+                if (divsteps >= maxDivsteps)
+                    return false;
 
-                if (Nat.Gte(len, u, v))
-                {
-                    Nat.SubFrom(len, v, u);
-                    Debug.Assert((u[0] & 1) == 0);
-                    ac += Nat.SubFrom(len, b, a) - bc;
-                    InversionStep(p, u, uvLen, a, ref ac);
-                    if (Nat.IsOne(len, u))
-                    {
-                        InversionResult(p, ac, a, z);
-                        return;
-                    }
-                }
-                else
+                divsteps += 30;
+
+                eta = Divsteps30Var(eta, F[0], G[0], t);
+                UpdateDE30(lenDE, D, E, t, m0Inv30x4, M);
+                UpdateFG30(lenFG, F, G, t);
+
+                int fn = F[lenFG - 1];
+                int gn = G[lenFG - 1];
+
+                int cond = (lenFG - 2) >> 31;
+                cond |= fn ^ (fn >> 31);
+                cond |= gn ^ (gn >> 31);
+
+                if (cond == 0)
                 {
-                    Nat.SubFrom(len, u, v);
-                    Debug.Assert((v[0] & 1) == 0);
-                    bc += Nat.SubFrom(len, a, b) - ac;
-                    InversionStep(p, v, uvLen, b, ref bc);
-                    if (Nat.IsOne(len, v))
-                    {
-                        InversionResult(p, bc, b, z);
-                        return;
-                    }
+                    F[lenFG - 2] |= fn << 30;
+                    G[lenFG - 2] |= gn << 30;
+                    --lenFG;
                 }
             }
+
+            int signF = F[lenFG - 1] >> 31;
+            Debug.Assert(-1 == signF | 0 == signF);
+
+            if (0 != signF)
+            {
+                Negate30(lenFG, F);
+                Negate30(lenDE, D);
+            }
+
+            if (!IsOne(lenFG, F))
+                return false;
+
+            Decode30(bits, D, 0, z, 0);
+
+            int signD = D[lenDE - 1] >> 31;
+            Debug.Assert(-1 == signD | 0 == signD);
+
+            if (signD < 0)
+            {
+                signD += (int)Nat.AddTo(len32, m, z);
+            }
+            Debug.Assert(0 == signD && !Nat.Gte(len32, z, m));
+
+            return true;
         }
 
         public static uint[] Random(uint[] p)
@@ -101,74 +217,346 @@ namespace Org.BouncyCastle.Math.Raw
             return s;
         }
 
-        public static void Add(uint[] p, uint[] x, uint[] y, uint[] z)
+        [Obsolete("Will be removed")]
+        public static void Subtract(uint[] p, uint[] x, uint[] y, uint[] z)
         {
             int len = p.Length;
-            uint c = Nat.Add(len, x, y, z);
+            int c = Nat.Sub(len, x, y, z);
             if (c != 0)
             {
-                Nat.SubFrom(len, p, z);
+                Nat.AddTo(len, p, z);
             }
         }
 
-        public static void Subtract(uint[] p, uint[] x, uint[] y, uint[] z)
+        private static void CNegate30(int len, int cond, int[] D)
         {
-            int len = p.Length;
-            int c = Nat.Sub(len, x, y, z);
-            if (c != 0)
+            Debug.Assert(len > 0);
+            Debug.Assert(D.Length >= len);
+
+            int last = len - 1;
+            long cd = 0L;
+
+            for (int i = 0; i < last; ++i)
             {
-                Nat.AddTo(len, p, z);
+                cd += (D[i] ^ cond) - cond;
+                D[i] = (int)cd & M30; cd >>= 30;
             }
+
+            cd += (D[last] ^ cond) - cond;
+            D[last] = (int)cd;
         }
 
-        private static void InversionResult(uint[] p, int ac, uint[] a, uint[] z)
+        private static void Decode30(int bits, int[] x, int xOff, uint[] z, int zOff)
         {
-            if (ac < 0)
+            Debug.Assert(bits > 0);
+
+            int avail = 0;
+            ulong data = 0L;
+
+            while (bits > 0)
             {
-                Nat.Add(p.Length, a, p, z);
+                while (avail < System.Math.Min(32, bits))
+                {
+                    data |= (ulong)x[xOff++] << avail;
+                    avail += 30;
+                }
+
+                z[zOff++] = (uint)data; data >>= 32;
+                avail -= 32;
+                bits -= 32;
             }
-            else
+        }
+
+        private static int Divsteps30(int eta, int f0, int g0, int[] t)
+        {
+            int u = 1, v = 0, q = 0, r = 1;
+            int f = f0, g = g0;
+
+            for (int i = 0; i < 30; ++i)
             {
-                Array.Copy(a, 0, z, 0, p.Length);
+                Debug.Assert((f & 1) == 1);
+                Debug.Assert((u * f0 + v * g0) == f << i);
+                Debug.Assert((q * f0 + r * g0) == g << i);
+
+                int p = -(g & 1);
+                int s = eta >> 31;
+
+                int c1 = p & s;
+                int c2 = p & ~s;
+
+                eta = (eta ^ c1) - (c1 + 1);
+
+                g += f & c2;
+                q += u & c2;
+                r += v & c2;
+
+                g -= f & c1;
+                q -= u & c1;
+                r -= v & c1;
+
+                f += g & c1;
+                u += q & c1;
+                v += r & c1;
+
+                g >>= 1;
+                u <<= 1;
+                v <<= 1;
             }
+
+            t[0] = u;
+            t[1] = v;
+            t[2] = q;
+            t[3] = r;
+
+            return eta;
         }
 
-        private static void InversionStep(uint[] p, uint[] u, int uLen, uint[] x, ref int xc)
+        private static int Divsteps30Var(int eta, int f0, int g0, int[] t)
         {
-            int len = p.Length;
-            int count = 0;
-            while (u[0] == 0)
+            int u = 1, v = 0, q = 0, r = 1;
+            int f = f0, g = g0, m, w, x, y, z;
+            int i = 30, limit, zeros;
+
+            for (;;)
             {
-                Nat.ShiftDownWord(uLen, u, 0);
-                count += 32;
+                // Use a sentinel bit to count zeros only up to i.
+                zeros = Integers.NumberOfTrailingZeros(g | (-1 << i));
+
+                g >>= zeros;
+                u <<= zeros;
+                v <<= zeros;
+                eta -= zeros;
+                i -= zeros;
+
+                if (i <= 0)
+                    break;
+
+                Debug.Assert((f & 1) == 1);
+                Debug.Assert((g & 1) == 1);
+                Debug.Assert((u * f0 + v * g0) == f << (30 - i));
+                Debug.Assert((q * f0 + r * g0) == g << (30 - i));
+
+                if (eta < 0)
+                {
+                    eta = -eta;
+                    x = f; f = g; g = -x;
+                    y = u; u = q; q = -y;
+                    z = v; v = r; r = -z;
+
+                    // Handle up to 6 divsteps at once, subject to eta and i.
+                    limit = (eta + 1) > i ? i : (eta + 1);
+                    m = (int)((uint.MaxValue >> (32 - limit)) & 63U);
+
+                    w = (f * g * (f * f - 2)) & m;
+                }
+                else
+                {
+                    // Handle up to 4 divsteps at once, subject to eta and i.
+                    limit = (eta + 1) > i ? i : (eta + 1);
+                    m = (int)((uint.MaxValue >> (32 - limit)) & 15U);
+
+                    w = f + (((f + 1) & 4) << 1);
+                    w = (-w * g) & m;
+                }
+
+                g += f * w;
+                q += u * w;
+                r += v * w;
+
+                Debug.Assert((g & m) == 0);
             }
 
+            t[0] = u;
+            t[1] = v;
+            t[2] = q;
+            t[3] = r;
+
+            return eta;
+        }
+
+        private static void Encode30(int bits, uint[] x, int xOff, int[] z, int zOff)
+        {
+            Debug.Assert(bits > 0);
+
+            int avail = 0;
+            ulong data = 0UL;
+
+            while (bits > 0)
             {
-                int zeroes = Integers.NumberOfTrailingZeros((int)u[0]);
-                if (zeroes > 0)
+                if (avail < System.Math.Min(30, bits))
                 {
-                    Nat.ShiftDownBits(uLen, u, zeroes, 0);
-                    count += zeroes;
+                    data |= (x[xOff++] & M32UL) << avail;
+                    avail += 32;
                 }
+
+                z[zOff++] = (int)data & M30; data >>= 30;
+                avail -= 30;
+                bits -= 30;
+            }
+        }
+
+        private static int EqualTo(int len, int[] x, int y)
+        {
+            int d = x[0] ^ y;
+            for (int i = 1; i < len; ++i)
+            {
+                d |= x[i];
+            }
+            d = (int)((uint)d >> 1) | (d & 1);
+            return (d - 1) >> 31;
+        }
+
+        private static int EqualToZero(int len, int[] x)
+        {
+            int d = 0;
+            for (int i = 0; i < len; ++i)
+            {
+                d |= x[i];
             }
+            d = (int)((uint)d >> 1) | (d & 1);
+            return (d - 1) >> 31;
+        }
+
+        private static int GetMaximumDivsteps(int bits)
+        {
+            return (49 * bits + (bits < 46 ? 80 : 47)) / 17;
+        }
 
-            for (int i = 0; i < count; ++i)
+        private static bool IsOne(int len, int[] x)
+        {
+            if (x[0] != 1)
             {
-                if ((x[0] & 1) != 0)
+                return false;
+            }
+            for (int i = 1; i < len; ++i)
+            {
+                if (x[i] != 0)
                 {
-                    if (xc < 0)
-                    {
-                        xc += (int)Nat.AddTo(len, p, x);
-                    }
-                    else
-                    {
-                        xc += Nat.SubFrom(len, p, x);
-                    }
+                    return false;
                 }
+            }
+            return true;
+        }
 
-                Debug.Assert(xc == 0 || xc == -1);
-                Nat.ShiftDownBit(len, x, (uint)xc);
+        private static bool IsZero(int len, int[] x)
+        {
+            if (x[0] != 0)
+            {
+                return false;
+            }
+            for (int i = 1; i < len; ++i)
+            {
+                if (x[i] != 0)
+                {
+                    return false;
+                }
+            }
+            return true;
+        }
+
+        private static void Negate30(int len, int[] D)
+        {
+            Debug.Assert(len > 0);
+            Debug.Assert(D.Length >= len);
+
+            int last = len - 1;
+            long cd = 0L;
+
+            for (int i = 0; i < last; ++i)
+            {
+                cd -= D[i];
+                D[i] = (int)cd & M30; cd >>= 30;
+            }
+
+            cd -= D[last];
+            D[last] = (int)cd;
+        }
+
+        private static void UpdateDE30(int len, int[] D, int[] E, int[] t, int m0Inv30x4, int[] M)
+        {
+            Debug.Assert(len > 0);
+            Debug.Assert(D.Length >= len);
+            Debug.Assert(E.Length >= len);
+            Debug.Assert(M.Length >= len);
+            Debug.Assert(m0Inv30x4 * M[0] == -1 << 2);
+
+            int u = t[0], v = t[1], q = t[2], r = t[3];
+            int di, ei, i, md, me;
+            long cd, ce;
+
+            di = D[0];
+            ei = E[0];
+
+            cd = (long)u * di + (long)v * ei;
+            ce = (long)q * di + (long)r * ei;
+
+            md = (m0Inv30x4 * (int)cd) >> 2;
+            me = (m0Inv30x4 * (int)ce) >> 2;
+
+            cd += (long)M[0] * md;
+            ce += (long)M[0] * me;
+
+            Debug.Assert(((int)cd & M30) == 0);
+            Debug.Assert(((int)ce & M30) == 0);
+
+            cd >>= 30;
+            ce >>= 30;
+
+            for (i = 1; i < len; ++i)
+            {
+                di = D[i];
+                ei = E[i];
+
+                cd += (long)u * di + (long)v * ei;
+                ce += (long)q * di + (long)r * ei;
+
+                cd += (long)M[i] * md;
+                ce += (long)M[i] * me;
+
+                D[i - 1] = (int)cd & M30; cd >>= 30;
+                E[i - 1] = (int)ce & M30; ce >>= 30;
+            }
+
+            D[len - 1] = (int)cd;
+            E[len - 1] = (int)ce;
+        }
+
+        private static void UpdateFG30(int len, int[] F, int[] G, int[] t)
+        {
+            Debug.Assert(len > 0);
+            Debug.Assert(F.Length >= len);
+            Debug.Assert(G.Length >= len);
+
+            int u = t[0], v = t[1], q = t[2], r = t[3];
+            int fi, gi, i;
+            long cf, cg;
+
+            fi = F[0];
+            gi = G[0];
+
+            cf = (long)u * fi + (long)v * gi;
+            cg = (long)q * fi + (long)r * gi;
+
+            Debug.Assert(((int)cf & M30) == 0);
+            Debug.Assert(((int)cg & M30) == 0);
+
+            cf >>= 30;
+            cg >>= 30;
+
+            for (i = 1; i < len; ++i)
+            {
+                fi = F[i];
+                gi = G[i];
+
+                cf += (long)u * fi + (long)v * gi;
+                cg += (long)q * fi + (long)r * gi;
+
+                F[i - 1] = (int)cf & M30; cf >>= 30;
+                G[i - 1] = (int)cg & M30; cg >>= 30;
             }
+
+            F[len - 1] = (int)cf;
+            G[len - 1] = (int)cg;
         }
     }
 }
diff --git a/crypto/src/math/raw/Nat.cs b/crypto/src/math/raw/Nat.cs
index 69942661f..9f2412580 100644
--- a/crypto/src/math/raw/Nat.cs
+++ b/crypto/src/math/raw/Nat.cs
@@ -447,6 +447,72 @@ namespace Org.BouncyCastle.Math.Raw
             return true;
         }
 
+        public static uint EqualTo(int len, uint[] x, uint y)
+        {
+            uint d = x[0] ^ y;
+            for (int i = 1; i < len; ++i)
+            {
+                d |= x[i];
+            }
+            d = (d >> 1) | (d & 1);
+            return (uint)(((int)d - 1) >> 31);
+        }
+
+        public static uint EqualTo(int len, uint[] x, int xOff, uint y)
+        {
+            uint d = x[xOff] ^ y;
+            for (int i = 1; i < len; ++i)
+            {
+                d |= x[xOff + i];
+            }
+            d = (d >> 1) | (d & 1);
+            return (uint)(((int)d - 1) >> 31);
+        }
+
+        public static uint EqualTo(int len, uint[] x, uint[] y)
+        {
+            uint d = 0;
+            for (int i = 0; i < len; ++i)
+            {
+                d |= x[i] ^ y[i];
+            }
+            d = (d >> 1) | (d & 1);
+            return (uint)(((int)d - 1) >> 31);
+        }
+
+        public static uint EqualTo(int len, uint[] x, int xOff, uint[] y, int yOff)
+        {
+            uint d = 0;
+            for (int i = 0; i < len; ++i)
+            {
+                d |= x[xOff + i] ^ y[yOff + i];
+            }
+            d = (d >> 1) | (d & 1);
+            return (uint)(((int)d - 1) >> 31);
+        }
+
+        public static uint EqualToZero(int len, uint[] x)
+        {
+            uint d = 0;
+            for (int i = 0; i < len; ++i)
+            {
+                d |= x[i];
+            }
+            d = (d >> 1) | (d & 1);
+            return (uint)(((int)d - 1) >> 31);
+        }
+
+        public static uint EqualToZero(int len, uint[] x, int xOff)
+        {
+            uint d = 0;
+            for (int i = 0; i < len; ++i)
+            {
+                d |= x[xOff + i];
+            }
+            d = (d >> 1) | (d & 1);
+            return (uint)(((int)d - 1) >> 31);
+        }
+
         public static uint[] FromBigInteger(int bits, BigInteger x)
         {
             if (x.SignValue < 0 || x.BitLength > bits)
diff --git a/crypto/src/util/BigIntegers.cs b/crypto/src/util/BigIntegers.cs
index bac5f12c0..d9c898676 100644
--- a/crypto/src/util/BigIntegers.cs
+++ b/crypto/src/util/BigIntegers.cs
@@ -1,6 +1,7 @@
 using System;
 
 using Org.BouncyCastle.Math;
+using Org.BouncyCastle.Math.Raw;
 using Org.BouncyCastle.Security;
 
 namespace Org.BouncyCastle.Utilities
@@ -10,6 +11,9 @@ namespace Org.BouncyCastle.Utilities
      */
     public abstract class BigIntegers
     {
+        public static readonly BigInteger Zero = BigInteger.Zero;
+        public static readonly BigInteger One = BigInteger.One;
+
         private const int MaxIterations = 1000;
 
         /**
@@ -131,6 +135,52 @@ namespace Org.BouncyCastle.Utilities
             return new BigInteger(max.Subtract(min).BitLength - 1, random).Add(min);
         }
 
+        public static BigInteger ModOddInverse(BigInteger M, BigInteger X)
+        {
+            if (!M.TestBit(0))
+                throw new ArgumentException("must be odd", "M");
+            if (M.SignValue != 1)
+                throw new ArithmeticException("BigInteger: modulus not positive");
+            if (X.SignValue < 0 || X.CompareTo(M) >= 0)
+            {
+                X = X.Mod(M);
+            }
+
+            int bits = M.BitLength;
+            uint[] m = Nat.FromBigInteger(bits, M);
+            uint[] x = Nat.FromBigInteger(bits, X);
+            int len = m.Length;
+            uint[] z = Nat.Create(len);
+            if (0 == Mod.ModOddInverse(m, x, z))
+                throw new ArithmeticException("BigInteger not invertible");
+            return Nat.ToBigInteger(len, z);
+        }
+
+        public static BigInteger ModOddInverseVar(BigInteger M, BigInteger X)
+        {
+            if (!M.TestBit(0))
+                throw new ArgumentException("must be odd", "M");
+            if (M.SignValue != 1)
+                throw new ArithmeticException("BigInteger: modulus not positive");
+            if (M.Equals(One))
+                return Zero;
+            if (X.SignValue < 0 || X.CompareTo(M) >= 0)
+            {
+                X = X.Mod(M);
+            }
+            if (X.Equals(One))
+                return One;
+
+            int bits = M.BitLength;
+            uint[] m = Nat.FromBigInteger(bits, M);
+            uint[] x = Nat.FromBigInteger(bits, X);
+            int len = m.Length;
+            uint[] z = Nat.Create(len);
+            if (!Mod.ModOddInverseVar(m, x, z))
+                throw new ArithmeticException("BigInteger not invertible");
+            return Nat.ToBigInteger(len, z);
+        }
+
         public static int GetUnsignedByteLength(BigInteger n)
         {
             return (n.BitLength + 7) / 8;
diff --git a/crypto/src/util/Integers.cs b/crypto/src/util/Integers.cs
index afb4b827f..7d98de586 100644
--- a/crypto/src/util/Integers.cs
+++ b/crypto/src/util/Integers.cs
@@ -4,6 +4,11 @@ namespace Org.BouncyCastle.Utilities
 {
     public abstract class Integers
     {
+        private static readonly byte[] DeBruijnTZ = {
+            0x00, 0x01, 0x02, 0x18, 0x03, 0x13, 0x06, 0x19, 0x16, 0x04, 0x14, 0x0A,
+            0x10, 0x07, 0x0C, 0x1A, 0x1F, 0x17, 0x12, 0x05, 0x15, 0x09, 0x0F, 0x0B,
+            0x1E, 0x11, 0x08, 0x0E, 0x1D, 0x0D, 0x1C, 0x1B };
+
         public static int NumberOfLeadingZeros(int i)
         {
             if (i <= 0)
@@ -21,16 +26,7 @@ namespace Org.BouncyCastle.Utilities
 
         public static int NumberOfTrailingZeros(int i)
         {
-            if (i == 0)
-                return 32;
-
-            int count = 0;
-            while ((i & 1) == 0)
-            {
-                i >>= 1;
-                ++count;
-            }
-            return count;
+            return DeBruijnTZ[(uint)((i & -i) * 0x04D7651F) >> 27];
         }
 
         public static int RotateLeft(int i, int distance)