Import gtank/ed25519#8 and refactor on top of it

This commit is contained in:
Filippo Valsorda 2019-01-26 22:20:45 -05:00
parent a3540ec35a
commit b0a75c0ab7
13 changed files with 420 additions and 462 deletions

162
fe.go
View File

@ -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)
}

View File

@ -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}
)

View File

@ -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 ← GC
field.FeSub(&F, &G, &C)
// H ← DB
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 ← GC
H.Sub(&D, &B) // H ← DB
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
}

View File

@ -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}
)

View File

@ -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())
}

View File

@ -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
}

View File

@ -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)

View File

@ -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

View File

@ -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
}

View File

@ -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)

View File

@ -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

View File

@ -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)

View File

@ -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)
}