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