From b0a75c0ab78b3e3ad7ab01ea7a8fdc89c022b348 Mon Sep 17 00:00:00 2001
From: Filippo Valsorda <hi@filippo.io>
Date: Sat, 26 Jan 2019 22:20:45 -0500
Subject: [PATCH] Import gtank/ed25519#8 and refactor on top of it

---
 fe.go                               | 162 +++++++---------
 internal/group/const.go             |  11 --
 internal/group/ge.go                | 233 +++++++++++------------
 internal/radix51/const.go           |  17 --
 internal/radix51/fe.go              | 277 ++++++++++++++--------------
 internal/radix51/fe_mul.go          |  15 +-
 internal/radix51/fe_mul_amd64.go    |   8 +-
 internal/radix51/fe_mul_amd64.s     |   4 +-
 internal/radix51/fe_square.go       |  15 +-
 internal/radix51/fe_square_amd64.go |   8 +-
 internal/radix51/fe_square_amd64.s  |   4 +-
 internal/radix51/fe_test.go         |  40 ++--
 ristretto255.go                     |  88 ++++-----
 13 files changed, 420 insertions(+), 462 deletions(-)
 delete mode 100644 internal/group/const.go
 delete mode 100644 internal/radix51/const.go

diff --git a/fe.go b/fe.go
index 812a159..c6cb6d0 100644
--- a/fe.go
+++ b/fe.go
@@ -8,126 +8,104 @@ package ristretto255
 import (
 	"math/big"
 
-	. "github.com/gtank/ristretto255/internal/radix51"
+	"github.com/gtank/ristretto255/internal/radix51"
 )
 
-// fePow22523 is from x/crypto/ed25519/internal/edwards25519.
-func fePow22523(out, z *FieldElement) {
-	var t0, t1, t2 FieldElement
-	var i int
+// fePow22523 sets out to z^((p-5)/8). TODO
+func fePow22523(out, z *radix51.FieldElement) *radix51.FieldElement {
+	// Refactored from golang.org/x/crypto/ed25519/internal/edwards25519.
 
-	FeSquare(&t0, z)
-	for i = 1; i < 1; i++ {
-		FeSquare(&t0, &t0)
+	var t0, t1, t2 radix51.FieldElement
+
+	t0.Square(z)
+	for i := 1; i < 1; i++ {
+		t0.Square(&t0)
 	}
-	FeSquare(&t1, &t0)
-	for i = 1; i < 2; i++ {
-		FeSquare(&t1, &t1)
+	t1.Square(&t0)
+	for i := 1; i < 2; i++ {
+		t1.Square(&t1)
 	}
-	FeMul(&t1, z, &t1)
-	FeMul(&t0, &t0, &t1)
-	FeSquare(&t0, &t0)
-	for i = 1; i < 1; i++ {
-		FeSquare(&t0, &t0)
+	t1.Mul(z, &t1)
+	t0.Mul(&t0, &t1)
+	t0.Square(&t0)
+	for i := 1; i < 1; i++ {
+		t0.Square(&t0)
 	}
-	FeMul(&t0, &t1, &t0)
-	FeSquare(&t1, &t0)
-	for i = 1; i < 5; i++ {
-		FeSquare(&t1, &t1)
+	t0.Mul(&t1, &t0)
+	t1.Square(&t0)
+	for i := 1; i < 5; i++ {
+		t1.Square(&t1)
 	}
-	FeMul(&t0, &t1, &t0)
-	FeSquare(&t1, &t0)
-	for i = 1; i < 10; i++ {
-		FeSquare(&t1, &t1)
+	t0.Mul(&t1, &t0)
+	t1.Square(&t0)
+	for i := 1; i < 10; i++ {
+		t1.Square(&t1)
 	}
-	FeMul(&t1, &t1, &t0)
-	FeSquare(&t2, &t1)
-	for i = 1; i < 20; i++ {
-		FeSquare(&t2, &t2)
+	t1.Mul(&t1, &t0)
+	t2.Square(&t1)
+	for i := 1; i < 20; i++ {
+		t2.Square(&t2)
 	}
-	FeMul(&t1, &t2, &t1)
-	FeSquare(&t1, &t1)
-	for i = 1; i < 10; i++ {
-		FeSquare(&t1, &t1)
+	t1.Mul(&t2, &t1)
+	t1.Square(&t1)
+	for i := 1; i < 10; i++ {
+		t1.Square(&t1)
 	}
-	FeMul(&t0, &t1, &t0)
-	FeSquare(&t1, &t0)
-	for i = 1; i < 50; i++ {
-		FeSquare(&t1, &t1)
+	t0.Mul(&t1, &t0)
+	t1.Square(&t0)
+	for i := 1; i < 50; i++ {
+		t1.Square(&t1)
 	}
-	FeMul(&t1, &t1, &t0)
-	FeSquare(&t2, &t1)
-	for i = 1; i < 100; i++ {
-		FeSquare(&t2, &t2)
+	t1.Mul(&t1, &t0)
+	t2.Square(&t1)
+	for i := 1; i < 100; i++ {
+		t2.Square(&t2)
 	}
-	FeMul(&t1, &t2, &t1)
-	FeSquare(&t1, &t1)
-	for i = 1; i < 50; i++ {
-		FeSquare(&t1, &t1)
+	t1.Mul(&t2, &t1)
+	t1.Square(&t1)
+	for i := 1; i < 50; i++ {
+		t1.Square(&t1)
 	}
-	FeMul(&t0, &t1, &t0)
-	FeSquare(&t0, &t0)
-	for i = 1; i < 2; i++ {
-		FeSquare(&t0, &t0)
+	t0.Mul(&t1, &t0)
+	t0.Square(&t0)
+	for i := 1; i < 2; i++ {
+		t0.Square(&t0)
 	}
-	FeMul(out, &t0, z)
+	return out.Mul(&t0, z)
 }
 
-func feSqrtRatio(out, u, v *FieldElement) int {
-	var a, b FieldElement
+// feSqrtRatio sets r to the square root of the ratio of u and v, according to
+// Section 3.1.3 of draft-hdevalence-cfrg-ristretto-00.
+func feSqrtRatio(r, u, v *radix51.FieldElement) (wasSquare int) {
+	var a, b radix51.FieldElement
 
-	// v^3, v^7
-	v3, v7 := &a, &b
-
-	FeSquare(v3, v)  // v^2 = v*v
-	FeMul(v3, v3, v) // v^3 = v^2 * v
-	FeSquare(v7, v3) // v^6 = v^3 * v^3
-	FeMul(v7, v7, v) // v^7 = v^6 * v
+	v3 := a.Mul(a.Square(v), v)  // v^3 = v^2 * v
+	v7 := b.Mul(b.Square(v3), v) // v^7 = (v^3)^2 * v
 
 	// r = (u * v3) * (u * v7)^((p-5)/8)
-	r := out
-	uv3, uv7 := v3, v7   // alias
-	FeMul(uv3, u, v3)    // (u * v3)
-	FeMul(uv7, u, v7)    // (u * v7)
-	fePow22523(uv7, uv7) // (u * v7) ^ ((q - 5)/8)
-	FeMul(r, uv3, uv7)
+	uv3 := a.Mul(u, v3) // (u * v3)
+	uv7 := b.Mul(u, v7) // (u * v7)
+	r.Mul(uv3, fePow22523(r, uv7))
 
-	// done with these ("freeing" a, b)
-	v3, v7, uv3, uv7 = nil, nil, nil, nil
+	check := a.Mul(v, a.Square(r)) // check = v * r^2
 
-	// check = v * r^2
-	check := &a
-	FeMul(check, r, r)     // r^2
-	FeMul(check, check, v) // v * r^2
+	uNeg := b.Neg(u)
+	correctSignSqrt := check.Equal(u)
+	flippedSignSqrt := check.Equal(uNeg)
+	flippedSignSqrtI := check.Equal(uNeg.Mul(uNeg, sqrtM1))
 
-	uneg := &b
-	FeNeg(uneg, u)
-	correct_sign_sqrt := FeEqual(check, u)
-	flipped_sign_sqrt := FeEqual(check, uneg)
-	FeMul(uneg, uneg, sqrtM1)
-	flipped_sign_sqrt_i := FeEqual(check, uneg)
-
-	// done with these ("freeing" a, b)
-	check, uneg = nil, nil
-
-	// r_prime = SQRT_M1 * r
+	rPrime := b.Mul(r, sqrtM1) // r_prime = SQRT_M1 * r
 	// r = CT_SELECT(r_prime IF flipped_sign_sqrt | flipped_sign_sqrt_i ELSE r)
-	r_prime := &a
-	FeMul(r_prime, r, sqrtM1)
-	FeSelect(r, r_prime, r, flipped_sign_sqrt|flipped_sign_sqrt_i)
+	r.Select(rPrime, r, flippedSignSqrt|flippedSignSqrtI)
 
-	FeAbs(r, r)
-	was_square := correct_sign_sqrt | flipped_sign_sqrt
-
-	return was_square
+	r.Abs(r) // Choose the nonnegative square root.
+	return correctSignSqrt | flippedSignSqrt
 }
 
-func fieldElementFromDecimal(s string) *FieldElement {
+func fieldElementFromDecimal(s string) *radix51.FieldElement {
 	n, ok := new(big.Int).SetString(s, 10)
 	if !ok {
 		panic("ristretto255: not a valid decimal: " + s)
 	}
-	var fe FieldElement
-	FeFromBig(&fe, n)
-	return &fe
+	return new(radix51.FieldElement).FromBig(n)
 }
diff --git a/internal/group/const.go b/internal/group/const.go
deleted file mode 100644
index f5f872b..0000000
--- a/internal/group/const.go
+++ /dev/null
@@ -1,11 +0,0 @@
-package group
-
-import "github.com/gtank/ristretto255/internal/radix51"
-
-var (
-	// d, a constant in the curve equation
-	D radix51.FieldElement = [5]uint64{929955233495203, 466365720129213, 1662059464998953, 2033849074728123, 1442794654840575}
-
-	// 2*d, used in addition formula
-	D2 radix51.FieldElement = [5]uint64{1859910466990425, 932731440258426, 1072319116312658, 1815898335770999, 633789495995903}
-)
diff --git a/internal/group/ge.go b/internal/group/ge.go
index 36094e2..2205922 100644
--- a/internal/group/ge.go
+++ b/internal/group/ge.go
@@ -1,13 +1,21 @@
-// Implements group logic for the Ed25519 curve.
+// Copyright (c) 2017 George Tankersley. All rights reserved.
+// Copyright (c) 2019 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 group mplements group logic for the Ed25519 curve.
 package group
 
 import (
 	"math/big"
 
-	field "github.com/gtank/ristretto255/internal/radix51"
+	"github.com/gtank/ristretto255/internal/radix51"
 )
 
+// D is a constant in the curve equation.
+var D = &radix51.FieldElement{929955233495203, 466365720129213,
+	1662059464998953, 2033849074728123, 1442794654840575}
+
 // From EFD https://hyperelliptic.org/EFD/g1p/auto-twisted-extended-1.html
 // An elliptic curve in twisted Edwards form has parameters a, d and coordinates
 // x, y satisfying the following equations:
@@ -24,18 +32,19 @@ import (
 // This representation was introduced in the HisilWongCarterDawson paper "Twisted
 // Edwards curves revisited" (Asiacrypt 2008).
 type ExtendedGroupElement struct {
-	X, Y, Z, T field.FieldElement
+	X, Y, Z, T radix51.FieldElement
 }
 
 // Converts (x,y) to (X:Y:T:Z) extended coordinates, or "P3" in ref10. As
 // described in "Twisted Edwards Curves Revisited", Hisil-Wong-Carter-Dawson
 // 2008, Section 3.1 (https://eprint.iacr.org/2008/522.pdf)
 // See also https://hyperelliptic.org/EFD/g1p/auto-twisted-extended-1.html#addition-add-2008-hwcd-3
-func (v *ExtendedGroupElement) FromAffine(x, y *big.Int) {
-	field.FeFromBig(&v.X, x)
-	field.FeFromBig(&v.Y, y)
-	field.FeMul(&v.T, &v.X, &v.Y)
-	field.FeOne(&v.Z)
+func (v *ExtendedGroupElement) FromAffine(x, y *big.Int) *ExtendedGroupElement {
+	v.X.FromBig(x)
+	v.Y.FromBig(y)
+	v.T.Mul(&v.X, &v.Y)
+	v.Z.One()
+	return v
 }
 
 // Extended coordinates are XYZT with x = X/Z, y = Y/Z, or the "P3"
@@ -43,57 +52,56 @@ func (v *ExtendedGroupElement) FromAffine(x, y *big.Int) {
 // from projective to affine. Per HWCD, it is safe to move from extended to
 // projective by simply ignoring T.
 func (v *ExtendedGroupElement) ToAffine() (*big.Int, *big.Int) {
-	var x, y, zinv field.FieldElement
+	var x, y, zinv radix51.FieldElement
 
-	field.FeInvert(&zinv, &v.Z)
-	field.FeMul(&x, &v.X, &zinv)
-	field.FeMul(&y, &v.Y, &zinv)
+	zinv.Invert(&v.Z)
+	x.Mul(&v.X, &zinv)
+	y.Mul(&v.Y, &zinv)
 
-	return field.FeToBig(&x), field.FeToBig(&y)
+	return x.ToBig(), y.ToBig()
 }
 
 // Per HWCD, it is safe to move from extended to projective by simply ignoring T.
-func (v *ExtendedGroupElement) ToProjective() *ProjectiveGroupElement {
-	var p ProjectiveGroupElement
-
-	field.FeCopy(&p.X, &v.X)
-	field.FeCopy(&p.Y, &v.Y)
-	field.FeCopy(&p.Z, &v.Z)
-
-	return &p
+func (v *ExtendedGroupElement) ToProjective(p *ProjectiveGroupElement) {
+	p.X.Set(&v.X)
+	p.Y.Set(&v.Y)
+	p.Z.Set(&v.Z)
 }
 
 func (v *ExtendedGroupElement) Zero() *ExtendedGroupElement {
-	field.FeZero(&v.X)
-	field.FeOne(&v.Y)
-	field.FeOne(&v.Z)
-	field.FeZero(&v.T)
+	v.X.Zero()
+	v.Y.One()
+	v.Z.One()
+	v.T.Zero()
 	return v
 }
 
+var twoD = &radix51.FieldElement{1859910466990425, 932731440258426,
+	1072319116312658, 1815898335770999, 633789495995903}
+
 // This is the same addition formula everyone uses, "add-2008-hwcd-3".
 // https://hyperelliptic.org/EFD/g1p/auto-twisted-extended-1.html#addition-add-2008-hwcd-3
 // TODO We know Z1=1 and Z2=1 here, so mmadd-2008-hwcd-3 (6M + 1S + 1*k + 9add) could apply
 func (v *ExtendedGroupElement) Add(p1, p2 *ExtendedGroupElement) *ExtendedGroupElement {
-	var tmp1, tmp2, A, B, C, D, E, F, G, H field.FieldElement
-	field.FeSub(&tmp1, &p1.Y, &p1.X) // tmp1 <-- Y1-X1
-	field.FeSub(&tmp2, &p2.Y, &p2.X) // tmp2 <-- Y2-X2
-	field.FeMul(&A, &tmp1, &tmp2)    // A <-- tmp1*tmp2 = (Y1-X1)*(Y2-X2)
-	field.FeAdd(&tmp1, &p1.Y, &p1.X) // tmp1 <-- Y1+X1
-	field.FeAdd(&tmp2, &p2.Y, &p2.X) // tmp2 <-- Y2+X2
-	field.FeMul(&B, &tmp1, &tmp2)    // B <-- tmp1*tmp2 = (Y1+X1)*(Y2+X2)
-	field.FeMul(&tmp1, &p1.T, &p2.T) // tmp1 <-- T1*T2
-	field.FeMul(&C, &tmp1, &D2)      // C <-- tmp1*2d = T1*2d*T2
-	field.FeMul(&tmp1, &p1.Z, &p2.Z) // tmp1 <-- Z1*Z2
-	field.FeAdd(&D, &tmp1, &tmp1)    // D <-- tmp1 + tmp1 = 2*Z1*Z2
-	field.FeSub(&E, &B, &A)          // E <-- B-A
-	field.FeSub(&F, &D, &C)          // F <-- D-C
-	field.FeAdd(&G, &D, &C)          // G <-- D+C
-	field.FeAdd(&H, &B, &A)          // H <-- B+A
-	field.FeMul(&v.X, &E, &F)        // X3 <-- E*F
-	field.FeMul(&v.Y, &G, &H)        // Y3 <-- G*H
-	field.FeMul(&v.T, &E, &H)        // T3 <-- E*H
-	field.FeMul(&v.Z, &F, &G)        // Z3 <-- F*G
+	var tmp1, tmp2, A, B, C, D, E, F, G, H radix51.FieldElement
+	tmp1.Sub(&p1.Y, &p1.X) // tmp1 <-- Y1-X1
+	tmp2.Sub(&p2.Y, &p2.X) // tmp2 <-- Y2-X2
+	A.Mul(&tmp1, &tmp2)    // A <-- tmp1*tmp2 = (Y1-X1)*(Y2-X2)
+	tmp1.Add(&p1.Y, &p1.X) // tmp1 <-- Y1+X1
+	tmp2.Add(&p2.Y, &p2.X) // tmp2 <-- Y2+X2
+	B.Mul(&tmp1, &tmp2)    // B <-- tmp1*tmp2 = (Y1+X1)*(Y2+X2)
+	tmp1.Mul(&p1.T, &p2.T) // tmp1 <-- T1*T2
+	C.Mul(&tmp1, twoD)     // C <-- tmp1*2d = T1*2*d*T2
+	tmp1.Mul(&p1.Z, &p2.Z) // tmp1 <-- Z1*Z2
+	D.Add(&tmp1, &tmp1)    // D <-- tmp1 + tmp1 = 2*Z1*Z2
+	E.Sub(&B, &A)          // E <-- B-A
+	F.Sub(&D, &C)          // F <-- D-C
+	G.Add(&D, &C)          // G <-- D+C
+	H.Add(&B, &A)          // H <-- B+A
+	v.X.Mul(&E, &F)        // X3 <-- E*F
+	v.Y.Mul(&G, &H)        // Y3 <-- G*H
+	v.T.Mul(&E, &H)        // T3 <-- E*H
+	v.Z.Mul(&F, &G)        // Z3 <-- F*G
 	return v
 }
 
@@ -120,46 +128,36 @@ func (v *ExtendedGroupElement) Add(p1, p2 *ExtendedGroupElement) *ExtendedGroupE
 // instead of another point in extended coordinates. I have implemented it
 // this way to see if more straightforward code is worth the (hopefully small)
 // performance tradeoff.
-func (v *ExtendedGroupElement) Double() *ExtendedGroupElement {
+func (v *ExtendedGroupElement) Double(u *ExtendedGroupElement) *ExtendedGroupElement {
 	// TODO: Convert to projective coordinates? Section 4.3 mixed doubling?
-	// TODO: make a decision about how these APIs work wrt chaining/smashing
-	// *v = *(v.ToProjective().Double().ToExtended())
-	// return v
 
-	var A, B, C, D, E, F, G, H field.FieldElement
+	var A, B, C, D, E, F, G, H radix51.FieldElement
 
 	// A ← X1^2, B ← Y1^2
-	field.FeSquare(&A, &v.X)
-	field.FeSquare(&B, &v.Y)
+	A.Square(&u.X)
+	B.Square(&u.Y)
 
 	// C ← 2*Z1^2
-	field.FeSquare(&C, &v.Z)
-	field.FeAdd(&C, &C, &C) // TODO should probably implement FeSquare2
+	C.Square(&u.Z)
+	C.Add(&C, &C) // TODO should probably implement FeSquare2
 
 	// D ← -1*A
-	field.FeNeg(&D, &A) // implemented as substraction
+	D.Neg(&A) // implemented as substraction
 
 	// E ← (X1+Y1)^2 − A − B
-	var t0 field.FieldElement
-	field.FeAdd(&t0, &v.X, &v.Y)
-	field.FeSquare(&t0, &t0)
-	field.FeSub(&E, &t0, &A)
-	field.FeSub(&E, &E, &B)
+	var t0 radix51.FieldElement
+	t0.Add(&u.X, &u.Y)
+	t0.Square(&t0)
+	E.Sub(&t0, &A)
+	E.Sub(&E, &B)
 
-	// G ← D+B
-	field.FeAdd(&G, &D, &B)
-	// F ← G−C
-	field.FeSub(&F, &G, &C)
-	// H ← D−B
-	field.FeSub(&H, &D, &B)
-	// X3 ← E*F
-	field.FeMul(&v.X, &E, &F)
-	// Y3 ← G*H
-	field.FeMul(&v.Y, &G, &H)
-	// T3 ← E*H
-	field.FeMul(&v.T, &E, &H)
-	// Z3 ← F*G
-	field.FeMul(&v.Z, &F, &G)
+	G.Add(&D, &B)   // G ← D+B
+	F.Sub(&G, &C)   // F ← G−C
+	H.Sub(&D, &B)   // H ← D−B
+	v.X.Mul(&E, &F) // X3 ← E*F
+	v.Y.Mul(&G, &H) // Y3 ← G*H
+	v.T.Mul(&E, &H) // T3 ← E*H
+	v.Z.Mul(&F, &G) // Z3 ← F*G
 
 	return v
 }
@@ -168,43 +166,40 @@ func (v *ExtendedGroupElement) Double() *ExtendedGroupElement {
 // representation in ref10. This representation has a cheaper doubling formula
 // than extended coordinates.
 type ProjectiveGroupElement struct {
-	X, Y, Z field.FieldElement
+	X, Y, Z radix51.FieldElement
 }
 
-func (v *ProjectiveGroupElement) FromAffine(x, y *big.Int) {
-	field.FeFromBig(&v.X, x)
-	field.FeFromBig(&v.Y, y)
-	field.FeOne(&v.Z)
+func (v *ProjectiveGroupElement) FromAffine(x, y *big.Int) *ProjectiveGroupElement {
+	v.X.FromBig(x)
+	v.Y.FromBig(y)
+	v.Z.Zero()
+	return v
 }
 
 func (v *ProjectiveGroupElement) ToAffine() (*big.Int, *big.Int) {
-	var x, y, zinv field.FieldElement
+	var x, y, zinv radix51.FieldElement
 
-	field.FeInvert(&zinv, &v.Z)
-	field.FeMul(&x, &v.X, &zinv)
-	field.FeMul(&y, &v.Y, &zinv)
+	zinv.Invert(&v.Z)
+	x.Mul(&v.X, &zinv)
+	y.Mul(&v.Y, &zinv)
 
-	return field.FeToBig(&x), field.FeToBig(&y)
+	return x.ToBig(), y.ToBig()
 }
 
 // HWCD Section 3: "Given (X : Y : Z) in [projective coordinates] passing to
 // [extended coordinates, (X : Y : T : Z)] can be performed in 3M+1S by computing
 // (XZ, YZ, XY, Z^2)"
-func (v *ProjectiveGroupElement) ToExtended() *ExtendedGroupElement {
-	var r ExtendedGroupElement
-
-	field.FeMul(&r.X, &v.X, &v.Z)
-	field.FeMul(&r.Y, &v.Y, &v.Z)
-	field.FeMul(&r.T, &v.X, &v.Y)
-	field.FeSquare(&r.Z, &v.Z)
-
-	return &r
+func (v *ProjectiveGroupElement) ToExtended(r *ExtendedGroupElement) {
+	r.X.Mul(&v.X, &v.Z)
+	r.Y.Mul(&v.Y, &v.Z)
+	r.T.Mul(&v.X, &v.Y)
+	r.Z.Square(&v.Z)
 }
 
 func (v *ProjectiveGroupElement) Zero() *ProjectiveGroupElement {
-	field.FeZero(&v.X)
-	field.FeOne(&v.Y)
-	field.FeOne(&v.Z)
+	v.X.Zero()
+	v.Y.One()
+	v.Z.One()
 	return v
 }
 
@@ -229,44 +224,30 @@ func (v *ProjectiveGroupElement) Zero() *ProjectiveGroupElement {
 // This assumption is one reason why this package is internal. For instance, it
 // will not hold throughout a Montgomery ladder, when we convert to projective
 // from possibly arbitrary extended coordinates.
-func (v *ProjectiveGroupElement) DoubleZ1() *ProjectiveGroupElement {
-	// TODO This function is inconsistent with the other ones in that it
-	// returns a copy rather than smashing the receiver. It doesn't matter
-	// because it is always called on ephemeral intermediate values, but should
-	// fix.
-	var p, q ProjectiveGroupElement
-	var t0, t1 field.FieldElement
+func (v *ProjectiveGroupElement) DoubleZ1(u *ProjectiveGroupElement) *ProjectiveGroupElement {
+	var B, C, D, E, F radix51.FieldElement
 
-	p = *v
+	if u.Z.Equal(radix51.Zero) != 1 {
+		panic("ed25519: DoubleZ1 called with Z != 1")
+	}
 
-	// C = X1^2, D = Y1^2
-	field.FeSquare(&t0, &p.X)
-	field.FeSquare(&t1, &p.Y)
-
-	// B = (X1+Y1)^2
-	field.FeAdd(&p.Z, &p.X, &p.Y) // Z is irrelevant but already allocated
-	field.FeSquare(&q.X, &p.Z)
-
-	// E = a*C where a = -1
-	field.FeNeg(&q.Z, &t0)
-
-	// F = E + D
-	field.FeAdd(&p.X, &q.Z, &t1)
+	B.Square(B.Add(&u.X, &u.Y)) // B = (X1+Y1)^2
+	C.Square(&u.X)              // C = X1^2
+	D.Square(&u.Y)              // D = Y1^2
+	E.Neg(&C)                   // E = a*C where a = -1
+	F.Add(&E, &D)               // F = E + D
 
 	// X3 = (B-C-D)*(F-2)
-	field.FeSub(&p.Y, &q.X, &t0)
-	field.FeSub(&p.Y, &p.Y, &t1)
-	field.FeSub(&p.Z, &p.X, &field.FieldTwo)
-	field.FeMul(&q.X, &p.Y, &p.Z)
+	v.Y.Sub(v.Y.Sub(&B, &C), &D)
+	v.X.Mul(&v.Y, v.X.Sub(&F, radix51.Two))
 
 	// Y3 = F*(E-D)
-	field.FeSub(&p.Y, &q.Z, &t1)
-	field.FeMul(&q.Y, &p.X, &p.Y)
+	v.Y.Mul(&F, v.Y.Sub(&E, &D))
 
 	// Z3 = F^2 - 2*F
-	field.FeSquare(&q.Z, &p.X)
-	field.FeSub(&q.Z, &q.Z, &p.X)
-	field.FeSub(&q.Z, &q.Z, &p.X)
+	v.Z.Square(&F)
+	v.Z.Sub(&v.Z, &F)
+	v.Z.Sub(&v.Z, &F)
 
-	return &q
+	return v
 }
diff --git a/internal/radix51/const.go b/internal/radix51/const.go
deleted file mode 100644
index b650ff8..0000000
--- a/internal/radix51/const.go
+++ /dev/null
@@ -1,17 +0,0 @@
-// Copyright (c) 2017 George Tankersley. All rights reserved.
-// Use of this source code is governed by a BSD-style
-// license that can be found in the LICENSE file.
-
-// Constants used in the implementation of GF(2^255-19) field arithmetic.
-package radix51
-
-const (
-	// The vaule 2^51-1, used in carry propagation
-	maskLow51Bits = uint64(1)<<51 - 1
-)
-
-var (
-	FieldZero FieldElement = [5]uint64{0, 0, 0, 0, 0}
-	FieldOne  FieldElement = [5]uint64{1, 0, 0, 0, 0}
-	FieldTwo  FieldElement = [5]uint64{2, 0, 0, 0, 0}
-)
diff --git a/internal/radix51/fe.go b/internal/radix51/fe.go
index 0756e87..1320839 100644
--- a/internal/radix51/fe.go
+++ b/internal/radix51/fe.go
@@ -15,90 +15,106 @@ import (
 
 // FieldElement represents an element of the field GF(2^255-19). An element t
 // represents the integer t[0] + t[1]*2^51 + t[2]*2^102 + t[3]*2^153 +
-// t[4]*2^204.
+// t[4]*2^204. The zero value is a valid zero element.
 type FieldElement [5]uint64
 
-func FeZero(v *FieldElement) {
+const (
+	// The vaule 2^51-1, used in carry propagation
+	maskLow51Bits = uint64(1)<<51 - 1
+)
+
+var (
+	Zero     = &FieldElement{0, 0, 0, 0, 0}
+	One      = &FieldElement{1, 0, 0, 0, 0}
+	Two      = &FieldElement{2, 0, 0, 0, 0}
+	MinusOne = new(FieldElement).Neg(One)
+)
+
+func (v *FieldElement) Zero() *FieldElement {
 	v[0] = 0
 	v[1] = 0
 	v[2] = 0
 	v[3] = 0
 	v[4] = 0
+	return v
 }
 
-func FeOne(v *FieldElement) {
+func (v *FieldElement) One() *FieldElement {
 	v[0] = 1
 	v[1] = 0
 	v[2] = 0
 	v[3] = 0
 	v[4] = 0
+	return v
 }
 
 // SetInt sets the receiving FieldElement to the specified small integer.
-func SetInt(v *FieldElement, x uint64) {
+func (v *FieldElement) SetInt(x uint64) *FieldElement {
 	v[0] = x
 	v[1] = 0
 	v[2] = 0
 	v[3] = 0
 	v[4] = 0
+	return v
 }
 
-func FeReduce(t, v *FieldElement) {
-	// Copy v
-	*t = *v
+func (v *FieldElement) Reduce(u *FieldElement) *FieldElement {
+	v.Set(u)
 
 	// Lev v = v[0] + v[1]*2^51 + v[2]*2^102 + v[3]*2^153 + v[4]*2^204
 	// Reduce each limb below 2^51, propagating carries.
-	t[1] += t[0] >> 51
-	t[0] = t[0] & maskLow51Bits
-	t[2] += t[1] >> 51
-	t[1] = t[1] & maskLow51Bits
-	t[3] += t[2] >> 51
-	t[2] = t[2] & maskLow51Bits
-	t[4] += t[3] >> 51
-	t[3] = t[3] & maskLow51Bits
-	t[0] += (t[4] >> 51) * 19
-	t[4] = t[4] & maskLow51Bits
+	v[1] += v[0] >> 51
+	v[0] = v[0] & maskLow51Bits
+	v[2] += v[1] >> 51
+	v[1] = v[1] & maskLow51Bits
+	v[3] += v[2] >> 51
+	v[2] = v[2] & maskLow51Bits
+	v[4] += v[3] >> 51
+	v[3] = v[3] & maskLow51Bits
+	v[0] += (v[4] >> 51) * 19
+	v[4] = v[4] & maskLow51Bits
 
-	// We now hate a field element t < 2^255, but need t <= 2^255-19
+	// We now hate a field element v < 2^255, but need v <= 2^255-19
 	// TODO Document why this works. It's the elaborate comment about r = h-pq etc etc.
 
 	// Get the carry bit
-	c := (t[0] + 19) >> 51
-	c = (t[1] + c) >> 51
-	c = (t[2] + c) >> 51
-	c = (t[3] + c) >> 51
-	c = (t[4] + c) >> 51
+	c := (v[0] + 19) >> 51
+	c = (v[1] + c) >> 51
+	c = (v[2] + c) >> 51
+	c = (v[3] + c) >> 51
+	c = (v[4] + c) >> 51
 
-	t[0] += 19 * c
+	v[0] += 19 * c
 
-	t[1] += t[0] >> 51
-	t[0] = t[0] & maskLow51Bits
-	t[2] += t[1] >> 51
-	t[1] = t[1] & maskLow51Bits
-	t[3] += t[2] >> 51
-	t[2] = t[2] & maskLow51Bits
-	t[4] += t[3] >> 51
-	t[3] = t[3] & maskLow51Bits
+	v[1] += v[0] >> 51
+	v[0] = v[0] & maskLow51Bits
+	v[2] += v[1] >> 51
+	v[1] = v[1] & maskLow51Bits
+	v[3] += v[2] >> 51
+	v[2] = v[2] & maskLow51Bits
+	v[4] += v[3] >> 51
+	v[3] = v[3] & maskLow51Bits
 	// no additional carry
-	t[4] = t[4] & maskLow51Bits
+	v[4] = v[4] & maskLow51Bits
+
+	return v
 }
 
-// FeAdd sets out = a + b. Long sequences of additions without reduction that
+// Add sets v = a + b. Long sequences of additions without reduction that
 // let coefficients grow larger than 54 bits would be a problem. Paper
 // cautions: "do not have such sequences of additions".
-func FeAdd(out, a, b *FieldElement) {
-	out[0] = a[0] + b[0]
-	out[1] = a[1] + b[1]
-	out[2] = a[2] + b[2]
-	out[3] = a[3] + b[3]
-	out[4] = a[4] + b[4]
+func (v *FieldElement) Add(a, b *FieldElement) *FieldElement {
+	v[0] = a[0] + b[0]
+	v[1] = a[1] + b[1]
+	v[2] = a[2] + b[2]
+	v[3] = a[3] + b[3]
+	v[4] = a[4] + b[4]
+	return v
 }
 
-// FeSub sets out = a - b
-func FeSub(out, a, b *FieldElement) {
-	var t FieldElement
-	t = *b
+// Sub sets v = a - b.
+func (v *FieldElement) Sub(a, b *FieldElement) *FieldElement {
+	t := *b
 
 	// Reduce each limb below 2^51, propagating carries. Ensures that results
 	// fit within the limbs. This would not be required for reduced input.
@@ -115,90 +131,91 @@ func FeSub(out, a, b *FieldElement) {
 
 	// This is slightly more complicated. Because we use unsigned coefficients, we
 	// first add a multiple of p and then subtract.
-	out[0] = (a[0] + 0xFFFFFFFFFFFDA) - t[0]
-	out[1] = (a[1] + 0xFFFFFFFFFFFFE) - t[1]
-	out[2] = (a[2] + 0xFFFFFFFFFFFFE) - t[2]
-	out[3] = (a[3] + 0xFFFFFFFFFFFFE) - t[3]
-	out[4] = (a[4] + 0xFFFFFFFFFFFFE) - t[4]
+	v[0] = (a[0] + 0xFFFFFFFFFFFDA) - t[0]
+	v[1] = (a[1] + 0xFFFFFFFFFFFFE) - t[1]
+	v[2] = (a[2] + 0xFFFFFFFFFFFFE) - t[2]
+	v[3] = (a[3] + 0xFFFFFFFFFFFFE) - t[3]
+	v[4] = (a[4] + 0xFFFFFFFFFFFFE) - t[4]
+
+	return v
 }
 
-// FeNeg sets out = -a
-func FeNeg(out, a *FieldElement) {
-	var t FieldElement
-	FeZero(&t)
-	FeSub(out, &t, a)
+// Neg sets v = -a.
+func (v *FieldElement) Neg(a *FieldElement) *FieldElement {
+	return v.Sub(Zero, a)
 }
 
-// FeInvert sets out = 1/z mod p by calculating z^(p-2), p-2 = 2^255 - 21.
-func FeInvert(out, z *FieldElement) {
+// Invert sets v = 1/z mod p by calculating z^(p-2), p-2 = 2^255 - 21.
+func (v *FieldElement) Invert(z *FieldElement) *FieldElement {
 	// Inversion is implemented as exponentiation with exponent p − 2. It uses the
 	// same sequence of 255 squarings and 11 multiplications as [Curve25519].
 	var z2, z9, z11, z2_5_0, z2_10_0, z2_20_0, z2_50_0, z2_100_0, t FieldElement
 
-	FeSquare(&z2, z)        // 2
-	FeSquare(&t, &z2)       // 4
-	FeSquare(&t, &t)        // 8
-	FeMul(&z9, &t, z)       // 9
-	FeMul(&z11, &z9, &z2)   // 11
-	FeSquare(&t, &z11)      // 22
-	FeMul(&z2_5_0, &t, &z9) // 2^5 - 2^0 = 31
+	z2.Square(z)        // 2
+	t.Square(&z2)       // 4
+	t.Square(&t)        // 8
+	z9.Mul(&t, z)       // 9
+	z11.Mul(&z9, &z2)   // 11
+	t.Square(&z11)      // 22
+	z2_5_0.Mul(&t, &z9) // 2^5 - 2^0 = 31
 
-	FeSquare(&t, &z2_5_0) // 2^6 - 2^1
+	t.Square(&z2_5_0) // 2^6 - 2^1
 	for i := 0; i < 4; i++ {
-		FeSquare(&t, &t) // 2^10 - 2^5
+		t.Square(&t) // 2^10 - 2^5
 	}
-	FeMul(&z2_10_0, &t, &z2_5_0) // 2^10 - 2^0
+	z2_10_0.Mul(&t, &z2_5_0) // 2^10 - 2^0
 
-	FeSquare(&t, &z2_10_0) // 2^11 - 2^1
+	t.Square(&z2_10_0) // 2^11 - 2^1
 	for i := 0; i < 9; i++ {
-		FeSquare(&t, &t) // 2^20 - 2^10
+		t.Square(&t) // 2^20 - 2^10
 	}
-	FeMul(&z2_20_0, &t, &z2_10_0) // 2^20 - 2^0
+	z2_20_0.Mul(&t, &z2_10_0) // 2^20 - 2^0
 
-	FeSquare(&t, &z2_20_0) // 2^21 - 2^1
+	t.Square(&z2_20_0) // 2^21 - 2^1
 	for i := 0; i < 19; i++ {
-		FeSquare(&t, &t) // 2^40 - 2^20
+		t.Square(&t) // 2^40 - 2^20
 	}
-	FeMul(&t, &t, &z2_20_0) // 2^40 - 2^0
+	t.Mul(&t, &z2_20_0) // 2^40 - 2^0
 
-	FeSquare(&t, &t) // 2^41 - 2^1
+	t.Square(&t) // 2^41 - 2^1
 	for i := 0; i < 9; i++ {
-		FeSquare(&t, &t) // 2^50 - 2^10
+		t.Square(&t) // 2^50 - 2^10
 	}
-	FeMul(&z2_50_0, &t, &z2_10_0) // 2^50 - 2^0
+	z2_50_0.Mul(&t, &z2_10_0) // 2^50 - 2^0
 
-	FeSquare(&t, &z2_50_0) // 2^51 - 2^1
+	t.Square(&z2_50_0) // 2^51 - 2^1
 	for i := 0; i < 49; i++ {
-		FeSquare(&t, &t) // 2^100 - 2^50
+		t.Square(&t) // 2^100 - 2^50
 	}
-	FeMul(&z2_100_0, &t, &z2_50_0) // 2^100 - 2^0
+	z2_100_0.Mul(&t, &z2_50_0) // 2^100 - 2^0
 
-	FeSquare(&t, &z2_100_0) // 2^101 - 2^1
+	t.Square(&z2_100_0) // 2^101 - 2^1
 	for i := 0; i < 99; i++ {
-		FeSquare(&t, &t) // 2^200 - 2^100
+		t.Square(&t) // 2^200 - 2^100
 	}
-	FeMul(&t, &t, &z2_100_0) // 2^200 - 2^0
+	t.Mul(&t, &z2_100_0) // 2^200 - 2^0
 
-	FeSquare(&t, &t) // 2^201 - 2^1
+	t.Square(&t) // 2^201 - 2^1
 	for i := 0; i < 49; i++ {
-		FeSquare(&t, &t) // 2^250 - 2^50
+		t.Square(&t) // 2^250 - 2^50
 	}
-	FeMul(&t, &t, &z2_50_0) // 2^250 - 2^0
+	t.Mul(&t, &z2_50_0) // 2^250 - 2^0
 
-	FeSquare(&t, &t) // 2^251 - 2^1
-	FeSquare(&t, &t) // 2^252 - 2^2
-	FeSquare(&t, &t) // 2^253 - 2^3
-	FeSquare(&t, &t) // 2^254 - 2^4
-	FeSquare(&t, &t) // 2^255 - 2^5
+	t.Square(&t) // 2^251 - 2^1
+	t.Square(&t) // 2^252 - 2^2
+	t.Square(&t) // 2^253 - 2^3
+	t.Square(&t) // 2^254 - 2^4
+	t.Square(&t) // 2^255 - 2^5
 
-	FeMul(out, &t, &z11) // 2^255 - 21
+	return v.Mul(&t, &z11) // 2^255 - 21
 }
 
-func FeCopy(out, in *FieldElement) {
-	copy(out[:], in[:])
+func (v *FieldElement) Set(a *FieldElement) *FieldElement {
+	*v = *a
+	return v
 }
 
-func FeFromBytes(v *FieldElement, x *[32]byte) {
+func (v *FieldElement) FromBytes(x *[32]byte) *FieldElement {
 	v[0] = uint64(x[0])
 	v[0] |= uint64(x[1]) << 8
 	v[0] |= uint64(x[2]) << 16
@@ -239,11 +256,12 @@ func FeFromBytes(v *FieldElement, x *[32]byte) {
 	v[4] |= uint64(x[29]) << 28
 	v[4] |= uint64(x[30]) << 36
 	v[4] |= uint64(x[31]&127) << 44
+
+	return v
 }
 
-func FeToBytes(r *[32]byte, v *FieldElement) {
-	var t FieldElement
-	FeReduce(&t, v)
+func (v *FieldElement) ToBytes(r *[32]byte) {
+	t := new(FieldElement).Reduce(v)
 
 	r[0] = byte(t[0] & 0xff)
 	r[1] = byte((t[0] >> 8) & 0xff)
@@ -287,7 +305,7 @@ func FeToBytes(r *[32]byte, v *FieldElement) {
 	r[31] = byte((t[4] >> 44))
 }
 
-func FeFromBig(h *FieldElement, num *big.Int) {
+func (v *FieldElement) FromBig(num *big.Int) *FieldElement {
 	var buf [32]byte
 
 	offset := 0
@@ -305,12 +323,12 @@ func FeFromBig(h *FieldElement, num *big.Int) {
 		}
 	}
 
-	FeFromBytes(h, &buf)
+	return v.FromBytes(&buf)
 }
 
-func FeToBig(h *FieldElement) *big.Int {
+func (v *FieldElement) ToBig() *big.Int {
 	var buf [32]byte
-	FeToBytes(&buf, h) // does a reduction
+	v.ToBytes(&buf) // does a reduction
 
 	numWords := 256 / bits.UintSize
 	words := make([]big.Word, numWords)
@@ -333,48 +351,39 @@ func FeToBig(h *FieldElement) *big.Int {
 	return out.SetBits(words)
 }
 
-// FeEqual returns 1 if a and b are equal, and 0 otherwise.
-func FeEqual(a, b *FieldElement) int {
-	var sa, sb [32]byte
-	FeToBytes(&sa, a)
-	FeToBytes(&sb, b)
-	return subtle.ConstantTimeCompare(sa[:], sb[:])
+// Equal returns 1 if v and u are equal, and 0 otherwise.
+func (v *FieldElement) Equal(u *FieldElement) int {
+	var sa, sv [32]byte
+	u.ToBytes(&sa)
+	v.ToBytes(&sv)
+	return subtle.ConstantTimeCompare(sa[:], sv[:])
 }
 
-// FeSelect sets out to v if cond == 1, and to u if cond == 0.
-// out, v and u are allowed to overlap.
-func FeSelect(out, v, u *FieldElement, cond int) {
-	b := uint64(cond) * 0xffffffffffffffff
-	out[0] = (b & v[0]) | (^b & u[0])
-	out[1] = (b & v[1]) | (^b & u[1])
-	out[2] = (b & v[2]) | (^b & u[2])
-	out[3] = (b & v[3]) | (^b & u[3])
-	out[4] = (b & v[4]) | (^b & u[4])
+// Select sets v to a if cond == 1, and to b if cond == 0.
+// v, a and b are allowed to overlap.
+func (v *FieldElement) Select(a, b *FieldElement, cond int) *FieldElement {
+	m := uint64(cond) * 0xffffffffffffffff
+	v[0] = (m & a[0]) | (^m & b[0])
+	v[1] = (m & a[1]) | (^m & b[1])
+	v[2] = (m & a[2]) | (^m & b[2])
+	v[3] = (m & a[3]) | (^m & b[3])
+	v[4] = (m & a[4]) | (^m & b[4])
+	return v
 }
 
-// FeCondNeg sets u to -u if cond == 1, and to u if cond == 0.
-func FeCondNeg(u *FieldElement, cond int) {
-	var neg FieldElement
-	FeNeg(&neg, u)
-
-	b := uint64(cond) * 0xffffffffffffffff
-	u[0] ^= b & (u[0] ^ neg[0])
-	u[1] ^= b & (u[1] ^ neg[1])
-	u[2] ^= b & (u[2] ^ neg[2])
-	u[3] ^= b & (u[3] ^ neg[3])
-	u[4] ^= b & (u[4] ^ neg[4])
+// CondNeg sets v to -u if cond == 1, and to u if cond == 0.
+func (v *FieldElement) CondNeg(u *FieldElement, cond int) *FieldElement {
+	return v.Select(v.Neg(u), u, cond)
 }
 
-// FeIsNegative returns 1 if u is negative, and 0 otherwise.
-func FeIsNegative(u *FieldElement) int {
+// IsNegative returns 1 if v is negative, and 0 otherwise.
+func (v *FieldElement) IsNegative() int {
 	var b [32]byte
-	FeToBytes(&b, u)
+	v.ToBytes(&b)
 	return int(b[0] & 1)
 }
 
-// FeAbs sets out to |u|. out and u are allowed to overlap.
-func FeAbs(out, u *FieldElement) {
-	var neg FieldElement
-	FeNeg(&neg, u)
-	FeSelect(out, &neg, u, FeIsNegative(u))
+// Abs sets v to |u|. v and u are allowed to overlap.
+func (v *FieldElement) Abs(u *FieldElement) *FieldElement {
+	return v.CondNeg(u, u.IsNegative())
 }
diff --git a/internal/radix51/fe_mul.go b/internal/radix51/fe_mul.go
index edc7bd2..08ffbb4 100644
--- a/internal/radix51/fe_mul.go
+++ b/internal/radix51/fe_mul.go
@@ -6,8 +6,8 @@
 
 package radix51
 
-// FeMul sets out = a * b
-func FeMul(out, x, y *FieldElement) {
+// Mul sets out = x * y.
+func (v *FieldElement) Mul(x, y *FieldElement) *FieldElement {
 	var x0, x1, x2, x3, x4 uint64
 	var y0, y1, y2, y3, y4 uint64
 
@@ -118,9 +118,10 @@ func FeMul(out, x, y *FieldElement) {
 	r00 += (r40 >> 51) * 19
 	r40 &= maskLow51Bits
 
-	out[0] = r00
-	out[1] = r10
-	out[2] = r20
-	out[3] = r30
-	out[4] = r40
+	v[0] = r00
+	v[1] = r10
+	v[2] = r20
+	v[3] = r30
+	v[4] = r40
+	return v
 }
diff --git a/internal/radix51/fe_mul_amd64.go b/internal/radix51/fe_mul_amd64.go
index b0d38be..7ad010f 100644
--- a/internal/radix51/fe_mul_amd64.go
+++ b/internal/radix51/fe_mul_amd64.go
@@ -6,5 +6,11 @@
 
 package radix51
 
+// Mul sets out = x * y.
+func (v *FieldElement) Mul(x, y *FieldElement) *FieldElement {
+	feMul(v, x, y)
+	return v
+}
+
 // go:noescape
-func FeMul(out, a, b *FieldElement)
+func feMul(out, a, b *FieldElement)
diff --git a/internal/radix51/fe_mul_amd64.s b/internal/radix51/fe_mul_amd64.s
index e246262..5a72401 100644
--- a/internal/radix51/fe_mul_amd64.s
+++ b/internal/radix51/fe_mul_amd64.s
@@ -7,8 +7,8 @@
 
 // +build amd64,!noasm
 
-// func FeMul(outp *uint64, xp *uint64, yp *uint64)
-TEXT ·FeMul(SB),$0-24
+// func feMul(outp *uint64, xp *uint64, yp *uint64)
+TEXT ·feMul(SB),$0-24
 	MOVQ outp+0(FP), DI
 	MOVQ xp+8(FP), BX
 	MOVQ yp+16(FP), CX
diff --git a/internal/radix51/fe_square.go b/internal/radix51/fe_square.go
index ec5bd62..0a2d4ac 100644
--- a/internal/radix51/fe_square.go
+++ b/internal/radix51/fe_square.go
@@ -6,8 +6,8 @@
 
 package radix51
 
-// FeSquare sets out = x*x
-func FeSquare(out, x *FieldElement) {
+// Square sets v = x * x.
+func (v *FieldElement) Square(x *FieldElement) *FieldElement {
 	// Squaring needs only 15 mul instructions. Some inputs are multiplied by 2;
 	// this is combined with multiplication by 19 where possible. The coefficient
 	// reduction after squaring is the same as for multiplication.
@@ -90,9 +90,10 @@ func FeSquare(out, x *FieldElement) {
 	r00 += (r40 >> 51) * 19
 	r40 &= maskLow51Bits
 
-	out[0] = r00
-	out[1] = r10
-	out[2] = r20
-	out[3] = r30
-	out[4] = r40
+	v[0] = r00
+	v[1] = r10
+	v[2] = r20
+	v[3] = r30
+	v[4] = r40
+	return v
 }
diff --git a/internal/radix51/fe_square_amd64.go b/internal/radix51/fe_square_amd64.go
index 8c38356..d91d669 100644
--- a/internal/radix51/fe_square_amd64.go
+++ b/internal/radix51/fe_square_amd64.go
@@ -6,5 +6,11 @@
 
 package radix51
 
+// Square sets v = x * x.
+func (v *FieldElement) Square(x *FieldElement) *FieldElement {
+	feSquare(v, x)
+	return v
+}
+
 // go:noescape
-func FeSquare(out, x *FieldElement)
+func feSquare(out, x *FieldElement)
diff --git a/internal/radix51/fe_square_amd64.s b/internal/radix51/fe_square_amd64.s
index fa01331..0103416 100644
--- a/internal/radix51/fe_square_amd64.s
+++ b/internal/radix51/fe_square_amd64.s
@@ -4,8 +4,8 @@
 
 // +build amd64,!noasm
 
-// func FeSquare(outp *uint64, xp *uint64)
-TEXT ·FeSquare(SB),4,$0-16
+// func feSquare(outp *uint64, xp *uint64)
+TEXT ·feSquare(SB),4,$0-16
     MOVQ outp+0(FP), DI
     MOVQ xp+8(FP), SI
 
diff --git a/internal/radix51/fe_test.go b/internal/radix51/fe_test.go
index 3450871..d7343c2 100644
--- a/internal/radix51/fe_test.go
+++ b/internal/radix51/fe_test.go
@@ -73,8 +73,8 @@ func TestFeFromBytesRoundTrip(t *testing.T) {
 	in = [32]byte{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
 		18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32}
 
-	FeFromBytes(&fe, &in)
-	FeToBytes(&out, &fe)
+	fe.FromBytes(&in)
+	fe.ToBytes(&out)
 
 	if !bytes.Equal(in[:], out[:]) {
 		t.Error("Bytes<>FE doesn't roundtrip")
@@ -87,8 +87,8 @@ func TestFeFromBytesRoundTrip(t *testing.T) {
 	fe[3] = 0x5e8fca9e0881c
 	fe[4] = 0x5c490f087d796
 
-	FeToBytes(&out, &fe)
-	FeFromBytes(&r, &out)
+	fe.ToBytes(&out)
+	r.FromBytes(&out)
 
 	for i := 0; i < len(fe); i++ {
 		if r[i] != fe[i] {
@@ -104,9 +104,9 @@ func TestSanity(t *testing.T) {
 	// var x2Go, x2sqGo FieldElement
 
 	x = [5]uint64{1, 1, 1, 1, 1}
-	FeMul(&x2, &x, &x)
+	x2.Mul(&x, &x)
 	// FeMulGo(&x2Go, &x, &x)
-	FeSquare(&x2sq, &x)
+	x2sq.Square(&x)
 	// FeSquareGo(&x2sqGo, &x)
 
 	// if !vartimeEqual(x2, x2Go) || !vartimeEqual(x2sq, x2sqGo) || !vartimeEqual(x2, x2sq) {
@@ -123,11 +123,11 @@ func TestSanity(t *testing.T) {
 	if err != nil {
 		t.Fatal(err)
 	}
-	FeFromBytes(&x, &bytes)
+	x.FromBytes(&bytes)
 
-	FeMul(&x2, &x, &x)
+	x2.Mul(&x, &x)
 	// FeMulGo(&x2Go, &x, &x)
-	FeSquare(&x2sq, &x)
+	x2sq.Square(&x)
 	// FeSquareGo(&x2sqGo, &x)
 
 	// if !vartimeEqual(x2, x2Go) || !vartimeEqual(x2sq, x2sqGo) || !vartimeEqual(x2, x2sq) {
@@ -152,13 +152,13 @@ func TestFeEqual(t *testing.T) {
 	var x FieldElement = [5]uint64{1, 1, 1, 1, 1}
 	var y FieldElement = [5]uint64{5, 4, 3, 2, 1}
 
-	eq := FeEqual(&x, &x)
-	if !eq {
+	eq := x.Equal(&x)
+	if eq != 1 {
 		t.Errorf("wrong about equality")
 	}
 
-	eq = FeEqual(&x, &y)
-	if eq {
+	eq = x.Equal(&y)
+	if eq != 0 {
 		t.Errorf("wrong about inequality")
 	}
 }
@@ -168,9 +168,9 @@ func TestFeInvert(t *testing.T) {
 	var one FieldElement = [5]uint64{1, 0, 0, 0, 0}
 	var xinv, r FieldElement
 
-	FeInvert(&xinv, &x)
-	FeMul(&r, &x, &xinv)
-	FeReduce(&r, &r)
+	xinv.Invert(&x)
+	r.Mul(&x, &xinv)
+	r.Reduce(&r)
 
 	if !vartimeEqual(one, r) {
 		t.Errorf("inversion identity failed, got: %x", r)
@@ -182,11 +182,11 @@ func TestFeInvert(t *testing.T) {
 	if err != nil {
 		t.Fatal(err)
 	}
-	FeFromBytes(&x, &bytes)
+	x.FromBytes(&bytes)
 
-	FeInvert(&xinv, &x)
-	FeMul(&r, &x, &xinv)
-	FeReduce(&r, &r)
+	xinv.Invert(&x)
+	r.Mul(&x, &xinv)
+	r.Reduce(&r)
 
 	if !vartimeEqual(one, r) {
 		t.Errorf("random inversion identity failed, got: %x for field element %x", r, x)
diff --git a/ristretto255.go b/ristretto255.go
index eff2c87..097b277 100644
--- a/ristretto255.go
+++ b/ristretto255.go
@@ -35,13 +35,13 @@ type Element struct {
 func (e *Element) Equal(ee *Element) int {
 	var f0, f1 radix51.FieldElement
 
-	radix51.FeMul(&f0, &e.r.X, &ee.r.Y) // x1 * y2
-	radix51.FeMul(&f1, &e.r.Y, &ee.r.X) // y1 * x2
-	out := radix51.FeEqual(&f0, &f1)
+	f0.Mul(&e.r.X, &ee.r.Y) // x1 * y2
+	f1.Mul(&e.r.Y, &ee.r.X) // y1 * x2
+	out := f0.Equal(&f1)
 
-	radix51.FeMul(&f0, &e.r.Y, &ee.r.Y) // y1 * y2
-	radix51.FeMul(&f1, &e.r.X, &ee.r.X) // x1 * x2
-	out = out | radix51.FeEqual(&f0, &f1)
+	f0.Mul(&e.r.Y, &ee.r.Y) // y1 * y2
+	f1.Mul(&e.r.X, &ee.r.X) // x1 * x2
+	out = out | f0.Equal(&f1)
 
 	return out
 }
@@ -58,71 +58,75 @@ func (e *Element) FromUniformBytes(b []byte) {
 	f := &radix51.FieldElement{}
 
 	copy(buf[:], b[:32])
-	radix51.FeFromBytes(f, &buf)
+	f.FromBytes(&buf)
 	p1 := &group.ExtendedGroupElement{}
 	mapToPoint(p1, f)
 
 	copy(buf[:], b[32:])
-	radix51.FeFromBytes(f, &buf)
+	f.FromBytes(&buf)
 	p2 := &group.ExtendedGroupElement{}
 	mapToPoint(p2, f)
 
 	e.r.Add(p1, p2)
 }
 
+// mapToPoint implements MAP from Section 3.2.4 of draft-hdevalence-cfrg-ristretto-00.
 func mapToPoint(out *group.ExtendedGroupElement, t *radix51.FieldElement) {
+	// r = SQRT_M1 * t^2
 	r := &radix51.FieldElement{}
-	radix51.FeSquare(r, t)
-	radix51.FeMul(r, r, sqrtM1)
-
-	one := &radix51.FieldElement{}
-	radix51.FeOne(one)
-	minusOne := &radix51.FieldElement{}
-	radix51.FeNeg(minusOne, one)
+	r.Mul(sqrtM1, r.Square(t))
 
+	// u = (r + 1) * ONE_MINUS_D_SQ
 	u := &radix51.FieldElement{}
-	radix51.FeAdd(u, r, one)
-	radix51.FeMul(u, u, oneMinusDSQ)
+	u.Mul(u.Add(r, radix51.One), oneMinusDSQ)
 
+	// c = -1
+	c := &radix51.FieldElement{}
+	c.Set(radix51.MinusOne)
+
+	// v = (c - r*D) * (r + D)
 	rPlusD := &radix51.FieldElement{}
-	radix51.FeAdd(rPlusD, r, &group.D)
+	rPlusD.Add(r, group.D)
 	v := &radix51.FieldElement{}
-	radix51.FeMul(v, r, &group.D)
-	radix51.FeSub(v, minusOne, v)
-	radix51.FeMul(v, v, rPlusD)
+	v.Mul(v.Sub(c, v.Mul(r, group.D)), rPlusD)
 
+	// (was_square, s) = SQRT_RATIO_M1(u, v)
 	s := &radix51.FieldElement{}
 	wasSquare := feSqrtRatio(s, u, v)
+
+	// s_prime = -CT_ABS(s*t)
 	sPrime := &radix51.FieldElement{}
-	radix51.FeMul(sPrime, s, t)
-	radix51.FeAbs(sPrime, sPrime)
-	radix51.FeNeg(sPrime, sPrime)
+	sPrime.Neg(sPrime.Abs(sPrime.Mul(s, t)))
 
-	c := &radix51.FieldElement{}
-	radix51.FeSelect(s, s, sPrime, wasSquare)
-	radix51.FeSelect(c, minusOne, r, wasSquare)
+	// s = CT_SELECT(s IF was_square ELSE s_prime)
+	s.Select(s, sPrime, wasSquare)
+	// c = CT_SELECT(c IF was_square ELSE r)
+	c.Select(c, r, wasSquare)
 
+	// N = c * (r - 1) * D_MINUS_ONE_SQ - v
 	N := &radix51.FieldElement{}
-	radix51.FeSub(N, r, one)
-	radix51.FeMul(N, N, c)
-	radix51.FeMul(N, N, dMinusOneSQ)
-	radix51.FeSub(N, N, v)
+	N.Mul(c, N.Sub(r, radix51.One))
+	N.Sub(N.Mul(N, dMinusOneSQ), v)
 
-	sSquare := &radix51.FieldElement{}
-	radix51.FeSquare(sSquare, s)
+	s2 := &radix51.FieldElement{}
+	s2.Square(s)
 
+	// w0 = 2 * s * v
 	w0 := &radix51.FieldElement{}
-	radix51.FeMul(w0, s, v)
-	radix51.FeAdd(w0, w0, w0)
+	w0.Add(w0, w0.Mul(s, v))
+	// w1 = N * SQRT_AD_MINUS_ONE
 	w1 := &radix51.FieldElement{}
-	radix51.FeMul(w1, N, sqrtADMinusOne)
+	w1.Mul(N, sqrtADMinusOne)
+	// w2 = 1 - s^2
 	w2 := &radix51.FieldElement{}
-	radix51.FeSub(w2, one, sSquare)
+	w2.Sub(radix51.One, s2)
+	// w3 = 1 + s^2
 	w3 := &radix51.FieldElement{}
-	radix51.FeAdd(w3, one, sSquare)
+	w3.Add(radix51.One, s2)
 
-	radix51.FeMul(&out.X, w0, w3)
-	radix51.FeMul(&out.Y, w2, w1)
-	radix51.FeMul(&out.Z, w1, w3)
-	radix51.FeMul(&out.T, w0, w2)
+	// return (w0*w3, w2*w1, w1*w3, w0*w2)
+	out.X.Mul(w0, w3)
+	out.Y.Mul(w2, w1)
+	out.Z.Mul(w1, w3)
+	out.T.Mul(w0, w2)
 }