Including some feedback received from PR #84.

ecc/goldilocks/point.go

@@ -29,9 +29,9 @@ func FromAffine(x, y *fp.Elt) (*Point, error) {
 	return P, nil
 }
 
-// isLessThan returns true if 0 <= x < y, and assumes that slices have the same length.
+// isLessThan returns true if 0 <= x < y, and assumes that slices are of the
+// same length and are interpreted in little-endian order.
 func isLessThan(x, y []byte) bool {
-
 	i := len(x) - 1
 	for i > 0 && x[i] == y[i] {
 		i--
@@ -142,6 +142,9 @@ func (P *Point) Double() { P.Add(P) }
 
 // Add sets P =P+Q..
 func (P *Point) Add(Q *Point) {
+	// This is formula (5) from "Twisted Edwards Curves Revisited" by
+	// Hisil H., Wong K.KH., Carter G., Dawson E. (2008)
+	// https://doi.org/10.1007/978-3-540-89255-7_20
 	x1, y1, z1, ta1, tb1 := &P.x, &P.y, &P.z, &P.ta, &P.tb
 	x2, y2, z2, ta2, tb2 := &Q.x, &Q.y, &Q.z, &Q.ta, &Q.tb
 	x3, y3, z3, E, H := &P.x, &P.y, &P.z, &P.ta, &P.tb

ecc/goldilocks/twistPoint.go

@@ -30,6 +30,9 @@ func (P *twistPoint) cneg(b uint) {
 
 // Double updates P with 2P.
 func (P *twistPoint) Double() {
+	// This is formula (7) from "Twisted Edwards Curves Revisited" by
+	// Hisil H., Wong K.KH., Carter G., Dawson E. (2008)
+	// https://doi.org/10.1007/978-3-540-89255-7_20
 	Px, Py, Pz, Pta, Ptb := &P.x, &P.y, &P.z, &P.ta, &P.tb
 	a, b, c, e, f, g, h := Px, Py, Pz, Pta, Px, Py, Ptb
 	fp.Add(e, Px, Py) // x+y
@@ -55,6 +58,9 @@ func (P *twistPoint) mixAddZ1(Q *preTwistPointAffine) {
 
 // coreAddition calculates P=P+Q for curves with A=-1
 func (P *twistPoint) coreAddition(Q *preTwistPointAffine) {
+	// This is the formula following (5) from "Twisted Edwards Curves Revisited" by
+	// Hisil H., Wong K.KH., Carter G., Dawson E. (2008)
+	// https://doi.org/10.1007/978-3-540-89255-7_20
 	Px, Py, Pz, Pta, Ptb := &P.x, &P.y, &P.z, &P.ta, &P.tb
 	addYX2, subYX2, dt2 := &Q.addYX, &Q.subYX, &Q.dt2
 	a, b, c, d, e, f, g, h := Px, Py, &fp.Elt{}, Pz, Pta, Px, Py, Ptb
@@ -120,10 +126,10 @@ func (P *preTwistPointProy) cmov(Q *preTwistPointProy, b uint) {
 
 // FromTwistPoint precomputes some coordinates of Q for mised addition.
 func (P *preTwistPointProy) FromTwistPoint(Q *twistPoint) {
-	fp.Add(&P.addYX, &Q.y, &Q.x)
-	fp.Sub(&P.subYX, &Q.y, &Q.x)
-	fp.Mul(&P.dt2, &Q.ta, &Q.tb)
-	fp.Mul(&P.dt2, &P.dt2, &paramDTwist)
-	fp.Add(&P.dt2, &P.dt2, &P.dt2)
-	fp.Add(&P.z2, &Q.z, &Q.z)
+	fp.Add(&P.addYX, &Q.y, &Q.x)         // addYX = X + Y
+	fp.Sub(&P.subYX, &Q.y, &Q.x)         // subYX = Y - X
+	fp.Mul(&P.dt2, &Q.ta, &Q.tb)         // T = ta*tb
+	fp.Mul(&P.dt2, &P.dt2, &paramDTwist) // D*T
+	fp.Add(&P.dt2, &P.dt2, &P.dt2)       // dt2 = 2*D*T
+	fp.Add(&P.z2, &Q.z, &Q.z)            // z2 = 2*Z
 }

math/fp448/fp.go

@@ -57,10 +57,15 @@ func Neg(z, x *Elt) { Sub(z, &p, x) }
 // Modp ensures that z is between [0,p-1].
 func Modp(z *Elt) { Sub(z, z, &p) }
 
-// InvSqrt calculates z = sqrt(x/y) iff x/y is a quadratic-residue, which is
-// indicated by returning isQR = true. Otherwise, when x/y is a quadratic
-// non-residue, z will have an undetermined value and isQR = false.
+// InvSqrt calculates z = sqrt(x/y) iff x/y is a quadratic-residue. If so,
+// isQR = true; otherwise, isQR = false, since x/y is a quadratic non-residue,
+// and z = sqrt(-x/y).
 func InvSqrt(z, x, y *Elt) (isQR bool) {
+	// First note that x^(2(k+1)) = x^(p-1)/2 * x = legendre(x) * x
+	// so that's x if x is a quadratic residue and -x otherwise.
+	// Next, y^(6k+3) = y^(4k+2) * y^(2k+1) = y^(p-1) * y^((p-1)/2) = legendre(y).
+	// So the z we compute satisfies z^2 y = x^(2(k+1)) y^(6k+3) = legendre(x)*legendre(y).
+	// Thus if x and y are quadratic residues, then z is indeed sqrt(x/y).
 	t0, t1 := &Elt{}, &Elt{}
 	Mul(t0, x, y)         // x*y
 	Sqr(t1, y)            // y^2
@@ -77,6 +82,7 @@ func InvSqrt(z, x, y *Elt) (isQR bool) {
 
 // Inv calculates z = 1/x mod p.
 func Inv(z, x *Elt) {
+	// Calculates z = x^(4k+1) = x^(p-3+1) = x^(p-2) = x^-1, where k = (p-3)/4.
 	t := &Elt{}
 	powPminus3div4(t, x) // t = x^k
 	Sqr(t, t)            // t = x^2k