Significantly optimize signature verification.

This commit essentially rewrites all of the primitives needed to perform
the arithmetic for ECDSA signature verification of secp256k1 signatures to
significantly speed it up.  Benchmarking has shown signature verification
is roughly 10 times faster with this commit over the previous.

In particular, it introduces a new field value which is used to perform the
modular field arithmetic using fixed-precision operations specifically
tailored for the secp256k1 prime.  The field also takes advantage of
special properties of the prime for significantly faster modular reduction
than is available through generic methods.

In addition, the curve point addition and doubling have been optimized
minimize the number of field multiplications in favor field squarings
since they are quite a bit faster.  They routines also now look for
certain assumptions such as z values of 1 or equivalent z values which
can be used to further reduce the number of multiplicaitons needed when
possible.

Note there are still quite a few more optimizations that could be done
such as using precomputation for ScalarBaseMult, making use of the
secp256k1 endomorphism, and using windowed NAF, however this work already
offers significant performance improvements.

For example, testing 10000 random signature verifications resulted in:
New btcec took 15.9821565s
Old btcec took 2m34.1016716s

Closes conformal/btcd#26.

btcec.go

@@ -6,6 +6,10 @@
 
 package btcec
 
+// References:
+//   [SECG]: Recommended Elliptic Curve Domain Parameters
+//     http://www.secg.org/download/aid-784/sec2-v2.pdf
+
 // This package operates, internally, on Jacobian coordinates. For a given
 // (x, y) position on the curve, the Jacobian coordinates are (x1, y1, z1)
 // where x = x1/z1² and y = y1/z1³. The greatest speedups come when the whole
@@ -22,6 +26,13 @@ import (
 //TODO: examine if we need to care about EC optimization as descibed here
 //      https://bitcointalk.org/index.php?topic=155054.0;all
 
+var (
+	// fieldOne is simple the integer 1 in field representation.  It is
+	// used to avoid needing to create it multiple times during the internal
+	// arithmetic.
+	fieldOne = new(fieldVal).SetInt(1)
+)
+
 // KoblitzCurve supports a koblitz curve implementation that fits the ECC Curve
 // interface from crypto/elliptic.
 type KoblitzCurve struct {
@@ -29,200 +40,584 @@ type KoblitzCurve struct {
 	q *big.Int
 }
 
-// Params returns the parameters fro the curve.
+// Params returns the parameters for the curve.
 func (curve *KoblitzCurve) Params() *elliptic.CurveParams {
 	return curve.CurveParams
 }
 
+// bigAffineToField takes an affine point (x, y) as big integers and converts
+// it to an affine point as field values.
+func (curve *KoblitzCurve) bigAffineToField(x, y *big.Int) (*fieldVal, *fieldVal) {
+	x3, y3 := new(fieldVal), new(fieldVal)
+	x3.SetByteSlice(x.Bytes())
+	y3.SetByteSlice(y.Bytes())
+
+	return x3, y3
+}
+
+// fieldJacobianToBigAffine takes a Jacobian point (x, y, z) as field values and
+// converts it to an affine point as big integers.
+func (curve *KoblitzCurve) fieldJacobianToBigAffine(x, y, z *fieldVal) (*big.Int, *big.Int) {
+	// Inversions are expensive and both point addition and point doubling
+	// are faster when working with points that have a z value of one.  So,
+	// if the point needs to be converted to affine, go ahead and normalize
+	// the point itself at the same time as the calculation is the same.
+	var zInv, tempZ fieldVal
+	zInv.Set(z).Inverse()   // zInv = Z^-1
+	tempZ.SquareVal(&zInv)  // tempZ = Z^-2
+	x.Mul(&tempZ)           // X = X/Z^2 (mag: 1)
+	y.Mul(tempZ.Mul(&zInv)) // Y = Y/Z^3 (mag: 1)
+	z.SetInt(1)             // Z = 1 (mag: 1)
+
+	// Normalize the x and y values.
+	x.Normalize()
+	y.Normalize()
+
+	// Convert the field values for the now affine point to big.Ints.
+	x3, y3 := new(big.Int), new(big.Int)
+	x3.SetBytes(x.Bytes()[:])
+	y3.SetBytes(y.Bytes()[:])
+	return x3, y3
+}
+
 // IsOnCurve returns boolean if the point (x,y) is on the curve.
 // Part of the elliptic.Curve interface. This function differs from the
 // crypto/elliptic algorithm since a = 0 not -3.
 func (curve *KoblitzCurve) IsOnCurve(x, y *big.Int) bool {
-	// y² = x³ + b
+	// Convert big ints to field values for faster arithmetic.
-	y2 := new(big.Int).Mul(y, y) //y²
+	fx, fy := curve.bigAffineToField(x, y)
-	y2.Mod(y2, curve.P)          //y²%P
 
-	x3 := new(big.Int).Mul(x, x) //x²
+	// Elliptic curve equation for secp256k1 is: y^2 = x^3 + 7
-	x3.Mul(x3, x)                //x³
+	y2 := new(fieldVal).SquareVal(fy).Normalize()
-
+	result := new(fieldVal).SquareVal(fx).Mul(fx).AddInt(7).Normalize()
-	x3.Add(x3, curve.B) //x³+B
+	return y2.Equals(result)
-	x3.Mod(x3, curve.P) //(x³+B)%P
-
-	return x3.Cmp(y2) == 0
 }
 
-// zForAffine returns a Jacobian Z value for the affine point (x, y). If x and
+// addZ1AndZ2EqualsOne adds two Jacobian points that are already known to have
-// y are zero, it assumes that they represent the point at infinity because (0,
+// z values of 1 and stores the result in (x3, y3, z3).  That is to say
-// 0) is not on the any of the curves handled here.
+// (x1, y1, 1) + (x2, y2, 1) = (x3, y3, z3).  It performs faster addition than
-func zForAffine(x, y *big.Int) *big.Int {
+// the generic add routine since less arithmetic is needed due to the ability to
-	z := new(big.Int)
+// avoid the z value multiplications.
-	if x.Sign() != 0 || y.Sign() != 0 {
+func (curve *KoblitzCurve) addZ1AndZ2EqualsOne(x1, y1, x2, y2, x3, y3, z3 *fieldVal) {
-		z.SetInt64(1)
+	// To compute the point addition efficiently, this implementation splits
-	}
+	// the equation into intermediate elements which are used to minimize
-	return z
+	// the number of field multiplications using the method shown at:
-}
+	// http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-mmadd-2007-bl
+	//
+	// In particular it performs the calculations using the following:
+	// H = X2-X1, HH = H^2, I = 4*HH, J = H*I, r = 2*(Y2-Y1), V = X1*I
+	// X3 = r^2-J-2*V, Y3 = r*(V-X3)-2*Y1*J, Z3 = 2*H
+	//
+	// This results in a cost of 4 field multiplications, 2 field squarings,
+	// 6 field additions, and 5 integer multiplications.
 
-// affineFromJacobian reverses the Jacobian transform. See the comment at the
+	// When the x coordinates are the same for two points on the curve, the
-// top of the file. If the point is ∞ it returns 0, 0.
+	// y coordinates either must be the same, in which case it is point
-func (curve *KoblitzCurve) affineFromJacobian(x, y, z *big.Int) (xOut, yOut *big.Int) {
+	// doubling, or they are opposite and the result is the point at
-	if z.Sign() == 0 {
+	// infinity per the group law for elliptic curve cryptography.
-		return new(big.Int), new(big.Int)
+	x1.Normalize()
-	}
+	y1.Normalize()
-
+	x2.Normalize()
-	zinv := new(big.Int).ModInverse(z, curve.P)
+	y2.Normalize()
-	zinvsq := new(big.Int).Mul(zinv, zinv)
+	if x1.Equals(x2) {
-
+		if y1.Equals(y2) {
-	xOut = new(big.Int).Mul(x, zinvsq)
+			// Since x1 == x2 and y1 == y2, point doubling must be
-	xOut.Mod(xOut, curve.P)
+			// done, otherwise the addition would end up dividing
-	zinvsq.Mul(zinvsq, zinv)
+			// by zero.
-	yOut = new(big.Int).Mul(y, zinvsq)
+			curve.doubleJacobian(x1, y1, fieldOne, x3, y3, z3)
-	yOut.Mod(yOut, curve.P)
 			return
 		}
 
-// Add returns the sum of (x1,y1 and (x2,y2). Part of the elliptic.Curve
+		// Since x1 == x2 and y1 == -y2, the sum is the point at
-// interface.
+		// infinity per the group law.
-func (curve *KoblitzCurve) Add(x1, y1, x2, y2 *big.Int) (*big.Int, *big.Int) {
+		x3.SetInt(0)
-	z1 := zForAffine(x1, y1)
+		y3.SetInt(0)
-	z2 := zForAffine(x2, y2)
+		z3.SetInt(0)
-	return curve.affineFromJacobian(curve.addJacobian(x1, y1, z1, x2, y2, z2))
+		return
 	}
 
-// addJacobian takes two points in Jacobian coordinates, (x1, y1, z1) and
+	// Calculate X3, Y3, and Z3 according to the intermediate elements
-// (x2, y2, z2) and returns their sum, also in Jacobian form.
+	// breakdown above.
-func (curve *KoblitzCurve) addJacobian(x1, y1, z1, x2, y2, z2 *big.Int) (*big.Int, *big.Int, *big.Int) {
+	var h, i, j, r, v fieldVal
-	// See http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-add-2007-bl
+	var negJ, neg2V, negX3 fieldVal
-	x3, y3, z3 := new(big.Int), new(big.Int), new(big.Int)
+	h.Set(x1).Negate(1).Add(x2)                // H = X2-X1 (mag: 3)
-	if z1.Sign() == 0 {
+	i.SquareVal(&h).MulInt(4)                  // I = 4*H^2 (mag: 4)
+	j.Mul2(&h, &i)                             // J = H*I (mag: 1)
+	r.Set(y1).Negate(1).Add(y2).MulInt(2)      // r = 2*(Y2-Y1) (mag: 6)
+	v.Mul2(x1, &i)                             // V = X1*I (mag: 1)
+	negJ.Set(&j).Negate(1)                     // negJ = -J (mag: 2)
+	neg2V.Set(&v).MulInt(2).Negate(2)          // neg2V = -(2*V) (mag: 3)
+	x3.Set(&r).Square().Add(&negJ).Add(&neg2V) // X3 = r^2-J-2*V (mag: 6)
+	negX3.Set(x3).Negate(6)                    // negX3 = -X3 (mag: 7)
+	j.Mul(y1).MulInt(2).Negate(2)              // J = -(2*Y1*J) (mag: 3)
+	y3.Set(&v).Add(&negX3).Mul(&r).Add(&j)     // Y3 = r*(V-X3)-2*Y1*J (mag: 4)
+	z3.Set(&h).MulInt(2)                       // Z3 = 2*H (mag: 6)
+
+	// Normalize the resulting field values to a magnitude of 1 as needed.
+	x3.Normalize()
+	y3.Normalize()
+	z3.Normalize()
+}
+
+// addZ1EqualsZ2 adds two Jacobian points that are already known to have the
+// same z value and stores the result in (x3, y3, z3).  That is to say
+// (x1, y1, z1) + (x2, y2, z1) = (x3, y3, z3).  It performs faster addition than
+// the generic add routine since less arithmetic is needed due to the known
+// equivalence.
+func (curve *KoblitzCurve) addZ1EqualsZ2(x1, y1, z1, x2, y2, x3, y3, z3 *fieldVal) {
+	// To compute the point addition efficiently, this implementation splits
+	// the equation into intermediate elements which are used to minimize
+	// the number of field multiplications using a slightly modified version
+	// of the method shown at:
+	// http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-mmadd-2007-bl
+	//
+	// In particular it performs the calculations using the following:
+	// A = X2-X1, B = A^2, C=Y2-Y1, D = C^2, E = X1*B, F = X2*B
+	// X3 = D-E-F, Y3 = C*(E-X3)-Y1*(F-E), Z3 = Z1*A
+	//
+	// This results in a cost of 5 field multiplications, 2 field squarings,
+	// 9 field additions, and 0 integer multiplications.
+
+	// When the x coordinates are the same for two points on the curve, the
+	// y coordinates either must be the same, in which case it is point
+	// doubling, or they are opposite and the result is the point at
+	// infinity per the group law for elliptic curve cryptography.
+	x1.Normalize()
+	y1.Normalize()
+	x2.Normalize()
+	y2.Normalize()
+	if x1.Equals(x2) {
+		if y1.Equals(y2) {
+			// Since x1 == x2 and y1 == y2, point doubling must be
+			// done, otherwise the addition would end up dividing
+			// by zero.
+			curve.doubleJacobian(x1, y1, z1, x3, y3, z3)
+			return
+		}
+
+		// Since x1 == x2 and y1 == -y2, the sum is the point at
+		// infinity per the group law.
+		x3.SetInt(0)
+		y3.SetInt(0)
+		z3.SetInt(0)
+		return
+	}
+
+	// Calculate X3, Y3, and Z3 according to the intermediate elements
+	// breakdown above.
+	var a, b, c, d, e, f fieldVal
+	var negX1, negY1, negE, negX3 fieldVal
+	negX1.Set(x1).Negate(1)                // negX1 = -X1 (mag: 2)
+	negY1.Set(y1).Negate(1)                // negY1 = -Y1 (mag: 2)
+	a.Set(&negX1).Add(x2)                  // A = X2-X1 (mag: 3)
+	b.SquareVal(&a)                        // B = A^2 (mag: 1)
+	c.Set(&negY1).Add(y2)                  // C = Y2-Y1 (mag: 3)
+	d.SquareVal(&c)                        // D = C^2 (mag: 1)
+	e.Mul2(x1, &b)                         // E = X1*B (mag: 1)
+	negE.Set(&e).Negate(1)                 // negE = -E (mag: 2)
+	f.Mul2(x2, &b)                         // F = X2*B (mag: 1)
+	x3.Add2(&e, &f).Negate(3).Add(&d)      // X3 = D-E-F (mag: 5)
+	negX3.Set(x3).Negate(5).Normalize()    // negX3 = -X3 (mag: 1)
+	y3.Set(y1).Mul(f.Add(&negE)).Negate(3) // Y3 = -(Y1*(F-E)) (mag: 4)
+	y3.Add(e.Add(&negX3).Mul(&c))          // Y3 = C*(E-X3)+Y3 (mag: 5)
+	z3.Mul2(z1, &a)                        // Z3 = Z1*A (mag: 1)
+
+	// Normalize the resulting field values to a magnitude of 1 as needed.
+	x3.Normalize()
+	y3.Normalize()
+}
+
+// addZ2EqualsOne adds two Jacobian points when the second point is already
+// known to have a z value of 1 (and the z value for the first point is not 1)
+// and stores the result in (x3, y3, z3).  That is to say (x1, y1, z1) +
+// (x2, y2, 1) = (x3, y3, z3).  It performs faster addition than the generic
+// add routine since less arithmetic is needed due to the ability to avoid
+// multiplications by the second point's z value.
+func (curve *KoblitzCurve) addZ2EqualsOne(x1, y1, z1, x2, y2, x3, y3, z3 *fieldVal) {
+	// To compute the point addition efficiently, this implementation splits
+	// the equation into intermediate elements which are used to minimize
+	// the number of field multiplications using the method shown at:
+	// http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-madd-2007-bl
+	//
+	// In particular it performs the calculations using the following:
+	// Z1Z1 = Z1^2, U2 = X2*Z1Z1, S2 = Y2*Z1*Z1Z1, H = U2-X1, HH = H^2,
+	// I = 4*HH, J = H*I, r = 2*(S2-Y1), V = X1*I
+	// X3 = r^2-J-2*V, Y3 = r*(V-X3)-2*Y1*J, Z3 = (Z1+H)^2-Z1Z1-HH
+	//
+	// This results in a cost of 7 field multiplications, 4 field squarings,
+	// 9 field additions, and 4 integer multiplications.
+
+	// When the x coordinates are the same for two points on the curve, the
+	// y coordinates either must be the same, in which case it is point
+	// doubling, or they are opposite and the result is the point at
+	// infinity per the group law for elliptic curve cryptography.  Since
+	// any number of Jacobian coordinates can represent the same affine
+	// point, the x and y values need to be converted to like terms.  Due to
+	// the assumption made for this function that the second point has a z
+	// value of 1 (z2=1), the first point is already "converted".
+	var z1z1, u2, s2 fieldVal
+	x1.Normalize()
+	y1.Normalize()
+	z1z1.SquareVal(z1)                        // Z1Z1 = Z1^2 (mag: 1)
+	u2.Set(x2).Mul(&z1z1).Normalize()         // U2 = X2*Z1Z1 (mag: 1)
+	s2.Set(y2).Mul(&z1z1).Mul(z1).Normalize() // S2 = Y2*Z1*Z1Z1 (mag: 1)
+	if x1.Equals(&u2) {
+		if y1.Equals(&s2) {
+			// Since x1 == x2 and y1 == y2, point doubling must be
+			// done, otherwise the addition would end up dividing
+			// by zero.
+			curve.doubleJacobian(x1, y1, z1, x3, y3, z3)
+			return
+		}
+
+		// Since x1 == x2 and y1 == -y2, the sum is the point at
+		// infinity per the group law.
+		x3.SetInt(0)
+		y3.SetInt(0)
+		z3.SetInt(0)
+		return
+	}
+
+	// Calculate X3, Y3, and Z3 according to the intermediate elements
+	// breakdown above.
+	var h, hh, i, j, r, rr, v fieldVal
+	var negX1, negY1, negX3 fieldVal
+	negX1.Set(x1).Negate(1)                // negX1 = -X1 (mag: 2)
+	h.Add2(&u2, &negX1)                    // H = U2-X1 (mag: 3)
+	hh.SquareVal(&h)                       // HH = H^2 (mag: 1)
+	i.Set(&hh).MulInt(4)                   // I = 4 * HH (mag: 4)
+	j.Mul2(&h, &i)                         // J = H*I (mag: 1)
+	negY1.Set(y1).Negate(1)                // negY1 = -Y1 (mag: 2)
+	r.Set(&s2).Add(&negY1).MulInt(2)       // r = 2*(S2-Y1) (mag: 6)
+	rr.SquareVal(&r)                       // rr = r^2 (mag: 1)
+	v.Mul2(x1, &i)                         // V = X1*I (mag: 1)
+	x3.Set(&v).MulInt(2).Add(&j).Negate(3) // X3 = -(J+2*V) (mag: 4)
+	x3.Add(&rr)                            // X3 = r^2+X3 (mag: 5)
+	negX3.Set(x3).Negate(5)                // negX3 = -X3 (mag: 6)
+	y3.Set(y1).Mul(&j).MulInt(2).Negate(2) // Y3 = -(2*Y1*J) (mag: 3)
+	y3.Add(v.Add(&negX3).Mul(&r))          // Y3 = r*(V-X3)+Y3 (mag: 4)
+	z3.Add2(z1, &h).Square()               // Z3 = (Z1+H)^2 (mag: 1)
+	z3.Add(z1z1.Add(&hh).Negate(2))        // Z3 = Z3-(Z1Z1+HH) (mag: 4)
+
+	// Normalize the resulting field values to a magnitude of 1 as needed.
+	x3.Normalize()
+	y3.Normalize()
+	z3.Normalize()
+}
+
+// addGeneric adds two Jacobian points (x1, y1, z1) and (x2, y2, z2) without any
+// assumptions about the z values of the two points and stores the result in
+// (x3, y3, z3).  That is to say (x1, y1, z1) + (x2, y2, z2) = (x3, y3, z3).  It
+// is the slowest of the add routines due to requiring the most arithmetic.
+func (curve *KoblitzCurve) addGeneric(x1, y1, z1, x2, y2, z2, x3, y3, z3 *fieldVal) {
+	// To compute the point addition efficiently, this implementation splits
+	// the equation into intermediate elements which are used to minimize
+	// the number of field multiplications using the method shown at:
+	// http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-add-2007-bl
+	//
+	// In particular it performs the calculations using the following:
+	// Z1Z1 = Z1^2, Z2Z2 = Z2^2, U1 = X1*Z2Z2, U2 = X2*Z1Z1, S1 = Y1*Z2*Z2Z2
+	// S2 = Y2*Z1*Z1Z1, H = U2-U1, I = (2*H)^2, J = H*I, r = 2*(S2-S1)
+	// V = U1*I
+	// X3 = r^2-J-2*V, Y3 = r*(V-X3)-2*S1*J, Z3 = ((Z1+Z2)^2-Z1Z1-Z2Z2)*H
+	//
+	// This results in a cost of 11 field multiplications, 5 field squarings,
+	// 9 field additions, and 4 integer multiplications.
+
+	// When the x coordinates are the same for two points on the curve, the
+	// y coordinates either must be the same, in which case it is point
+	// doubling, or they are opposite and the result is the point at
+	// infinity.  Since any number of Jacobian coordinates can represent the
+	// same affine point, the x and y values need to be converted to like
+	// terms.
+	var z1z1, z2z2, u1, u2, s1, s2 fieldVal
+	z1z1.SquareVal(z1)                        // Z1Z1 = Z1^2 (mag: 1)
+	z2z2.SquareVal(z2)                        // Z2Z2 = Z2^2 (mag: 1)
+	u1.Set(x1).Mul(&z2z2).Normalize()         // U1 = X1*Z2Z2 (mag: 1)
+	u2.Set(x2).Mul(&z1z1).Normalize()         // U2 = X2*Z1Z1 (mag: 1)
+	s1.Set(y1).Mul(&z2z2).Mul(z2).Normalize() // S1 = Y1*Z2*Z2Z2 (mag: 1)
+	s2.Set(y2).Mul(&z1z1).Mul(z1).Normalize() // S2 = Y2*Z1*Z1Z1 (mag: 1)
+	if u1.Equals(&u2) {
+		if s1.Equals(&s2) {
+			// Since x1 == x2 and y1 == y2, point doubling must be
+			// done, otherwise the addition would end up dividing
+			// by zero.
+			curve.doubleJacobian(x1, y1, z1, x3, y3, z3)
+			return
+		}
+
+		// Since x1 == x2 and y1 == -y2, the sum is the point at
+		// infinity per the group law.
+		x3.SetInt(0)
+		y3.SetInt(0)
+		z3.SetInt(0)
+		return
+	}
+
+	// Calculate X3, Y3, and Z3 according to the intermediate elements
+	// breakdown above.
+	var h, i, j, r, rr, v fieldVal
+	var negU1, negS1, negX3 fieldVal
+	negU1.Set(&u1).Negate(1)               // negU1 = -U1 (mag: 2)
+	h.Add2(&u2, &negU1)                    // H = U2-U1 (mag: 3)
+	i.Set(&h).MulInt(2).Square()           // I = (2*H)^2 (mag: 2)
+	j.Mul2(&h, &i)                         // J = H*I (mag: 1)
+	negS1.Set(&s1).Negate(1)               // negS1 = -S1 (mag: 2)
+	r.Set(&s2).Add(&negS1).MulInt(2)       // r = 2*(S2-S1) (mag: 6)
+	rr.SquareVal(&r)                       // rr = r^2 (mag: 1)
+	v.Mul2(&u1, &i)                        // V = U1*I (mag: 1)
+	x3.Set(&v).MulInt(2).Add(&j).Negate(3) // X3 = -(J+2*V) (mag: 4)
+	x3.Add(&rr)                            // X3 = r^2+X3 (mag: 5)
+	negX3.Set(x3).Negate(5)                // negX3 = -X3 (mag: 6)
+	y3.Mul2(&s1, &j).MulInt(2).Negate(2)   // Y3 = -(2*S1*J) (mag: 3)
+	y3.Add(v.Add(&negX3).Mul(&r))          // Y3 = r*(V-X3)+Y3 (mag: 4)
+	z3.Add2(z1, z2).Square()               // Z3 = (Z1+Z2)^2 (mag: 1)
+	z3.Add(z1z1.Add(&z2z2).Negate(2))      // Z3 = Z3-(Z1Z1+Z2Z2) (mag: 4)
+	z3.Mul(&h)                             // Z3 = Z3*H (mag: 1)
+
+	// Normalize the resulting field values to a magnitude of 1 as needed.
+	x3.Normalize()
+	y3.Normalize()
+}
+
+// addJacobian adds the passed Jacobian points (x1, y1, z1) and (x2, y2, z2)
+// together and stores the result in (x3, y3, z3).
+func (curve *KoblitzCurve) addJacobian(x1, y1, z1, x2, y2, z2, x3, y3, z3 *fieldVal) {
+	// A point at infinity is the identity according to the group law for
+	// elliptic curve cryptography.  Thus, ∞ + P = P and P + ∞ = P.
+	if (x1.IsZero() && y1.IsZero()) || z1.IsZero() {
 		x3.Set(x2)
 		y3.Set(y2)
 		z3.Set(z2)
-		return x3, y3, z3
+		return
 	}
-	if z2.Sign() == 0 {
+	if (x2.IsZero() && y2.IsZero()) || z2.IsZero() {
 		x3.Set(x1)
 		y3.Set(y1)
 		z3.Set(z1)
-		return x3, y3, z3
+		return
 	}
 
-	z1z1 := new(big.Int).Mul(z1, z1)
+	// Faster point addition can be achieved when certain assumptions are
-	z1z1.Mod(z1z1, curve.P)
+	// met.  For example, when both points have the same z value, arithmetic
-	z2z2 := new(big.Int).Mul(z2, z2)
+	// on the z values can be avoided.  This section thus checks for these
-	z2z2.Mod(z2z2, curve.P)
+	// conditions and calls an appropriate add function which is accelerated
-
+	// by using those assumptions.
-	u1 := new(big.Int).Mul(x1, z2z2)
+	z1.Normalize()
-	u1.Mod(u1, curve.P)
+	z2.Normalize()
-	u2 := new(big.Int).Mul(x2, z1z1)
+	isZ1One := z1.Equals(fieldOne)
-	u2.Mod(u2, curve.P)
+	isZ2One := z2.Equals(fieldOne)
-	h := new(big.Int).Sub(u2, u1)
+	switch {
-	xEqual := h.Sign() == 0
+	case isZ1One && isZ2One:
-	if h.Sign() == -1 {
+		curve.addZ1AndZ2EqualsOne(x1, y1, x2, y2, x3, y3, z3)
-		h.Add(h, curve.P)
+		return
+	case z1.Equals(z2):
+		curve.addZ1EqualsZ2(x1, y1, z1, x2, y2, x3, y3, z3)
+		return
+	case isZ2One:
+		curve.addZ2EqualsOne(x1, y1, z1, x2, y2, x3, y3, z3)
+		return
 	}
-	i := new(big.Int).Lsh(h, 1)
-	i.Mul(i, i)
-	j := new(big.Int).Mul(h, i)
 
-	s1 := new(big.Int).Mul(y1, z2)
+	// None of the above assumptions are true, so fall back to generic
-	s1.Mul(s1, z2z2)
+	// point addition.
-	s1.Mod(s1, curve.P)
+	curve.addGeneric(x1, y1, z1, x2, y2, z2, x3, y3, z3)
-	s2 := new(big.Int).Mul(y2, z1)
-	s2.Mul(s2, z1z1)
-	s2.Mod(s2, curve.P)
-	r := new(big.Int).Sub(s2, s1)
-	if r.Sign() == -1 {
-		r.Add(r, curve.P)
 }
-	yEqual := r.Sign() == 0
+
-	if xEqual && yEqual {
+// Add returns the sum of (x1,y1) and (x2,y2). Part of the elliptic.Curve
-		return curve.doubleJacobian(x1, y1, z1)
+// interface.
+func (curve *KoblitzCurve) Add(x1, y1, x2, y2 *big.Int) (*big.Int, *big.Int) {
+	// A point at infinity is the identity according to the group law for
+	// elliptic curve cryptography.  Thus, ∞ + P = P and P + ∞ = P.
+	if x1.Sign() == 0 && y1.Sign() == 0 {
+		return x2, y2
+	}
+	if x2.Sign() == 0 && y2.Sign() == 0 {
+		return x1, y1
 	}
-	r.Lsh(r, 1)
-	v := new(big.Int).Mul(u1, i)
 
-	x3.Set(r)
+	// Convert the affine coordinates from big integers to field values
-	x3.Mul(x3, x3)
+	// and do the point addition in Jacobian projective space.
-	x3.Sub(x3, j)
+	fx1, fy1 := curve.bigAffineToField(x1, y1)
-	x3.Sub(x3, v)
+	fx2, fy2 := curve.bigAffineToField(x2, y2)
-	x3.Sub(x3, v)
+	fx3, fy3, fz3 := new(fieldVal), new(fieldVal), new(fieldVal)
-	x3.Mod(x3, curve.P)
+	curve.addJacobian(fx1, fy1, fieldOne, fx2, fy2, fieldOne, fx3, fy3, fz3)
 
-	y3.Set(r)
+	// Convert the Jacobian coordinate field values back to affine big
-	v.Sub(v, x3)
+	// integers.
-	y3.Mul(y3, v)
+	return curve.fieldJacobianToBigAffine(fx3, fy3, fz3)
-	s1.Mul(s1, j)
+}
-	s1.Lsh(s1, 1)
-	y3.Sub(y3, s1)
-	y3.Mod(y3, curve.P)
 
-	z3.Add(z1, z2)
+// doubleZ1EqualsOne performs point doubling on the passed Jacobian point
-	z3.Mul(z3, z3)
+// when the point is already known to have a z value of 1 and stores
-	z3.Sub(z3, z1z1)
+// the result in (x3, y3, z3).  That is to say (x3, y3, z3) = 2*(x1, y1, 1).  It
-	z3.Sub(z3, z2z2)
+// performs faster point doubling than the generic routine since less arithmetic
-	z3.Mul(z3, h)
+// is needed due to the ability to avoid multiplication by the z value.
-	z3.Mod(z3, curve.P)
+func (curve *KoblitzCurve) doubleZ1EqualsOne(x1, y1, x3, y3, z3 *fieldVal) {
+	// This function uses the assumptions that z1 is 1, thus the point
+	// doubling formulas reduce to:
+	//
+	// X3 = (3*X1^2)^2 - 8*X1*Y1^2
+	// Y3 = (3*X1^2)*(4*X1*Y1^2 - X3) - 8*Y1^4
+	// Z3 = 2*Y1
+	//
+	// To compute the above efficiently, this implementation splits the
+	// equation into intermediate elements which are used to minimize the
+	// number of field multiplications in favor of field squarings which
+	// are roughly 35% faster than field multiplications with the current
+	// implementation at the time this was written.
+	//
+	// This uses a slightly modified version of the method shown at:
+	// http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-mdbl-2007-bl
+	//
+	// In particular it performs the calculations using the following:
+	// A = X1^2, B = Y1^2, C = B^2, D = 2*((X1+B)^2-A-C)
+	// E = 3*A, F = E^2, X3 = F-2*D, Y3 = E*(D-X3)-8*C
+	// Z3 = 2*Y1
+	//
+	// This results in a cost of 1 field multiplication, 5 field squarings,
+	// 6 field additions, and 5 integer multiplications.
+	var a, b, c, d, e, f fieldVal
+	z3.Set(y1).MulInt(2)                     // Z3 = 2*Y1 (mag: 2)
+	a.SquareVal(x1)                          // A = X1^2 (mag: 1)
+	b.SquareVal(y1)                          // B = Y1^2 (mag: 1)
+	c.SquareVal(&b)                          // C = B^2 (mag: 1)
+	b.Add(x1).Square()                       // B = (X1+B)^2 (mag: 1)
+	d.Set(&a).Add(&c).Negate(2)              // D = -(A+C) (mag: 3)
+	d.Add(&b).MulInt(2)                      // D = 2*(B+D)(mag: 8)
+	e.Set(&a).MulInt(3)                      // E = 3*A (mag: 3)
+	f.SquareVal(&e)                          // F = E^2 (mag: 1)
+	x3.Set(&d).MulInt(2).Negate(16)          // X3 = -(2*D) (mag: 17)
+	x3.Add(&f)                               // X3 = F+X3 (mag: 18)
+	f.Set(x3).Negate(18).Add(&d).Normalize() // F = D-X3 (mag: 1)
+	y3.Set(&c).MulInt(8).Negate(8)           // Y3 = -(8*C) (mag: 9)
+	y3.Add(f.Mul(&e))                        // Y3 = E*F+Y3 (mag: 10)
 
-	return x3, y3, z3
+	// Normalize the field values back to a magnitude of 1.
+	x3.Normalize()
+	y3.Normalize()
+	z3.Normalize()
+}
+
+// doubleGeneric performs point doubling on the passed Jacobian point without
+// any assumptions about the z value and stores the result in (x3, y3, z3).
+// That is to say (x3, y3, z3) = 2*(x1, y1, z1).  It is the slowest of the point
+// doubling routines due to requiring the most arithmetic.
+func (cuve *KoblitzCurve) doubleGeneric(x1, y1, z1, x3, y3, z3 *fieldVal) {
+	// Point doubling formula for Jacobian coordinates for the secp256k1
+	// curve:
+	// X3 = (3*X1^2)^2 - 8*X1*Y1^2
+	// Y3 = (3*X1^2)*(4*X1*Y1^2 - X3) - 8*Y1^4
+	// Z3 = 2*Y1*Z1
+	//
+	// To compute the above efficiently, this implementation splits the
+	// equation into intermediate elements which are used to minimize the
+	// number of field multiplications in favor of field squarings which
+	// are roughly 35% faster than field multiplications with the current
+	// implementation at the time this was written.
+	//
+	// This uses a slightly modified version of the method shown at:
+	// http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-dbl-2009-l
+	//
+	// In particular it performs the calculations using the following:
+	// A = X1^2, B = Y1^2, C = B^2, D = 2*((X1+B)^2-A-C)
+	// E = 3*A, F = E^2, X3 = F-2*D, Y3 = E*(D-X3)-8*C
+	// Z3 = 2*Y1*Z1
+	//
+	// This results in a cost of 1 field multiplication, 5 field squarings,
+	// 6 field additions, and 5 integer multiplications.
+	var a, b, c, d, e, f fieldVal
+	z3.Mul2(y1, z1).MulInt(2)                // Z3 = 2*Y1*Z1 (mag: 2)
+	a.SquareVal(x1)                          // A = X1^2 (mag: 1)
+	b.SquareVal(y1)                          // B = Y1^2 (mag: 1)
+	c.SquareVal(&b)                          // C = B^2 (mag: 1)
+	b.Add(x1).Square()                       // B = (X1+B)^2 (mag: 1)
+	d.Set(&a).Add(&c).Negate(2)              // D = -(A+C) (mag: 3)
+	d.Add(&b).MulInt(2)                      // D = 2*(B+D)(mag: 8)
+	e.Set(&a).MulInt(3)                      // E = 3*A (mag: 3)
+	f.SquareVal(&e)                          // F = E^2 (mag: 1)
+	x3.Set(&d).MulInt(2).Negate(16)          // X3 = -(2*D) (mag: 17)
+	x3.Add(&f)                               // X3 = F+X3 (mag: 18)
+	f.Set(x3).Negate(18).Add(&d).Normalize() // F = D-X3 (mag: 1)
+	y3.Set(&c).MulInt(8).Negate(8)           // Y3 = -(8*C) (mag: 9)
+	y3.Add(f.Mul(&e))                        // Y3 = E*F+Y3 (mag: 10)
+
+	// Normalize the field values back to a magnitude of 1.
+	x3.Normalize()
+	y3.Normalize()
+	z3.Normalize()
+}
+
+// doubleJacobian doubles the passed Jacobian point (x1, y1, z1) and stores the
+// result in (x3, y3, z3).
+func (curve *KoblitzCurve) doubleJacobian(x1, y1, z1, x3, y3, z3 *fieldVal) {
+	// Doubling a point at infinity is still infinity.
+	if y1.IsZero() || z1.IsZero() {
+		x3.SetInt(0)
+		y3.SetInt(0)
+		z3.SetInt(0)
+		return
+	}
+
+	// Slightly faster point doubling can be achieved when the z value is 1
+	// by avoiding the multiplication on the z value.  This section calls
+	// a point doubling function which is accelerated by using that
+	// assumption when possible.
+	if z1.Normalize().Equals(fieldOne) {
+		curve.doubleZ1EqualsOne(x1, y1, x3, y3, z3)
+		return
+	}
+
+	// Fall back to generic point doubling which works with arbitrary z
+	// values.
+	curve.doubleGeneric(x1, y1, z1, x3, y3, z3)
 }
 
 // Double returns 2*(x1,y1). Part of the elliptic.Curve interface.
 func (curve *KoblitzCurve) Double(x1, y1 *big.Int) (*big.Int, *big.Int) {
-	z1 := zForAffine(x1, y1)
+	if y1.Sign() == 0 {
-	return curve.affineFromJacobian(curve.doubleJacobian(x1, y1, z1))
+		return new(big.Int), new(big.Int)
 	}
 
-// doubleJacobian takes a point in Jacobian coordinates, (x, y, z), and
+	// Convert the affine coordinates from big integers to field values
-// returns its double, also in Jacobian form.
+	// and do the point doubling in Jacobian projective space.
-func (curve *KoblitzCurve) doubleJacobian(x, y, z *big.Int) (*big.Int, *big.Int, *big.Int) {
+	fx1, fy1 := curve.bigAffineToField(x1, y1)
-	// See http://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-dbl-2009-l
+	fx3, fy3, fz3 := new(fieldVal), new(fieldVal), new(fieldVal)
+	curve.doubleJacobian(fx1, fy1, fieldOne, fx3, fy3, fz3)
 
-	a := new(big.Int).Mul(x, x) //X1²
+	// Convert the Jacobian coordinate field values back to affine big
-	b := new(big.Int).Mul(y, y) //Y1²
+	// integers.
-	c := new(big.Int).Mul(b, b) //B²
+	return curve.fieldJacobianToBigAffine(fx3, fy3, fz3)
-
-	d := new(big.Int).Add(x, b) //X1+B
-	d.Mul(d, d)                 //(X1+B)²
-	d.Sub(d, a)                 //(X1+B)²-A
-	d.Sub(d, c)                 //(X1+B)²-A-C
-	d.Mul(d, big.NewInt(2))     //2*((X1+B)²-A-C)
-
-	e := new(big.Int).Mul(big.NewInt(3), a) //3*A
-	f := new(big.Int).Mul(e, e)             //E²
-
-	x3 := new(big.Int).Mul(big.NewInt(2), d) //2*D
-	x3.Sub(f, x3)                            //F-2*D
-	x3.Mod(x3, curve.P)
-
-	y3 := new(big.Int).Sub(d, x3)                  //D-X3
-	y3.Mul(e, y3)                                  //E*(D-X3)
-	y3.Sub(y3, new(big.Int).Mul(big.NewInt(8), c)) //E*(D-X3)-8*C
-	y3.Mod(y3, curve.P)
-
-	z3 := new(big.Int).Mul(y, z) //Y1*Z1
-	z3.Mul(big.NewInt(2), z3)    //3*Y1*Z1
-	z3.Mod(z3, curve.P)
-
-	return x3, y3, z3
 }
 
 // ScalarMult returns k*(Bx, By) where k is a big endian integer.
 // Part of the elliptic.Curve interface.
 func (curve *KoblitzCurve) ScalarMult(Bx, By *big.Int, k []byte) (*big.Int, *big.Int) {
-	Bz := new(big.Int).SetInt64(1)
+	// This uses the left to right binary method for point multiplication:
-	x, y, z := new(big.Int), new(big.Int), new(big.Int)
 
-	for _, byte := range k {
+	// Point Q = ∞ (point at infinity).
+	qx, qy, qz := new(fieldVal), new(fieldVal), new(fieldVal)
+
+	// Point P = the point to multiply the scalar with.
+	px, py := curve.bigAffineToField(Bx, By)
+	pz := fieldOne
+
+	// Double and add as necessary depending on the bits set in the scalar.
+	for _, byteVal := range k {
 		for bitNum := 0; bitNum < 8; bitNum++ {
-			x, y, z = curve.doubleJacobian(x, y, z)
+			// Q = 2*Q
-			if byte&0x80 == 0x80 {
+			curve.doubleJacobian(qx, qy, qz, qx, qy, qz)
-				x, y, z = curve.addJacobian(Bx, By, Bz, x, y, z)
+			if byteVal&0x80 == 0x80 {
+				// Q = Q + P
+				curve.addJacobian(qx, qy, qz, px, py, pz, qx,
+					qy, qz)
 			}
-			byte <<= 1
+			byteVal <<= 1
 		}
 	}
 
-	return curve.affineFromJacobian(x, y, z)
+	// Convert the Jacobian coordinate field values back to affine big.Ints.
+	return curve.fieldJacobianToBigAffine(qx, qy, qz)
 }
 
 // ScalarBaseMult returns k*G where G is the base point of the group and k is a
@@ -247,7 +642,7 @@ func initAll() {
 }
 
 func initS256() {
-	// See SEC 2 section 2.7.1
+	// See [SECG] section 2.7.1
 	secp256k1.CurveParams = new(elliptic.CurveParams)
 	secp256k1.P, _ = new(big.Int).SetString("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F", 16)
 	secp256k1.N, _ = new(big.Int).SetString("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141", 16)