crypto/elliptic: move P-521 group logic to internal/nistec

This abstracts the clunky and not constant time math/big elliptic.Curve
compatibility layer away from the pure fiat-backed group logic.

Change-Id: I3b7a7495034d0c569b21c442ae36958763b8b2d0
Reviewed-on: https://go-review.googlesource.com/c/go/+/320074
Trust: Filippo Valsorda <filippo@golang.org>
Run-TryBot: Filippo Valsorda <filippo@golang.org>
TryBot-Result: Go Bot <gobot@golang.org>
Reviewed-by: Julie Qiu <julie@golang.org>

7 files changed, 516 insertions, 250 deletions

diff --git a/src/crypto/elliptic/elliptic.go b/src/crypto/elliptic/elliptic.go
index f072960bfe..cdde0c4e60 100644
--- a/src/crypto/elliptic/elliptic.go
+++ b/src/crypto/elliptic/elliptic.go
@@ -21,9 +21,12 @@ import (
 
 // A Curve represents a short-form Weierstrass curve with a=-3.
 //
-// Note that the point at infinity (0, 0) is not considered on the curve, and
-// although it can be returned by Add, Double, ScalarMult, or ScalarBaseMult, it
-// can't be marshaled or unmarshaled, and IsOnCurve will return false for it.
+// The output of Add, Double, and ScalarMult when the input is not a point on
+// the curve is undefined.
+//
+// Note that the conventional point at infinity (0, 0) is not considered on the
+// curve, although it can be returned by Add, Double, ScalarMult, or
+// ScalarBaseMult (but not Unmarshal or UnmarshalCompressed).
 type Curve interface {
 	// Params returns the parameters for the curve.
 	Params() *CurveParams
diff --git a/src/crypto/elliptic/elliptic_test.go b/src/crypto/elliptic/elliptic_test.go
index 183861a54b..c9744b5a51 100644
--- a/src/crypto/elliptic/elliptic_test.go
+++ b/src/crypto/elliptic/elliptic_test.go
@@ -109,6 +109,15 @@ func testInfinity(t *testing.T, curve Curve) {
 	if curve.IsOnCurve(x, y) {
 		t.Errorf("IsOnCurve(∞) == true")
 	}
+
+	if xx, yy := Unmarshal(curve, Marshal(curve, x, y)); xx != nil || yy != nil {
+		t.Errorf("Unmarshal(Marshal(∞)) did not return an error")
+	}
+	// We don't test UnmarshalCompressed(MarshalCompressed(∞)) because there are
+	// two valid points with x = 0.
+	if xx, yy := Unmarshal(curve, []byte{0x00}); xx != nil || yy != nil {
+		t.Errorf("Unmarshal(∞) did not return an error")
+	}
 }
 
 func TestMarshal(t *testing.T) {
@@ -274,3 +283,29 @@ func BenchmarkScalarMult(b *testing.B) {
 		}
 	})
 }
+
+func BenchmarkMarshalUnmarshal(b *testing.B) {
+	benchmarkAllCurves(b, func(b *testing.B, curve Curve) {
+		_, x, y, _ := GenerateKey(curve, rand.Reader)
+		b.Run("Uncompressed", func(b *testing.B) {
+			b.ReportAllocs()
+			for i := 0; i < b.N; i++ {
+				buf := Marshal(curve, x, y)
+				xx, yy := Unmarshal(curve, buf)
+				if xx.Cmp(x) != 0 || yy.Cmp(y) != 0 {
+					b.Error("Unmarshal output different from Marshal input")
+				}
+			}
+		})
+		b.Run("Compressed", func(b *testing.B) {
+			b.ReportAllocs()
+			for i := 0; i < b.N; i++ {
+				buf := Marshal(curve, x, y)
+				xx, yy := Unmarshal(curve, buf)
+				if xx.Cmp(x) != 0 || yy.Cmp(y) != 0 {
+					b.Error("Unmarshal output different from Marshal input")
+				}
+			}
+		})
+	})
+}
diff --git a/src/crypto/elliptic/internal/fiat/p521.go b/src/crypto/elliptic/internal/fiat/p521.go
index dc677327e6..647c3f914f 100644
--- a/src/crypto/elliptic/internal/fiat/p521.go
+++ b/src/crypto/elliptic/internal/fiat/p521.go
@@ -53,28 +53,40 @@ func (e *P521Element) Set(t *P521Element) *P521Element {
 	return e
 }
 
-// Bytes returns the 66-byte little-endian encoding of e.
+// Bytes returns the 66-byte big-endian encoding of e.
 func (e *P521Element) Bytes() []byte {
-	// This function must be inlined to move the allocation to the parent and
-	// save it from escaping to the heap.
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
 	var out [66]byte
-	p521ToBytes(&out, &e.x)
+	return e.bytes(&out)
+}
+
+func (e *P521Element) bytes(out *[66]byte) []byte {
+	p521ToBytes(out, &e.x)
+	invertEndianness(out[:])
 	return out[:]
 }
 
-// SetBytes sets e = v, where v is a little-endian 66-byte encoding, and returns
+// SetBytes sets e = v, where v is a big-endian 66-byte encoding, and returns
 // e. If v is not 66 bytes or it encodes a value higher than 2^521 - 1, SetBytes
 // returns nil and an error, and e is unchanged.
 func (e *P521Element) SetBytes(v []byte) (*P521Element, error) {
-	if len(v) != 66 || v[65] > 1 {
+	if len(v) != 66 || v[0] > 1 {
 		return nil, errors.New("invalid P-521 field encoding")
 	}
 	var in [66]byte
 	copy(in[:], v)
+	invertEndianness(in[:])
 	p521FromBytes(&e.x, &in)
 	return e, nil
 }
 
+func invertEndianness(v []byte) {
+	for i := 0; i < len(v)/2; i++ {
+		v[i], v[len(v)-1-i] = v[len(v)-1-i], v[i]
+	}
+}
+
 // Add sets e = t1 + t2, and returns e.
 func (e *P521Element) Add(t1, t2 *P521Element) *P521Element {
 	p521Add(&e.x, &t1.x, &t2.x)
diff --git a/src/crypto/elliptic/internal/fiat/p521_test.go b/src/crypto/elliptic/internal/fiat/p521_test.go
index 661bde397e..2b374faa27 100644
--- a/src/crypto/elliptic/internal/fiat/p521_test.go
+++ b/src/crypto/elliptic/internal/fiat/p521_test.go
@@ -15,7 +15,7 @@ func p521Random(t *testing.T) *fiat.P521Element {
 	if _, err := rand.Read(buf); err != nil {
 		t.Fatal(err)
 	}
-	buf[65] &= 1
+	buf[0] &= 1
 	e, err := new(fiat.P521Element).SetBytes(buf)
 	if err != nil {
 		t.Fatal(err)
diff --git a/src/crypto/elliptic/internal/nistec/p521.go b/src/crypto/elliptic/internal/nistec/p521.go
new file mode 100644
index 0000000000..e5b4e46d4b
--- /dev/null
+++ b/src/crypto/elliptic/internal/nistec/p521.go
@@ -0,0 +1,306 @@
+// Copyright 2021 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package nistec implements the NIST P elliptic curves from FIPS 186-4.
+//
+// This package uses fiat-crypto for its backend field arithmetic (not math/big)
+// and exposes constant-time, heap allocation-free, byte slice-based safe APIs.
+// Group operations use modern and safe complete addition formulas. The point at
+// infinity is handled and encoded according to SEC 1, Version 2.0, and invalid
+// curve points can't be represented.
+package nistec
+
+import (
+	"crypto/elliptic/internal/fiat"
+	"crypto/subtle"
+	"errors"
+)
+
+var p521B, _ = new(fiat.P521Element).SetBytes([]byte{
+	0x00, 0x51, 0x95, 0x3e, 0xb9, 0x61, 0x8e, 0x1c, 0x9a, 0x1f, 0x92, 0x9a,
+	0x21, 0xa0, 0xb6, 0x85, 0x40, 0xee, 0xa2, 0xda, 0x72, 0x5b, 0x99, 0xb3,
+	0x15, 0xf3, 0xb8, 0xb4, 0x89, 0x91, 0x8e, 0xf1, 0x09, 0xe1, 0x56, 0x19,
+	0x39, 0x51, 0xec, 0x7e, 0x93, 0x7b, 0x16, 0x52, 0xc0, 0xbd, 0x3b, 0xb1,
+	0xbf, 0x07, 0x35, 0x73, 0xdf, 0x88, 0x3d, 0x2c, 0x34, 0xf1, 0xef, 0x45,
+	0x1f, 0xd4, 0x6b, 0x50, 0x3f, 0x00})
+
+var p521G, _ = NewP521Point().SetBytes([]byte{0x04,
+	0x00, 0xc6, 0x85, 0x8e, 0x06, 0xb7, 0x04, 0x04, 0xe9, 0xcd, 0x9e, 0x3e,
+	0xcb, 0x66, 0x23, 0x95, 0xb4, 0x42, 0x9c, 0x64, 0x81, 0x39, 0x05, 0x3f,
+	0xb5, 0x21, 0xf8, 0x28, 0xaf, 0x60, 0x6b, 0x4d, 0x3d, 0xba, 0xa1, 0x4b,
+	0x5e, 0x77, 0xef, 0xe7, 0x59, 0x28, 0xfe, 0x1d, 0xc1, 0x27, 0xa2, 0xff,
+	0xa8, 0xde, 0x33, 0x48, 0xb3, 0xc1, 0x85, 0x6a, 0x42, 0x9b, 0xf9, 0x7e,
+	0x7e, 0x31, 0xc2, 0xe5, 0xbd, 0x66, 0x01, 0x18, 0x39, 0x29, 0x6a, 0x78,
+	0x9a, 0x3b, 0xc0, 0x04, 0x5c, 0x8a, 0x5f, 0xb4, 0x2c, 0x7d, 0x1b, 0xd9,
+	0x98, 0xf5, 0x44, 0x49, 0x57, 0x9b, 0x44, 0x68, 0x17, 0xaf, 0xbd, 0x17,
+	0x27, 0x3e, 0x66, 0x2c, 0x97, 0xee, 0x72, 0x99, 0x5e, 0xf4, 0x26, 0x40,
+	0xc5, 0x50, 0xb9, 0x01, 0x3f, 0xad, 0x07, 0x61, 0x35, 0x3c, 0x70, 0x86,
+	0xa2, 0x72, 0xc2, 0x40, 0x88, 0xbe, 0x94, 0x76, 0x9f, 0xd1, 0x66, 0x50})
+
+const p521ElementLength = 66
+
+// P521Point is a P-521 point. The zero value is NOT valid.
+type P521Point struct {
+	// The point is represented in projective coordinates (X:Y:Z),
+	// where x = X/Z and y = Y/Z.
+	x, y, z *fiat.P521Element
+}
+
+// NewP521Point returns a new P521Point representing the point at infinity point.
+func NewP521Point() *P521Point {
+	return &P521Point{
+		x: new(fiat.P521Element),
+		y: new(fiat.P521Element).One(),
+		z: new(fiat.P521Element),
+	}
+}
+
+// NewP521Generator returns a new P521Point set to the canonical generator.
+func NewP521Generator() *P521Point {
+	return NewP521Point().Set(p521G)
+}
+
+// Set sets p = q and returns p.
+func (p *P521Point) Set(q *P521Point) *P521Point {
+	p.x.Set(q.x)
+	p.y.Set(q.y)
+	p.z.Set(q.z)
+	return p
+}
+
+// SetBytes sets p to the compressed, uncompressed, or infinity value encoded in
+// b, as specified in SEC 1, Version 2.0, Section 2.3.4. If the point is not on
+// the curve, it returns nil and an error, and the receiver is unchanged.
+// Otherwise, it returns p.
+func (p *P521Point) SetBytes(b []byte) (*P521Point, error) {
+	switch {
+	// Point at infinity.
+	case len(b) == 1 && b[0] == 0:
+		return p.Set(NewP521Point()), nil
+
+	// Uncompressed form.
+	case len(b) == 1+2*p521ElementLength && b[0] == 4:
+		x, err := new(fiat.P521Element).SetBytes(b[1 : 1+p521ElementLength])
+		if err != nil {
+			return nil, err
+		}
+		y, err := new(fiat.P521Element).SetBytes(b[1+p521ElementLength:])
+		if err != nil {
+			return nil, err
+		}
+		if err := p521CheckOnCurve(x, y); err != nil {
+			return nil, err
+		}
+		p.x.Set(x)
+		p.y.Set(y)
+		p.z.One()
+		return p, nil
+
+	// Compressed form
+	case len(b) == 1+p521ElementLength && b[0] == 0:
+		return nil, errors.New("unimplemented") // TODO(filippo)
+
+	default:
+		return nil, errors.New("invalid P521 point encoding")
+	}
+}
+
+func p521CheckOnCurve(x, y *fiat.P521Element) error {
+	// x³ - 3x + b.
+	x3 := new(fiat.P521Element).Square(x)
+	x3.Mul(x3, x)
+
+	threeX := new(fiat.P521Element).Add(x, x)
+	threeX.Add(threeX, x)
+
+	x3.Sub(x3, threeX)
+	x3.Add(x3, p521B)
+
+	// y² = x³ - 3x + b
+	y2 := new(fiat.P521Element).Square(y)
+
+	if x3.Equal(y2) != 1 {
+		return errors.New("P521 point not on curve")
+	}
+	return nil
+}
+
+// Bytes returns the uncompressed or infinity encoding of p, as specified in
+// SEC 1, Version 2.0, Section 2.3.3. Note that the encoding of the point at
+// infinity is shorter than all other encodings.
+func (p *P521Point) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [133]byte
+	return p.bytes(&out)
+}
+
+func (p *P521Point) bytes(out *[133]byte) []byte {
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new(fiat.P521Element).Invert(p.z)
+	xx := new(fiat.P521Element).Mul(p.x, zinv)
+	yy := new(fiat.P521Element).Mul(p.y, zinv)
+
+	buf := append(out[:0], 4)
+	buf = append(buf, xx.Bytes()...)
+	buf = append(buf, yy.Bytes()...)
+	return buf
+}
+
+// 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.
+
+	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
+	x3 := new(fiat.P521Element).Add(p2.y, p2.z) // X3 := Y2 + Z2
+	t4.Mul(t4, x3)                              // t4 := t4 * X3
+	x3.Add(t1, t2)                              // X3 := t1 + t2
+	t4.Sub(t4, x3)                              // t4 := t4 - X3
+	x3.Add(p1.x, p1.z)                          // X3 := X1 + Z1
+	y3 := new(fiat.P521Element).Add(p2.x, p2.z) // Y3 := X2 + Z2
+	x3.Mul(x3, y3)                              // X3 := X3 * Y3
+	y3.Add(t0, t2)                              // Y3 := t0 + t2
+	y3.Sub(x3, y3)                              // Y3 := X3 - Y3
+	z3 := new(fiat.P521Element).Mul(p521B, t2)  // Z3 := b * t2
+	x3.Sub(y3, z3)                              // X3 := Y3 - Z3
+	z3.Add(x3, x3)                              // Z3 := X3 + X3
+	x3.Add(x3, z3)                              // X3 := X3 + Z3
+	z3.Sub(t1, x3)                              // Z3 := t1 - X3
+	x3.Add(t1, x3)                              // X3 := t1 + X3
+	y3.Mul(p521B, y3)                           // Y3 := b * Y3
+	t1.Add(t2, t2)                              // t1 := t2 + t2
+	t2.Add(t1, t2)                              // t2 := t1 + t2
+	y3.Sub(y3, t2)                              // Y3 := Y3 - t2
+	y3.Sub(y3, t0)                              // Y3 := Y3 - t0
+	t1.Add(y3, y3)                              // t1 := Y3 + Y3
+	y3.Add(t1, y3)                              // 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, y3)                              // t1 := t4 * Y3
+	t2.Mul(t0, y3)                              // t2 := t0 * Y3
+	y3.Mul(x3, z3)                              // Y3 := X3 * Z3
+	y3.Add(y3, t2)                              // Y3 := Y3 + t2
+	x3.Mul(t3, x3)                              // X3 := t3 * X3
+	x3.Sub(x3, t1)                              // X3 := X3 - t1
+	z3.Mul(t4, z3)                              // Z3 := t4 * Z3
+	t1.Mul(t3, t0)                              // t1 := t3 * t0
+	z3.Add(z3, t1)                              // Z3 := Z3 + t1
+
+	q.x.Set(x3)
+	q.y.Set(y3)
+	q.z.Set(z3)
+	return q
+}
+
+// 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.
+
+	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
+	z3 := new(fiat.P521Element).Mul(p.x, p.z)  // Z3 := X * Z
+	z3.Add(z3, z3)                             // Z3 := Z3 + Z3
+	y3 := new(fiat.P521Element).Mul(p521B, t2) // Y3 := b * t2
+	y3.Sub(y3, z3)                             // Y3 := Y3 - Z3
+	x3 := new(fiat.P521Element).Add(y3, y3)    // X3 := Y3 + Y3
+	y3.Add(x3, y3)                             // Y3 := X3 + Y3
+	x3.Sub(t1, y3)                             // X3 := t1 - Y3
+	y3.Add(t1, y3)                             // Y3 := t1 + Y3
+	y3.Mul(x3, y3)                             // Y3 := X3 * Y3
+	x3.Mul(x3, t3)                             // X3 := X3 * t3
+	t3.Add(t2, t2)                             // t3 := t2 + t2
+	t2.Add(t2, t3)                             // t2 := t2 + t3
+	z3.Mul(p521B, z3)                          // Z3 := b * Z3
+	z3.Sub(z3, t2)                             // Z3 := Z3 - t2
+	z3.Sub(z3, t0)                             // Z3 := Z3 - t0
+	t3.Add(z3, z3)                             // t3 := Z3 + Z3
+	z3.Add(z3, 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, z3)                             // t0 := t0 * Z3
+	y3.Add(y3, t0)                             // Y3 := Y3 + t0
+	t0.Mul(p.y, p.z)                           // t0 := Y * Z
+	t0.Add(t0, t0)                             // t0 := t0 + t0
+	z3.Mul(t0, z3)                             // Z3 := t0 * Z3
+	x3.Sub(x3, z3)                             // X3 := X3 - Z3
+	z3.Mul(t0, t1)                             // Z3 := t0 * t1
+	z3.Add(z3, z3)                             // Z3 := Z3 + Z3
+	z3.Add(z3, z3)                             // Z3 := Z3 + Z3
+
+	q.x.Set(x3)
+	q.y.Set(y3)
+	q.z.Set(z3)
+	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
+}
+
+// ScalarMult sets p = scalar * q, and returns p.
+func (p *P521Point) ScalarMult(q *P521Point, scalar []byte) *P521Point {
+	// table holds the first 16 multiples of q. The explicit newP521Point calls
+	// get inlined, letting the allocations live on the stack.
+	var table = [16]*P521Point{
+		NewP521Point(), NewP521Point(), NewP521Point(), NewP521Point(),
+		NewP521Point(), NewP521Point(), NewP521Point(), NewP521Point(),
+		NewP521Point(), NewP521Point(), NewP521Point(), NewP521Point(),
+		NewP521Point(), NewP521Point(), NewP521Point(), NewP521Point(),
+	}
+	for i := 1; i < 16; i++ {
+		table[i].Add(table[i-1], q)
+	}
+
+	// Instead of doing the classic double-and-add chain, we do it with a
+	// four-bit window: we double four times, and then add [0-15]P.
+	t := NewP521Point()
+	p.Set(NewP521Point())
+	for _, byte := range scalar {
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+
+		for i := uint8(0); i < 16; i++ {
+			cond := subtle.ConstantTimeByteEq(byte>>4, i)
+			t.Select(table[i], t, cond)
+		}
+		p.Add(p, t)
+
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+
+		for i := uint8(0); i < 16; i++ {
+			cond := subtle.ConstantTimeByteEq(byte&0b1111, i)
+			t.Select(table[i], t, cond)
+		}
+		p.Add(p, t)
+	}
+
+	return p
+}
diff --git a/src/crypto/elliptic/internal/nistec/p521_test.go b/src/crypto/elliptic/internal/nistec/p521_test.go
new file mode 100644
index 0000000000..e62c1cbf29
--- /dev/null
+++ b/src/crypto/elliptic/internal/nistec/p521_test.go
@@ -0,0 +1,44 @@
+// Copyright 2021 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package nistec_test
+
+import (
+	"crypto/elliptic/internal/nistec"
+	"math/rand"
+	"os"
+	"strings"
+	"testing"
+)
+
+func TestP521Allocations(t *testing.T) {
+	if strings.HasSuffix(os.Getenv("GO_BUILDER_NAME"), "-noopt") {
+		t.Skip("skipping allocations test without relevant optimizations")
+	}
+	if allocs := testing.AllocsPerRun(100, func() {
+		p := nistec.NewP521Generator()
+		scalar := make([]byte, 66)
+		rand.Read(scalar)
+		p.ScalarMult(p, scalar)
+		out := p.Bytes()
+		if _, err := p.SetBytes(out); err != nil {
+			t.Fatal(err)
+		}
+	}); allocs > 0 {
+		t.Errorf("expected zero allocations, got %0.1f", allocs)
+	}
+}
+
+func BenchmarkScalarMult(b *testing.B) {
+	b.Run("P521", func(b *testing.B) {
+		scalar := make([]byte, 66)
+		rand.Read(scalar)
+		p := nistec.NewP521Generator()
+		b.ReportAllocs()
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			p.ScalarMult(p, scalar)
+		}
+	})
+}
diff --git a/src/crypto/elliptic/p521.go b/src/crypto/elliptic/p521.go
index 569a58c6f3..4cc5f86d6d 100644
--- a/src/crypto/elliptic/p521.go
+++ b/src/crypto/elliptic/p521.go
@@ -5,286 +5,152 @@
 package elliptic
 
 import (
-	"crypto/elliptic/internal/fiat"
-	"crypto/subtle"
+	"crypto/elliptic/internal/nistec"
+	"crypto/rand"
 	"math/big"
 )
 
+// p521Curve is a Curve implementation based on nistec.P521Point.
+//
+// It's a wrapper that exposes the big.Int-based Curve interface and encodes the
+// legacy idiosyncrasies it requires, such as invalid and infinity point
+// handling.
+//
+// To interact with the nistec package, points are encoded into and decoded from
+// properly formatted byte slices. All big.Int use is limited to this package.
+// Encoding and decoding is 1/1000th of the runtime of a scalar multiplication,
+// so the overhead is acceptable.
 type p521Curve struct {
-	*CurveParams
-	b *fiat.P521Element
+	params *CurveParams
 }
 
 var p521 p521Curve
-var p521Params *CurveParams
+var _ Curve = p521
 
 func initP521() {
-	// See FIPS 186-3, section D.2.5
-	p521.CurveParams = &CurveParams{Name: "P-521"}
-	p521.P, _ = new(big.Int).SetString("6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151", 10)
-	p521.N, _ = new(big.Int).SetString("6864797660130609714981900799081393217269435300143305409394463459185543183397655394245057746333217197532963996371363321113864768612440380340372808892707005449", 10)
-	p521.B, _ = new(big.Int).SetString("051953eb9618e1c9a1f929a21a0b68540eea2da725b99b315f3b8b489918ef109e156193951ec7e937b1652c0bd3bb1bf073573df883d2c34f1ef451fd46b503f00", 16)
-	p521.Gx, _ = new(big.Int).SetString("c6858e06b70404e9cd9e3ecb662395b4429c648139053fb521f828af606b4d3dbaa14b5e77efe75928fe1dc127a2ffa8de3348b3c1856a429bf97e7e31c2e5bd66", 16)
-	p521.Gy, _ = new(big.Int).SetString("11839296a789a3bc0045c8a5fb42c7d1bd998f54449579b446817afbd17273e662c97ee72995ef42640c550b9013fad0761353c7086a272c24088be94769fd16650", 16)
-	p521.BitSize = 521
-	p521.b = bigIntToFiatP521(p521.B)
+	p521.params = &CurveParams{
+		Name:    "P-521",
+		BitSize: 521,
+		// FIPS 186-4, section D.1.2.5
+		P: bigFromDecimal("68647976601306097149819007990813932172694353001433" +
+			"0540939446345918554318339765605212255964066145455497729631139148" +
+			"0858037121987999716643812574028291115057151"),
+		N: bigFromDecimal("68647976601306097149819007990813932172694353001433" +
+			"0540939446345918554318339765539424505774633321719753296399637136" +
+			"3321113864768612440380340372808892707005449"),
+		B: bigFromHex("0051953eb9618e1c9a1f929a21a0b68540eea2da725b99b315f3b8" +
+			"b489918ef109e156193951ec7e937b1652c0bd3bb1bf073573df883d2c34f1ef" +
+			"451fd46b503f00"),
+		Gx: bigFromHex("00c6858e06b70404e9cd9e3ecb662395b4429c648139053fb521f8" +
+			"28af606b4d3dbaa14b5e77efe75928fe1dc127a2ffa8de3348b3c1856a429bf9" +
+			"7e7e31c2e5bd66"),
+		Gy: bigFromHex("011839296a789a3bc0045c8a5fb42c7d1bd998f54449579b446817" +
+			"afbd17273e662c97ee72995ef42640c550b9013fad0761353c7086a272c24088" +
+			"be94769fd16650"),
+	}
 }
 
 func (curve p521Curve) Params() *CurveParams {
-	return curve.CurveParams
+	return curve.params
 }
 
 func (curve p521Curve) IsOnCurve(x, y *big.Int) bool {
-	x1 := bigIntToFiatP521(x)
-	y1 := bigIntToFiatP521(y)
-
-	// x³ - 3x + b.
-	x3 := new(fiat.P521Element).Square(x1)
-	x3.Mul(x3, x1)
-
-	threeX := new(fiat.P521Element).Add(x1, x1)
-	threeX.Add(threeX, x1)
-
-	x3.Sub(x3, threeX)
-	x3.Add(x3, curve.b)
-
-	// y² = x³ - 3x + b
-	y2 := new(fiat.P521Element).Square(y1)
-
-	return x3.Equal(y2) == 1
+	// IsOnCurve is documented to reject (0, 0), so we don't use
+	// p521PointFromAffine, but let SetBytes reject the invalid Marshal output.
+	_, err := nistec.NewP521Point().SetBytes(Marshal(curve, x, y))
+	return err == nil
 }
 
-// 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 p521PointFromAffine(x, y *big.Int) (p *nistec.P521Point, ok bool) {
+	// (0, 0) is by convention the point at infinity, which can't be represented
+	// in affine coordinates. Marshal incorrectly encodes it as an uncompressed
+	// point, which SetBytes correctly rejects. See Issue 37294.
+	if x.Sign() == 0 && y.Sign() == 0 {
+		return nistec.NewP521Point(), true
 	}
-}
-
-func fiatP521ToBigInt(x *fiat.P521Element) *big.Int {
-	xBytes := x.Bytes()
-	for i := range xBytes[:len(xBytes)/2] {
-		xBytes[i], xBytes[len(xBytes)-i-1] = xBytes[len(xBytes)-i-1], xBytes[i]
+	p, err := nistec.NewP521Point().SetBytes(Marshal(P521(), x, y))
+	if err != nil {
+		return nil, false
 	}
-	return new(big.Int).SetBytes(xBytes)
+	return p, true
 }
 
-// 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 {
+func p521PointToAffine(p *nistec.P521Point) (x, y *big.Int) {
+	out := p.Bytes()
+	if len(out) == 1 && out[0] == 0 {
+		// This is the correct encoding of the point at infinity, which
+		// Unmarshal does not support. See Issue 37294.
 		return new(big.Int), new(big.Int)
 	}
-
-	zinv := new(fiat.P521Element).Invert(p.z)
-	xx := new(fiat.P521Element).Mul(p.x, zinv)
-	yy := new(fiat.P521Element).Mul(p.y, zinv)
-
-	return fiatP521ToBigInt(xx), fiatP521ToBigInt(yy)
-}
-
-func bigIntToFiatP521(x *big.Int) *fiat.P521Element {
-	xBytes := new(big.Int).Mod(x, p521.P).FillBytes(make([]byte, 66))
-	for i := range xBytes[:len(xBytes)/2] {
-		xBytes[i], xBytes[len(xBytes)-i-1] = xBytes[len(xBytes)-i-1], xBytes[i]
+	x, y = Unmarshal(P521(), out)
+	if x == nil {
+		panic("crypto/elliptic: internal error: Unmarshal rejected a valid point encoding")
 	}
-	x1, err := new(fiat.P521Element).SetBytes(xBytes)
+	return x, y
+}
+
+// p521RandomPoint returns a random point on the curve. It's used when Add,
+// Double, or ScalarMult are fed a point not on the curve, which is undefined
+// behavior. Originally, we used to do the math on it anyway (which allows
+// invalid curve attacks) and relied on the caller and Unmarshal to avoid this
+// happening in the first place. Now, we just can't construct a nistec.P521Point
+// for an invalid pair of coordinates, because that API is safer. If we panic,
+// we risk introducing a DoS. If we return nil, we risk a panic. If we return
+// the input, ecdsa.Verify might fail open. The safest course seems to be to
+// return a valid, random point, which hopefully won't help the attacker.
+func p521RandomPoint() (x, y *big.Int) {
+	_, x, y, err := GenerateKey(P521(), rand.Reader)
 	if err != nil {
-		// The input is reduced modulo P and encoded in a fixed size bytes
-		// slice, this should be impossible.
-		panic("internal error: bigIntToFiatP521")
+		panic("crypto/elliptic: failed to generate random point")
 	}
-	return x1
+	return x, y
 }
 
-// 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 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 projective coordinates.
-	if x.Sign() == 0 && y.Sign() == 0 {
-		return newP521Point()
+func (curve p521Curve) Add(x1, y1, x2, y2 *big.Int) (*big.Int, *big.Int) {
+	p1, ok := p521PointFromAffine(x1, y1)
+	if !ok {
+		return p521RandomPoint()
 	}
-	return &p521Point{
-		x: bigIntToFiatP521(x),
-		y: bigIntToFiatP521(y),
-		z: new(fiat.P521Element).One(),
+	p2, ok := p521PointFromAffine(x2, y2)
+	if !ok {
+		return p521RandomPoint()
 	}
-}
-
-func (curve p521Curve) Add(x1, y1, x2, y2 *big.Int) (*big.Int, *big.Int) {
-	p1 := newP521PointFromAffine(x1, y1)
-	p2 := newP521PointFromAffine(x2, y2)
-	return p1.Add(p1, p2).Affine()
-}
-
-// 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.
-
-	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)
-	q.z.Set(z)
-	return q
+	return p521PointToAffine(p1.Add(p1, p2))
 }
 
 func (curve p521Curve) Double(x1, y1 *big.Int) (*big.Int, *big.Int) {
-	p := newP521PointFromAffine(x1, y1)
-	return p.Double(p).Affine()
-}
-
-// 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.
-
-	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.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
+	p, ok := p521PointFromAffine(x1, y1)
+	if !ok {
+		return p521RandomPoint()
+	}
+	return p521PointToAffine(p.Double(p))
 }
 
 func (curve p521Curve) ScalarMult(Bx, By *big.Int, scalar []byte) (*big.Int, *big.Int) {
-	B := newP521PointFromAffine(Bx, By)
-	p, t := newP521Point(), newP521Point()
-
-	// table holds the first 16 multiples of q. The explicit newP521Point calls
-	// get inlined, letting the allocations live on the stack.
-	var table = [16]*p521Point{
-		newP521Point(), newP521Point(), newP521Point(), newP521Point(),
-		newP521Point(), newP521Point(), newP521Point(), newP521Point(),
-		newP521Point(), newP521Point(), newP521Point(), newP521Point(),
-		newP521Point(), newP521Point(), newP521Point(), newP521Point(),
+	p, ok := p521PointFromAffine(Bx, By)
+	if !ok {
+		return p521RandomPoint()
 	}
-	for i := 1; i < 16; i++ {
-		table[i].Add(table[i-1], B)
-	}
-
-	// Instead of doing the classic double-and-add chain, we do it with a
-	// four-bit window: we double four times, and then add [0-15]P.
-	for _, byte := range scalar {
-		p.Double(p)
-		p.Double(p)
-		p.Double(p)
-		p.Double(p)
-
-		for i := uint8(0); i < 16; i++ {
-			cond := subtle.ConstantTimeByteEq(byte>>4, i)
-			t.Select(table[i], t, cond)
-		}
-		p.Add(p, t)
+	return p521PointToAffine(p.ScalarMult(p, scalar))
+}
 
-		p.Double(p)
-		p.Double(p)
-		p.Double(p)
-		p.Double(p)
+func (curve p521Curve) ScalarBaseMult(scalar []byte) (*big.Int, *big.Int) {
+	p := nistec.NewP521Generator()
+	return p521PointToAffine(p.ScalarMult(p, scalar))
+}
 
-		for i := uint8(0); i < 16; i++ {
-			cond := subtle.ConstantTimeByteEq(byte&0b1111, i)
-			t.Select(table[i], t, cond)
-		}
-		p.Add(p, t)
+func bigFromDecimal(s string) *big.Int {
+	b, ok := new(big.Int).SetString(s, 10)
+	if !ok {
+		panic("invalid encoding")
 	}
-
-	return p.Affine()
+	return b
 }
 
-func (curve p521Curve) ScalarBaseMult(k []byte) (*big.Int, *big.Int) {
-	return curve.ScalarMult(curve.Gx, curve.Gy, k)
+func bigFromHex(s string) *big.Int {
+	b, ok := new(big.Int).SetString(s, 16)
+	if !ok {
+		panic("invalid encoding")
+	}
+	return b
 }