Add initial nistec project files

LICENSE

@@ -0,0 +1,27 @@
+Copyright 2009 The Go Authors.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+   * Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+   * Redistributions in binary form must reproduce the above
+copyright notice, this list of conditions and the following disclaimer
+in the documentation and/or other materials provided with the
+distribution.
+   * Neither the name of Google LLC nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

README.md

@@ -0,0 +1,24 @@
+# sources.truenas.cloud/code/nistec
+
+```
+import "sources.truenas.cloud/code/nistec"
+```
+
+This package implements the NIST P elliptic curves, according to FIPS 186-4
+and SEC 1, Version 2.0, exposing the necessary APIs to build a wide array of
+higher-level primitives.
+
+It's an exported version of `crypto/internal/fips140/nistec` in the standard library,
+which powers `crypto/elliptic`, `crypto/ecdsa`, and `crypto/ecdh`.
+The git history has been preserved, and new upstream changes are applied periodically.
+
+This package uses fiat-crypto or specialized assembly and Go code 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 where possible. The point at infinity is
+handled and encoded according to SEC 1, Version 2.0, and invalid curve points
+can't be represented. This makes it particularly suitable to be used as a
+prime order group implementation.
+
+Use the `purego` build tag to exclude the assembly and rely entirely on formally
+verified fiat-crypto arithmetic and complete addition formulas.

_asm/go.mod

@@ -0,0 +1,11 @@
+module crypto/internal/fips140/nistec/_asm
+
+go 1.24
+
+require github.com/mmcloughlin/avo v0.6.0
+
+require (
+	golang.org/x/mod v0.20.0 // indirect
+	golang.org/x/sync v0.8.0 // indirect
+	golang.org/x/tools v0.24.0 // indirect
+)

_asm/go.sum

@@ -0,0 +1,8 @@
+github.com/mmcloughlin/avo v0.6.0 h1:QH6FU8SKoTLaVs80GA8TJuLNkUYl4VokHKlPhVDg4YY=
+github.com/mmcloughlin/avo v0.6.0/go.mod h1:8CoAGaCSYXtCPR+8y18Y9aB/kxb8JSS6FRI7mSkvD+8=
+golang.org/x/mod v0.20.0 h1:utOm6MM3R3dnawAiJgn0y+xvuYRsm1RKM/4giyfDgV0=
+golang.org/x/mod v0.20.0/go.mod h1:hTbmBsO62+eylJbnUtE2MGJUyE7QWk4xUqPFrRgJ+7c=
+golang.org/x/sync v0.8.0 h1:3NFvSEYkUoMifnESzZl15y791HH1qU2xm6eCJU5ZPXQ=
+golang.org/x/sync v0.8.0/go.mod h1:Czt+wKu1gCyEFDUtn0jG5QVvpJ6rzVqr5aXyt9drQfk=
+golang.org/x/tools v0.24.0 h1:J1shsA93PJUEVaUSaay7UXAyE8aimq3GW0pjlolpa24=
+golang.org/x/tools v0.24.0/go.mod h1:YhNqVBIfWHdzvTLs0d8LCuMhkKUgSUKldakyV7W/WDQ=

benchmark_test.go

@@ -0,0 +1,62 @@
+// Copyright 2024 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/rand"
+	"testing"
+
+	"sources.truenas.cloud/code/nistec"
+)
+
+func BenchmarkScalarMult(b *testing.B) {
+	b.Run("P224", func(b *testing.B) {
+		benchmarkScalarMult(b, nistec.NewP224Point().SetGenerator(), 28)
+	})
+	b.Run("P256", func(b *testing.B) {
+		benchmarkScalarMult(b, nistec.NewP256Point().SetGenerator(), 32)
+	})
+	b.Run("P384", func(b *testing.B) {
+		benchmarkScalarMult(b, nistec.NewP384Point().SetGenerator(), 48)
+	})
+	b.Run("P521", func(b *testing.B) {
+		benchmarkScalarMult(b, nistec.NewP521Point().SetGenerator(), 66)
+	})
+}
+
+func benchmarkScalarMult[P nistPoint[P]](b *testing.B, p P, scalarSize int) {
+	scalar := make([]byte, scalarSize)
+	rand.Read(scalar)
+	b.ReportAllocs()
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		p.ScalarMult(p, scalar)
+	}
+}
+
+func BenchmarkScalarBaseMult(b *testing.B) {
+	b.Run("P224", func(b *testing.B) {
+		benchmarkScalarBaseMult(b, nistec.NewP224Point().SetGenerator(), 28)
+	})
+	b.Run("P256", func(b *testing.B) {
+		benchmarkScalarBaseMult(b, nistec.NewP256Point().SetGenerator(), 32)
+	})
+	b.Run("P384", func(b *testing.B) {
+		benchmarkScalarBaseMult(b, nistec.NewP384Point().SetGenerator(), 48)
+	})
+	b.Run("P521", func(b *testing.B) {
+		benchmarkScalarBaseMult(b, nistec.NewP521Point().SetGenerator(), 66)
+	})
+}
+
+func benchmarkScalarBaseMult[P nistPoint[P]](b *testing.B, p P, scalarSize int) {
+	scalar := make([]byte, scalarSize)
+	rand.Read(scalar)
+	b.ReportAllocs()
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		p.ScalarBaseMult(scalar)
+	}
+}

extra.go

@@ -0,0 +1,85 @@
+// Copyright 2023 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
+
+import "sources.truenas.cloud/code/nistec/internal/fiat"
+
+// Negate sets p = -q and returns p.
+func (p *P224Point) Negate(q *P224Point) *P224Point {
+	p.x.Set(q.x)
+	p.y.Sub(new(fiat.P224Element), q.y)
+	p.z.Set(q.z)
+	return p
+}
+
+// Negate sets p = -q and returns p.
+func (p *P384Point) Negate(q *P384Point) *P384Point {
+	p.x.Set(q.x)
+	p.y.Sub(new(fiat.P384Element), q.y)
+	p.z.Set(q.z)
+	return p
+}
+
+// Negate sets p = -q and returns p.
+func (p *P521Point) Negate(q *P521Point) *P521Point {
+	p.x.Set(q.x)
+	p.y.Sub(new(fiat.P521Element), q.y)
+	p.z.Set(q.z)
+	return p
+}
+
+// IsInfinity returns 1 if p is the point-at-infinity, 0 otherwise.
+func (p *P224Point) IsInfinity() int {
+	return p.z.IsZero()
+}
+
+// IsInfinity returns 1 if p is the point-at-infinity, 0 otherwise.
+func (p *P384Point) IsInfinity() int {
+	return p.z.IsZero()
+}
+
+// IsInfinity returns 1 if p is the point-at-infinity, 0 otherwise.
+func (p *P521Point) IsInfinity() int {
+	return p.z.IsZero()
+}
+
+// Equal returns 1 if p and q represent the same point, 0 otherwise.
+func (p *P224Point) Equal(q *P224Point) int {
+	pinf := p.z.IsZero()
+	qinf := q.z.IsZero()
+	bothinf := pinf & qinf
+	noneinf := (1 - pinf) & (1 - qinf)
+	px := new(fiat.P224Element).Mul(p.x, q.z)
+	qx := new(fiat.P224Element).Mul(q.x, p.z)
+	py := new(fiat.P224Element).Mul(p.y, q.z)
+	qy := new(fiat.P224Element).Mul(q.y, p.z)
+	return bothinf | (noneinf & px.Equal(qx) & py.Equal(qy))
+}
+
+// Equal returns 1 if p and q represent the same point, 0 otherwise.
+func (p *P384Point) Equal(q *P384Point) int {
+	pinf := p.z.IsZero()
+	qinf := q.z.IsZero()
+	bothinf := pinf & qinf
+	noneinf := (1 - pinf) & (1 - qinf)
+	px := new(fiat.P384Element).Mul(p.x, q.z)
+	qx := new(fiat.P384Element).Mul(q.x, p.z)
+	py := new(fiat.P384Element).Mul(p.y, q.z)
+	qy := new(fiat.P384Element).Mul(q.y, p.z)
+	return bothinf | (noneinf & px.Equal(qx) & py.Equal(qy))
+}
+
+// Equal returns 1 if p and q represent the same point, 0 otherwise.
+func (p *P521Point) Equal(q *P521Point) int {
+	pinf := p.z.IsZero()
+	qinf := q.z.IsZero()
+	bothinf := pinf & qinf
+	noneinf := (1 - pinf) & (1 - qinf)
+	px := new(fiat.P521Element).Mul(p.x, q.z)
+	qx := new(fiat.P521Element).Mul(q.x, p.z)
+	py := new(fiat.P521Element).Mul(p.y, q.z)
+	qy := new(fiat.P521Element).Mul(q.y, p.z)
+	return bothinf | (noneinf & px.Equal(qx) & py.Equal(qy))
+}

extra_asm.go

@@ -0,0 +1,61 @@
+// Copyright 2023 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.
+
+//go:build !purego && (amd64 || arm64 || ppc64le || s390x)
+
+package nistec
+
+import "sources.truenas.cloud/code/nistec/internal/fiat"
+
+// Negate sets p = -q and returns p.
+func (p *P256Point) Negate(q *P256Point) *P256Point {
+	// fiat.P256Element is a little-endian Montgomery domain fully-reduced
+	// element, like p256Element, so they are actually interchangable.
+	qy := new(fiat.P256Element)
+	*qy.Bits() = q.y
+	py := new(fiat.P256Element).Sub(new(fiat.P256Element), qy)
+
+	p.x = q.x
+	p.y = *py.Bits()
+	p.z = q.z
+	return p
+}
+
+// IsInfinity returns 1 if p is the point-at-infinity, 0 otherwise.
+func (p *P256Point) IsInfinity() int {
+	return p.isInfinity()
+}
+
+// Equal returns 1 if p and q represent the same point, 0 otherwise.
+func (p *P256Point) Equal(q *P256Point) int {
+	pinf := p256Equal(&p.z, &p256Zero)
+	qinf := p256Equal(&q.z, &p256Zero)
+	bothinf := pinf & qinf
+	noneinf := (1 - pinf) & (1 - qinf)
+
+	// xp = Xp / Zp²
+	// yp = Yp / Zp³
+	// xq = Xq / Zq²
+	// yq = Yq / Zq³
+	// If Zp != 0 and Zq != 0, then:
+	//    xp == yp  <=>  Xp*Zq² == Xq*Zp²
+	//    xq == yq  <=>  Yp*Zq³ == Yq*Zp³
+	px := new(p256Element)
+	qx := new(p256Element)
+	py := new(p256Element)
+	qy := new(p256Element)
+	pz := new(p256Element)
+	qz := new(p256Element)
+	p256Sqr(pz, &p.z, 1)
+	p256Sqr(qz, &q.z, 1)
+	p256Mul(px, &p.x, qz)
+	p256Mul(qx, &q.x, pz)
+	samex := p256Equal(px, qx)
+	p256Mul(pz, pz, &p.z)
+	p256Mul(qz, qz, &q.z)
+	p256Mul(py, &p.y, qz)
+	p256Mul(qy, &q.y, pz)
+	samey := p256Equal(py, qy)
+	return bothinf | (noneinf & samex & samey)
+}

extra_noasm.go

@@ -0,0 +1,35 @@
+// Copyright 2023 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.
+
+//go:build purego || (!amd64 && !arm64 && !ppc64le && !s390x)
+
+package nistec
+
+import "sources.truenas.cloud/code/nistec/internal/fiat"
+
+// Negate sets p = -q and returns p.
+func (p *P256Point) Negate(q *P256Point) *P256Point {
+	p.x.Set(&q.x)
+	p.y.Sub(new(fiat.P256Element), &q.y)
+	p.z.Set(&q.z)
+	return p
+}
+
+// IsInfinity returns 1 if p is the point-at-infinity, 0 otherwise.
+func (p *P256Point) IsInfinity() int {
+	return p.z.IsZero()
+}
+
+// Equal returns 1 if p and q represent the same point, 0 otherwise.
+func (p *P256Point) Equal(q *P256Point) int {
+	pinf := p.z.IsZero()
+	qinf := q.z.IsZero()
+	bothinf := pinf & qinf
+	noneinf := (1 - pinf) & (1 - qinf)
+	px := new(fiat.P256Element).Mul(&p.x, &q.z)
+	qx := new(fiat.P256Element).Mul(&q.x, &p.z)
+	py := new(fiat.P256Element).Mul(&p.y, &q.z)
+	qy := new(fiat.P256Element).Mul(&q.y, &p.z)
+	return bothinf | (noneinf & px.Equal(qx) & py.Equal(qy))
+}

extra_norace_test.go

@@ -0,0 +1,5 @@
+//go:build !race
+
+package nistec_test
+
+const raceEnabled = false

extra_race_test.go

@@ -0,0 +1,5 @@
+//go:build race
+
+package nistec_test
+
+const raceEnabled = true

extra_test.go

@@ -0,0 +1,174 @@
+// Copyright 2023 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 (
+	"bytes"
+	"crypto/elliptic"
+	"math/big"
+	"testing"
+
+	"sources.truenas.cloud/code/nistec"
+)
+
+type nistPointExtra[T any] interface {
+	Bytes() []byte
+	SetGenerator() T
+	SetBytes([]byte) (T, error)
+	Set(T) T
+	Add(T, T) T
+	Double(T) T
+	Negate(T) T
+	ScalarMult(T, []byte) (T, error)
+	ScalarBaseMult([]byte) (T, error)
+	IsInfinity() int
+	Equal(T) int
+}
+
+func TestNegate(t *testing.T) {
+	t.Run("P224", func(t *testing.T) {
+		testNegate(t, nistec.NewP224Point, elliptic.P224())
+	})
+	t.Run("P256", func(t *testing.T) {
+		testNegate(t, nistec.NewP256Point, elliptic.P256())
+	})
+	t.Run("P384", func(t *testing.T) {
+		testNegate(t, nistec.NewP384Point, elliptic.P384())
+	})
+	t.Run("P521", func(t *testing.T) {
+		testNegate(t, nistec.NewP521Point, elliptic.P521())
+	})
+}
+
+func testNegate[P nistPointExtra[P]](t *testing.T, newPoint func() P, c elliptic.Curve) {
+	p := newPoint().SetGenerator()
+	p.Add(p, p) // make a test point with z != 1
+
+	p1 := newPoint().Negate(p)
+
+	minusOne := new(big.Int).Sub(c.Params().N, big.NewInt(1)).Bytes()
+	p2, err := newPoint().ScalarMult(p, minusOne)
+	fatalIfErr(t, err)
+	if !bytes.Equal(p1.Bytes(), p2.Bytes()) {
+		t.Error("-P != [-1]P")
+	}
+
+	p.Negate(p)
+	if !bytes.Equal(p.Bytes(), p1.Bytes()) {
+		t.Error("-P (aliasing) != -P")
+	}
+}
+
+func TestEqual(t *testing.T) {
+	t.Run("P224", func(t *testing.T) {
+		testEqual(t, nistec.NewP224Point, elliptic.P224())
+	})
+	t.Run("P256", func(t *testing.T) {
+		testEqual(t, nistec.NewP256Point, elliptic.P256())
+	})
+	t.Run("P384", func(t *testing.T) {
+		testEqual(t, nistec.NewP384Point, elliptic.P384())
+	})
+	t.Run("P521", func(t *testing.T) {
+		testEqual(t, nistec.NewP521Point, elliptic.P521())
+	})
+}
+
+func testEqual[P nistPointExtra[P]](t *testing.T, newPoint func() P, c elliptic.Curve) {
+	inf := newPoint()
+	g := newPoint().SetGenerator()
+	g.Add(g, g)
+	if inf.Equal(inf) != 1 {
+		t.Error("inf != inf")
+	}
+	if g.Equal(g) != 1 {
+		t.Error("G != G")
+	}
+	if inf.Equal(g) != 0 || g.Equal(inf) != 0 {
+		t.Error("G == inf")
+	}
+	p1 := newPoint().Set(g)
+	p2 := newPoint()
+	p3 := newPoint()
+	sc := make([]byte, (c.Params().N.BitLen()+7)>>3)
+	for i := 1; i < 30; i++ {
+		sc[len(sc)-1] = byte(i)
+		_, err := p2.ScalarMult(g, sc)
+		fatalIfErr(t, err)
+		sc[len(sc)-1] = byte(i * 2)
+		_, err = p3.ScalarBaseMult(sc)
+		fatalIfErr(t, err)
+		if p1.Equal(p2) != 1 || p2.Equal(p1) != 1 {
+			t.Error("n*G != n*G (1)")
+		}
+		if p1.Equal(p3) != 1 || p3.Equal(p1) != 1 {
+			t.Error("n*G != n*G (2)")
+		}
+		if p2.Equal(p3) != 1 || p3.Equal(p2) != 1 {
+			t.Error("n*G != n*G (3)")
+		}
+		p1.Add(p1, g)
+		if p1.Equal(g) != 0 || g.Equal(p1) != 0 {
+			t.Error("n*G == G")
+		}
+	}
+}
+
+func TestInfinity(t *testing.T) {
+	t.Run("P224", func(t *testing.T) {
+		testInfinity(t, nistec.NewP224Point, elliptic.P224())
+	})
+	t.Run("P256", func(t *testing.T) {
+		testInfinity(t, nistec.NewP256Point, elliptic.P256())
+	})
+	t.Run("P384", func(t *testing.T) {
+		testInfinity(t, nistec.NewP384Point, elliptic.P384())
+	})
+	t.Run("P521", func(t *testing.T) {
+		testInfinity(t, nistec.NewP521Point, elliptic.P521())
+	})
+}
+
+func testInfinity[P nistPointExtra[P]](t *testing.T, newPoint func() P, c elliptic.Curve) {
+	inf := newPoint()
+	g := newPoint().SetGenerator()
+	g.Add(g, g)
+	if inf.IsInfinity() != 1 {
+		t.Error("inf != inf")
+	}
+	if g.IsInfinity() != 0 {
+		t.Error("G == inf")
+	}
+	p1 := newPoint().Set(g)
+	p2 := newPoint()
+	sc := make([]byte, (c.Params().N.BitLen()+7)>>3)
+	for i := 1; i < 30; i++ {
+		sc[len(sc)-1] = byte(i)
+		_, err := p2.ScalarMult(g, sc)
+		fatalIfErr(t, err)
+		if p2.IsInfinity() != 0 {
+			t.Error("n*G == inf (1)")
+		}
+		p2.Negate(p2)
+		if p2.IsInfinity() != 0 {
+			t.Error("n*G == inf (2)")
+		}
+		p2.Add(p2, p1)
+		if p2.IsInfinity() != 1 {
+			t.Error("n*G - n*G != inf")
+		}
+		p1.Add(p1, g)
+		if p1.IsInfinity() != 0 {
+			t.Error("n*G == inf (3)")
+		}
+	}
+}
+
+func fatalIfErr(t *testing.T, err error) {
+	t.Helper()
+	if err != nil {
+		t.Fatal(err)
+	}
+}

generate.go

@@ -0,0 +1,628 @@
+// Copyright 2022 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.
+
+//go:build ignore
+
+package main
+
+// Running this generator requires addchain v0.4.0, which can be installed with
+//
+//   go install github.com/mmcloughlin/addchain/cmd/addchain@v0.4.0
+//
+
+import (
+	"bytes"
+	"crypto/elliptic"
+	"fmt"
+	"go/format"
+	"io"
+	"log"
+	"math/big"
+	"os"
+	"os/exec"
+	"strings"
+	"text/template"
+)
+
+var curves = []struct {
+	P       string
+	Element string
+	Params  *elliptic.CurveParams
+}{
+	{
+		P:       "P224",
+		Element: "fiat.P224Element",
+		Params:  elliptic.P224().Params(),
+	},
+	{
+		P:       "P384",
+		Element: "fiat.P384Element",
+		Params:  elliptic.P384().Params(),
+	},
+	{
+		P:       "P521",
+		Element: "fiat.P521Element",
+		Params:  elliptic.P521().Params(),
+	},
+}
+
+func main() {
+	t := template.Must(template.New("tmplNISTEC").Parse(tmplNISTEC))
+
+	tmplAddchainFile, err := os.CreateTemp("", "addchain-template")
+	if err != nil {
+		log.Fatal(err)
+	}
+	defer os.Remove(tmplAddchainFile.Name())
+	if _, err := io.WriteString(tmplAddchainFile, tmplAddchain); err != nil {
+		log.Fatal(err)
+	}
+	if err := tmplAddchainFile.Close(); err != nil {
+		log.Fatal(err)
+	}
+
+	for _, c := range curves {
+		p := strings.ToLower(c.P)
+		elementLen := (c.Params.BitSize + 7) / 8
+		B := fmt.Sprintf("%#v", c.Params.B.FillBytes(make([]byte, elementLen)))
+		Gx := fmt.Sprintf("%#v", c.Params.Gx.FillBytes(make([]byte, elementLen)))
+		Gy := fmt.Sprintf("%#v", c.Params.Gy.FillBytes(make([]byte, elementLen)))
+
+		log.Printf("Generating %s.go...", p)
+		f, err := os.Create(p + ".go")
+		if err != nil {
+			log.Fatal(err)
+		}
+		defer f.Close()
+		buf := &bytes.Buffer{}
+		if err := t.Execute(buf, map[string]interface{}{
+			"P": c.P, "p": p, "B": B, "Gx": Gx, "Gy": Gy,
+			"Element": c.Element, "ElementLen": elementLen,
+		}); err != nil {
+			log.Fatal(err)
+		}
+		out, err := format.Source(buf.Bytes())
+		if err != nil {
+			log.Fatal(err)
+		}
+		if _, err := f.Write(out); err != nil {
+			log.Fatal(err)
+		}
+
+		// If p = 3 mod 4, implement modular square root by exponentiation.
+		mod4 := new(big.Int).Mod(c.Params.P, big.NewInt(4))
+		if mod4.Cmp(big.NewInt(3)) != 0 {
+			continue
+		}
+
+		exp := new(big.Int).Add(c.Params.P, big.NewInt(1))
+		exp.Div(exp, big.NewInt(4))
+
+		tmp, err := os.CreateTemp("", "addchain-"+p)
+		if err != nil {
+			log.Fatal(err)
+		}
+		defer os.Remove(tmp.Name())
+		cmd := exec.Command("addchain", "search", fmt.Sprintf("%d", exp))
+		cmd.Stderr = os.Stderr
+		cmd.Stdout = tmp
+		if err := cmd.Run(); err != nil {
+			log.Fatal(err)
+		}
+		if err := tmp.Close(); err != nil {
+			log.Fatal(err)
+		}
+		cmd = exec.Command("addchain", "gen", "-tmpl", tmplAddchainFile.Name(), tmp.Name())
+		cmd.Stderr = os.Stderr
+		out, err = cmd.Output()
+		if err != nil {
+			log.Fatal(err)
+		}
+		out = bytes.Replace(out, []byte("Element"), []byte(c.Element), -1)
+		out = bytes.Replace(out, []byte("sqrtCandidate"), []byte(p+"SqrtCandidate"), -1)
+		out, err = format.Source(out)
+		if err != nil {
+			log.Fatal(err)
+		}
+		if _, err := f.Write(out); err != nil {
+			log.Fatal(err)
+		}
+	}
+}
+
+const tmplNISTEC = `// Copyright 2022 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.
+
+// Code generated by generate.go. DO NOT EDIT.
+
+package nistec
+
+import (
+	"crypto/subtle"
+	"errors"
+	"sync"
+
+	"filippo.io/nistec/internal/fiat"
+)
+
+// {{.p}}ElementLength is the length of an element of the base or scalar field,
+// which have the same bytes length for all NIST P curves.
+const {{.p}}ElementLength = {{ .ElementLen }}
+
+// {{.P}}Point is a {{.P}} point. The zero value is NOT valid.
+type {{.P}}Point struct {
+	// The point is represented in projective coordinates (X:Y:Z),
+	// where x = X/Z and y = Y/Z.
+	x, y, z *{{.Element}}
+}
+
+// New{{.P}}Point returns a new {{.P}}Point representing the point at infinity point.
+func New{{.P}}Point() *{{.P}}Point {
+	return &{{.P}}Point{
+		x: new({{.Element}}),
+		y: new({{.Element}}).One(),
+		z: new({{.Element}}),
+	}
+}
+
+// SetGenerator sets p to the canonical generator and returns p.
+func (p *{{.P}}Point) SetGenerator() *{{.P}}Point {
+	p.x.SetBytes({{.Gx}})
+	p.y.SetBytes({{.Gy}})
+	p.z.One()
+	return p
+}
+
+// Set sets p = q and returns p.
+func (p *{{.P}}Point) Set(q *{{.P}}Point) *{{.P}}Point {
+	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 *{{.P}}Point) SetBytes(b []byte) (*{{.P}}Point, error) {
+	switch {
+	// Point at infinity.
+	case len(b) == 1 && b[0] == 0:
+		return p.Set(New{{.P}}Point()), nil
+
+	// Uncompressed form.
+	case len(b) == 1+2*{{.p}}ElementLength && b[0] == 4:
+		x, err := new({{.Element}}).SetBytes(b[1 : 1+{{.p}}ElementLength])
+		if err != nil {
+			return nil, err
+		}
+		y, err := new({{.Element}}).SetBytes(b[1+{{.p}}ElementLength:])
+		if err != nil {
+			return nil, err
+		}
+		if err := {{.p}}CheckOnCurve(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+{{.p}}ElementLength && (b[0] == 2 || b[0] == 3):
+		x, err := new({{.Element}}).SetBytes(b[1:])
+		if err != nil {
+			return nil, err
+		}
+
+		// y² = x³ - 3x + b
+		y := {{.p}}Polynomial(new({{.Element}}), x)
+		if !{{.p}}Sqrt(y, y) {
+			return nil, errors.New("invalid {{.P}} compressed point encoding")
+		}
+
+		// Select the positive or negative root, as indicated by the least
+		// significant bit, based on the encoding type byte.
+		otherRoot := new({{.Element}})
+		otherRoot.Sub(otherRoot, y)
+		cond := y.Bytes()[{{.p}}ElementLength-1]&1 ^ b[0]&1
+		y.Select(otherRoot, y, int(cond))
+
+		p.x.Set(x)
+		p.y.Set(y)
+		p.z.One()
+		return p, nil
+
+	default:
+		return nil, errors.New("invalid {{.P}} point encoding")
+	}
+}
+
+
+var _{{.p}}B *{{.Element}}
+var _{{.p}}BOnce sync.Once
+
+func {{.p}}B() *{{.Element}} {
+	_{{.p}}BOnce.Do(func() {
+		_{{.p}}B, _ = new({{.Element}}).SetBytes({{.B}})
+	})
+	return _{{.p}}B
+}
+
+// {{.p}}Polynomial sets y2 to x³ - 3x + b, and returns y2.
+func {{.p}}Polynomial(y2, x *{{.Element}}) *{{.Element}} {
+	y2.Square(x)
+	y2.Mul(y2, x)
+
+	threeX := new({{.Element}}).Add(x, x)
+	threeX.Add(threeX, x)
+	y2.Sub(y2, threeX)
+
+	return y2.Add(y2, {{.p}}B())
+}
+
+func {{.p}}CheckOnCurve(x, y *{{.Element}}) error {
+	// y² = x³ - 3x + b
+	rhs := {{.p}}Polynomial(new({{.Element}}), x)
+	lhs := new({{.Element}}).Square(y)
+	if rhs.Equal(lhs) != 1 {
+		return errors.New("{{.P}} 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 *{{.P}}Point) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [1+2*{{.p}}ElementLength]byte
+	return p.bytes(&out)
+}
+
+func (p *{{.P}}Point) bytes(out *[1+2*{{.p}}ElementLength]byte) []byte {
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new({{.Element}}).Invert(p.z)
+	x := new({{.Element}}).Mul(p.x, zinv)
+	y := new({{.Element}}).Mul(p.y, zinv)
+
+	buf := append(out[:0], 4)
+	buf = append(buf, x.Bytes()...)
+	buf = append(buf, y.Bytes()...)
+	return buf
+}
+
+// BytesX returns the encoding of the x-coordinate of p, as specified in SEC 1,
+// Version 2.0, Section 2.3.5, or an error if p is the point at infinity.
+func (p *{{.P}}Point) BytesX() ([]byte, error) {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [{{.p}}ElementLength]byte
+	return p.bytesX(&out)
+}
+
+func (p *{{.P}}Point) bytesX(out *[{{.p}}ElementLength]byte) ([]byte, error) {
+	if p.z.IsZero() == 1 {
+		return nil, errors.New("{{.P}} point is the point at infinity")
+	}
+
+	zinv := new({{.Element}}).Invert(p.z)
+	x := new({{.Element}}).Mul(p.x, zinv)
+
+	return append(out[:0], x.Bytes()...), nil
+}
+
+// BytesCompressed returns the compressed 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 *{{.P}}Point) BytesCompressed() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [1 + {{.p}}ElementLength]byte
+	return p.bytesCompressed(&out)
+}
+
+func (p *{{.P}}Point) bytesCompressed(out *[1 + {{.p}}ElementLength]byte) []byte {
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new({{.Element}}).Invert(p.z)
+	x := new({{.Element}}).Mul(p.x, zinv)
+	y := new({{.Element}}).Mul(p.y, zinv)
+
+	// Encode the sign of the y coordinate (indicated by the least significant
+	// bit) as the encoding type (2 or 3).
+	buf := append(out[:0], 2)
+	buf[0] |= y.Bytes()[{{.p}}ElementLength-1] & 1
+	buf = append(buf, x.Bytes()...)
+	return buf
+}
+
+// Add sets q = p1 + p2, and returns q. The points may overlap.
+func (q *{{.P}}Point) Add(p1, p2 *{{.P}}Point) *{{.P}}Point {
+	// 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({{.Element}}).Mul(p1.x, p2.x)   // t0 := X1 * X2
+	t1 := new({{.Element}}).Mul(p1.y, p2.y)   // t1 := Y1 * Y2
+	t2 := new({{.Element}}).Mul(p1.z, p2.z)   // t2 := Z1 * Z2
+	t3 := new({{.Element}}).Add(p1.x, p1.y)   // t3 := X1 + Y1
+	t4 := new({{.Element}}).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({{.Element}}).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({{.Element}}).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({{.Element}}).Mul({{.p}}B(), 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({{.p}}B(), 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 *{{.P}}Point) Double(p *{{.P}}Point) *{{.P}}Point {
+	// 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({{.Element}}).Square(p.x)      // t0 := X ^ 2
+	t1 := new({{.Element}}).Square(p.y)      // t1 := Y ^ 2
+	t2 := new({{.Element}}).Square(p.z)      // t2 := Z ^ 2
+	t3 := new({{.Element}}).Mul(p.x, p.y)    // t3 := X * Y
+	t3.Add(t3, t3)                           // t3 := t3 + t3
+	z3 := new({{.Element}}).Mul(p.x, p.z)    // Z3 := X * Z
+	z3.Add(z3, z3)                           // Z3 := Z3 + Z3
+	y3 := new({{.Element}}).Mul({{.p}}B(), t2) // Y3 := b * t2
+	y3.Sub(y3, z3)                           // Y3 := Y3 - Z3
+	x3 := new({{.Element}}).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({{.p}}B(), 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 *{{.P}}Point) Select(p1, p2 *{{.P}}Point, cond int) *{{.P}}Point {
+	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
+}
+
+// A {{.p}}Table holds the first 15 multiples of a point at offset -1, so [1]P
+// is at table[0], [15]P is at table[14], and [0]P is implicitly the identity
+// point.
+type {{.p}}Table [15]*{{.P}}Point
+
+// Select selects the n-th multiple of the table base point into p. It works in
+// constant time by iterating over every entry of the table. n must be in [0, 15].
+func (table *{{.p}}Table) Select(p *{{.P}}Point, n uint8) {
+	if n >= 16 {
+		panic("nistec: internal error: {{.p}}Table called with out-of-bounds value")
+	}
+	p.Set(New{{.P}}Point())
+	for i := uint8(1); i < 16; i++ {
+		cond := subtle.ConstantTimeByteEq(i, n)
+		p.Select(table[i-1], p, cond)
+	}
+}
+
+// ScalarMult sets p = scalar * q, and returns p.
+func (p *{{.P}}Point) ScalarMult(q *{{.P}}Point, scalar []byte) (*{{.P}}Point, error) {
+	// Compute a {{.p}}Table for the base point q. The explicit New{{.P}}Point
+	// calls get inlined, letting the allocations live on the stack.
+	var table = {{.p}}Table{New{{.P}}Point(), New{{.P}}Point(), New{{.P}}Point(),
+		New{{.P}}Point(), New{{.P}}Point(), New{{.P}}Point(), New{{.P}}Point(),
+		New{{.P}}Point(), New{{.P}}Point(), New{{.P}}Point(), New{{.P}}Point(),
+		New{{.P}}Point(), New{{.P}}Point(), New{{.P}}Point(), New{{.P}}Point()}
+	table[0].Set(q)
+	for i := 1; i < 15; i += 2 {
+		table[i].Double(table[i/2])
+		table[i+1].Add(table[i], 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 := New{{.P}}Point()
+	p.Set(New{{.P}}Point())
+	for i, byte := range scalar {
+		// No need to double on the first iteration, as p is the identity at
+		// this point, and [N]∞ = ∞.
+		if i != 0 {
+			p.Double(p)
+			p.Double(p)
+			p.Double(p)
+			p.Double(p)
+		}
+
+		windowValue := byte >> 4
+		table.Select(t, windowValue)
+		p.Add(p, t)
+
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+
+		windowValue = byte & 0b1111
+		table.Select(t, windowValue)
+		p.Add(p, t)
+	}
+
+	return p, nil
+}
+
+var {{.p}}GeneratorTable *[{{.p}}ElementLength * 2]{{.p}}Table
+var {{.p}}GeneratorTableOnce sync.Once
+
+// generatorTable returns a sequence of {{.p}}Tables. The first table contains
+// multiples of G. Each successive table is the previous table doubled four
+// times.
+func (p *{{.P}}Point) generatorTable() *[{{.p}}ElementLength * 2]{{.p}}Table {
+	{{.p}}GeneratorTableOnce.Do(func() {
+		{{.p}}GeneratorTable = new([{{.p}}ElementLength * 2]{{.p}}Table)
+		base := New{{.P}}Point().SetGenerator()
+		for i := 0; i < {{.p}}ElementLength*2; i++ {
+			{{.p}}GeneratorTable[i][0] = New{{.P}}Point().Set(base)
+			for j := 1; j < 15; j++ {
+				{{.p}}GeneratorTable[i][j] = New{{.P}}Point().Add({{.p}}GeneratorTable[i][j-1], base)
+			}
+			base.Double(base)
+			base.Double(base)
+			base.Double(base)
+			base.Double(base)
+		}
+	})
+	return {{.p}}GeneratorTable
+}
+
+// ScalarBaseMult sets p = scalar * B, where B is the canonical generator, and
+// returns p.
+func (p *{{.P}}Point) ScalarBaseMult(scalar []byte) (*{{.P}}Point, error) {
+	if len(scalar) != {{.p}}ElementLength {
+		return nil, errors.New("invalid scalar length")
+	}
+	tables := p.generatorTable()
+
+	// This is also a scalar multiplication with a four-bit window like in
+	// ScalarMult, but in this case the doublings are precomputed. The value
+	// [windowValue]G added at iteration k would normally get doubled
+	// (totIterations-k)×4 times, but with a larger precomputation we can
+	// instead add [2^((totIterations-k)×4)][windowValue]G and avoid the
+	// doublings between iterations.
+	t := New{{.P}}Point()
+	p.Set(New{{.P}}Point())
+	tableIndex := len(tables) - 1
+	for _, byte := range scalar {
+		windowValue := byte >> 4
+		tables[tableIndex].Select(t, windowValue)
+		p.Add(p, t)
+		tableIndex--
+
+		windowValue = byte & 0b1111
+		tables[tableIndex].Select(t, windowValue)
+		p.Add(p, t)
+		tableIndex--
+	}
+
+	return p, nil
+}
+
+// {{.p}}Sqrt sets e to a square root of x. If x is not a square, {{.p}}Sqrt returns
+// false and e is unchanged. e and x can overlap.
+func {{.p}}Sqrt(e, x *{{ .Element }}) (isSquare bool) {
+	candidate := new({{ .Element }})
+	{{.p}}SqrtCandidate(candidate, x)
+	square := new({{ .Element }}).Square(candidate)
+	if square.Equal(x) != 1 {
+		return false
+	}
+	e.Set(candidate)
+	return true
+}
+`
+
+const tmplAddchain = `
+// sqrtCandidate sets z to a square root candidate for x. z and x must not overlap.
+func sqrtCandidate(z, x *Element) {
+	// Since p = 3 mod 4, exponentiation by (p + 1) / 4 yields a square root candidate.
+	//
+	// The sequence of {{ .Ops.Adds }} multiplications and {{ .Ops.Doubles }} squarings is derived from the
+	// following addition chain generated with {{ .Meta.Module }} {{ .Meta.ReleaseTag }}.
+	//
+	{{- range lines (format .Script) }}
+	//	{{ . }}
+	{{- end }}
+	//
+
+	{{- range .Program.Temporaries }}
+	var {{ . }} = new(Element)
+	{{- end }}
+	{{ range $i := .Program.Instructions -}}
+	{{- with add $i.Op }}
+	{{ $i.Output }}.Mul({{ .X }}, {{ .Y }})
+	{{- end -}}
+
+	{{- with double $i.Op }}
+	{{ $i.Output }}.Square({{ .X }})
+	{{- end -}}
+
+	{{- with shift $i.Op -}}
+	{{- $first := 0 -}}
+	{{- if ne $i.Output.Identifier .X.Identifier }}
+	{{ $i.Output }}.Square({{ .X }})
+	{{- $first = 1 -}}
+	{{- end }}
+	for s := {{ $first }}; s < {{ .S }}; s++ {
+		{{ $i.Output }}.Square({{ $i.Output }})
+	}
+	{{- end -}}
+	{{- end }}
+}
+`

go.mod

@@ -0,0 +1,5 @@
+module sources.truenas.cloud/code/nistec
+
+go 1.24.0
+
+require golang.org/x/sys v0.36.0

go.sum

@@ -0,0 +1,2 @@
+golang.org/x/sys v0.36.0 h1:KVRy2GtZBrk1cBYA7MKu5bEZFxQk4NIDV6RLVcC8o0k=
+golang.org/x/sys v0.36.0/go.mod h1:OgkHotnGiDImocRcuBABYBEXf8A9a87e/uXjp9XT3ks=

internal/byteorder/byteorder.go

@@ -0,0 +1,149 @@
+// Copyright 2024 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 byteorder provides functions for decoding and encoding
+// little and big endian integer types from/to byte slices.
+package byteorder
+
+func LEUint16(b []byte) uint16 {
+	_ = b[1] // bounds check hint to compiler; see golang.org/issue/14808
+	return uint16(b[0]) | uint16(b[1])<<8
+}
+
+func LEPutUint16(b []byte, v uint16) {
+	_ = b[1] // early bounds check to guarantee safety of writes below
+	b[0] = byte(v)
+	b[1] = byte(v >> 8)
+}
+
+func LEAppendUint16(b []byte, v uint16) []byte {
+	return append(b,
+		byte(v),
+		byte(v>>8),
+	)
+}
+
+func LEUint32(b []byte) uint32 {
+	_ = b[3] // bounds check hint to compiler; see golang.org/issue/14808
+	return uint32(b[0]) | uint32(b[1])<<8 | uint32(b[2])<<16 | uint32(b[3])<<24
+}
+
+func LEPutUint32(b []byte, v uint32) {
+	_ = b[3] // early bounds check to guarantee safety of writes below
+	b[0] = byte(v)
+	b[1] = byte(v >> 8)
+	b[2] = byte(v >> 16)
+	b[3] = byte(v >> 24)
+}
+
+func LEAppendUint32(b []byte, v uint32) []byte {
+	return append(b,
+		byte(v),
+		byte(v>>8),
+		byte(v>>16),
+		byte(v>>24),
+	)
+}
+
+func LEUint64(b []byte) uint64 {
+	_ = b[7] // bounds check hint to compiler; see golang.org/issue/14808
+	return uint64(b[0]) | uint64(b[1])<<8 | uint64(b[2])<<16 | uint64(b[3])<<24 |
+		uint64(b[4])<<32 | uint64(b[5])<<40 | uint64(b[6])<<48 | uint64(b[7])<<56
+}
+
+func LEPutUint64(b []byte, v uint64) {
+	_ = b[7] // early bounds check to guarantee safety of writes below
+	b[0] = byte(v)
+	b[1] = byte(v >> 8)
+	b[2] = byte(v >> 16)
+	b[3] = byte(v >> 24)
+	b[4] = byte(v >> 32)
+	b[5] = byte(v >> 40)
+	b[6] = byte(v >> 48)
+	b[7] = byte(v >> 56)
+}
+
+func LEAppendUint64(b []byte, v uint64) []byte {
+	return append(b,
+		byte(v),
+		byte(v>>8),
+		byte(v>>16),
+		byte(v>>24),
+		byte(v>>32),
+		byte(v>>40),
+		byte(v>>48),
+		byte(v>>56),
+	)
+}
+
+func BEUint16(b []byte) uint16 {
+	_ = b[1] // bounds check hint to compiler; see golang.org/issue/14808
+	return uint16(b[1]) | uint16(b[0])<<8
+}
+
+func BEPutUint16(b []byte, v uint16) {
+	_ = b[1] // early bounds check to guarantee safety of writes below
+	b[0] = byte(v >> 8)
+	b[1] = byte(v)
+}
+
+func BEAppendUint16(b []byte, v uint16) []byte {
+	return append(b,
+		byte(v>>8),
+		byte(v),
+	)
+}
+
+func BEUint32(b []byte) uint32 {
+	_ = b[3] // bounds check hint to compiler; see golang.org/issue/14808
+	return uint32(b[3]) | uint32(b[2])<<8 | uint32(b[1])<<16 | uint32(b[0])<<24
+}
+
+func BEPutUint32(b []byte, v uint32) {
+	_ = b[3] // early bounds check to guarantee safety of writes below
+	b[0] = byte(v >> 24)
+	b[1] = byte(v >> 16)
+	b[2] = byte(v >> 8)
+	b[3] = byte(v)
+}
+
+func BEAppendUint32(b []byte, v uint32) []byte {
+	return append(b,
+		byte(v>>24),
+		byte(v>>16),
+		byte(v>>8),
+		byte(v),
+	)
+}
+
+func BEUint64(b []byte) uint64 {
+	_ = b[7] // bounds check hint to compiler; see golang.org/issue/14808
+	return uint64(b[7]) | uint64(b[6])<<8 | uint64(b[5])<<16 | uint64(b[4])<<24 |
+		uint64(b[3])<<32 | uint64(b[2])<<40 | uint64(b[1])<<48 | uint64(b[0])<<56
+}
+
+func BEPutUint64(b []byte, v uint64) {
+	_ = b[7] // early bounds check to guarantee safety of writes below
+	b[0] = byte(v >> 56)
+	b[1] = byte(v >> 48)
+	b[2] = byte(v >> 40)
+	b[3] = byte(v >> 32)
+	b[4] = byte(v >> 24)
+	b[5] = byte(v >> 16)
+	b[6] = byte(v >> 8)
+	b[7] = byte(v)
+}
+
+func BEAppendUint64(b []byte, v uint64) []byte {
+	return append(b,
+		byte(v>>56),
+		byte(v>>48),
+		byte(v>>40),
+		byte(v>>32),
+		byte(v>>24),
+		byte(v>>16),
+		byte(v>>8),
+		byte(v),
+	)
+}

internal/fiat/Dockerfile

@@ -0,0 +1,12 @@
+# 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.
+
+FROM coqorg/coq:8.13.2
+
+RUN git clone https://github.com/mit-plv/fiat-crypto && cd fiat-crypto && \
+    git checkout 23d2dbc4ab897d14bde4404f70cd6991635f9c01 && \
+    git submodule update --init --recursive
+RUN cd fiat-crypto && eval $(opam env) && make -j4 standalone-ocaml SKIP_BEDROCK2=1
+
+ENV PATH /home/coq/fiat-crypto/src/ExtractionOCaml:$PATH

internal/fiat/README

@@ -0,0 +1,34 @@
+The code in this package was autogenerated by the fiat-crypto project
+at version v0.0.9 from a formally verified model, and by the addchain
+project at a recent tip version.
+
+    docker build -t fiat-crypto:v0.0.9 .
+    go install github.com/mmcloughlin/addchain/cmd/addchain@v0.3.1-0.20211027081849-6a7d3decbe08
+    go run generate.go
+
+fiat-crypto code comes under the following license.
+
+    Copyright (c) 2015-2020 The fiat-crypto Authors. All rights reserved.
+
+    Redistribution and use in source and binary forms, with or without
+    modification, are permitted provided that the following conditions are
+    met:
+
+        1. Redistributions of source code must retain the above copyright
+        notice, this list of conditions and the following disclaimer.
+
+    THIS SOFTWARE IS PROVIDED BY the fiat-crypto authors "AS IS"
+    AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
+    THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL Berkeley Software Design,
+    Inc. BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+The authors are listed at
+
+    https://github.com/mit-plv/fiat-crypto/blob/master/AUTHORS

internal/fiat/fiat_test.go

@@ -0,0 +1,65 @@
+// 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 fiat_test
+
+import (
+	"testing"
+
+	"sources.truenas.cloud/code/nistec/internal/fiat"
+)
+
+func BenchmarkMul(b *testing.B) {
+	b.Run("P224", func(b *testing.B) {
+		v := new(fiat.P224Element).One()
+		b.ReportAllocs()
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			v.Mul(v, v)
+		}
+	})
+	b.Run("P384", func(b *testing.B) {
+		v := new(fiat.P384Element).One()
+		b.ReportAllocs()
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			v.Mul(v, v)
+		}
+	})
+	b.Run("P521", func(b *testing.B) {
+		v := new(fiat.P521Element).One()
+		b.ReportAllocs()
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			v.Mul(v, v)
+		}
+	})
+}
+
+func BenchmarkSquare(b *testing.B) {
+	b.Run("P224", func(b *testing.B) {
+		v := new(fiat.P224Element).One()
+		b.ReportAllocs()
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			v.Square(v)
+		}
+	})
+	b.Run("P384", func(b *testing.B) {
+		v := new(fiat.P384Element).One()
+		b.ReportAllocs()
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			v.Square(v)
+		}
+	})
+	b.Run("P521", func(b *testing.B) {
+		v := new(fiat.P521Element).One()
+		b.ReportAllocs()
+		b.ResetTimer()
+		for i := 0; i < b.N; i++ {
+			v.Square(v)
+		}
+	})
+}

internal/fiat/generate.go

@@ -0,0 +1,325 @@
+// 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.
+
+//go:build ignore
+
+package main
+
+import (
+	"bytes"
+	"go/format"
+	"io"
+	"log"
+	"os"
+	"os/exec"
+	"text/template"
+)
+
+var curves = []struct {
+	Element  string
+	Prime    string
+	Prefix   string
+	FiatType string
+	BytesLen int
+}{
+	{
+		Element:  "P224Element",
+		Prime:    "2^224 - 2^96 + 1",
+		Prefix:   "p224",
+		FiatType: "[4]uint64",
+		BytesLen: 28,
+	},
+	// The P-256 fiat implementation is used only on 32-bit architectures, but
+	// the uint32 fiat code is for some reason slower than the uint64 one. That
+	// suggests there is a wide margin for improvement.
+	{
+		Element:  "P256Element",
+		Prime:    "2^256 - 2^224 + 2^192 + 2^96 - 1",
+		Prefix:   "p256",
+		FiatType: "[4]uint64",
+		BytesLen: 32,
+	},
+	{
+		Element:  "P384Element",
+		Prime:    "2^384 - 2^128 - 2^96 + 2^32 - 1",
+		Prefix:   "p384",
+		FiatType: "[6]uint64",
+		BytesLen: 48,
+	},
+	// Note that unsaturated_solinas would be about 2x faster than
+	// word_by_word_montgomery for P-521, but this curve is used rarely enough
+	// that it's not worth carrying unsaturated_solinas support for it.
+	{
+		Element:  "P521Element",
+		Prime:    "2^521 - 1",
+		Prefix:   "p521",
+		FiatType: "[9]uint64",
+		BytesLen: 66,
+	},
+}
+
+func main() {
+	t := template.Must(template.New("montgomery").Parse(tmplWrapper))
+
+	tmplAddchainFile, err := os.CreateTemp("", "addchain-template")
+	if err != nil {
+		log.Fatal(err)
+	}
+	defer os.Remove(tmplAddchainFile.Name())
+	if _, err := io.WriteString(tmplAddchainFile, tmplAddchain); err != nil {
+		log.Fatal(err)
+	}
+	if err := tmplAddchainFile.Close(); err != nil {
+		log.Fatal(err)
+	}
+
+	for _, c := range curves {
+		log.Printf("Generating %s.go...", c.Prefix)
+		f, err := os.Create(c.Prefix + ".go")
+		if err != nil {
+			log.Fatal(err)
+		}
+		if err := t.Execute(f, c); err != nil {
+			log.Fatal(err)
+		}
+		if err := f.Close(); err != nil {
+			log.Fatal(err)
+		}
+
+		log.Printf("Generating %s_fiat64.go...", c.Prefix)
+		cmd := exec.Command("docker", "run", "--rm", "--entrypoint", "word_by_word_montgomery",
+			"fiat-crypto:v0.0.9", "--lang", "Go", "--no-wide-int", "--cmovznz-by-mul",
+			"--relax-primitive-carry-to-bitwidth", "32,64", "--internal-static",
+			"--public-function-case", "camelCase", "--public-type-case", "camelCase",
+			"--private-function-case", "camelCase", "--private-type-case", "camelCase",
+			"--doc-text-before-function-name", "", "--doc-newline-before-package-declaration",
+			"--doc-prepend-header", "Code generated by Fiat Cryptography. DO NOT EDIT.",
+			"--package-name", "fiat", "--no-prefix-fiat", c.Prefix, "64", c.Prime,
+			"mul", "square", "add", "sub", "one", "from_montgomery", "to_montgomery",
+			"selectznz", "to_bytes", "from_bytes")
+		cmd.Stderr = os.Stderr
+		out, err := cmd.Output()
+		if err != nil {
+			log.Fatal(err)
+		}
+		out, err = format.Source(out)
+		if err != nil {
+			log.Fatal(err)
+		}
+		if err := os.WriteFile(c.Prefix+"_fiat64.go", out, 0644); err != nil {
+			log.Fatal(err)
+		}
+
+		log.Printf("Generating %s_invert.go...", c.Prefix)
+		f, err = os.CreateTemp("", "addchain-"+c.Prefix)
+		if err != nil {
+			log.Fatal(err)
+		}
+		defer os.Remove(f.Name())
+		cmd = exec.Command("addchain", "search", c.Prime+" - 2")
+		cmd.Stderr = os.Stderr
+		cmd.Stdout = f
+		if err := cmd.Run(); err != nil {
+			log.Fatal(err)
+		}
+		if err := f.Close(); err != nil {
+			log.Fatal(err)
+		}
+		cmd = exec.Command("addchain", "gen", "-tmpl", tmplAddchainFile.Name(), f.Name())
+		cmd.Stderr = os.Stderr
+		out, err = cmd.Output()
+		if err != nil {
+			log.Fatal(err)
+		}
+		out = bytes.Replace(out, []byte("Element"), []byte(c.Element), -1)
+		out, err = format.Source(out)
+		if err != nil {
+			log.Fatal(err)
+		}
+		if err := os.WriteFile(c.Prefix+"_invert.go", out, 0644); err != nil {
+			log.Fatal(err)
+		}
+	}
+}
+
+const tmplWrapper = `// 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.
+
+// Code generated by generate.go. DO NOT EDIT.
+
+package fiat
+
+import (
+	"crypto/subtle"
+	"errors"
+)
+
+// {{ .Element }} is an integer modulo {{ .Prime }}.
+//
+// The zero value is a valid zero element.
+type {{ .Element }} struct {
+	// Values are represented internally always in the Montgomery domain, and
+	// converted in Bytes and SetBytes.
+	x {{ .Prefix }}MontgomeryDomainFieldElement
+}
+
+const {{ .Prefix }}ElementLen = {{ .BytesLen }}
+
+type {{ .Prefix }}UntypedFieldElement = {{ .FiatType }}
+
+// One sets e = 1, and returns e.
+func (e *{{ .Element }}) One() *{{ .Element }} {
+	{{ .Prefix }}SetOne(&e.x)
+	return e
+}
+
+// Equal returns 1 if e == t, and zero otherwise.
+func (e *{{ .Element }}) Equal(t *{{ .Element }}) int {
+	eBytes := e.Bytes()
+	tBytes := t.Bytes()
+	return subtle.ConstantTimeCompare(eBytes, tBytes)
+}
+
+// IsZero returns 1 if e == 0, and zero otherwise.
+func (e *{{ .Element }}) IsZero() int {
+	zero := make([]byte, {{ .Prefix }}ElementLen)
+	eBytes := e.Bytes()
+	return subtle.ConstantTimeCompare(eBytes, zero)
+}
+
+// Set sets e = t, and returns e.
+func (e *{{ .Element }}) Set(t *{{ .Element }}) *{{ .Element }} {
+	e.x = t.x
+	return e
+}
+
+// Bytes returns the {{ .BytesLen }}-byte big-endian encoding of e.
+func (e *{{ .Element }}) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [{{ .Prefix }}ElementLen]byte
+	return e.bytes(&out)
+}
+
+func (e *{{ .Element }}) bytes(out *[{{ .Prefix }}ElementLen]byte) []byte {
+	var tmp {{ .Prefix }}NonMontgomeryDomainFieldElement
+	{{ .Prefix }}FromMontgomery(&tmp, &e.x)
+	{{ .Prefix }}ToBytes(out, (*{{ .Prefix }}UntypedFieldElement)(&tmp))
+	{{ .Prefix }}InvertEndianness(out[:])
+	return out[:]
+}
+
+// SetBytes sets e = v, where v is a big-endian {{ .BytesLen }}-byte encoding, and returns e.
+// If v is not {{ .BytesLen }} bytes or it encodes a value higher than {{ .Prime }},
+// SetBytes returns nil and an error, and e is unchanged.
+func (e *{{ .Element }}) SetBytes(v []byte) (*{{ .Element }}, error) {
+	if len(v) != {{ .Prefix }}ElementLen {
+		return nil, errors.New("invalid {{ .Element }} encoding")
+	}
+
+	// Check for non-canonical encodings (p + k, 2p + k, etc.) by comparing to
+	// the encoding of -1 mod p, so p - 1, the highest canonical encoding.
+	var minusOneEncoding = new({{ .Element }}).Sub(
+		new({{ .Element }}), new({{ .Element }}).One()).Bytes()
+	if subtle.ConstantTimeLessOrEqBytes(v, minusOneEncoding) == 0 {
+		return nil, errors.New("invalid {{ .Element }} encoding")
+	}
+
+	var in [{{ .Prefix }}ElementLen]byte
+	copy(in[:], v)
+	{{ .Prefix }}InvertEndianness(in[:])
+	var tmp {{ .Prefix }}NonMontgomeryDomainFieldElement
+	{{ .Prefix }}FromBytes((*{{ .Prefix }}UntypedFieldElement)(&tmp), &in)
+	{{ .Prefix }}ToMontgomery(&e.x, &tmp)
+	return e, nil
+}
+
+// Add sets e = t1 + t2, and returns e.
+func (e *{{ .Element }}) Add(t1, t2 *{{ .Element }}) *{{ .Element }} {
+	{{ .Prefix }}Add(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Sub sets e = t1 - t2, and returns e.
+func (e *{{ .Element }}) Sub(t1, t2 *{{ .Element }}) *{{ .Element }} {
+	{{ .Prefix }}Sub(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Mul sets e = t1 * t2, and returns e.
+func (e *{{ .Element }}) Mul(t1, t2 *{{ .Element }}) *{{ .Element }} {
+	{{ .Prefix }}Mul(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Square sets e = t * t, and returns e.
+func (e *{{ .Element }}) Square(t *{{ .Element }}) *{{ .Element }} {
+	{{ .Prefix }}Square(&e.x, &t.x)
+	return e
+}
+
+// Select sets v to a if cond == 1, and to b if cond == 0.
+func (v *{{ .Element }}) Select(a, b *{{ .Element }}, cond int) *{{ .Element }} {
+	{{ .Prefix }}Selectznz((*{{ .Prefix }}UntypedFieldElement)(&v.x), {{ .Prefix }}Uint1(cond),
+		(*{{ .Prefix }}UntypedFieldElement)(&b.x), (*{{ .Prefix }}UntypedFieldElement)(&a.x))
+	return v
+}
+
+func {{ .Prefix }}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]
+	}
+}
+`
+
+const tmplAddchain = `// 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.
+
+// Code generated by {{ .Meta.Name }}. DO NOT EDIT.
+
+package fiat
+
+// Invert sets e = 1/x, and returns e.
+//
+// If x == 0, Invert returns e = 0.
+func (e *Element) Invert(x *Element) *Element {
+	// Inversion is implemented as exponentiation with exponent p − 2.
+	// The sequence of {{ .Ops.Adds }} multiplications and {{ .Ops.Doubles }} squarings is derived from the
+	// following addition chain generated with {{ .Meta.Module }} {{ .Meta.ReleaseTag }}.
+	//
+	{{- range lines (format .Script) }}
+	//	{{ . }}
+	{{- end }}
+	//
+
+	var z = new(Element).Set(e)
+	{{- range .Program.Temporaries }}
+	var {{ . }} = new(Element)
+	{{- end }}
+	{{ range $i := .Program.Instructions -}}
+	{{- with add $i.Op }}
+	{{ $i.Output }}.Mul({{ .X }}, {{ .Y }})
+	{{- end -}}
+
+	{{- with double $i.Op }}
+	{{ $i.Output }}.Square({{ .X }})
+	{{- end -}}
+
+	{{- with shift $i.Op -}}
+	{{- $first := 0 -}}
+	{{- if ne $i.Output.Identifier .X.Identifier }}
+	{{ $i.Output }}.Square({{ .X }})
+	{{- $first = 1 -}}
+	{{- end }}
+	for s := {{ $first }}; s < {{ .S }}; s++ {
+		{{ $i.Output }}.Square({{ $i.Output }})
+	}
+	{{- end -}}
+	{{- end }}
+
+	return e.Set(z)
+}
+`

internal/fiat/p224.go

@@ -0,0 +1,130 @@
+// 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.
+
+// Code generated by generate.go. DO NOT EDIT.
+
+package fiat
+
+import (
+	"errors"
+
+	"sources.truenas.cloud/code/nistec/internal/subtle"
+)
+
+// P224Element is an integer modulo 2^224 - 2^96 + 1.
+//
+// The zero value is a valid zero element.
+type P224Element struct {
+	// Values are represented internally always in the Montgomery domain, and
+	// converted in Bytes and SetBytes.
+	x p224MontgomeryDomainFieldElement
+}
+
+const p224ElementLen = 28
+
+type p224UntypedFieldElement = [4]uint64
+
+// One sets e = 1, and returns e.
+func (e *P224Element) One() *P224Element {
+	p224SetOne(&e.x)
+	return e
+}
+
+// Equal returns 1 if e == t, and zero otherwise.
+func (e *P224Element) Equal(t *P224Element) int {
+	eBytes := e.Bytes()
+	tBytes := t.Bytes()
+	return subtle.ConstantTimeCompare(eBytes, tBytes)
+}
+
+// IsZero returns 1 if e == 0, and zero otherwise.
+func (e *P224Element) IsZero() int {
+	zero := make([]byte, p224ElementLen)
+	eBytes := e.Bytes()
+	return subtle.ConstantTimeCompare(eBytes, zero)
+}
+
+// Set sets e = t, and returns e.
+func (e *P224Element) Set(t *P224Element) *P224Element {
+	e.x = t.x
+	return e
+}
+
+// Bytes returns the 28-byte big-endian encoding of e.
+func (e *P224Element) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p224ElementLen]byte
+	return e.bytes(&out)
+}
+
+func (e *P224Element) bytes(out *[p224ElementLen]byte) []byte {
+	var tmp p224NonMontgomeryDomainFieldElement
+	p224FromMontgomery(&tmp, &e.x)
+	p224ToBytes(out, (*p224UntypedFieldElement)(&tmp))
+	p224InvertEndianness(out[:])
+	return out[:]
+}
+
+// SetBytes sets e = v, where v is a big-endian 28-byte encoding, and returns e.
+// If v is not 28 bytes or it encodes a value higher than 2^224 - 2^96 + 1,
+// SetBytes returns nil and an error, and e is unchanged.
+func (e *P224Element) SetBytes(v []byte) (*P224Element, error) {
+	if len(v) != p224ElementLen {
+		return nil, errors.New("invalid P224Element encoding")
+	}
+
+	// Check for non-canonical encodings (p + k, 2p + k, etc.) by comparing to
+	// the encoding of -1 mod p, so p - 1, the highest canonical encoding.
+	var minusOneEncoding = new(P224Element).Sub(
+		new(P224Element), new(P224Element).One()).Bytes()
+	if subtle.ConstantTimeLessOrEqBytes(v, minusOneEncoding) == 0 {
+		return nil, errors.New("invalid P224Element encoding")
+	}
+
+	var in [p224ElementLen]byte
+	copy(in[:], v)
+	p224InvertEndianness(in[:])
+	var tmp p224NonMontgomeryDomainFieldElement
+	p224FromBytes((*p224UntypedFieldElement)(&tmp), &in)
+	p224ToMontgomery(&e.x, &tmp)
+	return e, nil
+}
+
+// Add sets e = t1 + t2, and returns e.
+func (e *P224Element) Add(t1, t2 *P224Element) *P224Element {
+	p224Add(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Sub sets e = t1 - t2, and returns e.
+func (e *P224Element) Sub(t1, t2 *P224Element) *P224Element {
+	p224Sub(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Mul sets e = t1 * t2, and returns e.
+func (e *P224Element) Mul(t1, t2 *P224Element) *P224Element {
+	p224Mul(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Square sets e = t * t, and returns e.
+func (e *P224Element) Square(t *P224Element) *P224Element {
+	p224Square(&e.x, &t.x)
+	return e
+}
+
+// Select sets v to a if cond == 1, and to b if cond == 0.
+func (v *P224Element) Select(a, b *P224Element, cond int) *P224Element {
+	p224Selectznz((*p224UntypedFieldElement)(&v.x), p224Uint1(cond),
+		(*p224UntypedFieldElement)(&b.x), (*p224UntypedFieldElement)(&a.x))
+	return v
+}
+
+func p224InvertEndianness(v []byte) {
+	for i := 0; i < len(v)/2; i++ {
+		v[i], v[len(v)-1-i] = v[len(v)-1-i], v[i]
+	}
+}

internal/fiat/p224_invert.go

@@ -0,0 +1,87 @@
+// 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.
+
+// Code generated by addchain. DO NOT EDIT.
+
+package fiat
+
+// Invert sets e = 1/x, and returns e.
+//
+// If x == 0, Invert returns e = 0.
+func (e *P224Element) Invert(x *P224Element) *P224Element {
+	// Inversion is implemented as exponentiation with exponent p − 2.
+	// The sequence of 11 multiplications and 223 squarings is derived from the
+	// following addition chain generated with github.com/mmcloughlin/addchain v0.4.0.
+	//
+	//	_10     = 2*1
+	//	_11     = 1 + _10
+	//	_110    = 2*_11
+	//	_111    = 1 + _110
+	//	_111000 = _111 << 3
+	//	_111111 = _111 + _111000
+	//	x12     = _111111 << 6 + _111111
+	//	x14     = x12 << 2 + _11
+	//	x17     = x14 << 3 + _111
+	//	x31     = x17 << 14 + x14
+	//	x48     = x31 << 17 + x17
+	//	x96     = x48 << 48 + x48
+	//	x127    = x96 << 31 + x31
+	//	return    x127 << 97 + x96
+	//
+
+	var z = new(P224Element).Set(e)
+	var t0 = new(P224Element)
+	var t1 = new(P224Element)
+	var t2 = new(P224Element)
+
+	z.Square(x)
+	t0.Mul(x, z)
+	z.Square(t0)
+	z.Mul(x, z)
+	t1.Square(z)
+	for s := 1; s < 3; s++ {
+		t1.Square(t1)
+	}
+	t1.Mul(z, t1)
+	t2.Square(t1)
+	for s := 1; s < 6; s++ {
+		t2.Square(t2)
+	}
+	t1.Mul(t1, t2)
+	for s := 0; s < 2; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	t1.Square(t0)
+	for s := 1; s < 3; s++ {
+		t1.Square(t1)
+	}
+	z.Mul(z, t1)
+	t1.Square(z)
+	for s := 1; s < 14; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	t1.Square(t0)
+	for s := 1; s < 17; s++ {
+		t1.Square(t1)
+	}
+	z.Mul(z, t1)
+	t1.Square(z)
+	for s := 1; s < 48; s++ {
+		t1.Square(t1)
+	}
+	z.Mul(z, t1)
+	t1.Square(z)
+	for s := 1; s < 31; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	for s := 0; s < 97; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+
+	return e.Set(z)
+}

internal/fiat/p256.go

@@ -0,0 +1,130 @@
+// 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.
+
+// Code generated by generate.go. DO NOT EDIT.
+
+package fiat
+
+import (
+	"errors"
+
+	"sources.truenas.cloud/code/nistec/internal/subtle"
+)
+
+// P256Element is an integer modulo 2^256 - 2^224 + 2^192 + 2^96 - 1.
+//
+// The zero value is a valid zero element.
+type P256Element struct {
+	// Values are represented internally always in the Montgomery domain, and
+	// converted in Bytes and SetBytes.
+	x p256MontgomeryDomainFieldElement
+}
+
+const p256ElementLen = 32
+
+type p256UntypedFieldElement = [4]uint64
+
+// One sets e = 1, and returns e.
+func (e *P256Element) One() *P256Element {
+	p256SetOne(&e.x)
+	return e
+}
+
+// Equal returns 1 if e == t, and zero otherwise.
+func (e *P256Element) Equal(t *P256Element) int {
+	eBytes := e.Bytes()
+	tBytes := t.Bytes()
+	return subtle.ConstantTimeCompare(eBytes, tBytes)
+}
+
+// IsZero returns 1 if e == 0, and zero otherwise.
+func (e *P256Element) IsZero() int {
+	zero := make([]byte, p256ElementLen)
+	eBytes := e.Bytes()
+	return subtle.ConstantTimeCompare(eBytes, zero)
+}
+
+// Set sets e = t, and returns e.
+func (e *P256Element) Set(t *P256Element) *P256Element {
+	e.x = t.x
+	return e
+}
+
+// Bytes returns the 32-byte big-endian encoding of e.
+func (e *P256Element) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p256ElementLen]byte
+	return e.bytes(&out)
+}
+
+func (e *P256Element) bytes(out *[p256ElementLen]byte) []byte {
+	var tmp p256NonMontgomeryDomainFieldElement
+	p256FromMontgomery(&tmp, &e.x)
+	p256ToBytes(out, (*p256UntypedFieldElement)(&tmp))
+	p256InvertEndianness(out[:])
+	return out[:]
+}
+
+// SetBytes sets e = v, where v is a big-endian 32-byte encoding, and returns e.
+// If v is not 32 bytes or it encodes a value higher than 2^256 - 2^224 + 2^192 + 2^96 - 1,
+// SetBytes returns nil and an error, and e is unchanged.
+func (e *P256Element) SetBytes(v []byte) (*P256Element, error) {
+	if len(v) != p256ElementLen {
+		return nil, errors.New("invalid P256Element encoding")
+	}
+
+	// Check for non-canonical encodings (p + k, 2p + k, etc.) by comparing to
+	// the encoding of -1 mod p, so p - 1, the highest canonical encoding.
+	var minusOneEncoding = new(P256Element).Sub(
+		new(P256Element), new(P256Element).One()).Bytes()
+	if subtle.ConstantTimeLessOrEqBytes(v, minusOneEncoding) == 0 {
+		return nil, errors.New("invalid P256Element encoding")
+	}
+
+	var in [p256ElementLen]byte
+	copy(in[:], v)
+	p256InvertEndianness(in[:])
+	var tmp p256NonMontgomeryDomainFieldElement
+	p256FromBytes((*p256UntypedFieldElement)(&tmp), &in)
+	p256ToMontgomery(&e.x, &tmp)
+	return e, nil
+}
+
+// Add sets e = t1 + t2, and returns e.
+func (e *P256Element) Add(t1, t2 *P256Element) *P256Element {
+	p256Add(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Sub sets e = t1 - t2, and returns e.
+func (e *P256Element) Sub(t1, t2 *P256Element) *P256Element {
+	p256Sub(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Mul sets e = t1 * t2, and returns e.
+func (e *P256Element) Mul(t1, t2 *P256Element) *P256Element {
+	p256Mul(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Square sets e = t * t, and returns e.
+func (e *P256Element) Square(t *P256Element) *P256Element {
+	p256Square(&e.x, &t.x)
+	return e
+}
+
+// Select sets v to a if cond == 1, and to b if cond == 0.
+func (v *P256Element) Select(a, b *P256Element, cond int) *P256Element {
+	p256Selectznz((*p256UntypedFieldElement)(&v.x), p256Uint1(cond),
+		(*p256UntypedFieldElement)(&b.x), (*p256UntypedFieldElement)(&a.x))
+	return v
+}
+
+func p256InvertEndianness(v []byte) {
+	for i := 0; i < len(v)/2; i++ {
+		v[i], v[len(v)-1-i] = v[len(v)-1-i], v[i]
+	}
+}

internal/fiat/p256_extra.go

@@ -0,0 +1,12 @@
+// Copyright 2023 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 fiat
+
+// Bits returns a reference to the underlying little-endian fully-reduced
+// Montgomery representation of e. Handle with care.
+func (e *P256Element) Bits() *[4]uint64 {
+	var _ p256MontgomeryDomainFieldElement = e.x
+	return (*[4]uint64)(&e.x)
+}

internal/fiat/p256_invert.go

@@ -0,0 +1,84 @@
+// 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.
+
+// Code generated by addchain. DO NOT EDIT.
+
+package fiat
+
+// Invert sets e = 1/x, and returns e.
+//
+// If x == 0, Invert returns e = 0.
+func (e *P256Element) Invert(x *P256Element) *P256Element {
+	// Inversion is implemented as exponentiation with exponent p − 2.
+	// The sequence of 12 multiplications and 255 squarings is derived from the
+	// following addition chain generated with github.com/mmcloughlin/addchain v0.4.0.
+	//
+	//	_10     = 2*1
+	//	_11     = 1 + _10
+	//	_110    = 2*_11
+	//	_111    = 1 + _110
+	//	_111000 = _111 << 3
+	//	_111111 = _111 + _111000
+	//	x12     = _111111 << 6 + _111111
+	//	x15     = x12 << 3 + _111
+	//	x16     = 2*x15 + 1
+	//	x32     = x16 << 16 + x16
+	//	i53     = x32 << 15
+	//	x47     = x15 + i53
+	//	i263    = ((i53 << 17 + 1) << 143 + x47) << 47
+	//	return    (x47 + i263) << 2 + 1
+	//
+
+	var z = new(P256Element).Set(e)
+	var t0 = new(P256Element)
+	var t1 = new(P256Element)
+
+	z.Square(x)
+	z.Mul(x, z)
+	z.Square(z)
+	z.Mul(x, z)
+	t0.Square(z)
+	for s := 1; s < 3; s++ {
+		t0.Square(t0)
+	}
+	t0.Mul(z, t0)
+	t1.Square(t0)
+	for s := 1; s < 6; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	for s := 0; s < 3; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	t0.Square(z)
+	t0.Mul(x, t0)
+	t1.Square(t0)
+	for s := 1; s < 16; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	for s := 0; s < 15; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	for s := 0; s < 17; s++ {
+		t0.Square(t0)
+	}
+	t0.Mul(x, t0)
+	for s := 0; s < 143; s++ {
+		t0.Square(t0)
+	}
+	t0.Mul(z, t0)
+	for s := 0; s < 47; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	for s := 0; s < 2; s++ {
+		z.Square(z)
+	}
+	z.Mul(x, z)
+
+	return e.Set(z)
+}

internal/fiat/p384.go

@@ -0,0 +1,130 @@
+// 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.
+
+// Code generated by generate.go. DO NOT EDIT.
+
+package fiat
+
+import (
+	"errors"
+
+	"sources.truenas.cloud/code/nistec/internal/subtle"
+)
+
+// P384Element is an integer modulo 2^384 - 2^128 - 2^96 + 2^32 - 1.
+//
+// The zero value is a valid zero element.
+type P384Element struct {
+	// Values are represented internally always in the Montgomery domain, and
+	// converted in Bytes and SetBytes.
+	x p384MontgomeryDomainFieldElement
+}
+
+const p384ElementLen = 48
+
+type p384UntypedFieldElement = [6]uint64
+
+// One sets e = 1, and returns e.
+func (e *P384Element) One() *P384Element {
+	p384SetOne(&e.x)
+	return e
+}
+
+// Equal returns 1 if e == t, and zero otherwise.
+func (e *P384Element) Equal(t *P384Element) int {
+	eBytes := e.Bytes()
+	tBytes := t.Bytes()
+	return subtle.ConstantTimeCompare(eBytes, tBytes)
+}
+
+// IsZero returns 1 if e == 0, and zero otherwise.
+func (e *P384Element) IsZero() int {
+	zero := make([]byte, p384ElementLen)
+	eBytes := e.Bytes()
+	return subtle.ConstantTimeCompare(eBytes, zero)
+}
+
+// Set sets e = t, and returns e.
+func (e *P384Element) Set(t *P384Element) *P384Element {
+	e.x = t.x
+	return e
+}
+
+// Bytes returns the 48-byte big-endian encoding of e.
+func (e *P384Element) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p384ElementLen]byte
+	return e.bytes(&out)
+}
+
+func (e *P384Element) bytes(out *[p384ElementLen]byte) []byte {
+	var tmp p384NonMontgomeryDomainFieldElement
+	p384FromMontgomery(&tmp, &e.x)
+	p384ToBytes(out, (*p384UntypedFieldElement)(&tmp))
+	p384InvertEndianness(out[:])
+	return out[:]
+}
+
+// SetBytes sets e = v, where v is a big-endian 48-byte encoding, and returns e.
+// If v is not 48 bytes or it encodes a value higher than 2^384 - 2^128 - 2^96 + 2^32 - 1,
+// SetBytes returns nil and an error, and e is unchanged.
+func (e *P384Element) SetBytes(v []byte) (*P384Element, error) {
+	if len(v) != p384ElementLen {
+		return nil, errors.New("invalid P384Element encoding")
+	}
+
+	// Check for non-canonical encodings (p + k, 2p + k, etc.) by comparing to
+	// the encoding of -1 mod p, so p - 1, the highest canonical encoding.
+	var minusOneEncoding = new(P384Element).Sub(
+		new(P384Element), new(P384Element).One()).Bytes()
+	if subtle.ConstantTimeLessOrEqBytes(v, minusOneEncoding) == 0 {
+		return nil, errors.New("invalid P384Element encoding")
+	}
+
+	var in [p384ElementLen]byte
+	copy(in[:], v)
+	p384InvertEndianness(in[:])
+	var tmp p384NonMontgomeryDomainFieldElement
+	p384FromBytes((*p384UntypedFieldElement)(&tmp), &in)
+	p384ToMontgomery(&e.x, &tmp)
+	return e, nil
+}
+
+// Add sets e = t1 + t2, and returns e.
+func (e *P384Element) Add(t1, t2 *P384Element) *P384Element {
+	p384Add(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Sub sets e = t1 - t2, and returns e.
+func (e *P384Element) Sub(t1, t2 *P384Element) *P384Element {
+	p384Sub(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Mul sets e = t1 * t2, and returns e.
+func (e *P384Element) Mul(t1, t2 *P384Element) *P384Element {
+	p384Mul(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Square sets e = t * t, and returns e.
+func (e *P384Element) Square(t *P384Element) *P384Element {
+	p384Square(&e.x, &t.x)
+	return e
+}
+
+// Select sets v to a if cond == 1, and to b if cond == 0.
+func (v *P384Element) Select(a, b *P384Element, cond int) *P384Element {
+	p384Selectznz((*p384UntypedFieldElement)(&v.x), p384Uint1(cond),
+		(*p384UntypedFieldElement)(&b.x), (*p384UntypedFieldElement)(&a.x))
+	return v
+}
+
+func p384InvertEndianness(v []byte) {
+	for i := 0; i < len(v)/2; i++ {
+		v[i], v[len(v)-1-i] = v[len(v)-1-i], v[i]
+	}
+}

internal/fiat/p384_invert.go

@@ -0,0 +1,102 @@
+// 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.
+
+// Code generated by addchain. DO NOT EDIT.
+
+package fiat
+
+// Invert sets e = 1/x, and returns e.
+//
+// If x == 0, Invert returns e = 0.
+func (e *P384Element) Invert(x *P384Element) *P384Element {
+	// Inversion is implemented as exponentiation with exponent p − 2.
+	// The sequence of 15 multiplications and 383 squarings is derived from the
+	// following addition chain generated with github.com/mmcloughlin/addchain v0.4.0.
+	//
+	//	_10     = 2*1
+	//	_11     = 1 + _10
+	//	_110    = 2*_11
+	//	_111    = 1 + _110
+	//	_111000 = _111 << 3
+	//	_111111 = _111 + _111000
+	//	x12     = _111111 << 6 + _111111
+	//	x24     = x12 << 12 + x12
+	//	x30     = x24 << 6 + _111111
+	//	x31     = 2*x30 + 1
+	//	x32     = 2*x31 + 1
+	//	x63     = x32 << 31 + x31
+	//	x126    = x63 << 63 + x63
+	//	x252    = x126 << 126 + x126
+	//	x255    = x252 << 3 + _111
+	//	i397    = ((x255 << 33 + x32) << 94 + x30) << 2
+	//	return    1 + i397
+	//
+
+	var z = new(P384Element).Set(e)
+	var t0 = new(P384Element)
+	var t1 = new(P384Element)
+	var t2 = new(P384Element)
+	var t3 = new(P384Element)
+
+	z.Square(x)
+	z.Mul(x, z)
+	z.Square(z)
+	t1.Mul(x, z)
+	z.Square(t1)
+	for s := 1; s < 3; s++ {
+		z.Square(z)
+	}
+	z.Mul(t1, z)
+	t0.Square(z)
+	for s := 1; s < 6; s++ {
+		t0.Square(t0)
+	}
+	t0.Mul(z, t0)
+	t2.Square(t0)
+	for s := 1; s < 12; s++ {
+		t2.Square(t2)
+	}
+	t0.Mul(t0, t2)
+	for s := 0; s < 6; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	t0.Square(z)
+	t2.Mul(x, t0)
+	t0.Square(t2)
+	t0.Mul(x, t0)
+	t3.Square(t0)
+	for s := 1; s < 31; s++ {
+		t3.Square(t3)
+	}
+	t2.Mul(t2, t3)
+	t3.Square(t2)
+	for s := 1; s < 63; s++ {
+		t3.Square(t3)
+	}
+	t2.Mul(t2, t3)
+	t3.Square(t2)
+	for s := 1; s < 126; s++ {
+		t3.Square(t3)
+	}
+	t2.Mul(t2, t3)
+	for s := 0; s < 3; s++ {
+		t2.Square(t2)
+	}
+	t1.Mul(t1, t2)
+	for s := 0; s < 33; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	for s := 0; s < 94; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	for s := 0; s < 2; s++ {
+		z.Square(z)
+	}
+	z.Mul(x, z)
+
+	return e.Set(z)
+}

internal/fiat/p521.go

@@ -0,0 +1,130 @@
+// 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.
+
+// Code generated by generate.go. DO NOT EDIT.
+
+package fiat
+
+import (
+	"errors"
+
+	"sources.truenas.cloud/code/nistec/internal/subtle"
+)
+
+// P521Element is an integer modulo 2^521 - 1.
+//
+// The zero value is a valid zero element.
+type P521Element struct {
+	// Values are represented internally always in the Montgomery domain, and
+	// converted in Bytes and SetBytes.
+	x p521MontgomeryDomainFieldElement
+}
+
+const p521ElementLen = 66
+
+type p521UntypedFieldElement = [9]uint64
+
+// One sets e = 1, and returns e.
+func (e *P521Element) One() *P521Element {
+	p521SetOne(&e.x)
+	return e
+}
+
+// Equal returns 1 if e == t, and zero otherwise.
+func (e *P521Element) Equal(t *P521Element) int {
+	eBytes := e.Bytes()
+	tBytes := t.Bytes()
+	return subtle.ConstantTimeCompare(eBytes, tBytes)
+}
+
+// IsZero returns 1 if e == 0, and zero otherwise.
+func (e *P521Element) IsZero() int {
+	zero := make([]byte, p521ElementLen)
+	eBytes := e.Bytes()
+	return subtle.ConstantTimeCompare(eBytes, zero)
+}
+
+// Set sets e = t, and returns e.
+func (e *P521Element) Set(t *P521Element) *P521Element {
+	e.x = t.x
+	return e
+}
+
+// Bytes returns the 66-byte big-endian encoding of e.
+func (e *P521Element) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p521ElementLen]byte
+	return e.bytes(&out)
+}
+
+func (e *P521Element) bytes(out *[p521ElementLen]byte) []byte {
+	var tmp p521NonMontgomeryDomainFieldElement
+	p521FromMontgomery(&tmp, &e.x)
+	p521ToBytes(out, (*p521UntypedFieldElement)(&tmp))
+	p521InvertEndianness(out[:])
+	return out[:]
+}
+
+// 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) != p521ElementLen {
+		return nil, errors.New("invalid P521Element encoding")
+	}
+
+	// Check for non-canonical encodings (p + k, 2p + k, etc.) by comparing to
+	// the encoding of -1 mod p, so p - 1, the highest canonical encoding.
+	var minusOneEncoding = new(P521Element).Sub(
+		new(P521Element), new(P521Element).One()).Bytes()
+	if subtle.ConstantTimeLessOrEqBytes(v, minusOneEncoding) == 0 {
+		return nil, errors.New("invalid P521Element encoding")
+	}
+
+	var in [p521ElementLen]byte
+	copy(in[:], v)
+	p521InvertEndianness(in[:])
+	var tmp p521NonMontgomeryDomainFieldElement
+	p521FromBytes((*p521UntypedFieldElement)(&tmp), &in)
+	p521ToMontgomery(&e.x, &tmp)
+	return e, nil
+}
+
+// Add sets e = t1 + t2, and returns e.
+func (e *P521Element) Add(t1, t2 *P521Element) *P521Element {
+	p521Add(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Sub sets e = t1 - t2, and returns e.
+func (e *P521Element) Sub(t1, t2 *P521Element) *P521Element {
+	p521Sub(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Mul sets e = t1 * t2, and returns e.
+func (e *P521Element) Mul(t1, t2 *P521Element) *P521Element {
+	p521Mul(&e.x, &t1.x, &t2.x)
+	return e
+}
+
+// Square sets e = t * t, and returns e.
+func (e *P521Element) Square(t *P521Element) *P521Element {
+	p521Square(&e.x, &t.x)
+	return e
+}
+
+// Select sets v to a if cond == 1, and to b if cond == 0.
+func (v *P521Element) Select(a, b *P521Element, cond int) *P521Element {
+	p521Selectznz((*p521UntypedFieldElement)(&v.x), p521Uint1(cond),
+		(*p521UntypedFieldElement)(&b.x), (*p521UntypedFieldElement)(&a.x))
+	return v
+}
+
+func p521InvertEndianness(v []byte) {
+	for i := 0; i < len(v)/2; i++ {
+		v[i], v[len(v)-1-i] = v[len(v)-1-i], v[i]
+	}
+}

internal/fiat/p521_invert.go

@@ -0,0 +1,89 @@
+// 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.
+
+// Code generated by addchain. DO NOT EDIT.
+
+package fiat
+
+// Invert sets e = 1/x, and returns e.
+//
+// If x == 0, Invert returns e = 0.
+func (e *P521Element) Invert(x *P521Element) *P521Element {
+	// Inversion is implemented as exponentiation with exponent p − 2.
+	// The sequence of 13 multiplications and 520 squarings is derived from the
+	// following addition chain generated with github.com/mmcloughlin/addchain v0.4.0.
+	//
+	//	_10       = 2*1
+	//	_11       = 1 + _10
+	//	_1100     = _11 << 2
+	//	_1111     = _11 + _1100
+	//	_11110000 = _1111 << 4
+	//	_11111111 = _1111 + _11110000
+	//	x16       = _11111111 << 8 + _11111111
+	//	x32       = x16 << 16 + x16
+	//	x64       = x32 << 32 + x32
+	//	x65       = 2*x64 + 1
+	//	x129      = x65 << 64 + x64
+	//	x130      = 2*x129 + 1
+	//	x259      = x130 << 129 + x129
+	//	x260      = 2*x259 + 1
+	//	x519      = x260 << 259 + x259
+	//	return      x519 << 2 + 1
+	//
+
+	var z = new(P521Element).Set(e)
+	var t0 = new(P521Element)
+
+	z.Square(x)
+	z.Mul(x, z)
+	t0.Square(z)
+	for s := 1; s < 2; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	t0.Square(z)
+	for s := 1; s < 4; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	t0.Square(z)
+	for s := 1; s < 8; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	t0.Square(z)
+	for s := 1; s < 16; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	t0.Square(z)
+	for s := 1; s < 32; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	t0.Square(z)
+	t0.Mul(x, t0)
+	for s := 0; s < 64; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	t0.Square(z)
+	t0.Mul(x, t0)
+	for s := 0; s < 129; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	t0.Square(z)
+	t0.Mul(x, t0)
+	for s := 0; s < 259; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	for s := 0; s < 2; s++ {
+		z.Square(z)
+	}
+	z.Mul(x, z)
+
+	return e.Set(z)
+}

internal/subtle/subtle.go

@@ -0,0 +1,47 @@
+// Copyright 2009 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 subtle
+
+import (
+	"crypto/subtle"
+	"math/bits"
+
+	"sources.truenas.cloud/code/nistec/internal/byteorder"
+)
+
+func ConstantTimeCompare(x, y []byte) int {
+	return subtle.ConstantTimeCompare(x, y)
+}
+
+// ConstantTimeLessOrEqBytes returns 1 if x <= y and 0 otherwise. The comparison
+// is lexigraphical, or big-endian. The time taken is a function of the length of
+// the slices and is independent of the contents. If the lengths of x and y do not
+// match it returns 0 immediately.
+func ConstantTimeLessOrEqBytes(x, y []byte) int {
+	if len(x) != len(y) {
+		return 0
+	}
+
+	// Do a constant time subtraction chain y - x.
+	// If there is no borrow at the end, then x <= y.
+	var b uint64
+	for len(x) > 8 {
+		x0 := byteorder.BEUint64(x[len(x)-8:])
+		y0 := byteorder.BEUint64(y[len(y)-8:])
+		_, b = bits.Sub64(y0, x0, b)
+		x = x[:len(x)-8]
+		y = y[:len(y)-8]
+	}
+	if len(x) > 0 {
+		xb := make([]byte, 8)
+		yb := make([]byte, 8)
+		copy(xb[8-len(x):], x)
+		copy(yb[8-len(y):], y)
+		x0 := byteorder.BEUint64(xb)
+		y0 := byteorder.BEUint64(yb)
+		_, b = bits.Sub64(y0, x0, b)
+	}
+	return int(b ^ 1)
+}

nistec.go

@@ -0,0 +1,15 @@
+// Copyright 2022 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 elliptic curves from NIST SP 800-186.
+//
+// This package uses fiat-crypto or specialized assembly and Go code 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 where possible. 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
+
+//go:generate go run generate.go

nistec_test.go

@@ -0,0 +1,255 @@
+// 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 (
+	"bytes"
+	"crypto/elliptic"
+	"crypto/rand"
+	"fmt"
+	"math/big"
+	"testing"
+
+	"sources.truenas.cloud/code/nistec"
+)
+
+func TestNISTECAllocations(t *testing.T) {
+	if raceEnabled {
+		t.Skip("skipping allocation test in -race mode")
+	}
+	t.Run("P224", func(t *testing.T) {
+		if allocs := testing.AllocsPerRun(10, func() {
+			p := nistec.NewP224Point().SetGenerator()
+			scalar := make([]byte, 28)
+			rand.Read(scalar)
+			p.ScalarBaseMult(scalar)
+			p.ScalarMult(p, scalar)
+			out := p.Bytes()
+			if _, err := nistec.NewP224Point().SetBytes(out); err != nil {
+				t.Fatal(err)
+			}
+			out = p.BytesCompressed()
+			if _, err := p.SetBytes(out); err != nil {
+				t.Fatal(err)
+			}
+		}); allocs > 0 {
+			t.Errorf("expected zero allocations, got %0.1f", allocs)
+		}
+	})
+	t.Run("P256", func(t *testing.T) {
+		if allocs := testing.AllocsPerRun(10, func() {
+			p := nistec.NewP256Point().SetGenerator()
+			scalar := make([]byte, 32)
+			rand.Read(scalar)
+			p.ScalarBaseMult(scalar)
+			p.ScalarMult(p, scalar)
+			out := p.Bytes()
+			if _, err := nistec.NewP256Point().SetBytes(out); err != nil {
+				t.Fatal(err)
+			}
+			out = p.BytesCompressed()
+			if _, err := p.SetBytes(out); err != nil {
+				t.Fatal(err)
+			}
+		}); allocs > 0 {
+			t.Errorf("expected zero allocations, got %0.1f", allocs)
+		}
+	})
+	t.Run("P384", func(t *testing.T) {
+		if allocs := testing.AllocsPerRun(10, func() {
+			p := nistec.NewP384Point().SetGenerator()
+			scalar := make([]byte, 48)
+			rand.Read(scalar)
+			p.ScalarBaseMult(scalar)
+			p.ScalarMult(p, scalar)
+			out := p.Bytes()
+			if _, err := nistec.NewP384Point().SetBytes(out); err != nil {
+				t.Fatal(err)
+			}
+			out = p.BytesCompressed()
+			if _, err := p.SetBytes(out); err != nil {
+				t.Fatal(err)
+			}
+		}); allocs > 0 {
+			t.Errorf("expected zero allocations, got %0.1f", allocs)
+		}
+	})
+	t.Run("P521", func(t *testing.T) {
+		if allocs := testing.AllocsPerRun(10, func() {
+			p := nistec.NewP521Point().SetGenerator()
+			scalar := make([]byte, 66)
+			rand.Read(scalar)
+			p.ScalarBaseMult(scalar)
+			p.ScalarMult(p, scalar)
+			out := p.Bytes()
+			if _, err := nistec.NewP521Point().SetBytes(out); err != nil {
+				t.Fatal(err)
+			}
+			out = p.BytesCompressed()
+			if _, err := p.SetBytes(out); err != nil {
+				t.Fatal(err)
+			}
+		}); allocs > 0 {
+			t.Errorf("expected zero allocations, got %0.1f", allocs)
+		}
+	})
+}
+
+type nistPoint[T any] interface {
+	Bytes() []byte
+	SetGenerator() T
+	SetBytes([]byte) (T, error)
+	Add(T, T) T
+	Double(T) T
+	ScalarMult(T, []byte) (T, error)
+	ScalarBaseMult([]byte) (T, error)
+}
+
+func TestEquivalents(t *testing.T) {
+	t.Run("P224", func(t *testing.T) {
+		testEquivalents(t, nistec.NewP224Point, elliptic.P224())
+	})
+	t.Run("P256", func(t *testing.T) {
+		testEquivalents(t, nistec.NewP256Point, elliptic.P256())
+	})
+	t.Run("P384", func(t *testing.T) {
+		testEquivalents(t, nistec.NewP384Point, elliptic.P384())
+	})
+	t.Run("P521", func(t *testing.T) {
+		testEquivalents(t, nistec.NewP521Point, elliptic.P521())
+	})
+}
+
+func testEquivalents[P nistPoint[P]](t *testing.T, newPoint func() P, c elliptic.Curve) {
+	p := newPoint().SetGenerator()
+
+	elementSize := (c.Params().BitSize + 7) / 8
+	two := make([]byte, elementSize)
+	two[len(two)-1] = 2
+	nPlusTwo := make([]byte, elementSize)
+	new(big.Int).Add(c.Params().N, big.NewInt(2)).FillBytes(nPlusTwo)
+
+	p1 := newPoint().Double(p)
+	p2 := newPoint().Add(p, p)
+	p3, err := newPoint().ScalarMult(p, two)
+	fatalIfErr(t, err)
+	p4, err := newPoint().ScalarBaseMult(two)
+	fatalIfErr(t, err)
+	p5, err := newPoint().ScalarMult(p, nPlusTwo)
+	fatalIfErr(t, err)
+	p6, err := newPoint().ScalarBaseMult(nPlusTwo)
+	fatalIfErr(t, err)
+
+	if !bytes.Equal(p1.Bytes(), p2.Bytes()) {
+		t.Error("P+P != 2*P")
+	}
+	if !bytes.Equal(p1.Bytes(), p3.Bytes()) {
+		t.Error("P+P != [2]P")
+	}
+	if !bytes.Equal(p1.Bytes(), p4.Bytes()) {
+		t.Error("G+G != [2]G")
+	}
+	if !bytes.Equal(p1.Bytes(), p5.Bytes()) {
+		t.Error("P+P != [N+2]P")
+	}
+	if !bytes.Equal(p1.Bytes(), p6.Bytes()) {
+		t.Error("G+G != [N+2]G")
+	}
+}
+
+func TestScalarMult(t *testing.T) {
+	t.Run("P224", func(t *testing.T) {
+		testScalarMult(t, nistec.NewP224Point, elliptic.P224())
+	})
+	t.Run("P256", func(t *testing.T) {
+		testScalarMult(t, nistec.NewP256Point, elliptic.P256())
+	})
+	t.Run("P384", func(t *testing.T) {
+		testScalarMult(t, nistec.NewP384Point, elliptic.P384())
+	})
+	t.Run("P521", func(t *testing.T) {
+		testScalarMult(t, nistec.NewP521Point, elliptic.P521())
+	})
+}
+
+func testScalarMult[P nistPoint[P]](t *testing.T, newPoint func() P, c elliptic.Curve) {
+	G := newPoint().SetGenerator()
+	checkScalar := func(t *testing.T, scalar []byte) {
+		p1, err := newPoint().ScalarBaseMult(scalar)
+		fatalIfErr(t, err)
+		p2, err := newPoint().ScalarMult(G, scalar)
+		fatalIfErr(t, err)
+		if !bytes.Equal(p1.Bytes(), p2.Bytes()) {
+			t.Error("[k]G != ScalarBaseMult(k)")
+		}
+
+		expectInfinity := new(big.Int).Mod(new(big.Int).SetBytes(scalar), c.Params().N).Sign() == 0
+		if expectInfinity {
+			if !bytes.Equal(p1.Bytes(), newPoint().Bytes()) {
+				t.Error("ScalarBaseMult(k) != ∞")
+			}
+			if !bytes.Equal(p2.Bytes(), newPoint().Bytes()) {
+				t.Error("[k]G != ∞")
+			}
+		} else {
+			if bytes.Equal(p1.Bytes(), newPoint().Bytes()) {
+				t.Error("ScalarBaseMult(k) == ∞")
+			}
+			if bytes.Equal(p2.Bytes(), newPoint().Bytes()) {
+				t.Error("[k]G == ∞")
+			}
+		}
+
+		d := new(big.Int).SetBytes(scalar)
+		d.Sub(c.Params().N, d)
+		d.Mod(d, c.Params().N)
+		g1, err := newPoint().ScalarBaseMult(d.FillBytes(make([]byte, len(scalar))))
+		fatalIfErr(t, err)
+		g1.Add(g1, p1)
+		if !bytes.Equal(g1.Bytes(), newPoint().Bytes()) {
+			t.Error("[N - k]G + [k]G != ∞")
+		}
+	}
+
+	byteLen := len(c.Params().N.Bytes())
+	bitLen := c.Params().N.BitLen()
+	t.Run("0", func(t *testing.T) { checkScalar(t, make([]byte, byteLen)) })
+	t.Run("1", func(t *testing.T) {
+		checkScalar(t, big.NewInt(1).FillBytes(make([]byte, byteLen)))
+	})
+	t.Run("N-1", func(t *testing.T) {
+		checkScalar(t, new(big.Int).Sub(c.Params().N, big.NewInt(1)).Bytes())
+	})
+	t.Run("N", func(t *testing.T) { checkScalar(t, c.Params().N.Bytes()) })
+	t.Run("N+1", func(t *testing.T) {
+		checkScalar(t, new(big.Int).Add(c.Params().N, big.NewInt(1)).Bytes())
+	})
+	t.Run("all1s", func(t *testing.T) {
+		s := new(big.Int).Lsh(big.NewInt(1), uint(bitLen))
+		s.Sub(s, big.NewInt(1))
+		checkScalar(t, s.Bytes())
+	})
+	if testing.Short() {
+		return
+	}
+	for i := 0; i < bitLen; i++ {
+		t.Run(fmt.Sprintf("1<<%d", i), func(t *testing.T) {
+			s := new(big.Int).Lsh(big.NewInt(1), uint(i))
+			checkScalar(t, s.FillBytes(make([]byte, byteLen)))
+		})
+	}
+	for i := 0; i <= 64; i++ {
+		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
+			checkScalar(t, big.NewInt(int64(i)).FillBytes(make([]byte, byteLen)))
+		})
+	}
+	// Test N-64...N+64 since they risk overlapping with precomputed table values
+	// in the final additions.
+	for i := int64(-64); i <= 64; i++ {
+		t.Run(fmt.Sprintf("N%+d", i), func(t *testing.T) {
+			checkScalar(t, new(big.Int).Add(c.Params().N, big.NewInt(i)).Bytes())
+		})
+	}
+}

p224.go

@@ -0,0 +1,454 @@
+// Copyright 2022 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.
+
+// Code generated by generate.go. DO NOT EDIT.
+
+package nistec
+
+import (
+	"crypto/subtle"
+	"errors"
+	"sync"
+
+	"sources.truenas.cloud/code/nistec/internal/fiat"
+)
+
+// p224ElementLength is the length of an element of the base or scalar field,
+// which have the same bytes length for all NIST P curves.
+const p224ElementLength = 28
+
+// P224Point is a P224 point. The zero value is NOT valid.
+type P224Point struct {
+	// The point is represented in projective coordinates (X:Y:Z),
+	// where x = X/Z and y = Y/Z.
+	x, y, z *fiat.P224Element
+}
+
+// NewP224Point returns a new P224Point representing the point at infinity point.
+func NewP224Point() *P224Point {
+	return &P224Point{
+		x: new(fiat.P224Element),
+		y: new(fiat.P224Element).One(),
+		z: new(fiat.P224Element),
+	}
+}
+
+// SetGenerator sets p to the canonical generator and returns p.
+func (p *P224Point) SetGenerator() *P224Point {
+	p.x.SetBytes([]byte{0xb7, 0xe, 0xc, 0xbd, 0x6b, 0xb4, 0xbf, 0x7f, 0x32, 0x13, 0x90, 0xb9, 0x4a, 0x3, 0xc1, 0xd3, 0x56, 0xc2, 0x11, 0x22, 0x34, 0x32, 0x80, 0xd6, 0x11, 0x5c, 0x1d, 0x21})
+	p.y.SetBytes([]byte{0xbd, 0x37, 0x63, 0x88, 0xb5, 0xf7, 0x23, 0xfb, 0x4c, 0x22, 0xdf, 0xe6, 0xcd, 0x43, 0x75, 0xa0, 0x5a, 0x7, 0x47, 0x64, 0x44, 0xd5, 0x81, 0x99, 0x85, 0x0, 0x7e, 0x34})
+	p.z.One()
+	return p
+}
+
+// Set sets p = q and returns p.
+func (p *P224Point) Set(q *P224Point) *P224Point {
+	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 *P224Point) SetBytes(b []byte) (*P224Point, error) {
+	switch {
+	// Point at infinity.
+	case len(b) == 1 && b[0] == 0:
+		return p.Set(NewP224Point()), nil
+
+	// Uncompressed form.
+	case len(b) == 1+2*p224ElementLength && b[0] == 4:
+		x, err := new(fiat.P224Element).SetBytes(b[1 : 1+p224ElementLength])
+		if err != nil {
+			return nil, err
+		}
+		y, err := new(fiat.P224Element).SetBytes(b[1+p224ElementLength:])
+		if err != nil {
+			return nil, err
+		}
+		if err := p224CheckOnCurve(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+p224ElementLength && (b[0] == 2 || b[0] == 3):
+		x, err := new(fiat.P224Element).SetBytes(b[1:])
+		if err != nil {
+			return nil, err
+		}
+
+		// y² = x³ - 3x + b
+		y := p224Polynomial(new(fiat.P224Element), x)
+		if !p224Sqrt(y, y) {
+			return nil, errors.New("invalid P224 compressed point encoding")
+		}
+
+		// Select the positive or negative root, as indicated by the least
+		// significant bit, based on the encoding type byte.
+		otherRoot := new(fiat.P224Element)
+		otherRoot.Sub(otherRoot, y)
+		cond := y.Bytes()[p224ElementLength-1]&1 ^ b[0]&1
+		y.Select(otherRoot, y, int(cond))
+
+		p.x.Set(x)
+		p.y.Set(y)
+		p.z.One()
+		return p, nil
+
+	default:
+		return nil, errors.New("invalid P224 point encoding")
+	}
+}
+
+var _p224B *fiat.P224Element
+var _p224BOnce sync.Once
+
+func p224B() *fiat.P224Element {
+	_p224BOnce.Do(func() {
+		_p224B, _ = new(fiat.P224Element).SetBytes([]byte{0xb4, 0x5, 0xa, 0x85, 0xc, 0x4, 0xb3, 0xab, 0xf5, 0x41, 0x32, 0x56, 0x50, 0x44, 0xb0, 0xb7, 0xd7, 0xbf, 0xd8, 0xba, 0x27, 0xb, 0x39, 0x43, 0x23, 0x55, 0xff, 0xb4})
+	})
+	return _p224B
+}
+
+// p224Polynomial sets y2 to x³ - 3x + b, and returns y2.
+func p224Polynomial(y2, x *fiat.P224Element) *fiat.P224Element {
+	y2.Square(x)
+	y2.Mul(y2, x)
+
+	threeX := new(fiat.P224Element).Add(x, x)
+	threeX.Add(threeX, x)
+	y2.Sub(y2, threeX)
+
+	return y2.Add(y2, p224B())
+}
+
+func p224CheckOnCurve(x, y *fiat.P224Element) error {
+	// y² = x³ - 3x + b
+	rhs := p224Polynomial(new(fiat.P224Element), x)
+	lhs := new(fiat.P224Element).Square(y)
+	if rhs.Equal(lhs) != 1 {
+		return errors.New("P224 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 *P224Point) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [1 + 2*p224ElementLength]byte
+	return p.bytes(&out)
+}
+
+func (p *P224Point) bytes(out *[1 + 2*p224ElementLength]byte) []byte {
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new(fiat.P224Element).Invert(p.z)
+	x := new(fiat.P224Element).Mul(p.x, zinv)
+	y := new(fiat.P224Element).Mul(p.y, zinv)
+
+	buf := append(out[:0], 4)
+	buf = append(buf, x.Bytes()...)
+	buf = append(buf, y.Bytes()...)
+	return buf
+}
+
+// BytesX returns the encoding of the x-coordinate of p, as specified in SEC 1,
+// Version 2.0, Section 2.3.5, or an error if p is the point at infinity.
+func (p *P224Point) BytesX() ([]byte, error) {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p224ElementLength]byte
+	return p.bytesX(&out)
+}
+
+func (p *P224Point) bytesX(out *[p224ElementLength]byte) ([]byte, error) {
+	if p.z.IsZero() == 1 {
+		return nil, errors.New("P224 point is the point at infinity")
+	}
+
+	zinv := new(fiat.P224Element).Invert(p.z)
+	x := new(fiat.P224Element).Mul(p.x, zinv)
+
+	return append(out[:0], x.Bytes()...), nil
+}
+
+// BytesCompressed returns the compressed 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 *P224Point) BytesCompressed() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [1 + p224ElementLength]byte
+	return p.bytesCompressed(&out)
+}
+
+func (p *P224Point) bytesCompressed(out *[1 + p224ElementLength]byte) []byte {
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new(fiat.P224Element).Invert(p.z)
+	x := new(fiat.P224Element).Mul(p.x, zinv)
+	y := new(fiat.P224Element).Mul(p.y, zinv)
+
+	// Encode the sign of the y coordinate (indicated by the least significant
+	// bit) as the encoding type (2 or 3).
+	buf := append(out[:0], 2)
+	buf[0] |= y.Bytes()[p224ElementLength-1] & 1
+	buf = append(buf, x.Bytes()...)
+	return buf
+}
+
+// Add sets q = p1 + p2, and returns q. The points may overlap.
+func (q *P224Point) Add(p1, p2 *P224Point) *P224Point {
+	// 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.P224Element).Mul(p1.x, p2.x)  // t0 := X1 * X2
+	t1 := new(fiat.P224Element).Mul(p1.y, p2.y)  // t1 := Y1 * Y2
+	t2 := new(fiat.P224Element).Mul(p1.z, p2.z)  // t2 := Z1 * Z2
+	t3 := new(fiat.P224Element).Add(p1.x, p1.y)  // t3 := X1 + Y1
+	t4 := new(fiat.P224Element).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.P224Element).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.P224Element).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.P224Element).Mul(p224B(), 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(p224B(), 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 *P224Point) Double(p *P224Point) *P224Point {
+	// 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.P224Element).Square(p.x)      // t0 := X ^ 2
+	t1 := new(fiat.P224Element).Square(p.y)      // t1 := Y ^ 2
+	t2 := new(fiat.P224Element).Square(p.z)      // t2 := Z ^ 2
+	t3 := new(fiat.P224Element).Mul(p.x, p.y)    // t3 := X * Y
+	t3.Add(t3, t3)                               // t3 := t3 + t3
+	z3 := new(fiat.P224Element).Mul(p.x, p.z)    // Z3 := X * Z
+	z3.Add(z3, z3)                               // Z3 := Z3 + Z3
+	y3 := new(fiat.P224Element).Mul(p224B(), t2) // Y3 := b * t2
+	y3.Sub(y3, z3)                               // Y3 := Y3 - Z3
+	x3 := new(fiat.P224Element).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(p224B(), 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 *P224Point) Select(p1, p2 *P224Point, cond int) *P224Point {
+	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
+}
+
+// A p224Table holds the first 15 multiples of a point at offset -1, so [1]P
+// is at table[0], [15]P is at table[14], and [0]P is implicitly the identity
+// point.
+type p224Table [15]*P224Point
+
+// Select selects the n-th multiple of the table base point into p. It works in
+// constant time by iterating over every entry of the table. n must be in [0, 15].
+func (table *p224Table) Select(p *P224Point, n uint8) {
+	if n >= 16 {
+		panic("nistec: internal error: p224Table called with out-of-bounds value")
+	}
+	p.Set(NewP224Point())
+	for i := uint8(1); i < 16; i++ {
+		cond := subtle.ConstantTimeByteEq(i, n)
+		p.Select(table[i-1], p, cond)
+	}
+}
+
+// ScalarMult sets p = scalar * q, and returns p.
+func (p *P224Point) ScalarMult(q *P224Point, scalar []byte) (*P224Point, error) {
+	// Compute a p224Table for the base point q. The explicit NewP224Point
+	// calls get inlined, letting the allocations live on the stack.
+	var table = p224Table{NewP224Point(), NewP224Point(), NewP224Point(),
+		NewP224Point(), NewP224Point(), NewP224Point(), NewP224Point(),
+		NewP224Point(), NewP224Point(), NewP224Point(), NewP224Point(),
+		NewP224Point(), NewP224Point(), NewP224Point(), NewP224Point()}
+	table[0].Set(q)
+	for i := 1; i < 15; i += 2 {
+		table[i].Double(table[i/2])
+		table[i+1].Add(table[i], 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 := NewP224Point()
+	p.Set(NewP224Point())
+	for i, byte := range scalar {
+		// No need to double on the first iteration, as p is the identity at
+		// this point, and [N]∞ = ∞.
+		if i != 0 {
+			p.Double(p)
+			p.Double(p)
+			p.Double(p)
+			p.Double(p)
+		}
+
+		windowValue := byte >> 4
+		table.Select(t, windowValue)
+		p.Add(p, t)
+
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+
+		windowValue = byte & 0b1111
+		table.Select(t, windowValue)
+		p.Add(p, t)
+	}
+
+	return p, nil
+}
+
+var p224GeneratorTable *[p224ElementLength * 2]p224Table
+var p224GeneratorTableOnce sync.Once
+
+// generatorTable returns a sequence of p224Tables. The first table contains
+// multiples of G. Each successive table is the previous table doubled four
+// times.
+func (p *P224Point) generatorTable() *[p224ElementLength * 2]p224Table {
+	p224GeneratorTableOnce.Do(func() {
+		p224GeneratorTable = new([p224ElementLength * 2]p224Table)
+		base := NewP224Point().SetGenerator()
+		for i := 0; i < p224ElementLength*2; i++ {
+			p224GeneratorTable[i][0] = NewP224Point().Set(base)
+			for j := 1; j < 15; j++ {
+				p224GeneratorTable[i][j] = NewP224Point().Add(p224GeneratorTable[i][j-1], base)
+			}
+			base.Double(base)
+			base.Double(base)
+			base.Double(base)
+			base.Double(base)
+		}
+	})
+	return p224GeneratorTable
+}
+
+// ScalarBaseMult sets p = scalar * B, where B is the canonical generator, and
+// returns p.
+func (p *P224Point) ScalarBaseMult(scalar []byte) (*P224Point, error) {
+	if len(scalar) != p224ElementLength {
+		return nil, errors.New("invalid scalar length")
+	}
+	tables := p.generatorTable()
+
+	// This is also a scalar multiplication with a four-bit window like in
+	// ScalarMult, but in this case the doublings are precomputed. The value
+	// [windowValue]G added at iteration k would normally get doubled
+	// (totIterations-k)×4 times, but with a larger precomputation we can
+	// instead add [2^((totIterations-k)×4)][windowValue]G and avoid the
+	// doublings between iterations.
+	t := NewP224Point()
+	p.Set(NewP224Point())
+	tableIndex := len(tables) - 1
+	for _, byte := range scalar {
+		windowValue := byte >> 4
+		tables[tableIndex].Select(t, windowValue)
+		p.Add(p, t)
+		tableIndex--
+
+		windowValue = byte & 0b1111
+		tables[tableIndex].Select(t, windowValue)
+		p.Add(p, t)
+		tableIndex--
+	}
+
+	return p, nil
+}
+
+// p224Sqrt sets e to a square root of x. If x is not a square, p224Sqrt returns
+// false and e is unchanged. e and x can overlap.
+func p224Sqrt(e, x *fiat.P224Element) (isSquare bool) {
+	candidate := new(fiat.P224Element)
+	p224SqrtCandidate(candidate, x)
+	square := new(fiat.P224Element).Square(candidate)
+	if square.Equal(x) != 1 {
+		return false
+	}
+	e.Set(candidate)
+	return true
+}

p224_sqrt.go

@@ -0,0 +1,133 @@
+// Copyright 2022 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
+
+import (
+	"sync"
+
+	"sources.truenas.cloud/code/nistec/internal/fiat"
+)
+
+var p224GG *[96]fiat.P224Element
+var p224GGOnce sync.Once
+
+// p224SqrtCandidate sets r to a square root candidate for x. r and x must not overlap.
+func p224SqrtCandidate(r, x *fiat.P224Element) {
+	// Since p = 1 mod 4, we can't use the exponentiation by (p + 1) / 4 like
+	// for the other primes. Instead, implement a variation of Tonelli–Shanks.
+	// The constant-time implementation is adapted from Thomas Pornin's ecGFp5.
+	//
+	// https://github.com/pornin/ecgfp5/blob/82325b965/rust/src/field.rs#L337-L385
+
+	// p = q*2^n + 1 with q odd -> q = 2^128 - 1 and n = 96
+	// g^(2^n) = 1 -> g = 11 ^ q (where 11 is the smallest non-square)
+	// GG[j] = g^(2^j) for j = 0 to n-1
+
+	p224GGOnce.Do(func() {
+		p224GG = new([96]fiat.P224Element)
+		for i := range p224GG {
+			if i == 0 {
+				p224GG[i].SetBytes([]byte{0x6a, 0x0f, 0xec, 0x67,
+					0x85, 0x98, 0xa7, 0x92, 0x0c, 0x55, 0xb2, 0xd4,
+					0x0b, 0x2d, 0x6f, 0xfb, 0xbe, 0xa3, 0xd8, 0xce,
+					0xf3, 0xfb, 0x36, 0x32, 0xdc, 0x69, 0x1b, 0x74})
+			} else {
+				p224GG[i].Square(&p224GG[i-1])
+			}
+		}
+	})
+
+	// r <- x^((q+1)/2) = x^(2^127)
+	// v <- x^q = x^(2^128-1)
+
+	// Compute x^(2^127-1) first.
+	//
+	// The sequence of 10 multiplications and 126 squarings is derived from the
+	// following addition chain generated with github.com/mmcloughlin/addchain v0.4.0.
+	//
+	//	_10      = 2*1
+	//	_11      = 1 + _10
+	//	_110     = 2*_11
+	//	_111     = 1 + _110
+	//	_111000  = _111 << 3
+	//	_111111  = _111 + _111000
+	//	_1111110 = 2*_111111
+	//	_1111111 = 1 + _1111110
+	//	x12      = _1111110 << 5 + _111111
+	//	x24      = x12 << 12 + x12
+	//	i36      = x24 << 7
+	//	x31      = _1111111 + i36
+	//	x48      = i36 << 17 + x24
+	//	x96      = x48 << 48 + x48
+	//	return     x96 << 31 + x31
+	//
+	var t0 = new(fiat.P224Element)
+	var t1 = new(fiat.P224Element)
+
+	r.Square(x)
+	r.Mul(x, r)
+	r.Square(r)
+	r.Mul(x, r)
+	t0.Square(r)
+	for s := 1; s < 3; s++ {
+		t0.Square(t0)
+	}
+	t0.Mul(r, t0)
+	t1.Square(t0)
+	r.Mul(x, t1)
+	for s := 0; s < 5; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	t1.Square(t0)
+	for s := 1; s < 12; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	t1.Square(t0)
+	for s := 1; s < 7; s++ {
+		t1.Square(t1)
+	}
+	r.Mul(r, t1)
+	for s := 0; s < 17; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	t1.Square(t0)
+	for s := 1; s < 48; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	for s := 0; s < 31; s++ {
+		t0.Square(t0)
+	}
+	r.Mul(r, t0)
+
+	// v = x^(2^127-1)^2 * x
+	v := new(fiat.P224Element).Square(r)
+	v.Mul(v, x)
+
+	// r = x^(2^127-1) * x
+	r.Mul(r, x)
+
+	// for i = n-1 down to 1:
+	//     w = v^(2^(i-1))
+	//     if w == -1 then:
+	//         v <- v*GG[n-i]
+	//         r <- r*GG[n-i-1]
+
+	var p224MinusOne = new(fiat.P224Element).Sub(
+		new(fiat.P224Element), new(fiat.P224Element).One())
+
+	for i := 96 - 1; i >= 1; i-- {
+		w := new(fiat.P224Element).Set(v)
+		for j := 0; j < i-1; j++ {
+			w.Square(w)
+		}
+		cond := w.Equal(p224MinusOne)
+		v.Select(t0.Mul(v, &p224GG[96-i]), v, cond)
+		r.Select(t0.Mul(r, &p224GG[96-i-1]), r, cond)
+	}
+}

p256.go

@@ -0,0 +1,706 @@
+// Copyright 2022 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.
+
+//go:build (!amd64 && !arm64 && !ppc64le && !s390x) || purego
+
+package nistec
+
+import (
+	"crypto/subtle"
+	"errors"
+	"math/bits"
+	"sync"
+	"unsafe"
+
+	"golang.org/x/sys/cpu"
+	"sources.truenas.cloud/code/nistec/internal/byteorder"
+	"sources.truenas.cloud/code/nistec/internal/fiat"
+)
+
+// P256Point is a P-256 point. The zero value is NOT valid.
+type P256Point struct {
+	// The point is represented in projective coordinates (X:Y:Z), where x = X/Z
+	// and y = Y/Z. Infinity is (0:1:0).
+	//
+	// fiat.P256Element is a base field element in [0, P-1] in the Montgomery
+	// domain (with R 2²⁵⁶ and P 2²⁵⁶ - 2²²⁴ + 2¹⁹² + 2⁹⁶ - 1) as four limbs in
+	// little-endian order value.
+	x, y, z fiat.P256Element
+}
+
+// NewP256Point returns a new P256Point representing the point at infinity point.
+func NewP256Point() *P256Point {
+	p := &P256Point{}
+	p.y.One()
+	return p
+}
+
+// SetGenerator sets p to the canonical generator and returns p.
+func (p *P256Point) SetGenerator() *P256Point {
+	p.x.SetBytes([]byte{0x6b, 0x17, 0xd1, 0xf2, 0xe1, 0x2c, 0x42, 0x47, 0xf8, 0xbc, 0xe6, 0xe5, 0x63, 0xa4, 0x40, 0xf2, 0x77, 0x3, 0x7d, 0x81, 0x2d, 0xeb, 0x33, 0xa0, 0xf4, 0xa1, 0x39, 0x45, 0xd8, 0x98, 0xc2, 0x96})
+	p.y.SetBytes([]byte{0x4f, 0xe3, 0x42, 0xe2, 0xfe, 0x1a, 0x7f, 0x9b, 0x8e, 0xe7, 0xeb, 0x4a, 0x7c, 0xf, 0x9e, 0x16, 0x2b, 0xce, 0x33, 0x57, 0x6b, 0x31, 0x5e, 0xce, 0xcb, 0xb6, 0x40, 0x68, 0x37, 0xbf, 0x51, 0xf5})
+	p.z.One()
+	return p
+}
+
+// Set sets p = q and returns p.
+func (p *P256Point) Set(q *P256Point) *P256Point {
+	p.x.Set(&q.x)
+	p.y.Set(&q.y)
+	p.z.Set(&q.z)
+	return p
+}
+
+const p256ElementLength = 32
+const p256UncompressedLength = 1 + 2*p256ElementLength
+const p256CompressedLength = 1 + p256ElementLength
+
+// 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 *P256Point) SetBytes(b []byte) (*P256Point, error) {
+	switch {
+	// Point at infinity.
+	case len(b) == 1 && b[0] == 0:
+		return p.Set(NewP256Point()), nil
+
+	// Uncompressed form.
+	case len(b) == p256UncompressedLength && b[0] == 4:
+		x, err := new(fiat.P256Element).SetBytes(b[1 : 1+p256ElementLength])
+		if err != nil {
+			return nil, err
+		}
+		y, err := new(fiat.P256Element).SetBytes(b[1+p256ElementLength:])
+		if err != nil {
+			return nil, err
+		}
+		if err := p256CheckOnCurve(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) == p256CompressedLength && (b[0] == 2 || b[0] == 3):
+		x, err := new(fiat.P256Element).SetBytes(b[1:])
+		if err != nil {
+			return nil, err
+		}
+
+		// y² = x³ - 3x + b
+		y := p256Polynomial(new(fiat.P256Element), x)
+		if !p256Sqrt(y, y) {
+			return nil, errors.New("invalid P256 compressed point encoding")
+		}
+
+		// Select the positive or negative root, as indicated by the least
+		// significant bit, based on the encoding type byte.
+		otherRoot := new(fiat.P256Element)
+		otherRoot.Sub(otherRoot, y)
+		cond := y.Bytes()[p256ElementLength-1]&1 ^ b[0]&1
+		y.Select(otherRoot, y, int(cond))
+
+		p.x.Set(x)
+		p.y.Set(y)
+		p.z.One()
+		return p, nil
+
+	default:
+		return nil, errors.New("invalid P256 point encoding")
+	}
+}
+
+var _p256B *fiat.P256Element
+var _p256BOnce sync.Once
+
+func p256B() *fiat.P256Element {
+	_p256BOnce.Do(func() {
+		_p256B, _ = new(fiat.P256Element).SetBytes([]byte{0x5a, 0xc6, 0x35, 0xd8, 0xaa, 0x3a, 0x93, 0xe7, 0xb3, 0xeb, 0xbd, 0x55, 0x76, 0x98, 0x86, 0xbc, 0x65, 0x1d, 0x6, 0xb0, 0xcc, 0x53, 0xb0, 0xf6, 0x3b, 0xce, 0x3c, 0x3e, 0x27, 0xd2, 0x60, 0x4b})
+	})
+	return _p256B
+}
+
+// p256Polynomial sets y2 to x³ - 3x + b, and returns y2.
+func p256Polynomial(y2, x *fiat.P256Element) *fiat.P256Element {
+	y2.Square(x)
+	y2.Mul(y2, x)
+
+	threeX := new(fiat.P256Element).Add(x, x)
+	threeX.Add(threeX, x)
+	y2.Sub(y2, threeX)
+
+	return y2.Add(y2, p256B())
+}
+
+func p256CheckOnCurve(x, y *fiat.P256Element) error {
+	// y² = x³ - 3x + b
+	rhs := p256Polynomial(new(fiat.P256Element), x)
+	lhs := new(fiat.P256Element).Square(y)
+	if rhs.Equal(lhs) != 1 {
+		return errors.New("P256 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 *P256Point) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p256UncompressedLength]byte
+	return p.bytes(&out)
+}
+
+func (p *P256Point) bytes(out *[p256UncompressedLength]byte) []byte {
+	// The SEC 1 representation of the point at infinity is a single zero byte,
+	// and only infinity has z = 0.
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new(fiat.P256Element).Invert(&p.z)
+	x := new(fiat.P256Element).Mul(&p.x, zinv)
+	y := new(fiat.P256Element).Mul(&p.y, zinv)
+
+	buf := append(out[:0], 4)
+	buf = append(buf, x.Bytes()...)
+	buf = append(buf, y.Bytes()...)
+	return buf
+}
+
+// BytesX returns the encoding of the x-coordinate of p, as specified in SEC 1,
+// Version 2.0, Section 2.3.5, or an error if p is the point at infinity.
+func (p *P256Point) BytesX() ([]byte, error) {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p256ElementLength]byte
+	return p.bytesX(&out)
+}
+
+func (p *P256Point) bytesX(out *[p256ElementLength]byte) ([]byte, error) {
+	if p.z.IsZero() == 1 {
+		return nil, errors.New("P256 point is the point at infinity")
+	}
+
+	zinv := new(fiat.P256Element).Invert(&p.z)
+	x := new(fiat.P256Element).Mul(&p.x, zinv)
+
+	return append(out[:0], x.Bytes()...), nil
+}
+
+// BytesCompressed returns the compressed 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 *P256Point) BytesCompressed() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p256CompressedLength]byte
+	return p.bytesCompressed(&out)
+}
+
+func (p *P256Point) bytesCompressed(out *[p256CompressedLength]byte) []byte {
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new(fiat.P256Element).Invert(&p.z)
+	x := new(fiat.P256Element).Mul(&p.x, zinv)
+	y := new(fiat.P256Element).Mul(&p.y, zinv)
+
+	// Encode the sign of the y coordinate (indicated by the least significant
+	// bit) as the encoding type (2 or 3).
+	buf := append(out[:0], 2)
+	buf[0] |= y.Bytes()[p256ElementLength-1] & 1
+	buf = append(buf, x.Bytes()...)
+	return buf
+}
+
+// Add sets q = p1 + p2, and returns q. The points may overlap.
+func (q *P256Point) Add(p1, p2 *P256Point) *P256Point {
+	// 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.P256Element).Mul(&p1.x, &p2.x) // t0 := X1 * X2
+	t1 := new(fiat.P256Element).Mul(&p1.y, &p2.y) // t1 := Y1 * Y2
+	t2 := new(fiat.P256Element).Mul(&p1.z, &p2.z) // t2 := Z1 * Z2
+	t3 := new(fiat.P256Element).Add(&p1.x, &p1.y) // t3 := X1 + Y1
+	t4 := new(fiat.P256Element).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.P256Element).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.P256Element).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.P256Element).Mul(p256B(), 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(p256B(), 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 *P256Point) Double(p *P256Point) *P256Point {
+	// 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.P256Element).Square(&p.x)     // t0 := X ^ 2
+	t1 := new(fiat.P256Element).Square(&p.y)     // t1 := Y ^ 2
+	t2 := new(fiat.P256Element).Square(&p.z)     // t2 := Z ^ 2
+	t3 := new(fiat.P256Element).Mul(&p.x, &p.y)  // t3 := X * Y
+	t3.Add(t3, t3)                               // t3 := t3 + t3
+	z3 := new(fiat.P256Element).Mul(&p.x, &p.z)  // Z3 := X * Z
+	z3.Add(z3, z3)                               // Z3 := Z3 + Z3
+	y3 := new(fiat.P256Element).Mul(p256B(), t2) // Y3 := b * t2
+	y3.Sub(y3, z3)                               // Y3 := Y3 - Z3
+	x3 := new(fiat.P256Element).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(p256B(), 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
+}
+
+// p256AffinePoint is a point in affine coordinates (x, y). x and y are still
+// Montgomery domain elements. The point can't be the point at infinity.
+type p256AffinePoint struct {
+	x, y fiat.P256Element
+}
+
+func (p *p256AffinePoint) Projective() *P256Point {
+	pp := &P256Point{x: p.x, y: p.y}
+	pp.z.One()
+	return pp
+}
+
+// AddAffine sets q = p1 + p2, if infinity == 0, and to p1 if infinity == 1.
+// p2 can't be the point at infinity as it can't be represented in affine
+// coordinates, instead callers can set p2 to an arbitrary point and set
+// infinity to 1.
+func (q *P256Point) AddAffine(p1 *P256Point, p2 *p256AffinePoint, infinity int) *P256Point {
+	// Complete mixed addition formula for a = -3 from "Complete addition
+	// formulas for prime order elliptic curves"
+	// (https://eprint.iacr.org/2015/1060), Algorithm 5.
+
+	t0 := new(fiat.P256Element).Mul(&p1.x, &p2.x)   // t0 ← X1 · X2
+	t1 := new(fiat.P256Element).Mul(&p1.y, &p2.y)   // t1 ← Y1 · Y2
+	t3 := new(fiat.P256Element).Add(&p2.x, &p2.y)   // t3 ← X2 + Y2
+	t4 := new(fiat.P256Element).Add(&p1.x, &p1.y)   // t4 ← X1 + Y1
+	t3.Mul(t3, t4)                                  // t3 ← t3 · t4
+	t4.Add(t0, t1)                                  // t4 ← t0 + t1
+	t3.Sub(t3, t4)                                  // t3 ← t3 − t4
+	t4.Mul(&p2.y, &p1.z)                            // t4 ← Y2 · Z1
+	t4.Add(t4, &p1.y)                               // t4 ← t4 + Y1
+	y3 := new(fiat.P256Element).Mul(&p2.x, &p1.z)   // Y3 ← X2 · Z1
+	y3.Add(y3, &p1.x)                               // Y3 ← Y3 + X1
+	z3 := new(fiat.P256Element).Mul(p256B(), &p1.z) // Z3 ← b  · Z1
+	x3 := new(fiat.P256Element).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(p256B(), y3)                             // Y3 ← b  · Y3
+	t1.Add(&p1.z, &p1.z)                            // t1 ← Z1 + Z1
+	t2 := new(fiat.P256Element).Add(t1, &p1.z)      // t2 ← t1 + Z1
+	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.Select(&p1.x, x3, infinity)
+	q.y.Select(&p1.y, y3, infinity)
+	q.z.Select(&p1.z, z3, infinity)
+	return q
+}
+
+// Select sets q to p1 if cond == 1, and to p2 if cond == 0.
+func (q *P256Point) Select(p1, p2 *P256Point, cond int) *P256Point {
+	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
+}
+
+// p256OrdElement is a P-256 scalar field element in [0, ord(G)-1] in the
+// Montgomery domain (with R 2²⁵⁶) as four uint64 limbs in little-endian order.
+type p256OrdElement [4]uint64
+
+// SetBytes sets s to the big-endian value of x, reducing it as necessary.
+func (s *p256OrdElement) SetBytes(x []byte) (*p256OrdElement, error) {
+	if len(x) != 32 {
+		return nil, errors.New("invalid scalar length")
+	}
+
+	s[0] = byteorder.BEUint64(x[24:])
+	s[1] = byteorder.BEUint64(x[16:])
+	s[2] = byteorder.BEUint64(x[8:])
+	s[3] = byteorder.BEUint64(x[:])
+
+	// Ensure s is in the range [0, ord(G)-1]. Since 2 * ord(G) > 2²⁵⁶, we can
+	// just conditionally subtract ord(G), keeping the result if it doesn't
+	// underflow.
+	t0, b := bits.Sub64(s[0], 0xf3b9cac2fc632551, 0)
+	t1, b := bits.Sub64(s[1], 0xbce6faada7179e84, b)
+	t2, b := bits.Sub64(s[2], 0xffffffffffffffff, b)
+	t3, b := bits.Sub64(s[3], 0xffffffff00000000, b)
+	tMask := b - 1 // zero if subtraction underflowed
+	s[0] ^= (t0 ^ s[0]) & tMask
+	s[1] ^= (t1 ^ s[1]) & tMask
+	s[2] ^= (t2 ^ s[2]) & tMask
+	s[3] ^= (t3 ^ s[3]) & tMask
+
+	return s, nil
+}
+
+func (s *p256OrdElement) Bytes() []byte {
+	var out [32]byte
+	byteorder.BEPutUint64(out[24:], s[0])
+	byteorder.BEPutUint64(out[16:], s[1])
+	byteorder.BEPutUint64(out[8:], s[2])
+	byteorder.BEPutUint64(out[:], s[3])
+	return out[:]
+}
+
+// Rsh returns the 64 least significant bits of x >> n. n must be lower
+// than 256. The value of n leaks through timing side-channels.
+func (s *p256OrdElement) Rsh(n int) uint64 {
+	i := n / 64
+	n = n % 64
+	res := s[i] >> n
+	// Shift in the more significant limb, if present.
+	if i := i + 1; i < len(s) {
+		res |= s[i] << (64 - n)
+	}
+	return res
+}
+
+// p256Table is a table of the first 16 multiples of a point. Points are stored
+// at an index offset of -1 so [8]P is at index 7, P is at 0, and [16]P is at 15.
+// [0]P is the point at infinity and it's not stored.
+type p256Table [16]P256Point
+
+// Select selects the n-th multiple of the table base point into p. It works in
+// constant time. n must be in [0, 16]. If n is 0, p is set to the identity point.
+func (table *p256Table) Select(p *P256Point, n uint8) {
+	if n > 16 {
+		panic("nistec: internal error: p256Table called with out-of-bounds value")
+	}
+	p.Set(NewP256Point())
+	for i := uint8(1); i <= 16; i++ {
+		cond := subtle.ConstantTimeByteEq(i, n)
+		p.Select(&table[i-1], p, cond)
+	}
+}
+
+// Compute populates the table to the first 16 multiples of q.
+func (table *p256Table) Compute(q *P256Point) *p256Table {
+	table[0].Set(q)
+	for i := 1; i < 16; i += 2 {
+		table[i].Double(&table[i/2])
+		if i+1 < 16 {
+			table[i+1].Add(&table[i], q)
+		}
+	}
+	return table
+}
+
+func boothW5(in uint64) (uint8, int) {
+	s := ^((in >> 5) - 1)
+	d := (1 << 6) - in - 1
+	d = (d & s) | (in & (^s))
+	d = (d >> 1) + (d & 1)
+	return uint8(d), int(s & 1)
+}
+
+// ScalarMult sets r = scalar * q, where scalar is a 32-byte big endian value,
+// and returns r. If scalar is not 32 bytes long, ScalarMult returns an error
+// and the receiver is unchanged.
+func (p *P256Point) ScalarMult(q *P256Point, scalar []byte) (*P256Point, error) {
+	s, err := new(p256OrdElement).SetBytes(scalar)
+	if err != nil {
+		return nil, err
+	}
+
+	// Start scanning the window from the most significant bits. We move by
+	// 5 bits at a time and need to finish at -1, so -1 + 5 * 51 = 254.
+	index := 254
+
+	sel, sign := boothW5(s.Rsh(index))
+	// sign is always zero because the boothW5 input here is at
+	// most two bits long, so the top bit is never set.
+	_ = sign
+
+	// Neither Select nor Add have exceptions for the point at infinity /
+	// selector zero, so we don't need to check for it here or in the loop.
+	table := new(p256Table).Compute(q)
+	table.Select(p, sel)
+
+	t := NewP256Point()
+	for index >= 4 {
+		index -= 5
+
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+
+		if index >= 0 {
+			sel, sign = boothW5(s.Rsh(index) & 0b111111)
+		} else {
+			// Booth encoding considers a virtual zero bit at index -1,
+			// so we shift left the least significant limb.
+			wvalue := (s[0] << 1) & 0b111111
+			sel, sign = boothW5(wvalue)
+		}
+
+		table.Select(t, sel)
+		p256Negate(t, sign)
+		p.Add(p, t)
+	}
+
+	return p, nil
+}
+
+// Negate sets p to -p, if cond == 1, and to p if cond == 0.
+func p256Negate(p *P256Point, cond int) *P256Point {
+	negY := new(fiat.P256Element)
+	negY.Sub(negY, &p.y)
+	p.y.Select(negY, &p.y, cond)
+	return p
+}
+
+// p256AffineTable is a table of the first 32 multiples of a point. Points are
+// stored at an index offset of -1 like in p256Table, and [0]P is not stored.
+type p256AffineTable [32]p256AffinePoint
+
+// Select selects the n-th multiple of the table base point into p. It works in
+// constant time. n can be in [0, 32], but (unlike p256Table.Select) if n is 0,
+// p is set to an undefined value.
+func (table *p256AffineTable) Select(p *p256AffinePoint, n uint8) {
+	if n > 32 {
+		panic("nistec: internal error: p256AffineTable.Select called with out-of-bounds value")
+	}
+	for i := uint8(1); i <= 32; i++ {
+		cond := subtle.ConstantTimeByteEq(i, n)
+		p.x.Select(&table[i-1].x, &p.x, cond)
+		p.y.Select(&table[i-1].y, &p.y, cond)
+	}
+}
+
+// p256GeneratorTables is a series of precomputed multiples of G, the canonical
+// generator. The first p256AffineTable contains multiples of G. The second one
+// multiples of [2⁶]G, the third one of [2¹²]G, and so on, where each successive
+// table is the previous table doubled six times. Six is the width of the
+// sliding window used in ScalarBaseMult, and having each table already
+// pre-doubled lets us avoid the doublings between windows entirely. This table
+// aliases into p256PrecomputedEmbed.
+var p256GeneratorTables *[43]p256AffineTable
+
+func init() {
+	p256GeneratorTablesPtr := unsafe.Pointer(&p256PrecomputedEmbed)
+	if cpu.IsBigEndian {
+		var newTable [43 * 32 * 2 * 4]uint64
+		for i, x := range (*[43 * 32 * 2 * 4][8]byte)(p256GeneratorTablesPtr) {
+			newTable[i] = byteorder.LEUint64(x[:])
+		}
+		p256GeneratorTablesPtr = unsafe.Pointer(&newTable)
+	}
+	p256GeneratorTables = (*[43]p256AffineTable)(p256GeneratorTablesPtr)
+}
+
+func boothW6(in uint64) (uint8, int) {
+	s := ^((in >> 6) - 1)
+	d := (1 << 7) - in - 1
+	d = (d & s) | (in & (^s))
+	d = (d >> 1) + (d & 1)
+	return uint8(d), int(s & 1)
+}
+
+// ScalarBaseMult sets p = scalar * generator, where scalar is a 32-byte big
+// endian value, and returns r. If scalar is not 32 bytes long, ScalarBaseMult
+// returns an error and the receiver is unchanged.
+func (p *P256Point) ScalarBaseMult(scalar []byte) (*P256Point, error) {
+	// This function works like ScalarMult above, but the table is fixed and
+	// "pre-doubled" for each iteration, so instead of doubling we move to the
+	// next table at each iteration.
+
+	s, err := new(p256OrdElement).SetBytes(scalar)
+	if err != nil {
+		return nil, err
+	}
+
+	// Start scanning the window from the most significant bits. We move by
+	// 6 bits at a time and need to finish at -1, so -1 + 6 * 42 = 251.
+	index := 251
+
+	sel, sign := boothW6(s.Rsh(index))
+	// sign is always zero because the boothW6 input here is at
+	// most five bits long, so the top bit is never set.
+	_ = sign
+
+	t := &p256AffinePoint{}
+	table := &p256GeneratorTables[(index+1)/6]
+	table.Select(t, sel)
+
+	// Select's output is undefined if the selector is zero, when it should be
+	// the point at infinity (because infinity can't be represented in affine
+	// coordinates). Here we conditionally set p to the infinity if sel is zero.
+	// In the loop, that's handled by AddAffine.
+	selIsZero := subtle.ConstantTimeByteEq(sel, 0)
+	p.Select(NewP256Point(), t.Projective(), selIsZero)
+
+	for index >= 5 {
+		index -= 6
+
+		if index >= 0 {
+			sel, sign = boothW6(s.Rsh(index) & 0b1111111)
+		} else {
+			// Booth encoding considers a virtual zero bit at index -1,
+			// so we shift left the least significant limb.
+			wvalue := (s[0] << 1) & 0b1111111
+			sel, sign = boothW6(wvalue)
+		}
+
+		table := &p256GeneratorTables[(index+1)/6]
+		table.Select(t, sel)
+		t.Negate(sign)
+		selIsZero := subtle.ConstantTimeByteEq(sel, 0)
+		p.AddAffine(p, t, selIsZero)
+	}
+
+	return p, nil
+}
+
+// Negate sets p to -p, if cond == 1, and to p if cond == 0.
+func (p *p256AffinePoint) Negate(cond int) *p256AffinePoint {
+	negY := new(fiat.P256Element)
+	negY.Sub(negY, &p.y)
+	p.y.Select(negY, &p.y, cond)
+	return p
+}
+
+// p256Sqrt sets e to a square root of x. If x is not a square, p256Sqrt returns
+// false and e is unchanged. e and x can overlap.
+func p256Sqrt(e, x *fiat.P256Element) (isSquare bool) {
+	t0, t1 := new(fiat.P256Element), new(fiat.P256Element)
+
+	// Since p = 3 mod 4, exponentiation by (p + 1) / 4 yields a square root candidate.
+	//
+	// The sequence of 7 multiplications and 253 squarings is derived from the
+	// following addition chain generated with github.com/mmcloughlin/addchain v0.4.0.
+	//
+	//	_10       = 2*1
+	//	_11       = 1 + _10
+	//	_1100     = _11 << 2
+	//	_1111     = _11 + _1100
+	//	_11110000 = _1111 << 4
+	//	_11111111 = _1111 + _11110000
+	//	x16       = _11111111 << 8 + _11111111
+	//	x32       = x16 << 16 + x16
+	//	return      ((x32 << 32 + 1) << 96 + 1) << 94
+	//
+	p256Square(t0, x, 1)
+	t0.Mul(x, t0)
+	p256Square(t1, t0, 2)
+	t0.Mul(t0, t1)
+	p256Square(t1, t0, 4)
+	t0.Mul(t0, t1)
+	p256Square(t1, t0, 8)
+	t0.Mul(t0, t1)
+	p256Square(t1, t0, 16)
+	t0.Mul(t0, t1)
+	p256Square(t0, t0, 32)
+	t0.Mul(x, t0)
+	p256Square(t0, t0, 96)
+	t0.Mul(x, t0)
+	p256Square(t0, t0, 94)
+
+	// Check if the candidate t0 is indeed a square root of x.
+	t1.Square(t0)
+	if t1.Equal(x) != 1 {
+		return false
+	}
+	e.Set(t0)
+	return true
+}
+
+// p256Square sets e to the square of x, repeated n times > 1.
+func p256Square(e, x *fiat.P256Element, n int) {
+	e.Square(x)
+	for i := 1; i < n; i++ {
+		e.Square(e)
+	}
+}

p256_asm.go

@@ -0,0 +1,758 @@
+// Copyright 2015 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.
+
+// This file contains the Go wrapper for the constant-time, 64-bit assembly
+// implementation of P256. The optimizations performed here are described in
+// detail in:
+// S.Gueron and V.Krasnov, "Fast prime field elliptic-curve cryptography with
+//                          256-bit primes"
+// https://link.springer.com/article/10.1007%2Fs13389-014-0090-x
+// https://eprint.iacr.org/2013/816.pdf
+
+//go:build (amd64 || arm64 || ppc64le || s390x) && !purego
+
+package nistec
+
+import (
+	"errors"
+	"math/bits"
+	"runtime"
+	"unsafe"
+
+	"sources.truenas.cloud/code/nistec/internal/byteorder"
+)
+
+// p256Element is a P-256 base field element in [0, P-1] in the Montgomery
+// domain (with R 2²⁵⁶) as four limbs in little-endian order value.
+type p256Element [4]uint64
+
+// p256One is one in the Montgomery domain.
+var p256One = p256Element{0x0000000000000001, 0xffffffff00000000,
+	0xffffffffffffffff, 0x00000000fffffffe}
+
+var p256Zero = p256Element{}
+
+// p256P is 2²⁵⁶ - 2²²⁴ + 2¹⁹² + 2⁹⁶ - 1 in the Montgomery domain.
+var p256P = p256Element{0xffffffffffffffff, 0x00000000ffffffff,
+	0x0000000000000000, 0xffffffff00000001}
+
+// P256Point is a P-256 point. The zero value should not be assumed to be valid
+// (although it is in this implementation).
+type P256Point struct {
+	// (X:Y:Z) are Jacobian coordinates where x = X/Z² and y = Y/Z³. The point
+	// at infinity can be represented by any set of coordinates with Z = 0.
+	x, y, z p256Element
+}
+
+// NewP256Point returns a new P256Point representing the point at infinity.
+func NewP256Point() *P256Point {
+	return &P256Point{
+		x: p256One, y: p256One, z: p256Zero,
+	}
+}
+
+// SetGenerator sets p to the canonical generator and returns p.
+func (p *P256Point) SetGenerator() *P256Point {
+	p.x = p256Element{0x79e730d418a9143c, 0x75ba95fc5fedb601,
+		0x79fb732b77622510, 0x18905f76a53755c6}
+	p.y = p256Element{0xddf25357ce95560a, 0x8b4ab8e4ba19e45c,
+		0xd2e88688dd21f325, 0x8571ff1825885d85}
+	p.z = p256One
+	return p
+}
+
+// Set sets p = q and returns p.
+func (p *P256Point) Set(q *P256Point) *P256Point {
+	p.x, p.y, p.z = q.x, q.y, q.z
+	return p
+}
+
+const p256ElementLength = 32
+const p256UncompressedLength = 1 + 2*p256ElementLength
+const p256CompressedLength = 1 + p256ElementLength
+
+// 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 *P256Point) SetBytes(b []byte) (*P256Point, error) {
+	// p256Mul operates in the Montgomery domain with R = 2²⁵⁶ mod p. Thus rr
+	// here is R in the Montgomery domain, or R×R mod p. See comment in
+	// P256OrdInverse about how this is used.
+	rr := p256Element{0x0000000000000003, 0xfffffffbffffffff,
+		0xfffffffffffffffe, 0x00000004fffffffd}
+
+	switch {
+	// Point at infinity.
+	case len(b) == 1 && b[0] == 0:
+		return p.Set(NewP256Point()), nil
+
+	// Uncompressed form.
+	case len(b) == p256UncompressedLength && b[0] == 4:
+		var r P256Point
+		p256BigToLittle(&r.x, (*[32]byte)(b[1:33]))
+		p256BigToLittle(&r.y, (*[32]byte)(b[33:65]))
+		if p256LessThanP(&r.x) == 0 || p256LessThanP(&r.y) == 0 {
+			return nil, errors.New("invalid P256 element encoding")
+		}
+		p256Mul(&r.x, &r.x, &rr)
+		p256Mul(&r.y, &r.y, &rr)
+		if err := p256CheckOnCurve(&r.x, &r.y); err != nil {
+			return nil, err
+		}
+		r.z = p256One
+		return p.Set(&r), nil
+
+	// Compressed form.
+	case len(b) == p256CompressedLength && (b[0] == 2 || b[0] == 3):
+		var r P256Point
+		p256BigToLittle(&r.x, (*[32]byte)(b[1:33]))
+		if p256LessThanP(&r.x) == 0 {
+			return nil, errors.New("invalid P256 element encoding")
+		}
+		p256Mul(&r.x, &r.x, &rr)
+
+		// y² = x³ - 3x + b
+		p256Polynomial(&r.y, &r.x)
+		if !p256Sqrt(&r.y, &r.y) {
+			return nil, errors.New("invalid P256 compressed point encoding")
+		}
+
+		// Select the positive or negative root, as indicated by the least
+		// significant bit, based on the encoding type byte.
+		yy := new(p256Element)
+		p256FromMont(yy, &r.y)
+		cond := int(yy[0]&1) ^ int(b[0]&1)
+		p256NegCond(&r.y, cond)
+
+		r.z = p256One
+		return p.Set(&r), nil
+
+	default:
+		return nil, errors.New("invalid P256 point encoding")
+	}
+}
+
+// p256Polynomial sets y2 to x³ - 3x + b, and returns y2.
+func p256Polynomial(y2, x *p256Element) *p256Element {
+	x3 := new(p256Element)
+	p256Sqr(x3, x, 1)
+	p256Mul(x3, x3, x)
+
+	threeX := new(p256Element)
+	p256Add(threeX, x, x)
+	p256Add(threeX, threeX, x)
+	p256NegCond(threeX, 1)
+
+	p256B := &p256Element{0xd89cdf6229c4bddf, 0xacf005cd78843090,
+		0xe5a220abf7212ed6, 0xdc30061d04874834}
+
+	p256Add(x3, x3, threeX)
+	p256Add(x3, x3, p256B)
+
+	*y2 = *x3
+	return y2
+}
+
+func p256CheckOnCurve(x, y *p256Element) error {
+	// y² = x³ - 3x + b
+	rhs := p256Polynomial(new(p256Element), x)
+	lhs := new(p256Element)
+	p256Sqr(lhs, y, 1)
+	if p256Equal(lhs, rhs) != 1 {
+		return errors.New("P256 point not on curve")
+	}
+	return nil
+}
+
+// p256LessThanP returns 1 if x < p, and 0 otherwise. Note that a p256Element is
+// not allowed to be equal to or greater than p, so if this function returns 0
+// then x is invalid.
+func p256LessThanP(x *p256Element) int {
+	var b uint64
+	_, b = bits.Sub64(x[0], p256P[0], b)
+	_, b = bits.Sub64(x[1], p256P[1], b)
+	_, b = bits.Sub64(x[2], p256P[2], b)
+	_, b = bits.Sub64(x[3], p256P[3], b)
+	return int(b)
+}
+
+func p256BigToLittle(l *p256Element, b *[32]byte) {
+	bytesToLimbs((*[4]uint64)(l), b)
+}
+
+func bytesToLimbs(l *[4]uint64, b *[32]byte) {
+	l[0] = byteorder.BEUint64(b[24:])
+	l[1] = byteorder.BEUint64(b[16:])
+	l[2] = byteorder.BEUint64(b[8:])
+	l[3] = byteorder.BEUint64(b[:])
+}
+
+func p256LittleToBig(b *[32]byte, l *p256Element) {
+	limbsToBytes(b, (*[4]uint64)(l))
+}
+
+func limbsToBytes(b *[32]byte, l *[4]uint64) {
+	byteorder.BEPutUint64(b[24:], l[0])
+	byteorder.BEPutUint64(b[16:], l[1])
+	byteorder.BEPutUint64(b[8:], l[2])
+	byteorder.BEPutUint64(b[:], l[3])
+}
+
+// p256Add sets res = x + y.
+func p256Add(res, x, y *p256Element) {
+	var c, b uint64
+	t1 := make([]uint64, 4)
+	t1[0], c = bits.Add64(x[0], y[0], 0)
+	t1[1], c = bits.Add64(x[1], y[1], c)
+	t1[2], c = bits.Add64(x[2], y[2], c)
+	t1[3], c = bits.Add64(x[3], y[3], c)
+	t2 := make([]uint64, 4)
+	t2[0], b = bits.Sub64(t1[0], p256P[0], 0)
+	t2[1], b = bits.Sub64(t1[1], p256P[1], b)
+	t2[2], b = bits.Sub64(t1[2], p256P[2], b)
+	t2[3], b = bits.Sub64(t1[3], p256P[3], b)
+	// Three options:
+	//   - a+b < p
+	//     then c is 0, b is 1, and t1 is correct
+	//   - p <= a+b < 2^256
+	//     then c is 0, b is 0, and t2 is correct
+	//   - 2^256 <= a+b
+	//     then c is 1, b is 1, and t2 is correct
+	t2Mask := (c ^ b) - 1
+	res[0] = (t1[0] & ^t2Mask) | (t2[0] & t2Mask)
+	res[1] = (t1[1] & ^t2Mask) | (t2[1] & t2Mask)
+	res[2] = (t1[2] & ^t2Mask) | (t2[2] & t2Mask)
+	res[3] = (t1[3] & ^t2Mask) | (t2[3] & t2Mask)
+}
+
+// p256Sqrt sets e to a square root of x. If x is not a square, p256Sqrt returns
+// false and e is unchanged. e and x can overlap.
+func p256Sqrt(e, x *p256Element) (isSquare bool) {
+	t0, t1 := new(p256Element), new(p256Element)
+
+	// Since p = 3 mod 4, exponentiation by (p + 1) / 4 yields a square root candidate.
+	//
+	// The sequence of 7 multiplications and 253 squarings is derived from the
+	// following addition chain generated with github.com/mmcloughlin/addchain v0.4.0.
+	//
+	//	_10       = 2*1
+	//	_11       = 1 + _10
+	//	_1100     = _11 << 2
+	//	_1111     = _11 + _1100
+	//	_11110000 = _1111 << 4
+	//	_11111111 = _1111 + _11110000
+	//	x16       = _11111111 << 8 + _11111111
+	//	x32       = x16 << 16 + x16
+	//	return      ((x32 << 32 + 1) << 96 + 1) << 94
+	//
+	p256Sqr(t0, x, 1)
+	p256Mul(t0, x, t0)
+	p256Sqr(t1, t0, 2)
+	p256Mul(t0, t0, t1)
+	p256Sqr(t1, t0, 4)
+	p256Mul(t0, t0, t1)
+	p256Sqr(t1, t0, 8)
+	p256Mul(t0, t0, t1)
+	p256Sqr(t1, t0, 16)
+	p256Mul(t0, t0, t1)
+	p256Sqr(t0, t0, 32)
+	p256Mul(t0, x, t0)
+	p256Sqr(t0, t0, 96)
+	p256Mul(t0, x, t0)
+	p256Sqr(t0, t0, 94)
+
+	p256Sqr(t1, t0, 1)
+	if p256Equal(t1, x) != 1 {
+		return false
+	}
+	*e = *t0
+	return true
+}
+
+// The following assembly functions are implemented in p256_asm_*.s
+
+// Montgomery multiplication. Sets res = in1 * in2 * R⁻¹ mod p.
+//
+//go:noescape
+func p256Mul(res, in1, in2 *p256Element)
+
+// Montgomery square, repeated n times (n >= 1).
+//
+//go:noescape
+func p256Sqr(res, in *p256Element, n int)
+
+// Montgomery multiplication by R⁻¹, or 1 outside the domain.
+// Sets res = in * R⁻¹, bringing res out of the Montgomery domain.
+//
+//go:noescape
+func p256FromMont(res, in *p256Element)
+
+// If cond is not 0, sets val = -val mod p.
+//
+//go:noescape
+func p256NegCond(val *p256Element, cond int)
+
+// If cond is 0, sets res = b, otherwise sets res = a.
+//
+//go:noescape
+func p256MovCond(res, a, b *P256Point, cond int)
+
+// p256Table is a table of the first 16 multiples of a point. Points are stored
+// at an index offset of -1 so [8]P is at index 7, P is at 0, and [16]P is at 15.
+// [0]P is the point at infinity and it's not stored.
+type p256Table [16]P256Point
+
+// p256Select sets res to the point at index idx in the table.
+// idx must be in [0, 15]. It executes in constant time.
+//
+//go:noescape
+func p256Select(res *P256Point, table *p256Table, idx int)
+
+// p256AffinePoint is a point in affine coordinates (x, y). x and y are still
+// Montgomery domain elements. The point can't be the point at infinity.
+type p256AffinePoint struct {
+	x, y p256Element
+}
+
+// p256AffineTable is a table of the first 32 multiples of a point. Points are
+// stored at an index offset of -1 like in p256Table, and [0]P is not stored.
+type p256AffineTable [32]p256AffinePoint
+
+// p256Precomputed is a series of precomputed multiples of G, the canonical
+// generator. The first p256AffineTable contains multiples of G. The second one
+// multiples of [2⁶]G, the third one of [2¹²]G, and so on, where each successive
+// table is the previous table doubled six times. Six is the width of the
+// sliding window used in p256ScalarBaseMult, and having each table already
+// pre-doubled lets us avoid the doublings between windows entirely. This table
+// aliases into p256PrecomputedEmbed.
+var p256Precomputed *[43]p256AffineTable
+
+func init() {
+	p256PrecomputedPtr := unsafe.Pointer(&p256PrecomputedEmbed)
+	if runtime.GOARCH == "s390x" {
+		var newTable [43 * 32 * 2 * 4]uint64
+		for i, x := range (*[43 * 32 * 2 * 4][8]byte)(p256PrecomputedPtr) {
+			newTable[i] = byteorder.LEUint64(x[:])
+		}
+		p256PrecomputedPtr = unsafe.Pointer(&newTable)
+	}
+	p256Precomputed = (*[43]p256AffineTable)(p256PrecomputedPtr)
+}
+
+// p256SelectAffine sets res to the point at index idx in the table.
+// idx must be in [0, 31]. It executes in constant time.
+//
+//go:noescape
+func p256SelectAffine(res *p256AffinePoint, table *p256AffineTable, idx int)
+
+// Point addition with an affine point and constant time conditions.
+// If zero is 0, sets res = in2. If sel is 0, sets res = in1.
+// If sign is not 0, sets res = in1 + -in2. Otherwise, sets res = in1 + in2
+//
+//go:noescape
+func p256PointAddAffineAsm(res, in1 *P256Point, in2 *p256AffinePoint, sign, sel, zero int)
+
+// Point addition. Sets res = in1 + in2. Returns one if the two input points
+// were equal and zero otherwise. If in1 or in2 are the point at infinity, res
+// and the return value are undefined.
+//
+//go:noescape
+func p256PointAddAsm(res, in1, in2 *P256Point) int
+
+// Point doubling. Sets res = in + in. in can be the point at infinity.
+//
+//go:noescape
+func p256PointDoubleAsm(res, in *P256Point)
+
+// p256OrdElement is a P-256 scalar field element in [0, ord(G)-1] in the
+// Montgomery domain (with R 2²⁵⁶) as four uint64 limbs in little-endian order.
+type p256OrdElement [4]uint64
+
+// p256OrdReduce ensures s is in the range [0, ord(G)-1].
+func p256OrdReduce(s *p256OrdElement) {
+	// Since 2 * ord(G) > 2²⁵⁶, we can just conditionally subtract ord(G),
+	// keeping the result if it doesn't underflow.
+	t0, b := bits.Sub64(s[0], 0xf3b9cac2fc632551, 0)
+	t1, b := bits.Sub64(s[1], 0xbce6faada7179e84, b)
+	t2, b := bits.Sub64(s[2], 0xffffffffffffffff, b)
+	t3, b := bits.Sub64(s[3], 0xffffffff00000000, b)
+	tMask := b - 1 // zero if subtraction underflowed
+	s[0] ^= (t0 ^ s[0]) & tMask
+	s[1] ^= (t1 ^ s[1]) & tMask
+	s[2] ^= (t2 ^ s[2]) & tMask
+	s[3] ^= (t3 ^ s[3]) & tMask
+}
+
+func p256OrdLittleToBig(b *[32]byte, l *p256OrdElement) {
+	limbsToBytes(b, (*[4]uint64)(l))
+}
+
+func p256OrdBigToLittle(l *p256OrdElement, b *[32]byte) {
+	bytesToLimbs((*[4]uint64)(l), b)
+}
+
+// Add sets q = p1 + p2, and returns q. The points may overlap.
+func (q *P256Point) Add(r1, r2 *P256Point) *P256Point {
+	var sum, double P256Point
+	r1IsInfinity := r1.isInfinity()
+	r2IsInfinity := r2.isInfinity()
+	pointsEqual := p256PointAddAsm(&sum, r1, r2)
+	p256PointDoubleAsm(&double, r1)
+	p256MovCond(&sum, &double, &sum, pointsEqual)
+	p256MovCond(&sum, r1, &sum, r2IsInfinity)
+	p256MovCond(&sum, r2, &sum, r1IsInfinity)
+	return q.Set(&sum)
+}
+
+// Double sets q = p + p, and returns q. The points may overlap.
+func (q *P256Point) Double(p *P256Point) *P256Point {
+	var double P256Point
+	p256PointDoubleAsm(&double, p)
+	return q.Set(&double)
+}
+
+// ScalarBaseMult sets r = scalar * generator, where scalar is a 32-byte big
+// endian value, and returns r. If scalar is not 32 bytes long, ScalarBaseMult
+// returns an error and the receiver is unchanged.
+func (r *P256Point) ScalarBaseMult(scalar []byte) (*P256Point, error) {
+	if len(scalar) != 32 {
+		return nil, errors.New("invalid scalar length")
+	}
+	scalarReversed := new(p256OrdElement)
+	p256OrdBigToLittle(scalarReversed, (*[32]byte)(scalar))
+	p256OrdReduce(scalarReversed)
+
+	r.p256BaseMult(scalarReversed)
+	return r, nil
+}
+
+// ScalarMult sets r = scalar * q, where scalar is a 32-byte big endian value,
+// and returns r. If scalar is not 32 bytes long, ScalarBaseMult returns an
+// error and the receiver is unchanged.
+func (r *P256Point) ScalarMult(q *P256Point, scalar []byte) (*P256Point, error) {
+	if len(scalar) != 32 {
+		return nil, errors.New("invalid scalar length")
+	}
+	scalarReversed := new(p256OrdElement)
+	p256OrdBigToLittle(scalarReversed, (*[32]byte)(scalar))
+	p256OrdReduce(scalarReversed)
+
+	r.Set(q).p256ScalarMult(scalarReversed)
+	return r, nil
+}
+
+// uint64IsZero returns 1 if x is zero and zero otherwise.
+func uint64IsZero(x uint64) int {
+	x = ^x
+	x &= x >> 32
+	x &= x >> 16
+	x &= x >> 8
+	x &= x >> 4
+	x &= x >> 2
+	x &= x >> 1
+	return int(x & 1)
+}
+
+// p256Equal returns 1 if a and b are equal and 0 otherwise.
+func p256Equal(a, b *p256Element) int {
+	var acc uint64
+	for i := range a {
+		acc |= a[i] ^ b[i]
+	}
+	return uint64IsZero(acc)
+}
+
+// isInfinity returns 1 if p is the point at infinity and 0 otherwise.
+func (p *P256Point) isInfinity() int {
+	return p256Equal(&p.z, &p256Zero)
+}
+
+// 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 *P256Point) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p256UncompressedLength]byte
+	return p.bytes(&out)
+}
+
+func (p *P256Point) bytes(out *[p256UncompressedLength]byte) []byte {
+	// The proper representation of the point at infinity is a single zero byte.
+	if p.isInfinity() == 1 {
+		return append(out[:0], 0)
+	}
+
+	x, y := new(p256Element), new(p256Element)
+	p.affineFromMont(x, y)
+
+	out[0] = 4 // Uncompressed form.
+	p256LittleToBig((*[32]byte)(out[1:33]), x)
+	p256LittleToBig((*[32]byte)(out[33:65]), y)
+
+	return out[:]
+}
+
+// affineFromMont sets (x, y) to the affine coordinates of p, converted out of the
+// Montgomery domain.
+func (p *P256Point) affineFromMont(x, y *p256Element) {
+	p256Inverse(y, &p.z)
+	p256Sqr(x, y, 1)
+	p256Mul(y, y, x)
+
+	p256Mul(x, &p.x, x)
+	p256Mul(y, &p.y, y)
+
+	p256FromMont(x, x)
+	p256FromMont(y, y)
+}
+
+// BytesX returns the encoding of the x-coordinate of p, as specified in SEC 1,
+// Version 2.0, Section 2.3.5, or an error if p is the point at infinity.
+func (p *P256Point) BytesX() ([]byte, error) {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p256ElementLength]byte
+	return p.bytesX(&out)
+}
+
+func (p *P256Point) bytesX(out *[p256ElementLength]byte) ([]byte, error) {
+	if p.isInfinity() == 1 {
+		return nil, errors.New("P256 point is the point at infinity")
+	}
+
+	x := new(p256Element)
+	p256Inverse(x, &p.z)
+	p256Sqr(x, x, 1)
+	p256Mul(x, &p.x, x)
+	p256FromMont(x, x)
+	p256LittleToBig((*[32]byte)(out[:]), x)
+
+	return out[:], nil
+}
+
+// BytesCompressed returns the compressed 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 *P256Point) BytesCompressed() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p256CompressedLength]byte
+	return p.bytesCompressed(&out)
+}
+
+func (p *P256Point) bytesCompressed(out *[p256CompressedLength]byte) []byte {
+	if p.isInfinity() == 1 {
+		return append(out[:0], 0)
+	}
+
+	x, y := new(p256Element), new(p256Element)
+	p.affineFromMont(x, y)
+
+	out[0] = 2 | byte(y[0]&1)
+	p256LittleToBig((*[32]byte)(out[1:33]), x)
+
+	return out[:]
+}
+
+// Select sets q to p1 if cond == 1, and to p2 if cond == 0.
+func (q *P256Point) Select(p1, p2 *P256Point, cond int) *P256Point {
+	p256MovCond(q, p1, p2, cond)
+	return q
+}
+
+// p256Inverse sets out to in⁻¹ mod p. If in is zero, out will be zero.
+func p256Inverse(out, in *p256Element) {
+	// Inversion is calculated through exponentiation by p - 2, per Fermat's
+	// little theorem.
+	//
+	// The sequence of 12 multiplications and 255 squarings is derived from the
+	// following addition chain generated with github.com/mmcloughlin/addchain
+	// v0.4.0.
+	//
+	//  _10     = 2*1
+	//  _11     = 1 + _10
+	//  _110    = 2*_11
+	//  _111    = 1 + _110
+	//  _111000 = _111 << 3
+	//  _111111 = _111 + _111000
+	//  x12     = _111111 << 6 + _111111
+	//  x15     = x12 << 3 + _111
+	//  x16     = 2*x15 + 1
+	//  x32     = x16 << 16 + x16
+	//  i53     = x32 << 15
+	//  x47     = x15 + i53
+	//  i263    = ((i53 << 17 + 1) << 143 + x47) << 47
+	//  return    (x47 + i263) << 2 + 1
+	//
+	var z = new(p256Element)
+	var t0 = new(p256Element)
+	var t1 = new(p256Element)
+
+	p256Sqr(z, in, 1)
+	p256Mul(z, in, z)
+	p256Sqr(z, z, 1)
+	p256Mul(z, in, z)
+	p256Sqr(t0, z, 3)
+	p256Mul(t0, z, t0)
+	p256Sqr(t1, t0, 6)
+	p256Mul(t0, t0, t1)
+	p256Sqr(t0, t0, 3)
+	p256Mul(z, z, t0)
+	p256Sqr(t0, z, 1)
+	p256Mul(t0, in, t0)
+	p256Sqr(t1, t0, 16)
+	p256Mul(t0, t0, t1)
+	p256Sqr(t0, t0, 15)
+	p256Mul(z, z, t0)
+	p256Sqr(t0, t0, 17)
+	p256Mul(t0, in, t0)
+	p256Sqr(t0, t0, 143)
+	p256Mul(t0, z, t0)
+	p256Sqr(t0, t0, 47)
+	p256Mul(z, z, t0)
+	p256Sqr(z, z, 2)
+	p256Mul(out, in, z)
+}
+
+func boothW5(in uint) (int, int) {
+	var s uint = ^((in >> 5) - 1)
+	var d uint = (1 << 6) - in - 1
+	d = (d & s) | (in & (^s))
+	d = (d >> 1) + (d & 1)
+	return int(d), int(s & 1)
+}
+
+func boothW6(in uint) (int, int) {
+	var s uint = ^((in >> 6) - 1)
+	var d uint = (1 << 7) - in - 1
+	d = (d & s) | (in & (^s))
+	d = (d >> 1) + (d & 1)
+	return int(d), int(s & 1)
+}
+
+func (p *P256Point) p256BaseMult(scalar *p256OrdElement) {
+	var t0 p256AffinePoint
+
+	wvalue := (scalar[0] << 1) & 0x7f
+	sel, sign := boothW6(uint(wvalue))
+	p256SelectAffine(&t0, &p256Precomputed[0], sel)
+	p.x, p.y, p.z = t0.x, t0.y, p256One
+	p256NegCond(&p.y, sign)
+
+	index := uint(5)
+	zero := sel
+
+	for i := 1; i < 43; i++ {
+		if index < 192 {
+			wvalue = ((scalar[index/64] >> (index % 64)) + (scalar[index/64+1] << (64 - (index % 64)))) & 0x7f
+		} else {
+			wvalue = (scalar[index/64] >> (index % 64)) & 0x7f
+		}
+		index += 6
+		sel, sign = boothW6(uint(wvalue))
+		p256SelectAffine(&t0, &p256Precomputed[i], sel)
+		p256PointAddAffineAsm(p, p, &t0, sign, sel, zero)
+		zero |= sel
+	}
+
+	// If the whole scalar was zero, set to the point at infinity.
+	p256MovCond(p, p, NewP256Point(), zero)
+}
+
+func (p *P256Point) p256ScalarMult(scalar *p256OrdElement) {
+	// precomp is a table of precomputed points that stores powers of p
+	// from p^1 to p^16.
+	var precomp p256Table
+	var t0, t1, t2, t3 P256Point
+
+	// Prepare the table
+	precomp[0] = *p // 1
+
+	p256PointDoubleAsm(&t0, p)
+	p256PointDoubleAsm(&t1, &t0)
+	p256PointDoubleAsm(&t2, &t1)
+	p256PointDoubleAsm(&t3, &t2)
+	precomp[1] = t0  // 2
+	precomp[3] = t1  // 4
+	precomp[7] = t2  // 8
+	precomp[15] = t3 // 16
+
+	p256PointAddAsm(&t0, &t0, p)
+	p256PointAddAsm(&t1, &t1, p)
+	p256PointAddAsm(&t2, &t2, p)
+	precomp[2] = t0 // 3
+	precomp[4] = t1 // 5
+	precomp[8] = t2 // 9
+
+	p256PointDoubleAsm(&t0, &t0)
+	p256PointDoubleAsm(&t1, &t1)
+	precomp[5] = t0 // 6
+	precomp[9] = t1 // 10
+
+	p256PointAddAsm(&t2, &t0, p)
+	p256PointAddAsm(&t1, &t1, p)
+	precomp[6] = t2  // 7
+	precomp[10] = t1 // 11
+
+	p256PointDoubleAsm(&t0, &t0)
+	p256PointDoubleAsm(&t2, &t2)
+	precomp[11] = t0 // 12
+	precomp[13] = t2 // 14
+
+	p256PointAddAsm(&t0, &t0, p)
+	p256PointAddAsm(&t2, &t2, p)
+	precomp[12] = t0 // 13
+	precomp[14] = t2 // 15
+
+	// Start scanning the window from top bit
+	index := uint(254)
+	var sel, sign int
+
+	wvalue := (scalar[index/64] >> (index % 64)) & 0x3f
+	sel, _ = boothW5(uint(wvalue))
+
+	p256Select(p, &precomp, sel)
+	zero := sel
+
+	for index > 4 {
+		index -= 5
+		p256PointDoubleAsm(p, p)
+		p256PointDoubleAsm(p, p)
+		p256PointDoubleAsm(p, p)
+		p256PointDoubleAsm(p, p)
+		p256PointDoubleAsm(p, p)
+
+		if index < 192 {
+			wvalue = ((scalar[index/64] >> (index % 64)) + (scalar[index/64+1] << (64 - (index % 64)))) & 0x3f
+		} else {
+			wvalue = (scalar[index/64] >> (index % 64)) & 0x3f
+		}
+
+		sel, sign = boothW5(uint(wvalue))
+
+		p256Select(&t0, &precomp, sel)
+		p256NegCond(&t0.y, sign)
+		p256PointAddAsm(&t1, p, &t0)
+		p256MovCond(&t1, &t1, p, sel)
+		p256MovCond(p, &t1, &t0, zero)
+		zero |= sel
+	}
+
+	p256PointDoubleAsm(p, p)
+	p256PointDoubleAsm(p, p)
+	p256PointDoubleAsm(p, p)
+	p256PointDoubleAsm(p, p)
+	p256PointDoubleAsm(p, p)
+
+	wvalue = (scalar[0] << 1) & 0x3f
+	sel, sign = boothW5(uint(wvalue))
+
+	p256Select(&t0, &precomp, sel)
+	p256NegCond(&t0.y, sign)
+	p256PointAddAsm(&t1, p, &t0)
+	p256MovCond(&t1, &t1, p, sel)
+	p256MovCond(p, &t1, &t0, zero)
+}

p256_asm_test.go

@@ -0,0 +1,53 @@
+// Copyright 2024 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.
+
+//go:build (amd64 || arm64 || ppc64le || s390x) && !purego && linux
+
+package nistec
+
+import (
+	"syscall"
+	"testing"
+	"unsafe"
+)
+
+// Lightly adapted from the bytes test package. Allocate a pair of T one at the start of a page, another at the
+// end. Any access beyond or before the page boundary should cause a fault. This is linux specific.
+func dangerousObjs[T any](t *testing.T) (start *T, end *T) {
+	pagesize := syscall.Getpagesize()
+	b, err := syscall.Mmap(0, 0, 3*pagesize, syscall.PROT_READ|syscall.PROT_WRITE, syscall.MAP_ANONYMOUS|syscall.MAP_PRIVATE)
+	if err != nil {
+		t.Fatalf("mmap failed %s", err)
+	}
+	err = syscall.Mprotect(b[:pagesize], syscall.PROT_NONE)
+	if err != nil {
+		t.Fatalf("mprotect low failed %s\n", err)
+	}
+	err = syscall.Mprotect(b[2*pagesize:], syscall.PROT_NONE)
+	if err != nil {
+		t.Fatalf("mprotect high failed %s\n", err)
+	}
+	b = b[pagesize : 2*pagesize]
+	end = (*T)(unsafe.Pointer(&b[len(b)-(int)(unsafe.Sizeof(*end))]))
+	start = (*T)(unsafe.Pointer(&b[0]))
+	return start, end
+}
+
+func TestP256SelectAffinePageBoundary(t *testing.T) {
+	var out p256AffinePoint
+	begintp, endtp := dangerousObjs[p256AffineTable](t)
+	for i := 0; i < 31; i++ {
+		p256SelectAffine(&out, begintp, i)
+		p256SelectAffine(&out, endtp, i)
+	}
+}
+
+func TestP256SelectPageBoundary(t *testing.T) {
+	var out P256Point
+	begintp, endtp := dangerousObjs[p256Table](t)
+	for i := 0; i < 15; i++ {
+		p256Select(&out, begintp, i)
+		p256Select(&out, endtp, i)
+	}
+}

p256_ordinv.go

@@ -0,0 +1,102 @@
+// Copyright 2022 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.
+
+//go:build (amd64 || arm64) && !purego
+
+package nistec
+
+import "errors"
+
+// Montgomery multiplication modulo org(G). Sets res = in1 * in2 * R⁻¹.
+//
+//go:noescape
+func p256OrdMul(res, in1, in2 *p256OrdElement)
+
+// Montgomery square modulo org(G), repeated n times (n >= 1).
+//
+//go:noescape
+func p256OrdSqr(res, in *p256OrdElement, n int)
+
+func p256OrdInverse(k []byte) ([]byte, error) {
+	if len(k) != 32 {
+		return nil, errors.New("invalid scalar length")
+	}
+
+	x := new(p256OrdElement)
+	p256OrdBigToLittle(x, (*[32]byte)(k))
+	p256OrdReduce(x)
+
+	// Inversion is implemented as exponentiation by n - 2, per Fermat's little theorem.
+	//
+	// The sequence of 38 multiplications and 254 squarings is derived from
+	// https://briansmith.org/ecc-inversion-addition-chains-01#p256_scalar_inversion
+	_1 := new(p256OrdElement)
+	_11 := new(p256OrdElement)
+	_101 := new(p256OrdElement)
+	_111 := new(p256OrdElement)
+	_1111 := new(p256OrdElement)
+	_10101 := new(p256OrdElement)
+	_101111 := new(p256OrdElement)
+	t := new(p256OrdElement)
+
+	// This code operates in the Montgomery domain where R = 2²⁵⁶ mod n and n is
+	// the order of the scalar field. Elements in the Montgomery domain take the
+	// form a×R and p256OrdMul calculates (a × b × R⁻¹) mod n. RR is R in the
+	// domain, or R×R mod n, thus p256OrdMul(x, RR) gives x×R, i.e. converts x
+	// into the Montgomery domain.
+	RR := &p256OrdElement{0x83244c95be79eea2, 0x4699799c49bd6fa6,
+		0x2845b2392b6bec59, 0x66e12d94f3d95620}
+
+	p256OrdMul(_1, x, RR)      // _1
+	p256OrdSqr(x, _1, 1)       // _10
+	p256OrdMul(_11, x, _1)     // _11
+	p256OrdMul(_101, x, _11)   // _101
+	p256OrdMul(_111, x, _101)  // _111
+	p256OrdSqr(x, _101, 1)     // _1010
+	p256OrdMul(_1111, _101, x) // _1111
+
+	p256OrdSqr(t, x, 1)          // _10100
+	p256OrdMul(_10101, t, _1)    // _10101
+	p256OrdSqr(x, _10101, 1)     // _101010
+	p256OrdMul(_101111, _101, x) // _101111
+	p256OrdMul(x, _10101, x)     // _111111 = x6
+	p256OrdSqr(t, x, 2)          // _11111100
+	p256OrdMul(t, t, _11)        // _11111111 = x8
+	p256OrdSqr(x, t, 8)          // _ff00
+	p256OrdMul(x, x, t)          // _ffff = x16
+	p256OrdSqr(t, x, 16)         // _ffff0000
+	p256OrdMul(t, t, x)          // _ffffffff = x32
+
+	p256OrdSqr(x, t, 64)
+	p256OrdMul(x, x, t)
+	p256OrdSqr(x, x, 32)
+	p256OrdMul(x, x, t)
+
+	sqrs := []int{
+		6, 5, 4, 5, 5,
+		4, 3, 3, 5, 9,
+		6, 2, 5, 6, 5,
+		4, 5, 5, 3, 10,
+		2, 5, 5, 3, 7, 6}
+	muls := []*p256OrdElement{
+		_101111, _111, _11, _1111, _10101,
+		_101, _101, _101, _111, _101111,
+		_1111, _1, _1, _1111, _111,
+		_111, _111, _101, _11, _101111,
+		_11, _11, _11, _1, _10101, _1111}
+
+	for i, s := range sqrs {
+		p256OrdSqr(x, x, s)
+		p256OrdMul(x, x, muls[i])
+	}
+
+	// Montgomery multiplication by R⁻¹, or 1 outside the domain as R⁻¹×R = 1,
+	// converts a Montgomery value out of the domain.
+	one := &p256OrdElement{1}
+	p256OrdMul(x, x, one)
+
+	var xOut [32]byte
+	p256OrdLittleToBig(&xOut, x)
+	return xOut[:], nil
+}

p256_ordinv_export_test.go

@@ -0,0 +1,13 @@
+// Copyright 2022 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
+
+// This file exports the P256OrdInverse function so it's accessible during tests
+// from the unmodified p256_ordinv_test.go from the stdlib, but not as part
+// of the public API of filippo.io/nistec.
+
+func P256OrdInverse(k []byte) ([]byte, error) {
+	return p256OrdInverse(k)
+}

p256_ordinv_noasm.go

@@ -0,0 +1,13 @@
+// Copyright 2022 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.
+
+//go:build (!amd64 && !arm64) || purego
+
+package nistec
+
+import "errors"
+
+func p256OrdInverse(k []byte) ([]byte, error) {
+	return nil, errors.New("unimplemented")
+}

p256_ordinv_test.go

@@ -0,0 +1,95 @@
+// Copyright 2022 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.
+
+//go:build (amd64 || arm64) && !purego
+
+package nistec_test
+
+import (
+	"bytes"
+	"crypto/elliptic"
+	"math/big"
+	"testing"
+
+	"sources.truenas.cloud/code/nistec"
+)
+
+func TestP256OrdInverse(t *testing.T) {
+	N := elliptic.P256().Params().N
+
+	// inv(0) is expected to be 0.
+	zero := make([]byte, 32)
+	out, err := nistec.P256OrdInverse(zero)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !bytes.Equal(out, zero) {
+		t.Error("unexpected output for inv(0)")
+	}
+
+	// inv(N) is also 0 mod N.
+	input := make([]byte, 32)
+	N.FillBytes(input)
+	out, err = nistec.P256OrdInverse(input)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !bytes.Equal(out, zero) {
+		t.Error("unexpected output for inv(N)")
+	}
+	if !bytes.Equal(input, N.Bytes()) {
+		t.Error("input was modified")
+	}
+
+	// Check inv(1) and inv(N+1) against math/big
+	exp := new(big.Int).ModInverse(big.NewInt(1), N).FillBytes(make([]byte, 32))
+	big.NewInt(1).FillBytes(input)
+	out, err = nistec.P256OrdInverse(input)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !bytes.Equal(out, exp) {
+		t.Error("unexpected output for inv(1)")
+	}
+	new(big.Int).Add(N, big.NewInt(1)).FillBytes(input)
+	out, err = nistec.P256OrdInverse(input)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !bytes.Equal(out, exp) {
+		t.Error("unexpected output for inv(N+1)")
+	}
+
+	// Check inv(20) and inv(N+20) against math/big
+	exp = new(big.Int).ModInverse(big.NewInt(20), N).FillBytes(make([]byte, 32))
+	big.NewInt(20).FillBytes(input)
+	out, err = nistec.P256OrdInverse(input)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !bytes.Equal(out, exp) {
+		t.Error("unexpected output for inv(20)")
+	}
+	new(big.Int).Add(N, big.NewInt(20)).FillBytes(input)
+	out, err = nistec.P256OrdInverse(input)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !bytes.Equal(out, exp) {
+		t.Error("unexpected output for inv(N+20)")
+	}
+
+	// Check inv(2^256-1) against math/big
+	bigInput := new(big.Int).Lsh(big.NewInt(1), 256)
+	bigInput.Sub(bigInput, big.NewInt(1))
+	exp = new(big.Int).ModInverse(bigInput, N).FillBytes(make([]byte, 32))
+	bigInput.FillBytes(input)
+	out, err = nistec.P256OrdInverse(input)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !bytes.Equal(out, exp) {
+		t.Error("unexpected output for inv(2^256-1)")
+	}
+}

p256_test.go

@@ -0,0 +1,49 @@
+// 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.
+
+//go:build (!amd64 && !arm64 && !ppc64le && !s390x) || purego
+
+package nistec
+
+import (
+	"bytes"
+	"fmt"
+	"testing"
+
+	"sources.truenas.cloud/code/nistec/internal/fiat"
+)
+
+func TestP256PrecomputedTable(t *testing.T) {
+	base := NewP256Point().SetGenerator()
+
+	for i := 0; i < 43; i++ {
+		t.Run(fmt.Sprintf("table[%d]", i), func(t *testing.T) {
+			testP256AffineTable(t, base, &p256GeneratorTables[i])
+		})
+
+		for k := 0; k < 6; k++ {
+			base.Double(base)
+		}
+	}
+}
+
+func testP256AffineTable(t *testing.T, base *P256Point, table *p256AffineTable) {
+	p := NewP256Point()
+	zInv := new(fiat.P256Element)
+
+	for j := 0; j < 32; j++ {
+		p.Add(p, base)
+
+		// Convert p to affine coordinates.
+		zInv.Invert(&p.z)
+		p.x.Mul(&p.x, zInv)
+		p.y.Mul(&p.y, zInv)
+		p.z.One()
+
+		if !bytes.Equal(table[j].x.Bytes(), p.x.Bytes()) ||
+			!bytes.Equal(table[j].y.Bytes(), p.y.Bytes()) {
+			t.Fatalf("incorrect table entry at index %d", j)
+		}
+	}
+}

p384.go

@@ -0,0 +1,541 @@
+// Copyright 2022 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.
+
+// Code generated by generate.go. DO NOT EDIT.
+
+package nistec
+
+import (
+	"crypto/subtle"
+	"errors"
+	"sync"
+
+	"sources.truenas.cloud/code/nistec/internal/fiat"
+)
+
+// p384ElementLength is the length of an element of the base or scalar field,
+// which have the same bytes length for all NIST P curves.
+const p384ElementLength = 48
+
+// P384Point is a P384 point. The zero value is NOT valid.
+type P384Point struct {
+	// The point is represented in projective coordinates (X:Y:Z),
+	// where x = X/Z and y = Y/Z.
+	x, y, z *fiat.P384Element
+}
+
+// NewP384Point returns a new P384Point representing the point at infinity point.
+func NewP384Point() *P384Point {
+	return &P384Point{
+		x: new(fiat.P384Element),
+		y: new(fiat.P384Element).One(),
+		z: new(fiat.P384Element),
+	}
+}
+
+// SetGenerator sets p to the canonical generator and returns p.
+func (p *P384Point) SetGenerator() *P384Point {
+	p.x.SetBytes([]byte{0xaa, 0x87, 0xca, 0x22, 0xbe, 0x8b, 0x5, 0x37, 0x8e, 0xb1, 0xc7, 0x1e, 0xf3, 0x20, 0xad, 0x74, 0x6e, 0x1d, 0x3b, 0x62, 0x8b, 0xa7, 0x9b, 0x98, 0x59, 0xf7, 0x41, 0xe0, 0x82, 0x54, 0x2a, 0x38, 0x55, 0x2, 0xf2, 0x5d, 0xbf, 0x55, 0x29, 0x6c, 0x3a, 0x54, 0x5e, 0x38, 0x72, 0x76, 0xa, 0xb7})
+	p.y.SetBytes([]byte{0x36, 0x17, 0xde, 0x4a, 0x96, 0x26, 0x2c, 0x6f, 0x5d, 0x9e, 0x98, 0xbf, 0x92, 0x92, 0xdc, 0x29, 0xf8, 0xf4, 0x1d, 0xbd, 0x28, 0x9a, 0x14, 0x7c, 0xe9, 0xda, 0x31, 0x13, 0xb5, 0xf0, 0xb8, 0xc0, 0xa, 0x60, 0xb1, 0xce, 0x1d, 0x7e, 0x81, 0x9d, 0x7a, 0x43, 0x1d, 0x7c, 0x90, 0xea, 0xe, 0x5f})
+	p.z.One()
+	return p
+}
+
+// Set sets p = q and returns p.
+func (p *P384Point) Set(q *P384Point) *P384Point {
+	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 *P384Point) SetBytes(b []byte) (*P384Point, error) {
+	switch {
+	// Point at infinity.
+	case len(b) == 1 && b[0] == 0:
+		return p.Set(NewP384Point()), nil
+
+	// Uncompressed form.
+	case len(b) == 1+2*p384ElementLength && b[0] == 4:
+		x, err := new(fiat.P384Element).SetBytes(b[1 : 1+p384ElementLength])
+		if err != nil {
+			return nil, err
+		}
+		y, err := new(fiat.P384Element).SetBytes(b[1+p384ElementLength:])
+		if err != nil {
+			return nil, err
+		}
+		if err := p384CheckOnCurve(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+p384ElementLength && (b[0] == 2 || b[0] == 3):
+		x, err := new(fiat.P384Element).SetBytes(b[1:])
+		if err != nil {
+			return nil, err
+		}
+
+		// y² = x³ - 3x + b
+		y := p384Polynomial(new(fiat.P384Element), x)
+		if !p384Sqrt(y, y) {
+			return nil, errors.New("invalid P384 compressed point encoding")
+		}
+
+		// Select the positive or negative root, as indicated by the least
+		// significant bit, based on the encoding type byte.
+		otherRoot := new(fiat.P384Element)
+		otherRoot.Sub(otherRoot, y)
+		cond := y.Bytes()[p384ElementLength-1]&1 ^ b[0]&1
+		y.Select(otherRoot, y, int(cond))
+
+		p.x.Set(x)
+		p.y.Set(y)
+		p.z.One()
+		return p, nil
+
+	default:
+		return nil, errors.New("invalid P384 point encoding")
+	}
+}
+
+var _p384B *fiat.P384Element
+var _p384BOnce sync.Once
+
+func p384B() *fiat.P384Element {
+	_p384BOnce.Do(func() {
+		_p384B, _ = new(fiat.P384Element).SetBytes([]byte{0xb3, 0x31, 0x2f, 0xa7, 0xe2, 0x3e, 0xe7, 0xe4, 0x98, 0x8e, 0x5, 0x6b, 0xe3, 0xf8, 0x2d, 0x19, 0x18, 0x1d, 0x9c, 0x6e, 0xfe, 0x81, 0x41, 0x12, 0x3, 0x14, 0x8, 0x8f, 0x50, 0x13, 0x87, 0x5a, 0xc6, 0x56, 0x39, 0x8d, 0x8a, 0x2e, 0xd1, 0x9d, 0x2a, 0x85, 0xc8, 0xed, 0xd3, 0xec, 0x2a, 0xef})
+	})
+	return _p384B
+}
+
+// p384Polynomial sets y2 to x³ - 3x + b, and returns y2.
+func p384Polynomial(y2, x *fiat.P384Element) *fiat.P384Element {
+	y2.Square(x)
+	y2.Mul(y2, x)
+
+	threeX := new(fiat.P384Element).Add(x, x)
+	threeX.Add(threeX, x)
+	y2.Sub(y2, threeX)
+
+	return y2.Add(y2, p384B())
+}
+
+func p384CheckOnCurve(x, y *fiat.P384Element) error {
+	// y² = x³ - 3x + b
+	rhs := p384Polynomial(new(fiat.P384Element), x)
+	lhs := new(fiat.P384Element).Square(y)
+	if rhs.Equal(lhs) != 1 {
+		return errors.New("P384 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 *P384Point) Bytes() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [1 + 2*p384ElementLength]byte
+	return p.bytes(&out)
+}
+
+func (p *P384Point) bytes(out *[1 + 2*p384ElementLength]byte) []byte {
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new(fiat.P384Element).Invert(p.z)
+	x := new(fiat.P384Element).Mul(p.x, zinv)
+	y := new(fiat.P384Element).Mul(p.y, zinv)
+
+	buf := append(out[:0], 4)
+	buf = append(buf, x.Bytes()...)
+	buf = append(buf, y.Bytes()...)
+	return buf
+}
+
+// BytesX returns the encoding of the x-coordinate of p, as specified in SEC 1,
+// Version 2.0, Section 2.3.5, or an error if p is the point at infinity.
+func (p *P384Point) BytesX() ([]byte, error) {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p384ElementLength]byte
+	return p.bytesX(&out)
+}
+
+func (p *P384Point) bytesX(out *[p384ElementLength]byte) ([]byte, error) {
+	if p.z.IsZero() == 1 {
+		return nil, errors.New("P384 point is the point at infinity")
+	}
+
+	zinv := new(fiat.P384Element).Invert(p.z)
+	x := new(fiat.P384Element).Mul(p.x, zinv)
+
+	return append(out[:0], x.Bytes()...), nil
+}
+
+// BytesCompressed returns the compressed 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 *P384Point) BytesCompressed() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [1 + p384ElementLength]byte
+	return p.bytesCompressed(&out)
+}
+
+func (p *P384Point) bytesCompressed(out *[1 + p384ElementLength]byte) []byte {
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new(fiat.P384Element).Invert(p.z)
+	x := new(fiat.P384Element).Mul(p.x, zinv)
+	y := new(fiat.P384Element).Mul(p.y, zinv)
+
+	// Encode the sign of the y coordinate (indicated by the least significant
+	// bit) as the encoding type (2 or 3).
+	buf := append(out[:0], 2)
+	buf[0] |= y.Bytes()[p384ElementLength-1] & 1
+	buf = append(buf, x.Bytes()...)
+	return buf
+}
+
+// Add sets q = p1 + p2, and returns q. The points may overlap.
+func (q *P384Point) Add(p1, p2 *P384Point) *P384Point {
+	// 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.P384Element).Mul(p1.x, p2.x)  // t0 := X1 * X2
+	t1 := new(fiat.P384Element).Mul(p1.y, p2.y)  // t1 := Y1 * Y2
+	t2 := new(fiat.P384Element).Mul(p1.z, p2.z)  // t2 := Z1 * Z2
+	t3 := new(fiat.P384Element).Add(p1.x, p1.y)  // t3 := X1 + Y1
+	t4 := new(fiat.P384Element).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.P384Element).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.P384Element).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.P384Element).Mul(p384B(), 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(p384B(), 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 *P384Point) Double(p *P384Point) *P384Point {
+	// 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.P384Element).Square(p.x)      // t0 := X ^ 2
+	t1 := new(fiat.P384Element).Square(p.y)      // t1 := Y ^ 2
+	t2 := new(fiat.P384Element).Square(p.z)      // t2 := Z ^ 2
+	t3 := new(fiat.P384Element).Mul(p.x, p.y)    // t3 := X * Y
+	t3.Add(t3, t3)                               // t3 := t3 + t3
+	z3 := new(fiat.P384Element).Mul(p.x, p.z)    // Z3 := X * Z
+	z3.Add(z3, z3)                               // Z3 := Z3 + Z3
+	y3 := new(fiat.P384Element).Mul(p384B(), t2) // Y3 := b * t2
+	y3.Sub(y3, z3)                               // Y3 := Y3 - Z3
+	x3 := new(fiat.P384Element).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(p384B(), 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 *P384Point) Select(p1, p2 *P384Point, cond int) *P384Point {
+	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
+}
+
+// A p384Table holds the first 15 multiples of a point at offset -1, so [1]P
+// is at table[0], [15]P is at table[14], and [0]P is implicitly the identity
+// point.
+type p384Table [15]*P384Point
+
+// Select selects the n-th multiple of the table base point into p. It works in
+// constant time by iterating over every entry of the table. n must be in [0, 15].
+func (table *p384Table) Select(p *P384Point, n uint8) {
+	if n >= 16 {
+		panic("nistec: internal error: p384Table called with out-of-bounds value")
+	}
+	p.Set(NewP384Point())
+	for i := uint8(1); i < 16; i++ {
+		cond := subtle.ConstantTimeByteEq(i, n)
+		p.Select(table[i-1], p, cond)
+	}
+}
+
+// ScalarMult sets p = scalar * q, and returns p.
+func (p *P384Point) ScalarMult(q *P384Point, scalar []byte) (*P384Point, error) {
+	// Compute a p384Table for the base point q. The explicit NewP384Point
+	// calls get inlined, letting the allocations live on the stack.
+	var table = p384Table{NewP384Point(), NewP384Point(), NewP384Point(),
+		NewP384Point(), NewP384Point(), NewP384Point(), NewP384Point(),
+		NewP384Point(), NewP384Point(), NewP384Point(), NewP384Point(),
+		NewP384Point(), NewP384Point(), NewP384Point(), NewP384Point()}
+	table[0].Set(q)
+	for i := 1; i < 15; i += 2 {
+		table[i].Double(table[i/2])
+		table[i+1].Add(table[i], 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 := NewP384Point()
+	p.Set(NewP384Point())
+	for i, byte := range scalar {
+		// No need to double on the first iteration, as p is the identity at
+		// this point, and [N]∞ = ∞.
+		if i != 0 {
+			p.Double(p)
+			p.Double(p)
+			p.Double(p)
+			p.Double(p)
+		}
+
+		windowValue := byte >> 4
+		table.Select(t, windowValue)
+		p.Add(p, t)
+
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+
+		windowValue = byte & 0b1111
+		table.Select(t, windowValue)
+		p.Add(p, t)
+	}
+
+	return p, nil
+}
+
+var p384GeneratorTable *[p384ElementLength * 2]p384Table
+var p384GeneratorTableOnce sync.Once
+
+// generatorTable returns a sequence of p384Tables. The first table contains
+// multiples of G. Each successive table is the previous table doubled four
+// times.
+func (p *P384Point) generatorTable() *[p384ElementLength * 2]p384Table {
+	p384GeneratorTableOnce.Do(func() {
+		p384GeneratorTable = new([p384ElementLength * 2]p384Table)
+		base := NewP384Point().SetGenerator()
+		for i := 0; i < p384ElementLength*2; i++ {
+			p384GeneratorTable[i][0] = NewP384Point().Set(base)
+			for j := 1; j < 15; j++ {
+				p384GeneratorTable[i][j] = NewP384Point().Add(p384GeneratorTable[i][j-1], base)
+			}
+			base.Double(base)
+			base.Double(base)
+			base.Double(base)
+			base.Double(base)
+		}
+	})
+	return p384GeneratorTable
+}
+
+// ScalarBaseMult sets p = scalar * B, where B is the canonical generator, and
+// returns p.
+func (p *P384Point) ScalarBaseMult(scalar []byte) (*P384Point, error) {
+	if len(scalar) != p384ElementLength {
+		return nil, errors.New("invalid scalar length")
+	}
+	tables := p.generatorTable()
+
+	// This is also a scalar multiplication with a four-bit window like in
+	// ScalarMult, but in this case the doublings are precomputed. The value
+	// [windowValue]G added at iteration k would normally get doubled
+	// (totIterations-k)×4 times, but with a larger precomputation we can
+	// instead add [2^((totIterations-k)×4)][windowValue]G and avoid the
+	// doublings between iterations.
+	t := NewP384Point()
+	p.Set(NewP384Point())
+	tableIndex := len(tables) - 1
+	for _, byte := range scalar {
+		windowValue := byte >> 4
+		tables[tableIndex].Select(t, windowValue)
+		p.Add(p, t)
+		tableIndex--
+
+		windowValue = byte & 0b1111
+		tables[tableIndex].Select(t, windowValue)
+		p.Add(p, t)
+		tableIndex--
+	}
+
+	return p, nil
+}
+
+// p384Sqrt sets e to a square root of x. If x is not a square, p384Sqrt returns
+// false and e is unchanged. e and x can overlap.
+func p384Sqrt(e, x *fiat.P384Element) (isSquare bool) {
+	candidate := new(fiat.P384Element)
+	p384SqrtCandidate(candidate, x)
+	square := new(fiat.P384Element).Square(candidate)
+	if square.Equal(x) != 1 {
+		return false
+	}
+	e.Set(candidate)
+	return true
+}
+
+// p384SqrtCandidate sets z to a square root candidate for x. z and x must not overlap.
+func p384SqrtCandidate(z, x *fiat.P384Element) {
+	// Since p = 3 mod 4, exponentiation by (p + 1) / 4 yields a square root candidate.
+	//
+	// The sequence of 14 multiplications and 381 squarings is derived from the
+	// following addition chain generated with github.com/mmcloughlin/addchain v0.4.0.
+	//
+	//	_10      = 2*1
+	//	_11      = 1 + _10
+	//	_110     = 2*_11
+	//	_111     = 1 + _110
+	//	_111000  = _111 << 3
+	//	_111111  = _111 + _111000
+	//	_1111110 = 2*_111111
+	//	_1111111 = 1 + _1111110
+	//	x12      = _1111110 << 5 + _111111
+	//	x24      = x12 << 12 + x12
+	//	x31      = x24 << 7 + _1111111
+	//	x32      = 2*x31 + 1
+	//	x63      = x32 << 31 + x31
+	//	x126     = x63 << 63 + x63
+	//	x252     = x126 << 126 + x126
+	//	x255     = x252 << 3 + _111
+	//	return     ((x255 << 33 + x32) << 64 + 1) << 30
+	//
+	var t0 = new(fiat.P384Element)
+	var t1 = new(fiat.P384Element)
+	var t2 = new(fiat.P384Element)
+
+	z.Square(x)
+	z.Mul(x, z)
+	z.Square(z)
+	t0.Mul(x, z)
+	z.Square(t0)
+	for s := 1; s < 3; s++ {
+		z.Square(z)
+	}
+	t1.Mul(t0, z)
+	t2.Square(t1)
+	z.Mul(x, t2)
+	for s := 0; s < 5; s++ {
+		t2.Square(t2)
+	}
+	t1.Mul(t1, t2)
+	t2.Square(t1)
+	for s := 1; s < 12; s++ {
+		t2.Square(t2)
+	}
+	t1.Mul(t1, t2)
+	for s := 0; s < 7; s++ {
+		t1.Square(t1)
+	}
+	t1.Mul(z, t1)
+	z.Square(t1)
+	z.Mul(x, z)
+	t2.Square(z)
+	for s := 1; s < 31; s++ {
+		t2.Square(t2)
+	}
+	t1.Mul(t1, t2)
+	t2.Square(t1)
+	for s := 1; s < 63; s++ {
+		t2.Square(t2)
+	}
+	t1.Mul(t1, t2)
+	t2.Square(t1)
+	for s := 1; s < 126; s++ {
+		t2.Square(t2)
+	}
+	t1.Mul(t1, t2)
+	for s := 0; s < 3; s++ {
+		t1.Square(t1)
+	}
+	t0.Mul(t0, t1)
+	for s := 0; s < 33; s++ {
+		t0.Square(t0)
+	}
+	z.Mul(z, t0)
+	for s := 0; s < 64; s++ {
+		z.Square(z)
+	}
+	z.Mul(x, z)
+	for s := 0; s < 30; s++ {
+		z.Square(z)
+	}
+}

p521.go

@@ -0,0 +1,470 @@
+// Copyright 2022 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.
+
+// Code generated by generate.go. DO NOT EDIT.
+
+package nistec
+
+import (
+	"crypto/subtle"
+	"errors"
+	"sync"
+
+	"sources.truenas.cloud/code/nistec/internal/fiat"
+)
+
+// p521ElementLength is the length of an element of the base or scalar field,
+// which have the same bytes length for all NIST P curves.
+const p521ElementLength = 66
+
+// P521Point is a P521 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),
+	}
+}
+
+// SetGenerator sets p to the canonical generator and returns p.
+func (p *P521Point) SetGenerator() *P521Point {
+	p.x.SetBytes([]byte{0x0, 0xc6, 0x85, 0x8e, 0x6, 0xb7, 0x4, 0x4, 0xe9, 0xcd, 0x9e, 0x3e, 0xcb, 0x66, 0x23, 0x95, 0xb4, 0x42, 0x9c, 0x64, 0x81, 0x39, 0x5, 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})
+	p.y.SetBytes([]byte{0x1, 0x18, 0x39, 0x29, 0x6a, 0x78, 0x9a, 0x3b, 0xc0, 0x4, 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, 0x1, 0x3f, 0xad, 0x7, 0x61, 0x35, 0x3c, 0x70, 0x86, 0xa2, 0x72, 0xc2, 0x40, 0x88, 0xbe, 0x94, 0x76, 0x9f, 0xd1, 0x66, 0x50})
+	p.z.One()
+	return p
+}
+
+// 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] == 2 || b[0] == 3):
+		x, err := new(fiat.P521Element).SetBytes(b[1:])
+		if err != nil {
+			return nil, err
+		}
+
+		// y² = x³ - 3x + b
+		y := p521Polynomial(new(fiat.P521Element), x)
+		if !p521Sqrt(y, y) {
+			return nil, errors.New("invalid P521 compressed point encoding")
+		}
+
+		// Select the positive or negative root, as indicated by the least
+		// significant bit, based on the encoding type byte.
+		otherRoot := new(fiat.P521Element)
+		otherRoot.Sub(otherRoot, y)
+		cond := y.Bytes()[p521ElementLength-1]&1 ^ b[0]&1
+		y.Select(otherRoot, y, int(cond))
+
+		p.x.Set(x)
+		p.y.Set(y)
+		p.z.One()
+		return p, nil
+
+	default:
+		return nil, errors.New("invalid P521 point encoding")
+	}
+}
+
+var _p521B *fiat.P521Element
+var _p521BOnce sync.Once
+
+func p521B() *fiat.P521Element {
+	_p521BOnce.Do(func() {
+		_p521B, _ = new(fiat.P521Element).SetBytes([]byte{0x0, 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, 0x9, 0xe1, 0x56, 0x19, 0x39, 0x51, 0xec, 0x7e, 0x93, 0x7b, 0x16, 0x52, 0xc0, 0xbd, 0x3b, 0xb1, 0xbf, 0x7, 0x35, 0x73, 0xdf, 0x88, 0x3d, 0x2c, 0x34, 0xf1, 0xef, 0x45, 0x1f, 0xd4, 0x6b, 0x50, 0x3f, 0x0})
+	})
+	return _p521B
+}
+
+// p521Polynomial sets y2 to x³ - 3x + b, and returns y2.
+func p521Polynomial(y2, x *fiat.P521Element) *fiat.P521Element {
+	y2.Square(x)
+	y2.Mul(y2, x)
+
+	threeX := new(fiat.P521Element).Add(x, x)
+	threeX.Add(threeX, x)
+	y2.Sub(y2, threeX)
+
+	return y2.Add(y2, p521B())
+}
+
+func p521CheckOnCurve(x, y *fiat.P521Element) error {
+	// y² = x³ - 3x + b
+	rhs := p521Polynomial(new(fiat.P521Element), x)
+	lhs := new(fiat.P521Element).Square(y)
+	if rhs.Equal(lhs) != 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 [1 + 2*p521ElementLength]byte
+	return p.bytes(&out)
+}
+
+func (p *P521Point) bytes(out *[1 + 2*p521ElementLength]byte) []byte {
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new(fiat.P521Element).Invert(p.z)
+	x := new(fiat.P521Element).Mul(p.x, zinv)
+	y := new(fiat.P521Element).Mul(p.y, zinv)
+
+	buf := append(out[:0], 4)
+	buf = append(buf, x.Bytes()...)
+	buf = append(buf, y.Bytes()...)
+	return buf
+}
+
+// BytesX returns the encoding of the x-coordinate of p, as specified in SEC 1,
+// Version 2.0, Section 2.3.5, or an error if p is the point at infinity.
+func (p *P521Point) BytesX() ([]byte, error) {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [p521ElementLength]byte
+	return p.bytesX(&out)
+}
+
+func (p *P521Point) bytesX(out *[p521ElementLength]byte) ([]byte, error) {
+	if p.z.IsZero() == 1 {
+		return nil, errors.New("P521 point is the point at infinity")
+	}
+
+	zinv := new(fiat.P521Element).Invert(p.z)
+	x := new(fiat.P521Element).Mul(p.x, zinv)
+
+	return append(out[:0], x.Bytes()...), nil
+}
+
+// BytesCompressed returns the compressed 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) BytesCompressed() []byte {
+	// This function is outlined to make the allocations inline in the caller
+	// rather than happen on the heap.
+	var out [1 + p521ElementLength]byte
+	return p.bytesCompressed(&out)
+}
+
+func (p *P521Point) bytesCompressed(out *[1 + p521ElementLength]byte) []byte {
+	if p.z.IsZero() == 1 {
+		return append(out[:0], 0)
+	}
+
+	zinv := new(fiat.P521Element).Invert(p.z)
+	x := new(fiat.P521Element).Mul(p.x, zinv)
+	y := new(fiat.P521Element).Mul(p.y, zinv)
+
+	// Encode the sign of the y coordinate (indicated by the least significant
+	// bit) as the encoding type (2 or 3).
+	buf := append(out[:0], 2)
+	buf[0] |= y.Bytes()[p521ElementLength-1] & 1
+	buf = append(buf, x.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
+}
+
+// A p521Table holds the first 15 multiples of a point at offset -1, so [1]P
+// is at table[0], [15]P is at table[14], and [0]P is implicitly the identity
+// point.
+type p521Table [15]*P521Point
+
+// Select selects the n-th multiple of the table base point into p. It works in
+// constant time by iterating over every entry of the table. n must be in [0, 15].
+func (table *p521Table) Select(p *P521Point, n uint8) {
+	if n >= 16 {
+		panic("nistec: internal error: p521Table called with out-of-bounds value")
+	}
+	p.Set(NewP521Point())
+	for i := uint8(1); i < 16; i++ {
+		cond := subtle.ConstantTimeByteEq(i, n)
+		p.Select(table[i-1], p, cond)
+	}
+}
+
+// ScalarMult sets p = scalar * q, and returns p.
+func (p *P521Point) ScalarMult(q *P521Point, scalar []byte) (*P521Point, error) {
+	// Compute a p521Table for the base point q. The explicit NewP521Point
+	// calls get inlined, letting the allocations live on the stack.
+	var table = p521Table{NewP521Point(), NewP521Point(), NewP521Point(),
+		NewP521Point(), NewP521Point(), NewP521Point(), NewP521Point(),
+		NewP521Point(), NewP521Point(), NewP521Point(), NewP521Point(),
+		NewP521Point(), NewP521Point(), NewP521Point(), NewP521Point()}
+	table[0].Set(q)
+	for i := 1; i < 15; i += 2 {
+		table[i].Double(table[i/2])
+		table[i+1].Add(table[i], 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 i, byte := range scalar {
+		// No need to double on the first iteration, as p is the identity at
+		// this point, and [N]∞ = ∞.
+		if i != 0 {
+			p.Double(p)
+			p.Double(p)
+			p.Double(p)
+			p.Double(p)
+		}
+
+		windowValue := byte >> 4
+		table.Select(t, windowValue)
+		p.Add(p, t)
+
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+		p.Double(p)
+
+		windowValue = byte & 0b1111
+		table.Select(t, windowValue)
+		p.Add(p, t)
+	}
+
+	return p, nil
+}
+
+var p521GeneratorTable *[p521ElementLength * 2]p521Table
+var p521GeneratorTableOnce sync.Once
+
+// generatorTable returns a sequence of p521Tables. The first table contains
+// multiples of G. Each successive table is the previous table doubled four
+// times.
+func (p *P521Point) generatorTable() *[p521ElementLength * 2]p521Table {
+	p521GeneratorTableOnce.Do(func() {
+		p521GeneratorTable = new([p521ElementLength * 2]p521Table)
+		base := NewP521Point().SetGenerator()
+		for i := 0; i < p521ElementLength*2; i++ {
+			p521GeneratorTable[i][0] = NewP521Point().Set(base)
+			for j := 1; j < 15; j++ {
+				p521GeneratorTable[i][j] = NewP521Point().Add(p521GeneratorTable[i][j-1], base)
+			}
+			base.Double(base)
+			base.Double(base)
+			base.Double(base)
+			base.Double(base)
+		}
+	})
+	return p521GeneratorTable
+}
+
+// ScalarBaseMult sets p = scalar * B, where B is the canonical generator, and
+// returns p.
+func (p *P521Point) ScalarBaseMult(scalar []byte) (*P521Point, error) {
+	if len(scalar) != p521ElementLength {
+		return nil, errors.New("invalid scalar length")
+	}
+	tables := p.generatorTable()
+
+	// This is also a scalar multiplication with a four-bit window like in
+	// ScalarMult, but in this case the doublings are precomputed. The value
+	// [windowValue]G added at iteration k would normally get doubled
+	// (totIterations-k)×4 times, but with a larger precomputation we can
+	// instead add [2^((totIterations-k)×4)][windowValue]G and avoid the
+	// doublings between iterations.
+	t := NewP521Point()
+	p.Set(NewP521Point())
+	tableIndex := len(tables) - 1
+	for _, byte := range scalar {
+		windowValue := byte >> 4
+		tables[tableIndex].Select(t, windowValue)
+		p.Add(p, t)
+		tableIndex--
+
+		windowValue = byte & 0b1111
+		tables[tableIndex].Select(t, windowValue)
+		p.Add(p, t)
+		tableIndex--
+	}
+
+	return p, nil
+}
+
+// p521Sqrt sets e to a square root of x. If x is not a square, p521Sqrt returns
+// false and e is unchanged. e and x can overlap.
+func p521Sqrt(e, x *fiat.P521Element) (isSquare bool) {
+	candidate := new(fiat.P521Element)
+	p521SqrtCandidate(candidate, x)
+	square := new(fiat.P521Element).Square(candidate)
+	if square.Equal(x) != 1 {
+		return false
+	}
+	e.Set(candidate)
+	return true
+}
+
+// p521SqrtCandidate sets z to a square root candidate for x. z and x must not overlap.
+func p521SqrtCandidate(z, x *fiat.P521Element) {
+	// Since p = 3 mod 4, exponentiation by (p + 1) / 4 yields a square root candidate.
+	//
+	// The sequence of 0 multiplications and 519 squarings is derived from the
+	// following addition chain generated with github.com/mmcloughlin/addchain v0.4.0.
+	//
+	//	return  1 << 519
+	//
+
+	z.Square(x)
+	for s := 1; s < 519; s++ {
+		z.Square(z)
+	}
+}