diff --git a/src/crypto/elliptic/p521.go b/src/crypto/elliptic/p521.go
index 3d355943ec..463b9f4e39 100644
--- a/src/crypto/elliptic/p521.go
+++ b/src/crypto/elliptic/p521.go
@@ -11,6 +11,7 @@ import (
 
 type p521Curve struct {
 	*CurveParams
+	b *fiat.P521Element
 }
 
 var p521 p521Curve
@@ -25,6 +26,7 @@ func initP521() {
 	p521.Gx, _ = new(big.Int).SetString("c6858e06b70404e9cd9e3ecb662395b4429c648139053fb521f828af606b4d3dbaa14b5e77efe75928fe1dc127a2ffa8de3348b3c1856a429bf97e7e31c2e5bd66", 16)
 	p521.Gy, _ = new(big.Int).SetString("11839296a789a3bc0045c8a5fb42c7d1bd998f54449579b446817afbd17273e662c97ee72995ef42640c550b9013fad0761353c7086a272c24088be94769fd16650", 16)
 	p521.BitSize = 521
+	p521.b = bigIntToFiatP521(p521.B)
 }
 
 func (curve p521Curve) Params() *CurveParams {
@@ -34,7 +36,6 @@ func (curve p521Curve) Params() *CurveParams {
 func (curve p521Curve) IsOnCurve(x, y *big.Int) bool {
 	x1 := bigIntToFiatP521(x)
 	y1 := bigIntToFiatP521(y)
-	b := bigIntToFiatP521(curve.B) // TODO: precompute this value.
 
 	// x³ - 3x + b.
 	x3 := new(fiat.P521Element).Square(x1)
@@ -44,7 +45,7 @@ func (curve p521Curve) IsOnCurve(x, y *big.Int) bool {
 	threeX.Add(threeX, x1)
 
 	x3.Sub(x3, threeX)
-	x3.Add(x3, b)
+	x3.Add(x3, curve.b)
 
 	// y² = x³ - 3x + b
 	y2 := new(fiat.P521Element).Square(y1)
@@ -52,10 +53,20 @@ func (curve p521Curve) IsOnCurve(x, y *big.Int) bool {
 	return x3.Equal(y2) == 1
 }
 
+// p521Point is a P-521 point in projective coordinates, where x = X/Z, y = Y/Z.
 type p521Point struct {
 	x, y, z *fiat.P521Element
 }
 
+// newP521Point returns a new p521Point representing the identity point.
+func newP521Point() *p521Point {
+	return &p521Point{
+		x: new(fiat.P521Element),
+		y: new(fiat.P521Element).One(),
+		z: new(fiat.P521Element),
+	}
+}
+
 func fiatP521ToBigInt(x *fiat.P521Element) *big.Int {
 	xBytes := x.Bytes()
 	for i := range xBytes[:len(xBytes)/2] {
@@ -64,20 +75,16 @@ func fiatP521ToBigInt(x *fiat.P521Element) *big.Int {
 	return new(big.Int).SetBytes(xBytes)
 }
 
-// affineFromJacobian brings a point in Jacobian coordinates back to affine
-// coordinates, with (0, 0) representing infinity by convention. It also goes
-// back to big.Int values to match the exposed API.
-func (curve p521Curve) affineFromJacobian(p *p521Point) (x, y *big.Int) {
+// Affine returns p in affine coordinates, with (0, 0) representing infinity by
+// convention. It also goes back to big.Int values to match the exposed API.
+func (p *p521Point) Affine() (x, y *big.Int) {
 	if p.z.IsZero() == 1 {
 		return new(big.Int), new(big.Int)
 	}
 
 	zinv := new(fiat.P521Element).Invert(p.z)
-	zinvsq := new(fiat.P521Element).Mul(zinv, zinv)
-
-	xx := new(fiat.P521Element).Mul(p.x, zinvsq)
-	zinvsq.Mul(zinvsq, zinv)
-	yy := new(fiat.P521Element).Mul(p.y, zinvsq)
+	xx := new(fiat.P521Element).Mul(p.x, zinv)
+	yy := new(fiat.P521Element).Mul(p.y, zinv)
 
 	return fiatP521ToBigInt(xx), fiatP521ToBigInt(yy)
 }
@@ -96,18 +103,14 @@ func bigIntToFiatP521(x *big.Int) *fiat.P521Element {
 	return x1
 }
 
-// jacobianFromAffine converts (x, y) affine coordinates into (x, y, z) Jacobian
+// newP521PointFromAffine converts (x, y) affine coordinates into (X, Y, Z) projective
 // coordinates. It also converts from big.Int to fiat, which is necessarily a
 // messy and variable-time operation, which we can't avoid due to the exposed API.
-func (curve p521Curve) jacobianFromAffine(x, y *big.Int) *p521Point {
+func newP521PointFromAffine(x, y *big.Int) *p521Point {
 	// (0, 0) is by convention the point at infinity, which can't be represented
-	// in affine coordinates, but is (0, 0, 0) in Jacobian.
+	// in affine coordinates, but is (0, 0, 0) in projective coordinates.
 	if x.Sign() == 0 && y.Sign() == 0 {
-		return &p521Point{
-			x: new(fiat.P521Element),
-			y: new(fiat.P521Element),
-			z: new(fiat.P521Element),
-		}
+		return newP521Point()
 	}
 	return &p521Point{
 		x: bigIntToFiatP521(x),
@@ -117,65 +120,59 @@ func (curve p521Curve) jacobianFromAffine(x, y *big.Int) *p521Point {
 }
 
 func (curve p521Curve) Add(x1, y1, x2, y2 *big.Int) (*big.Int, *big.Int) {
-	p1 := curve.jacobianFromAffine(x1, y1)
-	p2 := curve.jacobianFromAffine(x2, y2)
-	return curve.affineFromJacobian(p1.addJacobian(p1, p2))
+	p1 := newP521PointFromAffine(x1, y1)
+	p2 := newP521PointFromAffine(x2, y2)
+	return p1.Add(p1, p2).Affine()
 }
 
-// addJacobian sets q = p1 + p2, and returns q. The points may overlap.
-func (q *p521Point) addJacobian(p1, p2 *p521Point) *p521Point {
-	// https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl
-	z1IsZero := p1.z.IsZero()
-	z2IsZero := p2.z.IsZero()
+// Add sets q = p1 + p2, and returns q. The points may overlap.
+func (q *p521Point) Add(p1, p2 *p521Point) *p521Point {
+	// Complete addition formula for a = -3 from "Complete addition formulas for
+	// prime order elliptic curves" (https://eprint.iacr.org/2015/1060), §A.2.
 
-	z1z1 := new(fiat.P521Element).Square(p1.z)
-	z2z2 := new(fiat.P521Element).Square(p2.z)
-
-	u1 := new(fiat.P521Element).Mul(p1.x, z2z2)
-	u2 := new(fiat.P521Element).Mul(p2.x, z1z1)
-	h := new(fiat.P521Element).Sub(u2, u1)
-	xEqual := h.IsZero() == 1
-	i := new(fiat.P521Element).Add(h, h)
-	i.Square(i)
-	j := new(fiat.P521Element).Mul(h, i)
-
-	s1 := new(fiat.P521Element).Mul(p1.y, p2.z)
-	s1.Mul(s1, z2z2)
-	s2 := new(fiat.P521Element).Mul(p2.y, p1.z)
-	s2.Mul(s2, z1z1)
-	r := new(fiat.P521Element).Sub(s2, s1)
-	yEqual := r.IsZero() == 1
-	if xEqual && yEqual && z1IsZero == 0 && z2IsZero == 0 {
-		return q.doubleJacobian(p1)
-	}
-	r.Add(r, r)
-	v := new(fiat.P521Element).Mul(u1, i)
-
-	x := new(fiat.P521Element).Set(r)
-	x.Square(x)
-	x.Sub(x, j)
-	x.Sub(x, v)
-	x.Sub(x, v)
-
-	y := new(fiat.P521Element).Set(r)
-	v.Sub(v, x)
-	y.Mul(y, v)
-	s1.Mul(s1, j)
-	s1.Add(s1, s1)
-	y.Sub(y, s1)
-
-	z := new(fiat.P521Element).Add(p1.z, p2.z)
-	z.Square(z)
-	z.Sub(z, z1z1)
-	z.Sub(z, z2z2)
-	z.Mul(z, h)
-
-	x.Select(p2.x, x, z1IsZero)
-	x.Select(p1.x, x, z2IsZero)
-	y.Select(p2.y, y, z1IsZero)
-	y.Select(p1.y, y, z2IsZero)
-	z.Select(p2.z, z, z1IsZero)
-	z.Select(p1.z, z, z2IsZero)
+	t0 := new(fiat.P521Element).Mul(p1.x, p2.x) // t0 := X1 * X2
+	t1 := new(fiat.P521Element).Mul(p1.y, p2.y) // t1 := Y1 * Y2
+	t2 := new(fiat.P521Element).Mul(p1.z, p2.z) // t2 := Z1 * Z2
+	t3 := new(fiat.P521Element).Add(p1.x, p1.y) // t3 := X1 + Y1
+	t4 := new(fiat.P521Element).Add(p2.x, p2.y) // t4 := X2 + Y2
+	t3.Mul(t3, t4)                              // t3 := t3 * t4
+	t4.Add(t0, t1)                              // t4 := t0 + t1
+	t3.Sub(t3, t4)                              // t3 := t3 - t4
+	t4.Add(p1.y, p1.z)                          // t4 := Y1 + Z1
+	x := new(fiat.P521Element).Add(p2.y, p2.z)  // X3 := Y2 + Z2
+	t4.Mul(t4, x)                               // t4 := t4 * X3
+	x.Add(t1, t2)                               // X3 := t1 + t2
+	t4.Sub(t4, x)                               // t4 := t4 - X3
+	x.Add(p1.x, p1.z)                           // X3 := X1 + Z1
+	y := new(fiat.P521Element).Add(p2.x, p2.z)  // Y3 := X2 + Z2
+	x.Mul(x, y)                                 // X3 := X3 * Y3
+	y.Add(t0, t2)                               // Y3 := t0 + t2
+	y.Sub(x, y)                                 // Y3 := X3 - Y3
+	z := new(fiat.P521Element).Mul(p521.b, t2)  // Z3 := b * t2
+	x.Sub(y, z)                                 // X3 := Y3 - Z3
+	z.Add(x, x)                                 // Z3 := X3 + X3
+	x.Add(x, z)                                 // X3 := X3 + Z3
+	z.Sub(t1, x)                                // Z3 := t1 - X3
+	x.Add(t1, x)                                // X3 := t1 + X3
+	y.Mul(p521.b, y)                            // Y3 := b * Y3
+	t1.Add(t2, t2)                              // t1 := t2 + t2
+	t2.Add(t1, t2)                              // t2 := t1 + t2
+	y.Sub(y, t2)                                // Y3 := Y3 - t2
+	y.Sub(y, t0)                                // Y3 := Y3 - t0
+	t1.Add(y, y)                                // t1 := Y3 + Y3
+	y.Add(t1, y)                                // Y3 := t1 + Y3
+	t1.Add(t0, t0)                              // t1 := t0 + t0
+	t0.Add(t1, t0)                              // t0 := t1 + t0
+	t0.Sub(t0, t2)                              // t0 := t0 - t2
+	t1.Mul(t4, y)                               // t1 := t4 * Y3
+	t2.Mul(t0, y)                               // t2 := t0 * Y3
+	y.Mul(x, z)                                 // Y3 := X3 * Z3
+	y.Add(y, t2)                                // Y3 := Y3 + t2
+	x.Mul(t3, x)                                // X3 := t3 * X3
+	x.Sub(x, t1)                                // X3 := X3 - t1
+	z.Mul(t4, z)                                // Z3 := t4 * Z3
+	t1.Mul(t3, t0)                              // t1 := t3 * t0
+	z.Add(z, t1)                                // Z3 := Z3 + t1
 
 	q.x.Set(x)
 	q.y.Set(y)
@@ -184,74 +181,78 @@ func (q *p521Point) addJacobian(p1, p2 *p521Point) *p521Point {
 }
 
 func (curve p521Curve) Double(x1, y1 *big.Int) (*big.Int, *big.Int) {
-	p := curve.jacobianFromAffine(x1, y1)
-	return curve.affineFromJacobian(p.doubleJacobian(p))
+	p := newP521PointFromAffine(x1, y1)
+	return p.Double(p).Affine()
 }
 
-// doubleJacobian sets q = p + p, and returns q. The points may overlap.
-func (q *p521Point) doubleJacobian(p *p521Point) *p521Point {
-	// https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2001-b
-	delta := new(fiat.P521Element).Square(p.z)
-	gamma := new(fiat.P521Element).Square(p.y)
-	alpha := new(fiat.P521Element).Sub(p.x, delta)
-	alpha2 := new(fiat.P521Element).Add(p.x, delta)
-	alpha.Mul(alpha, alpha2)
-	alpha2.Set(alpha)
-	alpha.Add(alpha, alpha)
-	alpha.Add(alpha, alpha2)
+// Double sets q = p + p, and returns q. The points may overlap.
+func (q *p521Point) Double(p *p521Point) *p521Point {
+	// Complete addition formula for a = -3 from "Complete addition formulas for
+	// prime order elliptic curves" (https://eprint.iacr.org/2015/1060), §A.2.
 
-	beta := alpha2.Mul(p.x, gamma)
+	t0 := new(fiat.P521Element).Square(p.x)    // t0 := X ^ 2
+	t1 := new(fiat.P521Element).Square(p.y)    // t1 := Y ^ 2
+	t2 := new(fiat.P521Element).Square(p.z)    // t2 := Z ^ 2
+	t3 := new(fiat.P521Element).Mul(p.x, p.y)  // t3 := X * Y
+	t3.Add(t3, t3)                             // t3 := t3 + t3
+	z := new(fiat.P521Element).Mul(p.x, p.z)   // Z3 := X * Z
+	z.Add(z, z)                                // Z3 := Z3 + Z3
+	y := new(fiat.P521Element).Mul(p521.b, t2) // Y3 := b * t2
+	y.Sub(y, z)                                // Y3 := Y3 - Z3
+	x := new(fiat.P521Element).Add(y, y)       // X3 := Y3 + Y3
+	y.Add(x, y)                                // Y3 := X3 + Y3
+	x.Sub(t1, y)                               // X3 := t1 - Y3
+	y.Add(t1, y)                               // Y3 := t1 + Y3
+	y.Mul(x, y)                                // Y3 := X3 * Y3
+	x.Mul(x, t3)                               // X3 := X3 * t3
+	t3.Add(t2, t2)                             // t3 := t2 + t2
+	t2.Add(t2, t3)                             // t2 := t2 + t3
+	z.Mul(p521.b, z)                           // Z3 := b * Z3
+	z.Sub(z, t2)                               // Z3 := Z3 - t2
+	z.Sub(z, t0)                               // Z3 := Z3 - t0
+	t3.Add(z, z)                               // t3 := Z3 + Z3
+	z.Add(z, t3)                               // Z3 := Z3 + t3
+	t3.Add(t0, t0)                             // t3 := t0 + t0
+	t0.Add(t3, t0)                             // t0 := t3 + t0
+	t0.Sub(t0, t2)                             // t0 := t0 - t2
+	t0.Mul(t0, z)                              // t0 := t0 * Z3
+	y.Add(y, t0)                               // Y3 := Y3 + t0
+	t0.Mul(p.y, p.z)                           // t0 := Y * Z
+	t0.Add(t0, t0)                             // t0 := t0 + t0
+	z.Mul(t0, z)                               // Z3 := t0 * Z3
+	x.Sub(x, z)                                // X3 := X3 - Z3
+	z.Mul(t0, t1)                              // Z3 := t0 * t1
+	z.Add(z, z)                                // Z3 := Z3 + Z3
+	z.Add(z, z)                                // Z3 := Z3 + Z3
 
-	q.x.Square(alpha)
-	beta8 := new(fiat.P521Element).Add(beta, beta)
-	beta8.Add(beta8, beta8)
-	beta8.Add(beta8, beta8)
-	q.x.Sub(q.x, beta8)
-
-	q.z.Add(p.y, p.z)
-	q.z.Square(q.z)
-	q.z.Sub(q.z, gamma)
-	q.z.Sub(q.z, delta)
-
-	beta.Add(beta, beta)
-	beta.Add(beta, beta)
-	beta.Sub(beta, q.x)
-	q.y.Mul(alpha, beta)
-
-	gamma.Square(gamma)
-	gamma.Add(gamma, gamma)
-	gamma.Add(gamma, gamma)
-	gamma.Add(gamma, gamma)
-
-	q.y.Sub(q.y, gamma)
+	q.x.Set(x)
+	q.y.Set(y)
+	q.z.Set(z)
+	return q
+}
 
+// Select sets q to p1 if cond == 1, and to p2 if cond == 0.
+func (q *p521Point) Select(p1, p2 *p521Point, cond int) *p521Point {
+	q.x.Select(p1.x, p2.x, cond)
+	q.y.Select(p1.y, p2.y, cond)
+	q.z.Select(p1.z, p2.z, cond)
 	return q
 }
 
 func (curve p521Curve) ScalarMult(Bx, By *big.Int, scalar []byte) (*big.Int, *big.Int) {
-	B := curve.jacobianFromAffine(Bx, By)
-	p, t := &p521Point{
-		x: new(fiat.P521Element),
-		y: new(fiat.P521Element),
-		z: new(fiat.P521Element),
-	}, &p521Point{
-		x: new(fiat.P521Element),
-		y: new(fiat.P521Element),
-		z: new(fiat.P521Element),
-	}
+	B := newP521PointFromAffine(Bx, By)
+	p, t := newP521Point(), newP521Point()
 
 	for _, byte := range scalar {
 		for bitNum := 0; bitNum < 8; bitNum++ {
-			p.doubleJacobian(p)
+			p.Double(p)
+			t.Add(p, B)
 			bit := (byte >> (7 - bitNum)) & 1
-			t.addJacobian(p, B)
-			p.x.Select(t.x, p.x, int(bit))
-			p.y.Select(t.y, p.y, int(bit))
-			p.z.Select(t.z, p.z, int(bit))
+			p.Select(t, p, int(bit))
 		}
 	}
 
-	return curve.affineFromJacobian(p)
+	return p.Affine()
 }
 
 func (curve p521Curve) ScalarBaseMult(k []byte) (*big.Int, *big.Int) {