Makes SIDH independent on underlying field arithmetic.

Makefile

@@ -5,7 +5,7 @@ GOPATH_LOCAL = $(PRJ_DIR)/build
 GOPATH_DIR   = github.com/cloudflare/p751sidh
 CSHAKE_PKG   ?= github.com/henrydcase/nobs/hash/sha3
 CPU_PKG      = golang.org/x/sys/cpu
-TARGETS      = p751toolbox sidh sike
+TARGETS      = p751 sidh sike
 GO           ?= go
 GOARCH       ?=
 OPTS_GCCGO   ?= -compiler gccgo -O2 -g
@@ -26,6 +26,7 @@ ifeq ($(V),1)
 	BENCH_OPTS += -gcflags=-m     # Show results from inlining
 endif
 
+all: test
 clean:
 	rm -rf $(GOPATH_LOCAL)
 	rm -rf coverage*.txt
@@ -33,6 +34,7 @@ clean:
 build_env:
 	GOPATH=$(GOPATH_LOCAL) $(GO) get $(CSHAKE_PKG) $(CPU_PKG)
 	mkdir -p $(GOPATH_LOCAL)/src/$(GOPATH_DIR)
+	cp -rf internal $(GOPATH_LOCAL)/src/$(GOPATH_DIR)
 	cp -rf etc $(GOPATH_LOCAL)/src/$(GOPATH_DIR)
 
 copy-target-%:
@@ -48,7 +50,7 @@ test-%: prep_targets
 	GOPATH=$(GOPATH_LOCAL) $(GO) test $(OPTS) $(GOPATH_DIR)/$*
 
 bench-%: prep_targets
-	cd $*; GOPATH=$(GOPATH_LOCAL) $(GO) test $(OPTS) $(BENCH_OPTS)
+	GOPATH=$(GOPATH_LOCAL) $(GO) test $(OPTS) $(GOPATH_DIR)/$* $(BENCH_OPTS)
 
 cover-%: prep_targets
 	GOPATH=$(GOPATH_LOCAL) $(GO) test \

internal/isogeny/curve_ops.go

@@ -0,0 +1,440 @@
+package internal
+
+type CurveOperations struct {
+    Params *SidhParams
+}
+
+// Computes j-invariant for a curve y2=x3+A/Cx+x with A,C in F_(p^2). Result
+// is returned in jBytes buffer, encoded in little-endian format. Caller
+// provided jBytes buffer has to be big enough to j-invariant value. In case
+// of SIDH, buffer size must be at least size of shared secret.
+// Implementation corresponds to Algorithm 9 from SIKE.
+func (c *CurveOperations) Jinvariant(cparams *ProjectiveCurveParameters, jBytes []byte) {
+    var j, t0, t1 Fp2Element
+
+    op := c.Params.Op
+    op.Square(&j, &cparams.A)  // j  = A^2
+    op.Square(&t1, &cparams.C) // t1 = C^2
+    op.Add(&t0, &t1, &t1)      // t0 = t1 + t1
+    op.Sub(&t0, &j, &t0)       // t0 = j - t0
+    op.Sub(&t0, &t0, &t1)      // t0 = t0 - t1
+    op.Sub(&j, &t0, &t1)       // t0 = t0 - t1
+    op.Square(&t1, &t1)        // t1 = t1^2
+    op.Mul(&j, &j, &t1)        // j = j * t1
+    op.Add(&t0, &t0, &t0)      // t0 = t0 + t0
+    op.Add(&t0, &t0, &t0)      // t0 = t0 + t0
+    op.Square(&t1, &t0)        // t1 = t0^2
+    op.Mul(&t0, &t0, &t1)      // t0 = t0 * t1
+    op.Add(&t0, &t0, &t0)      // t0 = t0 + t0
+    op.Add(&t0, &t0, &t0)      // t0 = t0 + t0
+    op.Inv(&j, &j)             // j  = 1/j
+    op.Mul(&j, &t0, &j)        // j  = t0 * j
+
+    c.Fp2ToBytes(jBytes, &j)
+}
+
+// Given affine points x(P), x(Q) and x(Q-P) in a extension field F_{p^2}, function
+// recorvers projective coordinate A of a curve. This is Algorithm 10 from SIKE.
+func (c *CurveOperations) RecoverCoordinateA(curve *ProjectiveCurveParameters, xp, xq, xr *Fp2Element) {
+    var t0, t1 Fp2Element
+
+    op := c.Params.Op
+    op.Add(&t1, xp, xq)                            // t1 = Xp + Xq
+    op.Mul(&t0, xp, xq)                            // t0 = Xp * Xq
+    op.Mul(&curve.A, xr, &t1)                      // A  = X(q-p) * t1
+    op.Add(&curve.A, &curve.A, &t0)                // A  = A + t0
+    op.Mul(&t0, &t0, xr)                           // t0 = t0 * X(q-p)
+    op.Sub(&curve.A, &curve.A, &c.Params.OneFp2)   // A  = A - 1
+    op.Add(&t0, &t0, &t0)                          // t0 = t0 + t0
+    op.Add(&t1, &t1, xr)                           // t1 = t1 + X(q-p)
+    op.Add(&t0, &t0, &t0)                          // t0 = t0 + t0
+    op.Square(&curve.A, &curve.A)                  // A  = A^2
+    op.Inv(&t0, &t0)                               // t0 = 1/t0
+    op.Mul(&curve.A, &curve.A, &t0)                // A  = A * t0
+    op.Sub(&curve.A, &curve.A, &t1)                // A  = A - t1
+}
+
+// Computes equivalence (A:C) ~ (A+2C : A-2C)
+func (c *CurveOperations) CalcCurveParamsEquiv3(cparams *ProjectiveCurveParameters) CurveCoefficientsEquiv {
+    var coef CurveCoefficientsEquiv
+    var c2 Fp2Element
+    var op = c.Params.Op
+
+    op.Add(&c2, &cparams.C, &cparams.C)
+    // A24p = A+2*C
+    op.Add(&coef.A, &cparams.A, &c2)
+    // A24m = A-2*C
+    op.Sub(&coef.C, &cparams.A, &c2)
+    return coef
+}
+
+// Computes equivalence (A:C) ~ (A+2C : 4C)
+func (c *CurveOperations) CalcCurveParamsEquiv4(cparams *ProjectiveCurveParameters) CurveCoefficientsEquiv {
+    var coefEq CurveCoefficientsEquiv
+    var op = c.Params.Op
+
+    op.Add(&coefEq.C, &cparams.C, &cparams.C)
+    // A24p = A+2C
+    op.Add(&coefEq.A, &cparams.A, &coefEq.C)
+    // C24 = 4*C
+    op.Add(&coefEq.C, &coefEq.C, &coefEq.C)
+    return coefEq
+}
+
+// Helper function for RightToLeftLadder(). Returns A+2C / 4.
+func (c *CurveOperations) CalcAplus2Over4(cparams *ProjectiveCurveParameters) (ret Fp2Element) {
+    var tmp Fp2Element
+    var op = c.Params.Op
+
+    // 2C
+    op.Add(&tmp, &cparams.C, &cparams.C)
+    // A+2C
+    op.Add(&ret, &cparams.A, &tmp)
+    // 1/4C
+    op.Add(&tmp, &tmp, &tmp)
+    op.Inv(&tmp, &tmp)
+    // A+2C/4C
+    op.Mul(&ret, &ret, &tmp)
+    return
+}
+
+// Recovers (A:C) curve parameters from projectively equivalent (A+2C:A-2C).
+func (c *CurveOperations) RecoverCurveCoefficients3(cparams *ProjectiveCurveParameters, coefEq *CurveCoefficientsEquiv) {
+    var op = c.Params.Op
+
+    op.Add(&cparams.A, &coefEq.A, &coefEq.C)
+    // cparams.A = 2*(A+2C+A-2C) = 4A
+    op.Add(&cparams.A, &cparams.A, &cparams.A)
+    // cparams.C = (A+2C-A+2C) = 4C
+    op.Sub(&cparams.C, &coefEq.A, &coefEq.C)
+    return
+}
+
+// Recovers (A:C) curve parameters from projectively equivalent (A+2C:4C).
+func (c *CurveOperations) RecoverCurveCoefficients4(cparams *ProjectiveCurveParameters, coefEq *CurveCoefficientsEquiv) {
+    var op = c.Params.Op
+    // cparams.C = (4C)*1/2=2C
+    op.Mul(&cparams.C, &coefEq.C, &c.Params.HalfFp2)
+    // cparams.A = A+2C - 2C = A
+    op.Sub(&cparams.A, &coefEq.A, &cparams.C)
+    // cparams.C = 2C * 1/2 = C
+    op.Mul(&cparams.C, &cparams.C, &c.Params.HalfFp2)
+    return
+}
+
+// Combined coordinate doubling and differential addition. Takes projective points
+// P,Q,Q-P and (A+2C)/4C curve E coefficient. Returns 2*P and P+Q calculated on E.
+// Function is used only by RightToLeftLadder. Corresponds to Algorithm 5 of SIKE
+func (c *CurveOperations) xDblAdd(P, Q, QmP *ProjectivePoint, a24 *Fp2Element) (dblP, PaQ ProjectivePoint) {
+    var t0, t1, t2 Fp2Element
+    var op = c.Params.Op
+
+    xQmP, zQmP := &QmP.X, &QmP.Z
+    xPaQ, zPaQ := &PaQ.X, &PaQ.Z
+    x2P, z2P := &dblP.X, &dblP.Z
+    xP, zP := &P.X, &P.Z
+    xQ, zQ := &Q.X, &Q.Z
+
+    op.Add(&t0, xP, zP)      // t0   = Xp+Zp
+    op.Sub(&t1, xP, zP)      // t1   = Xp-Zp
+    op.Square(x2P, &t0)      // 2P.X = t0^2
+    op.Sub(&t2, xQ, zQ)      // t2   = Xq-Zq
+    op.Add(xPaQ, xQ, zQ)     // Xp+q = Xq+Zq
+    op.Mul(&t0, &t0, &t2)    // t0   = t0 * t2
+    op.Mul(z2P, &t1, &t1)    // 2P.Z = t1 * t1
+    op.Mul(&t1, &t1, xPaQ)   // t1   = t1 * Xp+q
+    op.Sub(&t2, x2P, z2P)    // t2   = 2P.X - 2P.Z
+    op.Mul(x2P, x2P, z2P)    // 2P.X = 2P.X * 2P.Z
+    op.Mul(xPaQ, a24, &t2)   // Xp+q = A24 * t2
+    op.Sub(zPaQ, &t0, &t1)   // Zp+q = t0 - t1
+    op.Add(z2P, xPaQ, z2P)   // 2P.Z = Xp+q + 2P.Z
+    op.Add(xPaQ, &t0, &t1)   // Xp+q = t0 + t1
+    op.Mul(z2P, z2P, &t2)    // 2P.Z = 2P.Z * t2
+    op.Square(zPaQ, zPaQ)    // Zp+q = Zp+q ^ 2
+    op.Square(xPaQ, xPaQ)    // Xp+q = Xp+q ^ 2
+    op.Mul(zPaQ, xQmP, zPaQ) // Zp+q = Xq-p * Zp+q
+    op.Mul(xPaQ, zQmP, xPaQ) // Xp+q = Zq-p * Xp+q
+    return
+}
+
+// Given the curve parameters, xP = x(P), computes xP = x([2^k]P)
+// Safe to overlap xP, x2P.
+func (c *CurveOperations) Pow2k(xP *ProjectivePoint, params *CurveCoefficientsEquiv, k uint32) {
+    var t0, t1 Fp2Element
+    var op = c.Params.Op
+
+    x, z := &xP.X, &xP.Z
+    for i := uint32(0); i < k; i++ {
+        op.Sub(&t0, x, z)           // t0  = Xp - Zp
+        op.Add(&t1, x, z)           // t1  = Xp + Zp
+        op.Square(&t0, &t0)         // t0  = t0 ^ 2
+        op.Square(&t1, &t1)         // t1  = t1 ^ 2
+        op.Mul(z, &params.C, &t0)   // Z2p = C24 * t0
+        op.Mul(x, z, &t1)           // X2p = Z2p * t1
+        op.Sub(&t1, &t1, &t0)       // t1  = t1 - t0
+        op.Mul(&t0, &params.A, &t1) // t0  = A24+ * t1
+        op.Add(z, z, &t0)           // Z2p = Z2p + t0
+        op.Mul(z, z, &t1)           // Zp  = Z2p * t1
+    }
+}
+
+// Given the curve parameters, xP = x(P), and k >= 0, compute xP = x([3^k]P).
+//
+// Safe to overlap xP, xR.
+func (c *CurveOperations) Pow3k(xP *ProjectivePoint, params *CurveCoefficientsEquiv, k uint32) {
+    var t0, t1, t2, t3, t4, t5, t6 Fp2Element
+    var op = c.Params.Op
+
+    x, z := &xP.X, &xP.Z
+    for i := uint32(0); i < k; i++ {
+        op.Sub(&t0, x, z)           // t0  = Xp - Zp
+        op.Square(&t2, &t0)         // t2  = t0^2
+        op.Add(&t1, x, z)           // t1  = Xp + Zp
+        op.Square(&t3, &t1)         // t3  = t1^2
+        op.Add(&t4, &t1, &t0)       // t4  = t1 + t0
+        op.Sub(&t0, &t1, &t0)       // t0  = t1 - t0
+        op.Square(&t1, &t4)         // t1  = t4^2
+        op.Sub(&t1, &t1, &t3)       // t1  = t1 - t3
+        op.Sub(&t1, &t1, &t2)       // t1  = t1 - t2
+        op.Mul(&t5, &t3, &params.A) // t5  = t3 * A24+
+        op.Mul(&t3, &t3, &t5)       // t3  = t5 * t3
+        op.Mul(&t6, &t2, &params.C) // t6  = t2 * A24-
+        op.Mul(&t2, &t2, &t6)       // t2  = t2 * t6
+        op.Sub(&t3, &t2, &t3)       // t3  = t2 - t3
+        op.Sub(&t2, &t5, &t6)       // t2  = t5 - t6
+        op.Mul(&t1, &t2, &t1)       // t1  = t2 * t1
+        op.Add(&t2, &t3, &t1)       // t2  = t3 + t1
+        op.Square(&t2, &t2)         // t2  = t2^2
+        op.Mul(x, &t2, &t4)        // X3p = t2 * t4
+        op.Sub(&t1, &t3, &t1)       // t1  = t3 - t1
+        op.Square(&t1, &t1)         // t1  = t1^2
+        op.Mul(z, &t1, &t0)        // Z3p = t1 * t0
+    }
+}
+
+// Set (y1, y2, y3)  = (1/x1, 1/x2, 1/x3).
+//
+// All xi, yi must be distinct.
+func (c *CurveOperations) Fp2Batch3Inv(x1, x2, x3, y1, y2, y3 *Fp2Element) {
+    var x1x2, t Fp2Element
+    var op = c.Params.Op
+
+    op.Mul(&x1x2, x1, x2)     // x1*x2
+    op.Mul(&t, &x1x2, x3)     // 1/(x1*x2*x3)
+    op.Inv(&t, &t)
+    op.Mul(y1, &t, x2)        // 1/x1
+    op.Mul(y1, y1, x3)
+    op.Mul(y2, &t, x1)        // 1/x2
+    op.Mul(y2, y2, x3)
+    op.Mul(y3, &t, &x1x2)     // 1/x3
+}
+
+// ScalarMul3Pt is a right-to-left point multiplication that given the
+// x-coordinate of P, Q and P-Q calculates the x-coordinate of R=Q+[scalar]P.
+// nbits must be smaller or equal to len(scalar).
+func (c *CurveOperations) ScalarMul3Pt(cparams *ProjectiveCurveParameters, P, Q, PmQ *ProjectivePoint, nbits uint, scalar []uint8) ProjectivePoint {
+    var R0, R2, R1 ProjectivePoint
+    var op = c.Params.Op
+    aPlus2Over4 := c.CalcAplus2Over4(cparams)
+    R1 = *P
+    R2 = *PmQ
+    R0 = *Q
+
+    // Iterate over the bits of the scalar, bottom to top
+    prevBit := uint8(0)
+    for i := uint(0); i < nbits; i++ {
+        bit := (scalar[i>>3] >> (i & 7) & 1)
+        swap := prevBit ^ bit
+        prevBit = bit
+        op.CondSwap(&R1.X, &R1.Z, &R2.X, &R2.Z, swap)
+        R0, R2 = c.xDblAdd(&R0, &R2, &R1, &aPlus2Over4)
+    }
+    op.CondSwap(&R1.X, &R1.Z, &R2.X, &R2.Z, prevBit)
+    return R1
+}
+
+// Convert the input to wire format.
+//
+// The output byte slice must be at least 2*bytelen(p) bytes long.
+func (c *CurveOperations) Fp2ToBytes(output []byte, fp2 *Fp2Element) {
+    if len(output) < 2*c.Params.Bytelen {
+        panic("output byte slice too short")
+    }
+    var a Fp2Element
+    c.Params.Op.FromMontgomery(fp2,&a)
+
+    // convert to bytes in little endian form
+    for i := 0; i < c.Params.Bytelen; i++ {
+        // set i = j*8 + k
+        fp2 := i / 8
+        k := uint64(i % 8)
+        output[i] = byte(a.A[fp2] >> (8 * k))
+        output[i+c.Params.Bytelen] = byte(a.B[fp2] >> (8 * k))
+    }
+}
+
+// Read 2*bytelen(p) bytes into the given ExtensionFieldElement.
+//
+// It is an error to call this function if the input byte slice is less than 2*bytelen(p) bytes long.
+func (c *CurveOperations) Fp2FromBytes(fp2 *Fp2Element, input []byte) {
+    if len(input) < 2*c.Params.Bytelen {
+        panic("input byte slice too short")
+    }
+
+    for i:=0; i<c.Params.Bytelen; i++ {
+        j := i / 8
+        k := uint64(i % 8)
+        fp2.A[j] |= uint64(input[i]) << (8 * k)
+        fp2.B[j] |= uint64(input[i+c.Params.Bytelen]) << (8 * k)
+    }
+    c.Params.Op.ToMontgomery(fp2)
+}
+
+/* -------------------------------------------------------------------------
+   Mechnisms used for isogeny calculations
+   -------------------------------------------------------------------------*/
+
+// Constructs isogeny3 objects
+func Newisogeny3(op FieldOps) Isogeny {
+    return &isogeny3{Field: op}
+}
+
+// Constructs isogeny4 objects
+func Newisogeny4(op FieldOps) Isogeny {
+    return &isogeny4{isogeny3: isogeny3{Field: op}}
+}
+
+// Given a three-torsion point p = x(PB) on the curve E_(A:C), construct the
+// three-isogeny phi : E_(A:C) -> E_(A:C)/<P_3> = E_(A':C').
+//
+// Input: (XP_3: ZP_3), where P_3 has exact order 3 on E_A/C
+// Output: * Curve coordinates (A' + 2C', A' - 2C') corresponding to E_A'/C' = A_E/C/<P3>
+//         * Isogeny phi with constants in F_p^2
+func (phi *isogeny3) GenerateCurve(p *ProjectivePoint) CurveCoefficientsEquiv {
+    var t0, t1, t2, t3, t4 Fp2Element
+    var coefEq CurveCoefficientsEquiv
+    var K1, K2 = &phi.K1, &phi.K2
+
+    op := phi.Field
+    op.Sub(K1, &p.X, &p.Z)           // K1 = XP3 - ZP3
+    op.Square(&t0, K1)                // t0 = K1^2
+    op.Add(K2, &p.X, &p.Z)           // K2 = XP3 + ZP3
+    op.Square(&t1, K2)                // t1 = K2^2
+    op.Add(&t2, &t0, &t1)             // t2 = t0 + t1
+    op.Add(&t3, K1, K2)               // t3 = K1 + K2
+    op.Square(&t3, &t3)               // t3 = t3^2
+    op.Sub(&t3, &t3, &t2)             // t3 = t3 - t2
+    op.Add(&t2, &t1, &t3)             // t2 = t1 + t3
+    op.Add(&t3, &t3, &t0)             // t3 = t3 + t0
+    op.Add(&t4, &t3, &t0)             // t4 = t3 + t0
+    op.Add(&t4, &t4, &t4)             // t4 = t4 + t4
+    op.Add(&t4, &t1, &t4)             // t4 = t1 + t4
+    op.Mul(&coefEq.C, &t2, &t4)       // A24m = t2 * t4
+    op.Add(&t4, &t1, &t2)             // t4 = t1 + t2
+    op.Add(&t4, &t4, &t4)             // t4 = t4 + t4
+    op.Add(&t4, &t0, &t4)             // t4 = t0 + t4
+    op.Mul(&t4, &t3, &t4)             // t4 = t3 * t4
+    op.Sub(&t0, &t4, &coefEq.C)       // t0 = t4 - A24m
+    op.Add(&coefEq.A, &coefEq.C, &t0) // A24p = A24m + t0
+    return coefEq
+}
+
+// Given a 3-isogeny phi and a point pB = x(PB), compute x(QB), the x-coordinate
+// of the image QB = phi(PB) of PB under phi : E_(A:C) -> E_(A':C').
+//
+// The output xQ = x(Q) is then a point on the curve E_(A':C'); the curve
+// parameters are returned by the GenerateCurve function used to construct phi.
+func (phi *isogeny3) EvaluatePoint(p *ProjectivePoint) ProjectivePoint {
+    var t0, t1, t2 Fp2Element
+    var q ProjectivePoint
+    var K1, K2 = &phi.K1, &phi.K2
+    var px, pz = &p.X, &p.Z
+
+    op := phi.Field
+    op.Add(&t0, px, pz)   // t0 = XQ + ZQ
+    op.Sub(&t1, px, pz)   // t1 = XQ - ZQ
+    op.Mul(&t0, K1, &t0)  // t2 = K1 * t0
+    op.Mul(&t1, K2, &t1)  // t1 = K2 * t1
+    op.Add(&t2, &t0, &t1) // t2 = t0 + t1
+    op.Sub(&t0, &t1, &t0) // t0 = t1 - t0
+    op.Square(&t2, &t2)   // t2 = t2 ^ 2
+    op.Square(&t0, &t0)   // t0 = t0 ^ 2
+    op.Mul(&q.X, px, &t2) // XQ'= XQ * t2
+    op.Mul(&q.Z, pz, &t0) // ZQ'= ZQ * t0
+    return q
+}
+
+// Given a four-torsion point p = x(PB) on the curve E_(A:C), construct the
+// four-isogeny phi : E_(A:C) -> E_(A:C)/<P_4> = E_(A':C').
+//
+// Input: (XP_4: ZP_4), where P_4 has exact order 4 on E_A/C
+// Output: * Curve coordinates (A' + 2C', 4C') corresponding to E_A'/C' = A_E/C/<P4>
+//         * Isogeny phi with constants in F_p^2
+func (phi *isogeny4) GenerateCurve(p *ProjectivePoint) CurveCoefficientsEquiv {
+    var coefEq CurveCoefficientsEquiv
+    var xp4, zp4 = &p.X, &p.Z
+    var K1, K2, K3 = &phi.K1, &phi.K2, &phi.K3
+
+    op := phi.Field
+    op.Sub(K2, xp4, zp4)
+    op.Add(K3, xp4, zp4)
+    op.Square(K1, zp4)
+    op.Add(K1, K1, K1)
+    op.Square(&coefEq.C, K1)
+    op.Add(K1, K1, K1)
+    op.Square(&coefEq.A, xp4)
+    op.Add(&coefEq.A, &coefEq.A, &coefEq.A)
+    op.Square(&coefEq.A, &coefEq.A)
+    return coefEq
+}
+
+// Given a 4-isogeny phi and a point xP = x(P), compute x(Q), the x-coordinate
+// of the image Q = phi(P) of P under phi : E_(A:C) -> E_(A':C').
+//
+// Input: Isogeny returned by GenerateCurve and point q=(Qx,Qz) from E0_A/C
+// Output: Corresponding point q from E1_A'/C', where E1 is 4-isogenous to E0
+func (phi *isogeny4) EvaluatePoint(p *ProjectivePoint) ProjectivePoint {
+    var t0, t1 Fp2Element
+    var q = *p
+    var xq, zq = &q.X, &q.Z
+    var K1, K2, K3 = &phi.K1, &phi.K2, &phi.K3
+
+    op := phi.Field
+    op.Add(&t0, xq, zq)
+    op.Sub(&t1, xq, zq)
+    op.Mul(xq, &t0, K2)
+    op.Mul(zq, &t1, K3)
+    op.Mul(&t0, &t0, &t1)
+    op.Mul(&t0, &t0, K1)
+    op.Add(&t1, xq, zq)
+    op.Sub(zq, xq, zq)
+    op.Square(&t1, &t1)
+    op.Square(zq, zq)
+    op.Add(xq, &t0, &t1)
+    op.Sub(&t0, zq, &t0)
+    op.Mul(xq, xq, &t1)
+    op.Mul(zq, zq, &t0)
+    return q
+}
+
+/* -------------------------------------------------------------------------
+   Utils
+   -------------------------------------------------------------------------*/
+func (point *ProjectivePoint) ToAffine(c *CurveOperations) *Fp2Element {
+    var affine_x Fp2Element
+    c.Params.Op.Inv(&affine_x, &point.Z)
+    c.Params.Op.Mul(&affine_x, &affine_x, &point.X)
+    return &affine_x
+}
+
+// Cleans data in fp
+func (fp *Fp2Element) Zeroize() {
+    // Zeroizing in 2 seperated loops tells compiler to
+    // use fast runtime.memclr()
+    for i:=range fp.A {
+        fp.A[i] = 0
+    }
+    for i:=range fp.B {
+        fp.B[i] = 0
+    }
+}

internal/isogeny/types.go

@@ -0,0 +1,146 @@
+package internal
+
+const (
+    FP_MAX_WORDS = 12 // Currently p751.NumWords
+)
+
+// I keep it bool in order to be able to apply logical NOT
+type KeyVariant uint
+type PrimeFieldId uint
+
+// Representation of an element of the base field F_p.
+//
+// No particular meaning is assigned to the representation -- it could represent
+// an element in Montgomery form, or not.  Tracking the meaning of the field
+// element is left to higher types.
+type FpElement [FP_MAX_WORDS]uint64
+// Represents an intermediate product of two elements of the base field F_p.
+type FpElementX2 [2 * FP_MAX_WORDS]uint64
+// Represents an element of the extended field Fp^2 = Fp(x+i)
+type Fp2Element struct {
+    A FpElement
+    B FpElement
+}
+
+type DomainParams struct {
+    // P, Q and R=P-Q base points
+    Affine_P, Affine_Q, Affine_R Fp2Element
+    // Max size of secret key for x-torsion group
+    SecretBitLen uint
+    // MaskBytes
+    MaskBytes []byte
+    // Size of a compuatation strategy for x-torsion group
+    IsogenyStrategy []uint32
+}
+
+type SidhParams struct {
+    Id PrimeFieldId
+    // Bytelen of P
+    Bytelen int
+    // The secret key size, in bytes.
+    SecretKeySize int
+    // The public key size, in bytes.
+    PublicKeySize int
+    // The shared secret size, in bytes.
+    SharedSecretSize uint
+    // 2- and 3-torsion group parameter definitions
+    A, B DomainParams
+    // Sample rate to obtain a value in [0,3^238]
+    SampleRate uint
+    // Precomputed identity element in the Fp2 in Montgomery domain
+    OneFp2 Fp2Element
+    // Precomputed 1/2 in the Fp2 in Montgomery domain
+    HalfFp2 Fp2Element
+    // Length of SIKE secret message. Must be one of {24,32,40},
+    // depending on size of prime field used (see [SIKE], 1.4 and 5.1)
+    MsgLen uint
+    // Length of SIKE ephemeral KEM key (see [SIKE], 1.4 and 5.1)
+    KemSize uint
+    // Access to field arithmetic
+    Op FieldOps
+}
+
+// Interface for working with isogenies.
+type Isogeny interface {
+    // Given a torsion point on a curve computes isogenous curve.
+    // Returns curve coefficients (A:C), so that E_(A/C) = E_(A/C)/<P>,
+    // where P is a provided projective point. Sets also isogeny constants
+    // that are needed for isogeny evaluation.
+    GenerateCurve(*ProjectivePoint) CurveCoefficientsEquiv
+    // Evaluates isogeny at caller provided point. Requires isogeny curve constants
+    // to be earlier computed by GenerateCurve.
+    EvaluatePoint(*ProjectivePoint) ProjectivePoint
+}
+
+// Stores curve projective parameters equivalent to A/C. Meaning of the
+// values depends on the context. When working with isogenies over
+// subgroup that are powers of:
+// * three then  (A:C) ~ (A+2C:A-2C)
+// * four then   (A:C) ~ (A+2C:  4C)
+// See Appendix A of SIKE for more details
+type CurveCoefficientsEquiv struct {
+    A Fp2Element
+    C Fp2Element
+}
+
+// A point on the projective line P^1(F_{p^2}).
+//
+// This represents a point on the Kummer line of a Montgomery curve.  The
+// curve is specified by a ProjectiveCurveParameters struct.
+type ProjectivePoint struct {
+    X Fp2Element
+    Z Fp2Element
+}
+
+// A point on the projective line P^1(F_{p^2}).
+//
+// This is used to work projectively with the curve coefficients.
+type ProjectiveCurveParameters struct {
+    A Fp2Element
+    C Fp2Element
+}
+
+// Stores Isogeny 3 curve constants
+type isogeny3 struct {
+    Field FieldOps
+    K1 Fp2Element
+    K2 Fp2Element
+}
+
+// Stores Isogeny 4 curve constants
+type isogeny4 struct {
+    isogeny3
+    K3 Fp2Element
+}
+
+type FieldOps interface {
+    // Set res = lhs + rhs.
+    //
+    // Allowed to overlap lhs or rhs with res.
+    Add(res, lhs, rhs *Fp2Element)
+
+    // Set res = lhs - rhs.
+    //
+    // Allowed to overlap lhs or rhs with res.
+    Sub(res, lhs, rhs *Fp2Element)
+
+    // Set res = lhs * rhs.
+    //
+    // Allowed to overlap lhs or rhs with res.
+    Mul(res, lhs, rhs *Fp2Element)
+    // Set res = x * x
+    //
+    // Allowed to overlap res with x.
+    Square(res, x *Fp2Element)
+    // Set res = 1/x
+    //
+    // Allowed to overlap res with x.
+    Inv(res, x *Fp2Element)
+    // If choice = 1u8, set (x,y) = (y,x). If choice = 0u8, set (x,y) = (x,y).
+    CondSwap(xPx, xPz, xQx, xQz *Fp2Element, choice uint8)
+    // Converts Fp2Element to Montgomery domain (x*R mod p)
+    ToMontgomery(x *Fp2Element)
+    // Converts 'a' in montgomery domain to element from Fp2Element
+    // and stores it in 'x'
+    FromMontgomery(x *Fp2Element, a *Fp2Element)
+}

p751toolbox/field_decl.go → p751/arith_decl.go

@@ -1,59 +1,60 @@
 // +build amd64,!noasm
 
-package p751toolbox
+package p751
 
 import (
 	"golang.org/x/sys/cpu"
+	. "github.com/cloudflare/p751sidh/internal/isogeny"
 )
 
 // If choice = 0, leave x,y unchanged. If choice = 1, set x,y = y,x.
 // If choice is neither 0 nor 1 then behaviour is undefined.
 // This function executes in constant time.
 //go:noescape
-func fp751ConditionalSwap(x, y *Fp751Element, choice uint8)
+func fp751ConditionalSwap(x, y *FpElement, choice uint8)
 
 // Compute z = x + y (mod p).
 //go:noescape
-func fp751AddReduced(z, x, y *Fp751Element)
+func fp751AddReduced(z, x, y *FpElement)
 
 // Compute z = x - y (mod p).
 //go:noescape
-func fp751SubReduced(z, x, y *Fp751Element)
+func fp751SubReduced(z, x, y *FpElement)
 
 // Compute z = x + y, without reducing mod p.
 //go:noescape
-func fp751AddLazy(z, x, y *Fp751Element)
+func fp751AddLazy(z, x, y *FpElement)
 
 // Compute z = x + y, without reducing mod p.
 //go:noescape
-func fp751X2AddLazy(z, x, y *fp751X2)
+func fp751X2AddLazy(z, x, y *FpElementX2)
 
 // Compute z = x - y, without reducing mod p.
 //go:noescape
-func fp751X2SubLazy(z, x, y *fp751X2)
+func fp751X2SubLazy(z, x, y *FpElementX2)
 
 // Compute z = x * y.
 //go:noescape
-func fp751Mul(z *fp751X2, x, y *Fp751Element)
+func fp751Mul(z *FpElementX2, x, y *FpElement)
 
 // Function pointer that should point to one of the
 // fp751MontgomeryReduce implementations below.
 // When set, it performs Montgomery reduction: set z = x R^{-1} (mod 2*p).
 // It may destroy the input value.
-var fp751MontgomeryReduce func(z *Fp751Element, x *fp751X2)
+var fp751MontgomeryReduce func(z *FpElement, x *FpElementX2)
 
 //go:noescape
-func fp751MontgomeryReduceBMI2ADX(z *Fp751Element, x *fp751X2)
+func fp751MontgomeryReduceBMI2ADX(z *FpElement, x *FpElementX2)
 
 //go:noescape
-func fp751MontgomeryReduceBMI2(z *Fp751Element, x *fp751X2)
+func fp751MontgomeryReduceBMI2(z *FpElement, x *FpElementX2)
 
 //go:noescape
-func fp751MontgomeryReduceFallback(z *Fp751Element, x *fp751X2)
+func fp751MontgomeryReduceFallback(z *FpElement, x *FpElementX2)
 
 // Reduce a field element in [0, 2*p) to one in [0,p).
 //go:noescape
-func fp751StrongReduce(x *Fp751Element)
+func fp751StrongReduce(x *FpElement)
 
 // On initialization, set the fp751MontgomeryReduce function pointer to the
 // fastest implementation depending on CPU capabilities.

p751toolbox/field_generic.go → p751/arith_generic.go

@@ -1,6 +1,10 @@
 // +build noasm arm64 arm
 
-package p751toolbox
+package p751
+
+import (
+	. "github.com/cloudflare/p751sidh/internal/isogeny"
+)
 
 // helper used for uint128 representation
 type uint128 struct {
@@ -66,40 +70,40 @@ func mul64(a, b uint64) (res uint128) {
 }
 
 // Compute z = x + y (mod p).
-func fp751AddReduced(z, x, y *Fp751Element) {
+func fp751AddReduced(z, x, y *FpElement) {
 	var carry uint64
 
 	// z=x+y % p751
-	for i := 0; i < fp751NumWords; i++ {
+	for i := 0; i < NumWords; i++ {
 		z[i], carry = addc64(carry, x[i], y[i])
 	}
 
 	// z = z - p751x2
 	carry = 0
-	for i := 0; i < fp751NumWords; i++ {
+	for i := 0; i < NumWords; i++ {
 		z[i], carry = subc64(carry, z[i], p751x2[i])
 	}
 
 	// z = z + p751x2
 	mask := uint64(0 - carry)
 	carry = 0
-	for i := 0; i < fp751NumWords; i++ {
+	for i := 0; i < NumWords; i++ {
 		z[i], carry = addc64(carry, z[i], p751x2[i]&mask)
 	}
 }
 
 // Compute z = x - y (mod p).
-func fp751SubReduced(z, x, y *Fp751Element) {
+func fp751SubReduced(z, x, y *FpElement) {
 	var borrow uint64
 
-	for i := 0; i < fp751NumWords; i++ {
+	for i := 0; i < NumWords; i++ {
 		z[i], borrow = subc64(borrow, x[i], y[i])
 	}
 
 	mask := uint64(0 - borrow)
 	borrow = 0
 
-	for i := 0; i < fp751NumWords; i++ {
+	for i := 0; i < NumWords; i++ {
 		z[i], borrow = addc64(borrow, z[i], p751x2[i]&mask)
 	}
 }
@@ -109,7 +113,7 @@ func fp751SubReduced(z, x, y *Fp751Element) {
 // For details see "Hackers Delight, 2.20"
 //
 // Implementation doesn't actually depend on a prime field.
-func fp751ConditionalSwap(x, y *Fp751Element, mask uint8) {
+func fp751ConditionalSwap(x, y *FpElement, mask uint8) {
 	var tmp, mask64 uint64
 
 	mask64 = 0 - uint64(mask)
@@ -122,14 +126,14 @@ func fp751ConditionalSwap(x, y *Fp751Element, mask uint8) {
 
 // Perform Montgomery reduction: set z = x R^{-1} (mod 2*p)
 // with R=2^768. Destroys the input value.
-func fp751MontgomeryReduce(z *Fp751Element, x *fp751X2) {
+func fp751MontgomeryReduce(z *FpElement, x *FpElementX2) {
 	var carry, t, u, v uint64
 	var uv uint128
 	var count int
 
 	count = 5 // number of 0 digits in the least significat part of p751 + 1
 
-	for i := 0; i < fp751NumWords; i++ {
+	for i := 0; i < NumWords; i++ {
 		for j := 0; j < i; j++ {
 			if j < (i - count + 1) {
 				uv = mul64(z[j], p751p1[i-j])
@@ -148,12 +152,12 @@ func fp751MontgomeryReduce(z *Fp751Element, x *fp751X2) {
 		t = 0
 	}
 
-	for i := fp751NumWords; i < 2*fp751NumWords-1; i++ {
+	for i := NumWords; i < 2*NumWords-1; i++ {
 		if count > 0 {
 			count--
 		}
-		for j := i - fp751NumWords + 1; j < fp751NumWords; j++ {
-			if j < (fp751NumWords - count) {
+		for j := i - NumWords + 1; j < NumWords; j++ {
+			if j < (NumWords - count) {
 				uv = mul64(z[j], p751p1[i-j])
 				v, carry = addc64(0, uv.L, v)
 				u, carry = addc64(carry, uv.H, u)
@@ -164,22 +168,22 @@ func fp751MontgomeryReduce(z *Fp751Element, x *fp751X2) {
 		u, carry = addc64(carry, u, 0)
 
 		t += carry
-		z[i-fp751NumWords] = v
+		z[i-NumWords] = v
 		v = u
 		u = t
 		t = 0
 	}
-	v, carry = addc64(0, v, x[2*fp751NumWords-1])
-	z[fp751NumWords-1] = v
+	v, carry = addc64(0, v, x[2*NumWords-1])
+	z[NumWords-1] = v
 }
 
 // Compute z = x * y.
-func fp751Mul(z *fp751X2, x, y *Fp751Element) {
+func fp751Mul(z *FpElementX2, x, y *FpElement) {
 	var u, v, t uint64
 	var carry uint64
 	var uv uint128
 
-	for i := uint64(0); i < fp751NumWords; i++ {
+	for i := uint64(0); i < NumWords; i++ {
 		for j := uint64(0); j <= i; j++ {
 			uv = mul64(x[j], y[i-j])
 			v, carry = addc64(0, uv.L, v)
@@ -192,8 +196,8 @@ func fp751Mul(z *fp751X2, x, y *Fp751Element) {
 		t = 0
 	}
 
-	for i := fp751NumWords; i < (2*fp751NumWords)-1; i++ {
-		for j := i - fp751NumWords + 1; j < fp751NumWords; j++ {
+	for i := NumWords; i < (2*NumWords)-1; i++ {
+		for j := i - NumWords + 1; j < NumWords; j++ {
 			uv = mul64(x[j], y[i-j])
 			v, carry = addc64(0, uv.L, v)
 			u, carry = addc64(carry, uv.H, u)
@@ -204,42 +208,42 @@ func fp751Mul(z *fp751X2, x, y *Fp751Element) {
 		u = t
 		t = 0
 	}
-	z[2*fp751NumWords-1] = v
+	z[2*NumWords-1] = v
 }
 
 // Compute z = x + y, without reducing mod p.
-func fp751AddLazy(z, x, y *Fp751Element) {
+func fp751AddLazy(z, x, y *FpElement) {
 	var carry uint64
-	for i := 0; i < fp751NumWords; i++ {
+	for i := 0; i < NumWords; i++ {
 		z[i], carry = addc64(carry, x[i], y[i])
 	}
 }
 
 // Compute z = x + y, without reducing mod p.
-func fp751X2AddLazy(z, x, y *fp751X2) {
+func fp751X2AddLazy(z, x, y *FpElementX2) {
 	var carry uint64
-	for i := 0; i < 2*fp751NumWords; i++ {
+	for i := 0; i < 2*NumWords; i++ {
 		z[i], carry = addc64(carry, x[i], y[i])
 	}
 }
 
 // Reduce a field element in [0, 2*p) to one in [0,p).
-func fp751StrongReduce(x *Fp751Element) {
+func fp751StrongReduce(x *FpElement) {
 	var borrow, mask uint64
-	for i := 0; i < fp751NumWords; i++ {
+	for i := 0; i < NumWords; i++ {
 		x[i], borrow = subc64(borrow, x[i], p751[i])
 	}
 
 	// Sets all bits if borrow = 1
 	mask = 0 - borrow
 	borrow = 0
-	for i := 0; i < fp751NumWords; i++ {
+	for i := 0; i < NumWords; i++ {
 		x[i], borrow = addc64(borrow, x[i], p751[i]&mask)
 	}
 }
 
 // Compute z = x - y, without reducing mod p.
-func fp751X2SubLazy(z, x, y *fp751X2) {
+func fp751X2SubLazy(z, x, y *FpElementX2) {
 	var borrow, mask uint64
 	for i := 0; i < len(z); i++ {
 		z[i], borrow = subc64(borrow, x[i], y[i])
@@ -248,7 +252,7 @@ func fp751X2SubLazy(z, x, y *fp751X2) {
 	// Sets all bits if borrow = 1
 	mask = 0 - borrow
 	borrow = 0
-	for i := fp751NumWords; i < len(z); i++ {
-		z[i], borrow = addc64(borrow, z[i], p751[i-fp751NumWords]&mask)
+	for i := NumWords; i < len(z); i++ {
+		z[i], borrow = addc64(borrow, z[i], p751[i-NumWords]&mask)
 	}
 }

p751/arith_test.go

@@ -0,0 +1,322 @@
+package p751
+
+import (
+	"math/big"
+	"testing"
+    "testing/quick"
+	. "github.com/cloudflare/p751sidh/internal/isogeny"
+)
+
+func TestPrimeFieldElementToBigInt(t *testing.T) {
+	// Chosen so that p < xR < 2p
+	x := primeFieldElement{A: FpElement{
+		1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 140737488355328,
+	}}
+	// Computed using Sage:
+	// sage: p = 2^372 * 3^239 - 1
+	// sage: R = 2^768
+	// sage: from_radix_64 = lambda xs: sum((xi * (2**64)**i for i,xi in enumerate(xs)))
+	// sage: xR = from_radix_64([1]*11 + [2^47])
+	// sage: assert(p < xR)
+	// sage: assert(xR < 2*p)
+	// sage: (xR / R) % p
+	xBig, _ := new(big.Int).SetString("4469946751055876387821312289373600189787971305258234719850789711074696941114031433609871105823930699680637820852699269802003300352597419024286385747737509380032982821081644521634652750355306547718505685107272222083450567982240", 10)
+	if xBig.Cmp(toBigInt(&x.A)) != 0 {
+		t.Error("Expected", xBig, "found", toBigInt(&x.A))
+	}
+}
+
+//------------------------------------------------------------------------------
+// Extended Field
+//------------------------------------------------------------------------------
+
+func TestOneFp2ToBytes(t *testing.T) {
+	var x = P751_OneFp2
+	var xBytes [188]byte
+
+	kCurveOps.Fp2ToBytes(xBytes[:], &x)
+	if xBytes[0] != 1 {
+		t.Error("Expected 1, got", xBytes[0])
+	}
+	for i := 1; i < 188; i++ {
+		if xBytes[i] != 0 {
+			t.Error("Expected 0, got", xBytes[0])
+		}
+	}
+}
+
+func TestFp2ElementToBytesRoundTrip(t *testing.T) {
+	roundTrips := func(x GeneratedTestParams) bool {
+		var xBytes [188]byte
+		var xPrime Fp2Element
+		kCurveOps.Fp2ToBytes(xBytes[:], &x.ExtElem)
+		kCurveOps.Fp2FromBytes(&xPrime, xBytes[:])
+
+		return VartimeEqFp2(&xPrime, &x.ExtElem)
+	}
+
+	if err := quick.Check(roundTrips, quickCheckConfig); err != nil {
+		t.Error(err)
+	}
+}
+
+func TestFp2ElementMulDistributesOverAdd(t *testing.T) {
+	mulDistributesOverAdd := func(x, y, z GeneratedTestParams) bool {
+		// Compute t1 = (x+y)*z
+		t1 := new(Fp2Element)
+		kFieldOps.Add(t1, &x.ExtElem, &y.ExtElem)
+		kFieldOps.Mul(t1, t1, &z.ExtElem)
+
+		// Compute t2 = x*z + y*z
+		t2 := new(Fp2Element)
+		t3 := new(Fp2Element)
+		kFieldOps.Mul(t2, &x.ExtElem, &z.ExtElem)
+		kFieldOps.Mul(t3, &y.ExtElem, &z.ExtElem)
+		kFieldOps.Add(t2, t2, t3)
+
+		return VartimeEqFp2(t1, t2)
+	}
+
+	if err := quick.Check(mulDistributesOverAdd, quickCheckConfig); err != nil {
+		t.Error(err)
+	}
+}
+
+func TestFp2ElementMulIsAssociative(t *testing.T) {
+	isAssociative := func(x, y, z GeneratedTestParams) bool {
+		// Compute t1 = (x*y)*z
+		t1 := new(Fp2Element)
+		kFieldOps.Mul(t1, &x.ExtElem, &y.ExtElem)
+		kFieldOps.Mul(t1, t1, &z.ExtElem)
+
+		// Compute t2 = (y*z)*x
+		t2 := new(Fp2Element)
+		kFieldOps.Mul(t2, &y.ExtElem, &z.ExtElem)
+		kFieldOps.Mul(t2, t2, &x.ExtElem)
+
+		return VartimeEqFp2(t1, t2)
+	}
+
+	if err := quick.Check(isAssociative, quickCheckConfig); err != nil {
+		t.Error(err)
+	}
+}
+
+func TestFp2ElementSquareMatchesMul(t *testing.T) {
+	sqrMatchesMul := func(x GeneratedTestParams) bool {
+		// Compute t1 = (x*x)
+		t1 := new(Fp2Element)
+		kFieldOps.Mul(t1, &x.ExtElem, &x.ExtElem)
+
+		// Compute t2 = x^2
+		t2 := new(Fp2Element)
+		kFieldOps.Square(t2, &x.ExtElem)
+
+		return VartimeEqFp2(t1, t2)
+	}
+
+	if err := quick.Check(sqrMatchesMul, quickCheckConfig); err != nil {
+		t.Error(err)
+	}
+}
+
+func TestFp2ElementInv(t *testing.T) {
+	inverseIsCorrect := func(x GeneratedTestParams) bool {
+		z := new(Fp2Element)
+		kFieldOps.Inv(z, &x.ExtElem)
+
+		// Now z = (1/x), so (z * x) * x == x
+		kFieldOps.Mul(z, z, &x.ExtElem)
+		kFieldOps.Mul(z, z, &x.ExtElem)
+
+		return VartimeEqFp2(z, &x.ExtElem)
+	}
+
+	// This is more expensive; run fewer tests
+	var quickCheckConfig = &quick.Config{MaxCount: (1 << (8 + quickCheckScaleFactor))}
+	if err := quick.Check(inverseIsCorrect, quickCheckConfig); err != nil {
+		t.Error(err)
+	}
+}
+
+func TestFp2ElementBatch3Inv(t *testing.T) {
+	batchInverseIsCorrect := func(x1, x2, x3 GeneratedTestParams) bool {
+		var x1Inv, x2Inv, x3Inv Fp2Element
+		kFieldOps.Inv(&x1Inv, &x1.ExtElem)
+		kFieldOps.Inv(&x2Inv, &x2.ExtElem)
+		kFieldOps.Inv(&x3Inv, &x3.ExtElem)
+
+		var y1, y2, y3 Fp2Element
+		kCurveOps.Fp2Batch3Inv(&x1.ExtElem, &x2.ExtElem, &x3.ExtElem, &y1, &y2, &y3)
+
+		return (VartimeEqFp2(&x1Inv,&y1) && VartimeEqFp2(&x2Inv,&y2) && VartimeEqFp2(&x3Inv,&y3))
+	}
+
+	// This is more expensive; run fewer tests
+	var quickCheckConfig = &quick.Config{MaxCount: (1 << (5 + quickCheckScaleFactor))}
+	if err := quick.Check(batchInverseIsCorrect, quickCheckConfig); err != nil {
+		t.Error(err)
+	}
+}
+
+//------------------------------------------------------------------------------
+// Prime Field
+//------------------------------------------------------------------------------
+
+func TestPrimeFieldElementMulVersusBigInt(t *testing.T) {
+	mulMatchesBigInt := func(x, y primeFieldElement) bool {
+		z := new(primeFieldElement)
+		z.Mul(&x, &y)
+
+		check := new(big.Int)
+		check.Mul(toBigInt(&x.A), toBigInt(&y.A))
+		check.Mod(check, cln16prime)
+
+		return check.Cmp(toBigInt(&z.A)) == 0
+	}
+
+	if err := quick.Check(mulMatchesBigInt, quickCheckConfig); err != nil {
+		t.Error(err)
+	}
+}
+
+func TestPrimeFieldElementP34VersusBigInt(t *testing.T) {
+	var p34, _ = new(big.Int).SetString("2588679435442326313244442059466701330356847411387267792529047419763669735170619711625720724140266678406138302904710050596300977994130638598261040117192787954244176710019728333589599932738193731745058771712747875468166412894207", 10)
+	p34MatchesBigInt := func(x primeFieldElement) bool {
+		z := new(primeFieldElement)
+		z.P34(&x)
+
+		check := toBigInt(&x.A)
+		check.Exp(check, p34, cln16prime)
+
+		return check.Cmp(toBigInt(&z.A)) == 0
+	}
+
+	// This is more expensive; run fewer tests
+	var quickCheckConfig = &quick.Config{MaxCount: (1 << (8 + quickCheckScaleFactor))}
+	if err := quick.Check(p34MatchesBigInt, quickCheckConfig); err != nil {
+		t.Error(err)
+	}
+}
+
+func BenchmarkFp2ElementMul(b *testing.B) {
+	z := &Fp2Element{A: bench_x, B: bench_y}
+	w := new(Fp2Element)
+
+	for n := 0; n < b.N; n++ {
+		kFieldOps.Mul(w, z, z)
+	}
+}
+
+func BenchmarkFp2ElementInv(b *testing.B) {
+	z := &Fp2Element{A: bench_x, B: bench_y}
+	w := new(Fp2Element)
+
+	for n := 0; n < b.N; n++ {
+		kFieldOps.Inv(w, z)
+	}
+}
+
+func BenchmarkFp2ElementSquare(b *testing.B) {
+	z := &Fp2Element{A: bench_x, B: bench_y}
+	w := new(Fp2Element)
+
+	for n := 0; n < b.N; n++ {
+		kFieldOps.Square(w, z)
+	}
+}
+
+func BenchmarkFp2ElementAdd(b *testing.B) {
+	z := &Fp2Element{A: bench_x, B: bench_y}
+	w := new(Fp2Element)
+
+	for n := 0; n < b.N; n++ {
+		kFieldOps.Add(w, z, z)
+	}
+}
+
+func BenchmarkFp2ElementSub(b *testing.B) {
+	z := &Fp2Element{A: bench_x, B: bench_y}
+	w := new(Fp2Element)
+
+	for n := 0; n < b.N; n++ {
+		kFieldOps.Sub(w, z, z)
+	}
+}
+
+func BenchmarkPrimeFieldElementMul(b *testing.B) {
+	z := &primeFieldElement{A: bench_x}
+	w := new(primeFieldElement)
+
+	for n := 0; n < b.N; n++ {
+		w.Mul(z, z)
+	}
+}
+
+// --- field operation functions
+
+func BenchmarkFp751Multiply(b *testing.B) {
+	for n := 0; n < b.N; n++ {
+		fp751Mul(&benchmarkFpElementX2, &bench_x, &bench_y)
+	}
+}
+
+func BenchmarkFp751MontgomeryReduce(b *testing.B) {
+	z := bench_z
+
+	// This benchmark actually computes garbage, because
+	// fp751MontgomeryReduce mangles its input, but since it's
+	// constant-time that shouldn't matter for the benchmarks.
+	for n := 0; n < b.N; n++ {
+		fp751MontgomeryReduce(&benchmarkFpElement, &z)
+	}
+}
+
+func BenchmarkFp751AddReduced(b *testing.B) {
+	for n := 0; n < b.N; n++ {
+		fp751AddReduced(&benchmarkFpElement, &bench_x, &bench_y)
+	}
+}
+
+func BenchmarkFp751SubReduced(b *testing.B) {
+	for n := 0; n < b.N; n++ {
+		fp751SubReduced(&benchmarkFpElement, &bench_x, &bench_y)
+	}
+}
+
+func BenchmarkFp751ConditionalSwap(b *testing.B) {
+	x, y := bench_x, bench_y
+	for n := 0; n < b.N; n++ {
+		fp751ConditionalSwap(&x, &y, 1)
+		fp751ConditionalSwap(&x, &y, 0)
+	}
+}
+
+func BenchmarkFp751StrongReduce(b *testing.B) {
+	x := bench_x
+	for n := 0; n < b.N; n++ {
+		fp751StrongReduce(&x)
+	}
+}
+
+func BenchmarkFp751AddLazy(b *testing.B) {
+	var z FpElement
+	x, y := bench_x, bench_y
+	for n := 0; n < b.N; n++ {
+		fp751AddLazy(&z, &x, &y)
+	}
+}
+
+func BenchmarkFp751X2AddLazy(b *testing.B) {
+	x, y, z := bench_z, bench_z, bench_z
+	for n := 0; n < b.N; n++ {
+		fp751X2AddLazy(&x, &y, &z)
+	}
+}
+
+func BenchmarkFp751X2SubLazy(b *testing.B) {
+	x, y, z := bench_z, bench_z, bench_z
+	for n := 0; n < b.N; n++ {
+		fp751X2SubLazy(&x, &y, &z)
+	}
+}

p751toolbox/consts.go → p751/consts.go

@@ -1,4 +1,6 @@
-package p751toolbox
+package p751
+
+import . "github.com/cloudflare/p751sidh/internal/isogeny"
 
 const (
 	// The secret key size, in bytes. Secret key is actually different for
@@ -21,21 +23,25 @@ const (
 	P751_MaskBobByte = 0x03
 	// Sample rate to obtain a value in [0,3^238]
 	P751_SampleRate = 102
+	// P751 bytelen ceil(751/8)
+	P751_Bytelen = 94
 	// Size of a compuatation strategy for 2-torsion group
 	strategySizeA = 185
 	// Size of a compuatation strategy for 3-torsion group
 	strategySizeB = 238
+	// Number of 64-bit limbs used to store Fp element
+	NumWords = 12
 )
 
 // The x-coordinate of PA
-var P751_affine_PA = ExtensionFieldElement{
-	A: Fp751Element{
+var P751_affine_PA = Fp2Element{
+	A: FpElement{
 		0xC2FC08CEAB50AD8B, 0x1D7D710F55E457B1, 0xE8738D92953DCD6E,
 		0xBAA7EBEE8A3418AA, 0xC9A288345F03F46F, 0xC8D18D167CFE2616,
 		0x02043761F6B1C045, 0xAA1975E13180E7E9, 0x9E13D3FDC6690DE6,
 		0x3A024640A3A3BB4F, 0x4E5AD44E6ACBBDAE, 0x0000544BEB561DAD,
 	},
-	B: Fp751Element{
+	B: FpElement{
 		0xE6CC41D21582E411, 0x07C2ECB7C5DF400A, 0xE8E34B521432AEC4,
 		0x50761E2AB085167D, 0x032CFBCAA6094B3C, 0x6C522F5FDF9DDD71,
 		0x1319217DC3A1887D, 0xDC4FB25803353A86, 0x362C8D7B63A6AB09,
@@ -44,14 +50,14 @@ var P751_affine_PA = ExtensionFieldElement{
 }
 
 // The x-coordinate of QA
-var P751_affine_QA = ExtensionFieldElement{
-	A: Fp751Element{
+var P751_affine_QA = Fp2Element{
+	A: FpElement{
 		0xD56FE52627914862, 0x1FAD60DC96B5BAEA, 0x01E137D0BF07AB91,
 		0x404D3E9252161964, 0x3C5385E4CD09A337, 0x4476426769E4AF73,
 		0x9790C6DB989DFE33, 0xE06E1C04D2AA8B5E, 0x38C08185EDEA73B9,
 		0xAA41F678A4396CA6, 0x92B9259B2229E9A0, 0x00002F9326818BE0,
 	},
-	B: Fp751Element{
+	B: FpElement{
 		0x0000000000000000, 0x0000000000000000, 0x0000000000000000,
 		0x0000000000000000, 0x0000000000000000, 0x0000000000000000,
 		0x0000000000000000, 0x0000000000000000, 0x0000000000000000,
@@ -60,14 +66,14 @@ var P751_affine_QA = ExtensionFieldElement{
 }
 
 // The x-coordinate of RA = PA-QA
-var P751_affine_RA = ExtensionFieldElement{
-	A: Fp751Element{
+var P751_affine_RA = Fp2Element{
+	A: FpElement{
 		0x0BB84441DFFD19B3, 0x84B4DEA99B48C18E, 0x692DE648AD313805,
 		0xE6D72761B6DFAEE0, 0x223975C672C3058D, 0xA0FDE0C3CBA26FDC,
 		0xA5326132A922A3CA, 0xCA5E7F5D5EA96FA4, 0x127C7EFE33FFA8C6,
 		0x4749B1567E2A23C4, 0x2B7DF5B4AF413BFA, 0x0000656595B9623C,
 	},
-	B: Fp751Element{
+	B: FpElement{
 		0xED78C17F1EC71BE8, 0xF824D6DF753859B1, 0x33A10839B2A8529F,
 		0xFC03E9E25FDEA796, 0xC4708A8054DF1762, 0x4034F2EC034C6467,
 		0xABFB70FBF06ECC79, 0xDABE96636EC108B7, 0x49CBCFB090605FD3,
@@ -76,14 +82,14 @@ var P751_affine_RA = ExtensionFieldElement{
 }
 
 // The x-coordinate of PB
-var P751_affine_PB = ExtensionFieldElement{
-	A: Fp751Element{
+var P751_affine_PB = Fp2Element{
+	A: FpElement{
 		0xCFB6D71EF867AB0B, 0x4A5FDD76E9A45C76, 0x38B1EE69194B1F03,
 		0xF6E7B18A7761F3F0, 0xFCF01A486A52C84C, 0xCBE2F63F5AA75466,
 		0x6487BCE837B5E4D6, 0x7747F5A8C622E9B8, 0x4CBFE1E4EE6AEBBA,
 		0x8A8616A13FA91512, 0x53DB980E1579E0A5, 0x000058FEBFF3BE69,
 	},
-	B: Fp751Element{
+	B: FpElement{
 		0xA492034E7C075CC3, 0x677BAF00B04AA430, 0x3AAE0C9A755C94C8,
 		0x1DC4B064E9EBB08B, 0x3684EDD04E826C66, 0x9BAA6CB661F01B22,
 		0x20285A00AD2EFE35, 0xDCE95ABD0497065F, 0x16C7FBB3778E3794,
@@ -92,14 +98,14 @@ var P751_affine_PB = ExtensionFieldElement{
 }
 
 // The x-coordinate of QB
-var P751_affine_QB = ExtensionFieldElement{
-	A: Fp751Element{
+var P751_affine_QB = Fp2Element{
+	A: FpElement{
 		0xF1A8C9ED7B96C4AB, 0x299429DA5178486E, 0xEF4926F20CD5C2F4,
 		0x683B2E2858B4716A, 0xDDA2FBCC3CAC3EEB, 0xEC055F9F3A600460,
 		0xD5A5A17A58C3848B, 0x4652D836F42EAED5, 0x2F2E71ED78B3A3B3,
 		0xA771C057180ADD1D, 0xC780A5D2D835F512, 0x0000114EA3B55AC1,
 	},
-	B: Fp751Element{
+	B: FpElement{
 		0x0000000000000000, 0x0000000000000000, 0x0000000000000000,
 		0x0000000000000000, 0x0000000000000000, 0x0000000000000000,
 		0x0000000000000000, 0x0000000000000000, 0x0000000000000000,
@@ -108,14 +114,14 @@ var P751_affine_QB = ExtensionFieldElement{
 }
 
 // The x-coordinate of RB = PB - QB
-var P751_affine_RB = ExtensionFieldElement{
-	A: Fp751Element{
+var P751_affine_RB = Fp2Element{
+	A: FpElement{
 		0x1C0D6733769D0F31, 0xF084C3086E2659D1, 0xE23D5DA27BCBD133,
 		0xF38EC9A8D5864025, 0x6426DC781B3B645B, 0x4B24E8E3C9FB03EE,
 		0x6432792F9D2CEA30, 0x7CC8E8B1AE76E857, 0x7F32BFB626BB8963,
 		0xB9F05995B48D7B74, 0x4D71200A7D67E042, 0x0000228457AF0637,
 	},
-	B: Fp751Element{
+	B: FpElement{
 		0x4AE37E7D8F72BD95, 0xDD2D504B3E993488, 0x5D14E7FA1ECB3C3E,
 		0x127610CEB75D6350, 0x255B4B4CAC446B11, 0x9EA12336C1F70CAF,
 		0x79FA68A2147BC2F8, 0x11E895CFDADBBC49, 0xE4B9D3C4D6356C18,
@@ -171,3 +177,27 @@ var P751_BobIsogenyStrategy = [strategySizeB]uint32{
 	0x01, 0x01, 0x01, 0x01, 0x02, 0x01, 0x01, 0x01, 0x09, 0x05,
 	0x03, 0x02, 0x01, 0x01, 0x01, 0x01, 0x02, 0x01, 0x01, 0x01,
 	0x04, 0x02, 0x01, 0x01, 0x01, 0x02, 0x01, 0x01}
+
+var P751_OneFp2 = Fp2Element{
+    A: FpElement{
+    	0x249ad, 0x0, 0x0, 0x0, 0x0, 0x8310000000000000, 0x5527b1e4375c6c66, 0x697797bf3f4f24d0, 0xc89db7b2ac5c4e2e, 0x4ca4b439d2076956, 0x10f7926c7512c7e9, 0x2d5b24bce5e2},
+    B: FpElement{0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0},
+}
+
+var P751_HalfFp2 = Fp2Element{
+    A: FpElement{
+        0x00000000000124D6, 0x0000000000000000, 0x0000000000000000,
+        0x0000000000000000, 0x0000000000000000, 0xB8E0000000000000,
+        0x9C8A2434C0AA7287, 0xA206996CA9A378A3, 0x6876280D41A41B52,
+        0xE903B49F175CE04F, 0x0F8511860666D227, 0x00004EA07CFF6E7F},
+    B: FpElement{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+}
+
+// (2^768)^2 mod p
+// This can't be a constant because Go doesn't allow array constants, so try
+// not to modify it.
+var montgomeryRsq = FpElement{
+	2535603850726686808, 15780896088201250090, 6788776303855402382,
+	17585428585582356230, 5274503137951975249, 2266259624764636289,
+	11695651972693921304, 13072885652150159301, 4908312795585420432,
+	6229583484603254826, 488927695601805643, 72213483953973}

p751/curve_test.go

@@ -0,0 +1,192 @@
+package p751
+
+import (
+	"bytes"
+	"testing"
+	"testing/quick"
+    . "github.com/cloudflare/p751sidh/internal/isogeny"
+)
+
+func TestOne(t *testing.T) {
+	var tmp Fp2Element
+
+	kFieldOps.Mul(&tmp, &P751_OneFp2, &affine_xP)
+	if !VartimeEqFp2(&tmp, &affine_xP) {
+		t.Error("Not equal 1")
+	}
+}
+
+// This test is here only to ensure that ScalarMult helper works correctly
+func TestScalarMultVersusSage(t *testing.T) {
+	var xP ProjectivePoint
+
+	xP = ProjectivePoint{X: affine_xP, Z: P751_OneFp2}
+	xP = ScalarMult(&curve, &xP, mScalarBytes[:]) // = x([m]P)
+	affine_xQ := xP.ToAffine(kCurveOps)
+	if !VartimeEqFp2(&affine_xaP, affine_xQ) {
+		t.Error("\nExpected\n", affine_xaP, "\nfound\n", affine_xQ)
+	}
+}
+
+func Test_jInvariant(t *testing.T) {
+	var curve = ProjectiveCurveParameters{A: curve_A, C: curve_C}
+	var jbufRes [P751_SharedSecretSize]byte
+	var jbufExp [P751_SharedSecretSize]byte
+	// Computed using Sage
+	// j = 3674553797500778604587777859668542828244523188705960771798425843588160903687122861541242595678107095655647237100722594066610650373491179241544334443939077738732728884873568393760629500307797547379838602108296735640313894560419*i + 3127495302417548295242630557836520229396092255080675419212556702820583041296798857582303163183558315662015469648040494128968509467224910895884358424271180055990446576645240058960358037224785786494172548090318531038910933793845
+	var known_j = Fp2Element{
+		A: FpElement{0xc7a8921c1fb23993, 0xa20aea321327620b, 0xf1caa17ed9676fa8, 0x61b780e6b1a04037, 0x47784af4c24acc7a, 0x83926e2e300b9adf, 0xcd891d56fae5b66, 0x49b66985beb733bc, 0xd4bcd2a473d518f, 0xe242239991abe224, 0xa8af5b20f98672f8, 0x139e4d4e4d98},
+		B: FpElement{0xb5b52a21f81f359, 0x715e3a865db6d920, 0x9bac2f9d8911978b, 0xef14acd8ac4c1e3d, 0xe81aacd90cfb09c8, 0xaf898288de4a09d9, 0xb85a7fb88c5c4601, 0x2c37c3f1dd303387, 0x7ad3277fe332367c, 0xd4cbee7f25a8e6f8, 0x36eacbe979eaeffa, 0x59eb5a13ac33},
+	}
+	kCurveOps.Jinvariant(&curve, jbufRes[:])
+	kCurveOps.Fp2ToBytes(jbufExp[:], &known_j)
+
+	if !bytes.Equal(jbufRes[:], jbufExp[:]) {
+		t.Error("Computed incorrect j-invariant: found\n", jbufRes, "\nexpected\n", jbufExp)
+	}
+}
+
+func TestProjectivePointVartimeEq(t *testing.T) {
+	var xP ProjectivePoint
+
+	xP = ProjectivePoint{X: affine_xP, Z: P751_OneFp2}
+	xQ := xP
+	// Scale xQ, which results in the same projective point
+	kFieldOps.Mul(&xQ.X, &xQ.X, &curve_A)
+	kFieldOps.Mul(&xQ.Z, &xQ.Z, &curve_A)
+	if !VartimeEqProjFp2(&xP, &xQ) {
+		t.Error("Expected the scaled point to be equal to the original")
+	}
+}
+
+func TestPointDoubleVersusSage(t *testing.T) {
+	var curve = ProjectiveCurveParameters{A: curve_A, C: curve_C}
+	var params = kCurveOps.CalcCurveParamsEquiv4(&curve)
+	var xP ProjectivePoint
+
+	xP = ProjectivePoint{X: affine_xP, Z: P751_OneFp2}
+	kCurveOps.Pow2k(&xP, &params, 1)
+	affine_xQ := xP.ToAffine(kCurveOps)
+	if !VartimeEqFp2(affine_xQ, &affine_xP2) {
+		t.Error("\nExpected\n", affine_xP2, "\nfound\n", affine_xQ)
+	}
+}
+
+func TestPointMul4VersusSage(t *testing.T) {
+	var params = kCurveOps.CalcCurveParamsEquiv4(&curve)
+	var xP ProjectivePoint
+
+	xP = ProjectivePoint{X: affine_xP, Z: P751_OneFp2}
+	kCurveOps.Pow2k(&xP, &params, 2)
+	affine_xQ := xP.ToAffine(kCurveOps)
+	if !VartimeEqFp2(affine_xQ, &affine_xP4) {
+		t.Error("\nExpected\n", affine_xP4, "\nfound\n", affine_xQ)
+	}
+}
+
+func TestPointMul9VersusSage(t *testing.T) {
+	var params = kCurveOps.CalcCurveParamsEquiv3(&curve)
+	var xP ProjectivePoint
+
+	xP = ProjectivePoint{X: affine_xP, Z: P751_OneFp2}
+	kCurveOps.Pow3k(&xP, &params, 2)
+	affine_xQ := xP.ToAffine(kCurveOps)
+	if !VartimeEqFp2(affine_xQ, &affine_xP9) {
+		t.Error("\nExpected\n", affine_xP9, "\nfound\n", affine_xQ)
+	}
+}
+
+func TestPointPow2kVersusScalarMult(t *testing.T) {
+	var xP, xQ, xR ProjectivePoint
+	var params = kCurveOps.CalcCurveParamsEquiv4(&curve)
+
+	xP = ProjectivePoint{X: affine_xP, Z: P751_OneFp2}
+	xQ = xP
+	kCurveOps.Pow2k(&xQ, &params, 5)
+	xR = ScalarMult(&curve, &xP, []byte{32})
+	affine_xQ := xQ.ToAffine(kCurveOps) // = x([32]P)
+	affine_xR := xR.ToAffine(kCurveOps) // = x([32]P)
+
+	if !VartimeEqFp2(affine_xQ, affine_xR) {
+		t.Error("\nExpected\n", affine_xQ, "\nfound\n", affine_xR)
+	}
+}
+
+func TestRecoverCoordinateA(t *testing.T) {
+	var cparam ProjectiveCurveParameters
+	// Vectors generated with SIKE reference implementation
+	var a = Fp2Element{
+		A: FpElement{0x9331D9C5AAF59EA4, 0xB32B702BE4046931, 0xCEBB333912ED4D34, 0x5628CE37CD29C7A2, 0x0BEAC5ED48B7F58E, 0x1FB9D3E281D65B07, 0x9C0CFACC1E195662, 0xAE4BCE0F6B70F7D9, 0x59E4E63D43FE71A0, 0xEF7CE57560CC8615, 0xE44A8FB7901E74E8, 0x000069D13C8366D1},
+		B: FpElement{0xF6DA1070279AB966, 0xA78FB0CE7268C762, 0x19B40F044A57ABFA, 0x7AC8EE6160C0C233, 0x93D4993442947072, 0x757D2B3FA4E44860, 0x073A920F8C4D5257, 0x2031F1B054734037, 0xDEFAA1D2406555CD, 0x26F9C70E1496BE3D, 0x5B3F335A0A4D0976, 0x000013628B2E9C59}}
+	var affine_xP = Fp2Element{
+		A: FpElement{0xea6b2d1e2aebb250, 0x35d0b205dc4f6386, 0xb198e93cb1830b8d, 0x3b5b456b496ddcc6, 0x5be3f0d41132c260, 0xce5f188807516a00, 0x54f3e7469ea8866d, 0x33809ef47f36286, 0x6fa45f83eabe1edb, 0x1b3391ae5d19fd86, 0x1e66daf48584af3f, 0xb430c14aaa87},
+		B: FpElement{0x97b41ebc61dcb2ad, 0x80ead31cb932f641, 0x40a940099948b642, 0x2a22fd16cdc7fe84, 0xaabf35b17579667f, 0x76c1d0139feb4032, 0x71467e1e7b1949be, 0x678ca8dadd0d6d81, 0x14445daea9064c66, 0x92d161eab4fa4691, 0x8dfbb01b6b238d36, 0x2e3718434e4e}}
+	var affine_xQ = Fp2Element{
+		A: FpElement{0xb055cf0ca1943439, 0xa9ff5de2fa6c69ed, 0x4f2761f934e5730a, 0x61a1dcaa1f94aa4b, 0xce3c8fadfd058543, 0xeac432aaa6701b8e, 0x8491d523093aea8b, 0xba273f9bd92b9b7f, 0xd8f59fd34439bb5a, 0xdc0350261c1fe600, 0x99375ab1eb151311, 0x14d175bbdbc5},
+		B: FpElement{0xffb0ef8c2111a107, 0x55ceca3825991829, 0xdbf8a1ccc075d34b, 0xb8e9187bd85d8494, 0x670aa2d5c34a03b0, 0xef9fe2ed2b064953, 0xc911f5311d645aee, 0xf4411f409e410507, 0x934a0a852d03e1a8, 0xe6274e67ae1ad544, 0x9f4bc563c69a87bc, 0x6f316019681e}}
+	var affine_xQmP = Fp2Element{
+		A: FpElement{0x6ffb44306a153779, 0xc0ffef21f2f918f3, 0x196c46d35d77f778, 0x4a73f80452edcfe6, 0x9b00836bce61c67f, 0x387879418d84219e, 0x20700cf9fc1ec5d1, 0x1dfe2356ec64155e, 0xf8b9e33038256b1c, 0xd2aaf2e14bada0f0, 0xb33b226e79a4e313, 0x6be576fad4e5},
+		B: FpElement{0x7db5dbc88e00de34, 0x75cc8cb9f8b6e11e, 0x8c8001c04ebc52ac, 0x67ef6c981a0b5a94, 0xc3654fbe73230738, 0xc6a46ee82983ceca, 0xed1aa61a27ef49f0, 0x17fe5a13b0858fe0, 0x9ae0ca945a4c6b3c, 0x234104a218ad8878, 0xa619627166104394, 0x556a01ff2e7e}}
+
+	cparam.C = P751_OneFp2
+	kCurveOps.RecoverCoordinateA(&cparam, &affine_xP, &affine_xQ, &affine_xQmP)
+
+	// Check A is correct
+	if !VartimeEqFp2(&cparam.A, &a) {
+		t.Error("\nExpected\n", a, "\nfound\n", cparam.A)
+	}
+
+	// Check C is not changed
+	if !VartimeEqFp2(&cparam.C, &P751_OneFp2) {
+		t.Error("\nExpected\n", cparam.C, "\nfound\n", P751_OneFp2)
+	}
+}
+
+func TestR2LVersusSage(t *testing.T) {
+	var xR ProjectivePoint
+
+	sageAffine_xR := Fp2Element{
+		A: FpElement{0x729465ba800d4fd5, 0x9398015b59e514a1, 0x1a59dd6be76c748e, 0x1a7db94eb28dd55c, 0x444686e680b1b8ec, 0xcc3d4ace2a2454ff, 0x51d3dab4ec95a419, 0xc3b0f33594acac6a, 0x9598a74e7fd44f8a, 0x4fbf8c638f1c2e37, 0x844e347033052f51, 0x6cd6de3eafcf},
+		B: FpElement{0x85da145412d73430, 0xd83c0e3b66eb3232, 0xd08ff2d453ec1369, 0xa64aaacfdb395b13, 0xe9cba211a20e806e, 0xa4f80b175d937cfc, 0x556ce5c64b1f7937, 0xb59b39ea2b3fdf7a, 0xc2526b869a4196b3, 0x8dad90bca9371750, 0xdfb4a30c9d9147a2, 0x346d2130629b}}
+	xR = kCurveOps.ScalarMul3Pt(&curve, &threePointLadderInputs[0], &threePointLadderInputs[1], &threePointLadderInputs[2], uint(len(mScalarBytes)*8), mScalarBytes[:])
+	affine_xR := xR.ToAffine(kCurveOps)
+
+	if !VartimeEqFp2(affine_xR, &sageAffine_xR) {
+		t.Error("\nExpected\n", sageAffine_xR, "\nfound\n", affine_xR)
+	}
+}
+
+func TestPointTripleVersusAddDouble(t *testing.T) {
+	tripleEqualsAddDouble := func(params GeneratedTestParams) bool {
+		var P2, P3, P2plusP ProjectivePoint
+
+		eqivParams4 := kCurveOps.CalcCurveParamsEquiv4(&params.Cparam)
+		eqivParams3 := kCurveOps.CalcCurveParamsEquiv3(&params.Cparam)
+		P2 = params.Point
+		P3 = params.Point
+		kCurveOps.Pow2k(&P2, &eqivParams4, 1)	// = x([2]P)
+		kCurveOps.Pow3k(&P3, &eqivParams3, 1)	// = x([3]P)
+		P2plusP = AddProjFp2(&P2, &params.Point, &params.Point)	// = x([2]P + P)
+		return VartimeEqProjFp2(&P3, &P2plusP)
+	}
+
+	if err := quick.Check(tripleEqualsAddDouble, quickCheckConfig); err != nil {
+		t.Error(err)
+	}
+}
+
+func BenchmarkThreePointLadder379BitScalar(b *testing.B) {
+	var mScalarBytes = [...]uint8{84, 222, 146, 63, 85, 18, 173, 162, 167, 38, 10, 8, 143, 176, 93, 228, 247, 128, 50, 128, 205, 42, 15, 137, 119, 67, 43, 3, 61, 91, 237, 24, 235, 12, 53, 96, 186, 164, 232, 223, 197, 224, 64, 109, 137, 63, 246, 4}
+
+	for n := 0; n < b.N; n++ {
+		kCurveOps.ScalarMul3Pt(&curve, &threePointLadderInputs[0], &threePointLadderInputs[1], &threePointLadderInputs[2], uint(len(mScalarBytes)*8), mScalarBytes[:])
+	}
+}
+
+func BenchmarkR2L379BitScalar(b *testing.B) {
+	var mScalarBytes = [...]uint8{84, 222, 146, 63, 85, 18, 173, 162, 167, 38, 10, 8, 143, 176, 93, 228, 247, 128, 50, 128, 205, 42, 15, 137, 119, 67, 43, 3, 61, 91, 237, 24, 235, 12, 53, 96, 186, 164, 232, 223, 197, 224, 64, 109, 137, 63, 246, 4}
+
+	for n := 0; n < b.N; n++ {
+		kCurveOps.ScalarMul3Pt(&curve, &threePointLadderInputs[0], &threePointLadderInputs[1], &threePointLadderInputs[2], uint(len(mScalarBytes)*8), mScalarBytes[:])
+	}
+}

p751toolbox/field_amd64_test.go → p751/field_amd64_test.go

@@ -1,11 +1,12 @@
 // +build amd64,!noasm
 
-package p751toolbox
+package p751
 
 import (
 	"golang.org/x/sys/cpu"
 	"testing"
 	"testing/quick"
+    . "github.com/cloudflare/p751sidh/internal/isogeny"
 )
 
 func TestFp751MontgomeryReduce(t *testing.T) {
@@ -23,9 +24,9 @@ func TestFp751MontgomeryReduce(t *testing.T) {
 
 	// Also check that the BMI2 implementation produces the same results
 	// as the fallback implementation.
-	compareMontgomeryReduce := func(x, y PrimeFieldElement) bool {
-		var z, zbackup fp751X2
-		var zred1, zred2 Fp751Element
+	compareMontgomeryReduce := func(x, y primeFieldElement) bool {
+		var z, zbackup FpElementX2
+		var zred1, zred2 FpElement
 
 		fp751Mul(&z, &x.A, &y.A)
 		zbackup = z
@@ -52,9 +53,9 @@ func TestFp751MontgomeryReduce(t *testing.T) {
 	// Check that the BMI2ADX implementation produces the same results as
 	// the BMI2 implementation. By transitivity, it should also produce the
 	// same results as the fallback implementation.
-	compareMontgomeryReduce = func(x, y PrimeFieldElement) bool {
-		var z, zbackup fp751X2
-		var zred1, zred2 Fp751Element
+	compareMontgomeryReduce = func(x, y primeFieldElement) bool {
+		var z, zbackup FpElementX2
+		var zred1, zred2 FpElement
 
 		fp751Mul(&z, &x.A, &y.A)
 		zbackup = z

p751/field_ops.go

@@ -0,0 +1,273 @@
+package p751
+
+import . "github.com/cloudflare/p751sidh/internal/isogeny"
+
+// 2*p751
+var (
+	// p751
+	p751 = FpElement{
+		0xffffffffffffffff, 0xffffffffffffffff, 0xffffffffffffffff,
+		0xffffffffffffffff, 0xffffffffffffffff, 0xeeafffffffffffff,
+		0xe3ec968549f878a8, 0xda959b1a13f7cc76, 0x084e9867d6ebe876,
+		0x8562b5045cb25748, 0x0e12909f97badc66, 0x00006fe5d541f71c}
+	// p751 + 1
+	p751p1 = FpElement{
+		0x0000000000000000, 0x0000000000000000, 0x0000000000000000,
+		0x0000000000000000, 0x0000000000000000, 0xeeb0000000000000,
+		0xe3ec968549f878a8, 0xda959b1a13f7cc76, 0x084e9867d6ebe876,
+		0x8562b5045cb25748, 0x0e12909f97badc66, 0x00006fe5d541f71c}
+	// 2*p751
+	p751x2 = FpElement{
+		0xFFFFFFFFFFFFFFFE, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF,
+		0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xDD5FFFFFFFFFFFFF,
+		0xC7D92D0A93F0F151, 0xB52B363427EF98ED, 0x109D30CFADD7D0ED,
+		0x0AC56A08B964AE90, 0x1C25213F2F75B8CD, 0x0000DFCBAA83EE38}
+)
+
+//------------------------------------------------------------------------------
+// Implementtaion of FieldOperations
+//------------------------------------------------------------------------------
+
+// Implements FieldOps
+type fp751Ops struct{}
+
+func FieldOperations() FieldOps {
+	return &fp751Ops{}
+}
+
+func (fp751Ops) Add(dest, lhs, rhs *Fp2Element) {
+	fp751AddReduced(&dest.A, &lhs.A, &rhs.A)
+	fp751AddReduced(&dest.B, &lhs.B, &rhs.B)
+}
+
+func (fp751Ops) Sub(dest, lhs, rhs *Fp2Element) {
+	fp751SubReduced(&dest.A, &lhs.A, &rhs.A)
+	fp751SubReduced(&dest.B, &lhs.B, &rhs.B)
+}
+
+func (fp751Ops) Mul(dest, lhs, rhs *Fp2Element) {
+	// Let (a,b,c,d) = (lhs.a,lhs.b,rhs.a,rhs.b).
+	a := &lhs.A
+	b := &lhs.B
+	c := &rhs.A
+	d := &rhs.B
+
+	// We want to compute
+	//
+	// (a + bi)*(c + di) = (a*c - b*d) + (a*d + b*c)i
+	//
+	// Use Karatsuba's trick: note that
+	//
+	// (b - a)*(c - d) = (b*c + a*d) - a*c - b*d
+	//
+	// so (a*d + b*c) = (b-a)*(c-d) + a*c + b*d.
+
+	var ac, bd FpElementX2
+	fp751Mul(&ac, a, c) // = a*c*R*R
+	fp751Mul(&bd, b, d) // = b*d*R*R
+
+	var b_minus_a, c_minus_d FpElement
+	fp751SubReduced(&b_minus_a, b, a) // = (b-a)*R
+	fp751SubReduced(&c_minus_d, c, d) // = (c-d)*R
+
+	var ad_plus_bc FpElementX2
+	fp751Mul(&ad_plus_bc, &b_minus_a, &c_minus_d) // = (b-a)*(c-d)*R*R
+	fp751X2AddLazy(&ad_plus_bc, &ad_plus_bc, &ac) // = ((b-a)*(c-d) + a*c)*R*R
+	fp751X2AddLazy(&ad_plus_bc, &ad_plus_bc, &bd) // = ((b-a)*(c-d) + a*c + b*d)*R*R
+
+	fp751MontgomeryReduce(&dest.B, &ad_plus_bc) // = (a*d + b*c)*R mod p
+
+	var ac_minus_bd FpElementX2
+	fp751X2SubLazy(&ac_minus_bd, &ac, &bd)       // = (a*c - b*d)*R*R
+	fp751MontgomeryReduce(&dest.A, &ac_minus_bd) // = (a*c - b*d)*R mod p
+}
+
+func (fp751Ops) Square(dest, x *Fp2Element) {
+	a := &x.A
+	b := &x.B
+
+	// We want to compute
+	//
+	// (a + bi)*(a + bi) = (a^2 - b^2) + 2abi.
+
+	var a2, a_plus_b, a_minus_b FpElement
+	fp751AddReduced(&a2, a, a)        // = a*R + a*R = 2*a*R
+	fp751AddReduced(&a_plus_b, a, b)  // = a*R + b*R = (a+b)*R
+	fp751SubReduced(&a_minus_b, a, b) // = a*R - b*R = (a-b)*R
+
+	var asq_minus_bsq, ab2 FpElementX2
+	fp751Mul(&asq_minus_bsq, &a_plus_b, &a_minus_b) // = (a+b)*(a-b)*R*R = (a^2 - b^2)*R*R
+	fp751Mul(&ab2, &a2, b)                          // = 2*a*b*R*R
+
+	fp751MontgomeryReduce(&dest.A, &asq_minus_bsq) // = (a^2 - b^2)*R mod p
+	fp751MontgomeryReduce(&dest.B, &ab2)           // = 2*a*b*R mod p
+}
+
+// Set dest = 1/x
+//
+// Allowed to overlap dest with x.
+//
+// Returns dest to allow chaining operations.
+func (fp751Ops) Inv(dest, x *Fp2Element) {
+	a := &x.A
+	b := &x.B
+
+	// We want to compute
+	//
+	//    1          1     (a - bi)	    (a - bi)
+	// -------- = -------- -------- = -----------
+	// (a + bi)   (a + bi) (a - bi)   (a^2 + b^2)
+	//
+	// Letting c = 1/(a^2 + b^2), this is
+	//
+	// 1/(a+bi) = a*c - b*ci.
+
+	var asq_plus_bsq primeFieldElement
+	var asq, bsq FpElementX2
+	fp751Mul(&asq, a, a)                         // = a*a*R*R
+	fp751Mul(&bsq, b, b)                         // = b*b*R*R
+	fp751X2AddLazy(&asq, &asq, &bsq)             // = (a^2 + b^2)*R*R
+	fp751MontgomeryReduce(&asq_plus_bsq.A, &asq) // = (a^2 + b^2)*R mod p
+	// Now asq_plus_bsq = a^2 + b^2
+
+	// Invert asq_plus_bsq
+	inv := asq_plus_bsq
+	inv.Mul(&asq_plus_bsq, &asq_plus_bsq)
+	inv.P34(&inv)
+	inv.Mul(&inv, &inv)
+	inv.Mul(&inv, &asq_plus_bsq)
+
+	var ac FpElementX2
+	fp751Mul(&ac, a, &inv.A)
+	fp751MontgomeryReduce(&dest.A, &ac)
+
+	var minus_b FpElement
+	fp751SubReduced(&minus_b, &minus_b, b)
+	var minus_bc FpElementX2
+	fp751Mul(&minus_bc, &minus_b, &inv.A)
+	fp751MontgomeryReduce(&dest.B, &minus_bc)
+}
+
+// In case choice == 1, performs following swap in constant time:
+// 	xPx <-> xQx
+//	xPz <-> xQz
+// Otherwise returns xPx, xPz, xQx, xQz unchanged
+func (fp751Ops) CondSwap(xPx, xPz, xQx, xQz *Fp2Element, choice uint8) {
+	fp751ConditionalSwap(&xPx.A, &xQx.A, choice)
+	fp751ConditionalSwap(&xPx.B, &xQx.B, choice)
+	fp751ConditionalSwap(&xPz.A, &xQz.A, choice)
+	fp751ConditionalSwap(&xPz.B, &xQz.B, choice)
+}
+
+// Converts values in x.A and x.B to Montgomery domain
+// x.A = x.A * R mod p
+// x.B = x.B * R mod p
+func (fp751Ops) ToMontgomery(x *Fp2Element) {
+	var aRR FpElementX2
+
+	// convert to montgomery domain
+	fp751Mul(&aRR, &x.A, &montgomeryRsq)   // = a*R*R
+	fp751MontgomeryReduce(&x.A, &aRR)      // = a*R mod p
+	fp751Mul(&aRR, &x.B, &montgomeryRsq)
+	fp751MontgomeryReduce(&x.B, &aRR)
+}
+
+// Converts values in x.A and x.B from Montgomery domain
+// a = x.A mod p
+// b = x.B mod p
+//
+// After returning from the call x is not modified.
+func (fp751Ops) FromMontgomery(x *Fp2Element, out *Fp2Element) {
+	var aR FpElementX2
+
+	// convert from montgomery domain
+	copy(aR[:], x.A[:])
+	fp751MontgomeryReduce(&out.A, &aR)    // = a mod p in [0, 2p)
+	fp751StrongReduce(&out.A)             // = a mod p in [0, p)
+	for i:=range(aR) {
+		aR[i] = 0
+	}
+	copy(aR[:], x.B[:])
+	fp751MontgomeryReduce(&out.B, &aR)
+	fp751StrongReduce(&out.B)
+}
+
+//------------------------------------------------------------------------------
+// Prime Field
+//------------------------------------------------------------------------------
+
+// Represents an element of the prime field F_p in Montgomery domain
+type primeFieldElement struct {
+	// The value `A`is represented by `aR mod p`.
+	A FpElement
+}
+
+// Set dest = lhs * rhs.
+//
+// Allowed to overlap lhs or rhs with dest.
+//
+// Returns dest to allow chaining operations.
+func (dest *primeFieldElement) Mul(lhs, rhs *primeFieldElement) *primeFieldElement {
+	a := &lhs.A // = a*R
+	b := &rhs.A // = b*R
+
+	var ab FpElementX2
+	fp751Mul(&ab, a, b)                 // = a*b*R*R
+	fp751MontgomeryReduce(&dest.A, &ab) // = a*b*R mod p
+
+	return dest
+}
+
+// Set dest = x^(2^k), for k >= 1, by repeated squarings.
+//
+// Allowed to overlap x with dest.
+//
+// Returns dest to allow chaining operations.
+func (dest *primeFieldElement) Pow2k(x *primeFieldElement, k uint8) *primeFieldElement {
+	dest.Mul(x, x)
+	for i := uint8(1); i < k; i++ {
+		dest.Mul(dest, dest)
+	}
+
+	return dest
+}
+
+// Set dest = x^((p-3)/4).  If x is square, this is 1/sqrt(x).
+//
+// Allowed to overlap x with dest.
+//
+// Returns dest to allow chaining operations.
+func (dest *primeFieldElement) P34(x *primeFieldElement) *primeFieldElement {
+	// Sliding-window strategy computed with Sage, awk, sed, and tr.
+	//
+	// This performs sum(powStrategy) = 744 squarings and len(mulStrategy)
+	// = 137 multiplications, in addition to 1 squaring and 15
+	// multiplications to build a lookup table.
+	//
+	// In total this is 745 squarings, 152 multiplications.  Since squaring
+	// is not implemented for the prime field, this is 897 multiplications
+	// in total.
+	powStrategy := [137]uint8{5, 7, 6, 2, 10, 4, 6, 9, 8, 5, 9, 4, 7, 5, 5, 4, 8, 3, 9, 5, 5, 4, 10, 4, 6, 6, 6, 5, 8, 9, 3, 4, 9, 4, 5, 6, 6, 2, 9, 4, 5, 5, 5, 7, 7, 9, 4, 6, 4, 8, 5, 8, 6, 6, 2, 9, 7, 4, 8, 8, 8, 4, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 2}
+	mulStrategy := [137]uint8{31, 23, 21, 1, 31, 7, 7, 7, 9, 9, 19, 15, 23, 23, 11, 7, 25, 5, 21, 17, 11, 5, 17, 7, 11, 9, 23, 9, 1, 19, 5, 3, 25, 15, 11, 29, 31, 1, 29, 11, 13, 9, 11, 27, 13, 19, 15, 31, 3, 29, 23, 31, 25, 11, 1, 21, 19, 15, 15, 21, 29, 13, 23, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 3}
+	initialMul := uint8(27)
+
+	// Build a lookup table of odd multiples of x.
+	lookup := [16]primeFieldElement{}
+	xx := &primeFieldElement{}
+	xx.Mul(x, x) // Set xx = x^2
+	lookup[0] = *x
+	for i := 1; i < 16; i++ {
+		lookup[i].Mul(&lookup[i-1], xx)
+	}
+	// Now lookup = {x, x^3, x^5, ... }
+	// so that lookup[i] = x^{2*i + 1}
+	// so that lookup[k/2] = x^k, for odd k
+
+	*dest = lookup[initialMul/2]
+	for i := uint8(0); i < 137; i++ {
+		dest.Pow2k(dest, powStrategy[i])
+		dest.Mul(dest, &lookup[mulStrategy[i]/2])
+	}
+
+	return dest
+}

p751/isogeny_test.go

@@ -0,0 +1,129 @@
+package p751
+
+import (
+	"testing"
+	. "github.com/cloudflare/p751sidh/internal/isogeny"
+)
+
+func TestFourIsogenyVersusSage(t *testing.T) {
+	var xR, xP4, resPhiXr, expPhiXr ProjectivePoint
+	var phi = Newisogeny4(kFieldOps)
+
+	// sage: p = 2^372 * 3^239 - 1; Fp = GF(p)
+	// sage: R.<x> = Fp[]
+	// sage: Fp2 = Fp.extension(x^2 + 1, 'i')
+	// sage: i = Fp2.gen()
+	// sage: E0Fp = EllipticCurve(Fp, [0,0,0,1,0])
+	// sage: E0Fp2 = EllipticCurve(Fp2, [0,0,0,1,0])
+	// sage: x_PA = 11
+	// sage: y_PA = -Fp(11^3 + 11).sqrt()
+	// sage: x_PB = 6
+	// sage: y_PB = -Fp(6^3 + 6).sqrt()
+	// sage: P_A = 3^239 * E0Fp((x_PA,y_PA))
+	// sage: P_B = 2^372 * E0Fp((x_PB,y_PB))
+	// sage: def tau(P):
+	// ....:     return E0Fp2( (-P.xy()[0], i*P.xy()[1]))
+	// ....:
+	// sage: m_B = 3*randint(0,3^238)
+	// sage: m_A = 2*randint(0,2^371)
+	// sage: R_A = E0Fp2(P_A) + m_A*tau(P_A)
+	// sage: def y_recover(x, a):
+	// ....:     return (x**3 + a*x**2 + x).sqrt()
+	// ....:
+	// sage: first_4_torsion_point = E0Fp2(1, y_recover(Fp2(1),0))
+	// sage: sage_first_4_isogeny = E0Fp2.isogeny(first_4_torsion_point)
+	// sage: a = Fp2(0)
+	// sage: E1A = EllipticCurve(Fp2, [0,(2*(a+6))/(a-2),0,1,0])
+	// sage: sage_isomorphism = sage_first_4_isogeny.codomain().isomorphism_to(E1A)
+	// sage: isogenized_R_A = sage_isomorphism(sage_first_4_isogeny(R_A))
+	// sage: P_4 = (2**(372-4))*isogenized_R_A
+	// sage: P_4._order = 4 #otherwise falls back to generic group methods for order
+	// sage: X4, Z4 = P_4.xy()[0], 1
+	// sage: phi4 = EllipticCurveIsogeny(E1A, P_4, None, 4)
+	// sage: E2A_sage = phi4.codomain() # not in monty form
+	// sage: Aprime, Cprime = 2*(2*X4^4 - Z4^4), Z4^4
+	// sage: E2A = EllipticCurve(Fp2, [0,Aprime/Cprime,0,1,0])
+	// sage: sage_iso = E2A_sage.isomorphism_to(E2A)
+	// sage: isogenized2_R_A = sage_iso(phi4(isogenized_R_A))
+
+	xP4 = ProjectivePoint{
+		X: Fp2Element{
+			A: FpElement{0x2afd75a913f3d5e7, 0x2918fba06f88c9ab, 0xa4ac4dc7cb526f05, 0x2d19e9391a607300, 0x7a79e2b34091b54, 0x3ad809dcb42f1792, 0xd46179328bd6402a, 0x1afa73541e2c4f3f, 0xf602d73ace9bdbd8, 0xd77ac58f6bab7004, 0x4689d97f6793b3b3, 0x4f26b00e42b7},
+			B: FpElement{0x6cdf918dafdcb890, 0x666f273cc29cfae2, 0xad00fcd31ba618e2, 0x5fbcf62bef2f6a33, 0xf408bb88318e5098, 0x84ab97849453d175, 0x501bbfcdcfb8e1ac, 0xf2370098e6b5542c, 0xc7dc73f5f0f6bd32, 0xdd76dcd86729d1cf, 0xca22c905029996e4, 0x5cf4a9373de3}},
+		Z: P751_OneFp2}
+	xR = ProjectivePoint{
+		X: Fp2Element{
+			A: FpElement{0xff99e76f78da1e05, 0xdaa36bd2bb8d97c4, 0xb4328cee0a409daf, 0xc28b099980c5da3f, 0xf2d7cd15cfebb852, 0x1935103dded6cdef, 0xade81528de1429c3, 0x6775b0fa90a64319, 0x25f89817ee52485d, 0x706e2d00848e697, 0xc4958ec4216d65c0, 0xc519681417f},
+			B: FpElement{0x742fe7dde60e1fb9, 0x801a3c78466a456b, 0xa9f945b786f48c35, 0x20ce89e1b144348f, 0xf633970b7776217e, 0x4c6077a9b38976e5, 0x34a513fc766c7825, 0xacccba359b9cd65, 0xd0ca8383f0fd0125, 0x77350437196287a, 0x9fe1ad7706d4ea21, 0x4d26129ee42d}},
+		Z: P751_OneFp2}
+	expPhiXr = ProjectivePoint{
+		X: Fp2Element{
+			A: FpElement{0x111efd8bd0b7a01e, 0x6ab75a4f3789ca9b, 0x939dbe518564cac4, 0xf9eeaba1601d0434, 0x8d41f8ba6edac998, 0xfcd2557efe9aa170, 0xb3c3549c098b7844, 0x52874fef6f81127c, 0xb2b9ac82aa518bb3, 0xee70820230520a86, 0xd4012b7f5efb184a, 0x573e4536329b},
+			B: FpElement{0xa99952281e932902, 0x569a89a571f2c7b1, 0x6150143846ba3f6b, 0x11fd204441e91430, 0x7f469bd55c9b07b, 0xb72db8b9de35b161, 0x455a9a37a940512a, 0xb0cff7670abaf906, 0x18c785b7583375fe, 0x603ab9ca403c9148, 0xab54ba3a6e6c62c1, 0x2726d7d57c4f}},
+		Z: P751_OneFp2}
+
+	phi.GenerateCurve(&xP4)
+	resPhiXr = phi.EvaluatePoint(&xR)
+	if !VartimeEqProjFp2(&expPhiXr, &resPhiXr) {
+		t.Error("\nExpected\n", expPhiXr.ToAffine(kCurveOps), "\nfound\n", resPhiXr.ToAffine(kCurveOps))
+	}
+}
+
+func TestThreeIsogenyVersusSage(t *testing.T) {
+	var xR, xP3, resPhiXr, expPhiXr ProjectivePoint
+	var phi = Newisogeny3(kFieldOps)
+
+	// sage: %colors Linux
+	// sage: p = 2^372 * 3^239 - 1; Fp = GF(p)
+	// sage: R.<x> = Fp[]
+	// sage: Fp2 = Fp.extension(x^2 + 1, 'i')
+	// sage: i = Fp2.gen()
+	// sage: E0Fp = EllipticCurve(Fp, [0,0,0,1,0])
+	// sage: E0Fp2 = EllipticCurve(Fp2, [0,0,0,1,0])
+	// sage: x_PA = 11
+	// sage: y_PA = -Fp(11^3 + 11).sqrt()
+	// sage: x_PB = 6
+	// sage: y_PB = -Fp(6^3 + 6).sqrt()
+	// sage: P_A = 3^239 * E0Fp((x_PA,y_PA))
+	// sage: P_B = 2^372 * E0Fp((x_PB,y_PB))
+	// sage: def tau(P):
+	// ....:     return E0Fp2( (-P.xy()[0], i*P.xy()[1]))
+	// ....:
+	// sage: m_B = 3*randint(0,3^238)
+	// sage: R_B = E0Fp2(P_B) + m_B*tau(P_B)
+	// sage: P_3 = (3^238)*R_B
+	// sage: def three_isog(P_3, P):
+	// ....:     X3, Z3 = P_3.xy()[0], 1
+	// ....:     XP, ZP = P.xy()[0], 1
+	// ....:     x = (XP*(X3*XP - Z3*ZP)^2)/(ZP*(Z3*XP - X3*ZP)^2)
+	// ....:     A3, C3 = (Z3^4 + 9*X3^2*(2*Z3^2 - 3*X3^2)), 4*X3*Z3^3
+	// ....:     cod = EllipticCurve(Fp2, [0,A3/C3,0,1,0])
+	// ....:     return cod.lift_x(x)
+	// ....:
+	// sage: isogenized_R_B = three_isog(P_3, R_B)
+
+	xR = ProjectivePoint{
+		X: Fp2Element{
+			A: FpElement{0xbd0737ed5cc9a3d7, 0x45ae6d476517c101, 0x6f228e9e7364fdb2, 0xbba4871225b3dbd, 0x6299ccd2e5da1a07, 0x38488fe4af5f2d0e, 0xec23cae5a86e980c, 0x26c804ba3f1edffa, 0xfbbed81932df60e5, 0x7e00e9d182ae9187, 0xc7654abb66d05f4b, 0x262d0567237b},
+			B: FpElement{0x3a3b5b6ad0b2ac33, 0x246602b5179127d3, 0x502ae0e9ad65077d, 0x10a3a37237e1bf70, 0x4a1ab9294dd05610, 0xb0f3adac30fe1fa6, 0x341995267faf70cb, 0xa14dd94d39cf4ec1, 0xce4b7527d1bf5568, 0xe0410423ed45c7e4, 0x38011809b6425686, 0x28f52472ebed}},
+		Z: P751_OneFp2}
+
+	xP3 = ProjectivePoint{
+		X: Fp2Element{
+			A: FpElement{0x7bb7a4a07b0788dc, 0xdc36a3f6607b21b0, 0x4750e18ee74cf2f0, 0x464e319d0b7ab806, 0xc25aa44c04f758ff, 0x392e8521a46e0a68, 0xfc4e76b63eff37df, 0x1f3566d892e67dd8, 0xf8d2eb0f73295e65, 0x457b13ebc470bccb, 0xfda1cc9efef5be33, 0x5dbf3d92cc02},
+			B: FpElement{0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0}},
+		Z: P751_OneFp2}
+
+	expPhiXr = ProjectivePoint{
+		X: Fp2Element{
+			A: FpElement{0x286db7d75913c5b1, 0xcb2049ad50189220, 0xccee90ef765fa9f4, 0x65e52ce2730e7d88, 0xa6b6b553bd0d06e7, 0xb561ecec14591590, 0x17b7a66d8c64d959, 0x77778cecbe1461e, 0x9405c9c0c41a57ce, 0x8f6b4847e8ca7d3d, 0xf625eb987b366937, 0x421b3590e345},
+			B: FpElement{0x566b893803e7d8d6, 0xe8c71a04d527e696, 0x5a1d8f87bf5eb51, 0x42ae08ae098724f, 0x4ee3d7c7af40ca2e, 0xd9f9ab9067bb10a7, 0xecd53d69edd6328c, 0xa581e9202dea107d, 0x8bcdfb6c8ecf9257, 0xe7cbbc2e5cbcf2af, 0x5f031a8701f0e53e, 0x18312d93e3cb}},
+		Z: P751_OneFp2}
+
+	phi.GenerateCurve(&xP3)
+	resPhiXr = phi.EvaluatePoint(&xR)
+
+	if !VartimeEqProjFp2(&expPhiXr, &resPhiXr) {
+		t.Error("\nExpected\n", expPhiXr.ToAffine(kCurveOps), "\nfound\n", resPhiXr.ToAffine(kCurveOps))
+	}
+}

p751/utils_test.go

@@ -0,0 +1,341 @@
+package p751
+
+// Tools used for testing and debugging
+
+import (
+	"fmt"
+    "math/big"
+    "math/rand"
+    "reflect"
+    "testing/quick"
+    . "github.com/cloudflare/p751sidh/internal/isogeny"
+)
+
+/* -------------------------------------------------------------------------
+   Underlying field configuration
+   -------------------------------------------------------------------------*/
+var (
+    kFieldOps = FieldOperations()
+    kParams   = &SidhParams{
+        Op: kFieldOps,
+        OneFp2: P751_OneFp2,
+        HalfFp2: P751_HalfFp2,
+        Bytelen: P751_Bytelen,
+    }
+    kCurveOps = &CurveOperations{Params: kParams}
+)
+
+/* -------------------------------------------------------------------------
+   Configure testing/quick
+   -------------------------------------------------------------------------*/
+var (
+    quickCheckScaleFactor = uint8(3)
+    quickCheckConfig = &quick.Config{MaxCount: (1 << (12 + quickCheckScaleFactor))}
+)
+
+/* -------------------------------------------------------------------------
+   Structure used by tests
+   -------------------------------------------------------------------------*/
+type GeneratedTestParams struct {
+    Point   ProjectivePoint
+    Cparam  ProjectiveCurveParameters
+    ExtElem Fp2Element
+}
+
+/* -------------------------------------------------------------------------
+  Test values
+   -------------------------------------------------------------------------*/
+
+// A = 4385300808024233870220415655826946795549183378139271271040522089756750951667981765872679172832050962894122367066234419550072004266298327417513857609747116903999863022476533671840646615759860564818837299058134292387429068536219*i + 1408083354499944307008104531475821995920666351413327060806684084512082259107262519686546161682384352696826343970108773343853651664489352092568012759783386151707999371397181344707721407830640876552312524779901115054295865393760
+var (
+    curve_A = Fp2Element{
+        A: FpElement{0x8319eb18ca2c435e, 0x3a93beae72cd0267, 0x5e465e1f72fd5a84, 0x8617fa4150aa7272, 0x887da24799d62a13, 0xb079b31b3c7667fe, 0xc4661b150fa14f2e, 0xd4d2b2967bc6efd6, 0x854215a8b7239003, 0x61c5302ccba656c2, 0xf93194a27d6f97a2, 0x1ed9532bca75},
+        B: FpElement{0xb6f541040e8c7db6, 0x99403e7365342e15, 0x457e9cee7c29cced, 0x8ece72dc073b1d67, 0x6e73cef17ad28d28, 0x7aed836ca317472, 0x89e1de9454263b54, 0x745329277aa0071b, 0xf623dfc73bc86b9b, 0xb8e3c1d8a9245882, 0x6ad0b3d317770bec, 0x5b406e8d502b}}
+
+    // C = 933177602672972392833143808100058748100491911694554386487433154761658932801917030685312352302083870852688835968069519091048283111836766101703759957146191882367397129269726925521881467635358356591977198680477382414690421049768*i + 9088894745865170214288643088620446862479558967886622582768682946704447519087179261631044546285104919696820250567182021319063155067584445633834024992188567423889559216759336548208016316396859149888322907914724065641454773776307
+    curve_C = Fp2Element{
+        A: FpElement{0x4fb2358bbf723107, 0x3a791521ac79e240, 0x283e24ef7c4c922f, 0xc89baa1205e33cc, 0x3031be81cff6fee1, 0xaf7a494a2f6a95c4, 0x248d251eaac83a1d, 0xc122fca1e2550c88, 0xbc0451b11b6cfd3d, 0x9c0a114ab046222c, 0x43b957b32f21f6ea, 0x5b9c87fa61de},
+        B: FpElement{0xacf142afaac15ec6, 0xfd1322a504a071d5, 0x56bb205e10f6c5c6, 0xe204d2849a97b9bd, 0x40b0122202fe7f2e, 0xecf72c6fafacf2cb, 0x45dfc681f869f60a, 0x11814c9aff4af66c, 0x9278b0c4eea54fe7, 0x9a633d5baf7f2e2e, 0x69a329e6f1a05112, 0x1d874ace23e4}}
+
+
+    // x(P) = 8172151271761071554796221948801462094972242987811852753144865524899433583596839357223411088919388342364651632180452081960511516040935428737829624206426287774255114241789158000915683252363913079335550843837650671094705509470594*i + 9326574858039944121604015439381720195556183422719505497448541073272720545047742235526963773359004021838961919129020087515274115525812121436661025030481584576474033630899768377131534320053412545346268645085054880212827284581557
+    affine_xP = Fp2Element{
+        A: FpElement{0xe8d05f30aac47247, 0x576ec00c55441de7, 0xbf1a8ec5fe558518, 0xd77cb17f77515881, 0x8e9852837ee73ec4, 0x8159634ad4f44a6b, 0x2e4eb5533a798c5, 0x9be8c4354d5bc849, 0xf47dc61806496b84, 0x25d0e130295120e0, 0xdbef54095f8139e3, 0x5a724f20862c},
+        B: FpElement{0x3ca30d7623602e30, 0xfb281eddf45f07b7, 0xd2bf62d5901a45bc, 0xc67c9baf86306dd2, 0x4e2bd93093f538ca, 0xcfd92075c25b9cbe, 0xceafe9a3095bcbab, 0x7d928ad380c85414, 0x37c5f38b2afdc095, 0x75325899a7b779f4, 0xf130568249f20fdd, 0x178f264767d1}}
+
+    // x([2]P) = 1476586462090705633631615225226507185986710728845281579274759750260315746890216330325246185232948298241128541272709769576682305216876843626191069809810990267291824247158062860010264352034514805065784938198193493333201179504845*i + 3623708673253635214546781153561465284135688791018117615357700171724097420944592557655719832228709144190233454198555848137097153934561706150196041331832421059972652530564323645509890008896574678228045006354394485640545367112224
+    affine_xP2 = Fp2Element{
+        A: FpElement{0x2a77afa8576ce979, 0xab1360e69b0aeba0, 0xd79e3e3cbffad660, 0x5fd0175aa10f106b, 0x1800ebafce9fbdbc, 0x228fc9142bdd6166, 0x867cf907314e34c3, 0xa58d18c94c13c31c, 0x699a5bc78b11499f, 0xa29fc29a01f7ccf1, 0x6c69c0c5347eebce, 0x38ecee0cc57},
+        B: FpElement{0x43607fd5f4837da0, 0x560bad4ce27f8f4a, 0x2164927f8495b4dd, 0x621103fdb831a997, 0xad740c4eea7db2db, 0x2cde0442205096cd, 0x2af51a70ede8324e, 0x41a4e680b9f3466, 0x5481f74660b8f476, 0xfcb2f3e656ff4d18, 0x42e3ce0837171acc, 0x44238c30530c}}
+
+    // x([3]P) = 9351941061182433396254169746041546943662317734130813745868897924918150043217746763025923323891372857734564353401396667570940585840576256269386471444236630417779544535291208627646172485976486155620044292287052393847140181703665*i + 9010417309438761934687053906541862978676948345305618417255296028956221117900864204687119686555681136336037659036201780543527957809743092793196559099050594959988453765829339642265399496041485088089691808244290286521100323250273
+    affine_xP3 = Fp2Element{
+        A: FpElement{0x2096e3f23feca947, 0xf36f635aa4ad8634, 0xdae3b1c6983c5e9a, 0xe08df6c262cb74b4, 0xd2ca4edc37452d3d, 0xfb5f3fe42f500c79, 0x73740aa3abc2b21f, 0xd535fd869f914cca, 0x4a558466823fb67f, 0x3e50a7a0e3bfc715, 0xf43c6da9183a132f, 0x61aca1e1b8b9},
+        B: FpElement{0x1e54ec26ea5077bd, 0x61380572d8769f9a, 0xc615170684f59818, 0x6309c3b93e84ef6e, 0x33c74b1318c3fcd0, 0xfe8d7956835afb14, 0x2d5a7b55423c1ecc, 0x869db67edfafea68, 0x1292632394f0a628, 0x10bba48225bfd141, 0x6466c28b408daba, 0x63cacfdb7c43}}
+
+    // x([2^2]P) = 441719501189485559222919502512761433931671682884872259563221427434901842337947564993718830905758163254463901652874331063768876314142359813382575876106725244985607032091781306919778265250690045578695338669105227100119314831452*i + 6961734028200975729170216310486458180126343885294922940439352055937945948015840788921225114530454649744697857047401608073256634790353321931728699534700109268264491160589480994022419317695690866764726967221310990488404411684053
+    affine_xP4 = Fp2Element{
+        A: FpElement{0x6f9dbe4c39175153, 0xf2fec757eb99e88, 0x43d7361a93733d91, 0x3abd10ed19c85a3d, 0xc4de9ab9c5ef7181, 0x53e375901684c900, 0x68ffc3e7d71c41ff, 0x47adab62c8d942fe, 0x226a33fd6fbb381d, 0x87ef4c8fdd83309a, 0xaca1cf44c5fa8799, 0x6cbae86c755f},
+        B: FpElement{0x4c80c37fe68282a7, 0xbd8b9d7248bf553a, 0x1fb0e8e74d5e1762, 0xb63fa0e4e5f91482, 0xc675ab8a45a1439, 0xdfa6772deace7820, 0xf0d813d71d9a9255, 0x53a1a58c634534bd, 0x4ebfc6485fdfd888, 0x6991fe4358bcf169, 0xc0547bdaca85b6fd, 0xf461548d632}}
+
+    // x([3^2]P) = 3957171963425208493644602380039721164492341594850197356580248639045894821895524981729970650520936632013218950972842867220898274664982599375786979902471523505057611521217523103474682939638645404445093536997296151472632038973463*i + 1357869545269286021642168835877253886774707209614159162748874474269328421720121175566245719916322684751967981171882659798149072149161259103020057556362998810229937432814792024248155991141511691087135859252304684633946087474060
+    affine_xP9 = Fp2Element{
+        A: FpElement{0x7c0daa0f04ded4e0, 0x52dc4f883d85e065, 0x91afbdc2c1714d0b, 0xb7b3db8e658cfeba, 0x43d4e72a692882f3, 0x535c56d83753da30, 0xc8a58724433cbf5d, 0x351153c0a5e74219, 0x2c81827d19f93dd5, 0x26ef8aca3370ea1a, 0x1cf939a6dd225dec, 0x3403cb28ad41},
+        B: FpElement{0x93e7bc373a9ff7b, 0x57b8cc47635ebc0f, 0x92eab55689106cf3, 0x93643111d421f24c, 0x1c58b519506f6b7a, 0xebd409fb998faa13, 0x5c86ed799d09d80e, 0xd9a1d764d6363562, 0xf95e87f92fb0c4cc, 0x6b2bbaf5632a5609, 0x2d9b6a809dfaff7f, 0x29c0460348b}}
+
+    // m = 96550223052359874398280314003345143371473380422728857598463622014420884224892
+    mScalarBytes = [...]uint8{0x7c, 0x7b, 0x95, 0xfa, 0xb4, 0x75, 0x6c, 0x48, 0x8c, 0x17, 0x55, 0xb4, 0x49, 0xf5, 0x1e, 0xa3, 0xb, 0x31, 0xf0, 0xa4, 0xa6, 0x81, 0xad, 0x94, 0x51, 0x11, 0xe7, 0xf5, 0x5b, 0x7d, 0x75, 0xd5}
+
+    // x([m]P) = 7893578558852400052689739833699289348717964559651707250677393044951777272628231794999463214496545377542328262828965953246725804301238040891993859185944339366910592967840967752138115122568615081881937109746463885908097382992642*i + 8293895847098220389503562888233557012043261770526854885191188476280014204211818299871679993460086974249554528517413590157845430186202704783785316202196966198176323445986064452630594623103149383929503089342736311904030571524837
+    affine_xaP = Fp2Element{
+        A: FpElement{0x2112f3c7d7f938bb, 0x704a677f0a4df08f, 0x825370e31fb4ef00, 0xddbf79b7469f902, 0x27640c899ea739fd, 0xfb7b8b19f244108e, 0x546a6679dd3baebc, 0xe9f0ecf398d5265f, 0x223d2b350e75e461, 0x84b322a0b6aff016, 0xfabe426f539f8b39, 0x4507a0604f50},
+        B: FpElement{0xac77737e5618a5fe, 0xf91c0e08c436ca52, 0xd124037bc323533c, 0xc9a772bf52c58b63, 0x3b30c8f38ef6af4d, 0xb9eed160e134f36e, 0x24e3836393b25017, 0xc828be1b11baf1d9, 0x7b7dab585df50e93, 0x1ca3852c618bd8e0, 0x4efa73bcb359fa00, 0x50b6a923c2d4}}
+
+    // Inputs for testing 3-point-ladder
+    threePointLadderInputs = []ProjectivePoint{
+        // x(P)
+        ProjectivePoint{
+            X: Fp2Element{
+                A: FpElement{0xe8d05f30aac47247, 0x576ec00c55441de7, 0xbf1a8ec5fe558518, 0xd77cb17f77515881, 0x8e9852837ee73ec4, 0x8159634ad4f44a6b, 0x2e4eb5533a798c5, 0x9be8c4354d5bc849, 0xf47dc61806496b84, 0x25d0e130295120e0, 0xdbef54095f8139e3, 0x5a724f20862c},
+                B: FpElement{0x3ca30d7623602e30, 0xfb281eddf45f07b7, 0xd2bf62d5901a45bc, 0xc67c9baf86306dd2, 0x4e2bd93093f538ca, 0xcfd92075c25b9cbe, 0xceafe9a3095bcbab, 0x7d928ad380c85414, 0x37c5f38b2afdc095, 0x75325899a7b779f4, 0xf130568249f20fdd, 0x178f264767d1}},
+            Z: P751_OneFp2,
+        },
+        // x(Q)
+        ProjectivePoint{
+            X: Fp2Element{
+                A: FpElement{0x2b71a2a93ad1e10e, 0xf0b9842a92cfb333, 0xae17373615a27f5c, 0x3039239f428330c4, 0xa0c4b735ed7dcf98, 0x6e359771ddf6af6a, 0xe986e4cac4584651, 0x8233a2b622d5518, 0xbfd67bf5f06b818b, 0xdffe38d0f5b966a6, 0xa86b36a3272ee00a, 0x193e2ea4f68f},
+                B: FpElement{0x5a0f396459d9d998, 0x479f42250b1b7dda, 0x4016b57e2a15bf75, 0xc59f915203fa3749, 0xd5f90257399cf8da, 0x1fb2dadfd86dcef4, 0x600f20e6429021dc, 0x17e347d380c57581, 0xc1b0d5fa8fe3e440, 0xbcf035330ac20e8, 0x50c2eb5f6a4f03e6, 0x86b7c4571}},
+            Z: P751_OneFp2,
+        },
+        // x(P-Q)
+        ProjectivePoint{
+            X: Fp2Element{
+                A: FpElement{0x4aafa9f378f7b5ff, 0x1172a683aa8eee0, 0xea518d8cbec2c1de, 0xe191bcbb63674557, 0x97bc19637b259011, 0xdbeae5c9f4a2e454, 0x78f64d1b72a42f95, 0xe71cb4ea7e181e54, 0xe4169d4c48543994, 0x6198c2286a98730f, 0xd21d675bbab1afa5, 0x2e7269fce391},
+                B: FpElement{0x23355783ce1d0450, 0x683164cf4ce3d93f, 0xae6d1c4d25970fd8, 0x7807007fb80b48cf, 0xa005a62ec2bbb8a2, 0x6b5649bd016004cb, 0xbb1a13fa1330176b, 0xbf38e51087660461, 0xe577fddc5dd7b930, 0x5f38116f56947cd3, 0x3124f30b98c36fde, 0x4ca9b6e6db37}},
+            Z: P751_OneFp2,
+        },
+    }
+    curve = ProjectiveCurveParameters{A: curve_A, C: curve_C}
+    cln16prime, _ = new(big.Int).SetString("10354717741769305252977768237866805321427389645549071170116189679054678940682478846502882896561066713624553211618840202385203911976522554393044160468771151816976706840078913334358399730952774926980235086850991501872665651576831", 10)
+)
+
+/* -------------------------------------------------------------------------
+   Values used by benchmarking tools
+   -------------------------------------------------------------------------*/
+
+// Package-level storage for this field element is intended to deter
+// compiler optimizations.
+var (
+    benchmarkFpElement FpElement
+    benchmarkFpElementX2 FpElementX2
+    bench_x = FpElement{17026702066521327207, 5108203422050077993, 10225396685796065916, 11153620995215874678, 6531160855165088358, 15302925148404145445, 1248821577836769963, 9789766903037985294, 7493111552032041328, 10838999828319306046, 18103257655515297935, 27403304611634}
+    bench_y = FpElement{4227467157325093378, 10699492810770426363, 13500940151395637365, 12966403950118934952, 16517692605450415877, 13647111148905630666, 14223628886152717087, 7167843152346903316, 15855377759596736571, 4300673881383687338, 6635288001920617779, 30486099554235}
+    bench_z = FpElementX2{1595347748594595712, 10854920567160033970, 16877102267020034574, 12435724995376660096, 3757940912203224231, 8251999420280413600, 3648859773438820227, 17622716832674727914, 11029567000887241528, 11216190007549447055, 17606662790980286987, 4720707159513626555, 12887743598335030915, 14954645239176589309, 14178817688915225254, 1191346797768989683, 12629157932334713723, 6348851952904485603, 16444232588597434895, 7809979927681678066, 14642637672942531613, 3092657597757640067, 10160361564485285723, 240071237}
+)
+
+// Helpers
+func (primeElement primeFieldElement) String() string {
+   b := toBigInt(&primeElement.A)
+   return fmt.Sprintf("%X", b.String())
+}
+
+/* -------------------------------------------------------------------------
+   Helpers
+   -------------------------------------------------------------------------*/
+
+// Given xP = x(P), xQ = x(Q), and xPmQ = x(P-Q), compute xR = x(P+Q).
+//
+// Returns xR to allow chaining.  Safe to overlap xP, xQ, xR.
+func AddProjFp2(xP, xQ, xPmQ *ProjectivePoint) ProjectivePoint {
+    // Algorithm 1 of Costello-Smith.
+    var v0, v1, v2, v3, v4 Fp2Element
+    var xR ProjectivePoint
+    kFieldOps.Add(&v0, &xP.X, &xP.Z)              // X_P + Z_P
+    kFieldOps.Sub(&v1, &xQ.X, &xQ.Z)
+    kFieldOps.Mul(&v1, &v1, &v0)                    // (X_Q - Z_Q)(X_P + Z_P)
+    kFieldOps.Sub(&v0, &xP.X, &xP.Z)              // X_P - Z_P
+    kFieldOps.Add(&v2, &xQ.X, &xQ.Z)
+    kFieldOps.Mul(&v2, &v2, &v0)                    // (X_Q + Z_Q)(X_P - Z_P)
+    kFieldOps.Add(&v3, &v1, &v2)
+    kFieldOps.Square(&v3, &v3)                      // 4(X_Q X_P - Z_Q Z_P)^2
+    kFieldOps.Sub(&v4, &v1, &v2)
+    kFieldOps.Square(&v4, &v4)                      // 4(X_Q Z_P - Z_Q X_P)^2
+    kFieldOps.Mul(&v0, &xPmQ.Z, &v3)              // 4X_{P-Q}(X_Q X_P - Z_Q Z_P)^2
+    kFieldOps.Mul(&xR.Z, &xPmQ.X, &v4)            // 4Z_{P-Q}(X_Q Z_P - Z_Q X_P)^2
+    xR.X = v0
+    return xR
+}
+
+// Given xP = x(P) and cached curve parameters Aplus2C = A + 2*C, C4 = 4*C,
+// compute xQ = x([2]P).
+//
+// Returns xQ to allow chaining.  Safe to overlap xP, xQ.
+func DoubleProjFp2(xP *ProjectivePoint, Aplus2C, C4 *Fp2Element) ProjectivePoint {
+    // Algorithm 2 of Costello-Smith, amended to work with projective curve coefficients.
+    var v1, v2, v3, xz4 Fp2Element
+    var xQ ProjectivePoint
+    kFieldOps.Add(&v1, &xP.X, &xP.Z)               // (X+Z)^2
+    kFieldOps.Square(&v1, &v1)
+    kFieldOps.Sub(&v2, &xP.X, &xP.Z)               // (X-Z)^2
+    kFieldOps.Square(&v2, &v2)
+    kFieldOps.Sub(&xz4, &v1, &v2)                // 4XZ = (X+Z)^2 - (X-Z)^2
+    kFieldOps.Mul(&v2, &v2, C4)                  // 4C(X-Z)^2
+    kFieldOps.Mul(&xQ.X, &v1, &v2)               // 4C(X+Z)^2(X-Z)^2
+    kFieldOps.Mul(&v3, &xz4, Aplus2C)            // 4XZ(A + 2C)
+    kFieldOps.Add(&v3, &v3, &v2)                 // 4XZ(A + 2C) + 4C(X-Z)^2
+    kFieldOps.Mul(&xQ.Z, &v3, &xz4)              // (4XZ(A + 2C) + 4C(X-Z)^2)4XZ
+    // Now (xQ.x : xQ.z)
+    //   = (4C(X+Z)^2(X-Z)^2 : (4XZ(A + 2C) + 4C(X-Z)^2)4XZ )
+    //   = ((X+Z)^2(X-Z)^2 : (4XZ((A + 2C)/4C) + (X-Z)^2)4XZ )
+    //   = ((X+Z)^2(X-Z)^2 : (4XZ((a + 2)/4) + (X-Z)^2)4XZ )
+    return xQ
+}
+
+// Given x(P) and a scalar m in little-endian bytes, compute x([m]P) using the
+// Montgomery ladder.  This is described in Algorithm 8 of Costello-Smith.
+//
+// This function's execution time is dependent only on the byte-length of the
+// input scalar.  All scalars of the same input length execute in uniform time.
+// The scalar can be padded with zero bytes to ensure a uniform length.
+//
+// Safe to overlap the source with the destination.
+func ScalarMult(curve *ProjectiveCurveParameters, xP *ProjectivePoint, scalar []uint8) ProjectivePoint {
+    var x0, x1, tmp ProjectivePoint
+    var Aplus2C, C4 Fp2Element
+
+    kFieldOps.Add(&Aplus2C, &curve.C, &curve.C) // = 2*C
+    kFieldOps.Add(&C4, &Aplus2C, &Aplus2C)      // = 4*C
+    kFieldOps.Add(&Aplus2C, &Aplus2C, &curve.A) // = 2*C + A
+
+    x0.X = P751_OneFp2
+    x1 = *xP
+
+    // Iterate over the bits of the scalar, top to bottom
+    prevBit := uint8(0)
+    for i := len(scalar) - 1; i >= 0; i-- {
+        scalarByte := scalar[i]
+        for j := 7; j >= 0; j-- {
+            bit := (scalarByte >> uint(j)) & 0x1
+            kCurveOps.Params.Op.CondSwap(&x0.X, &x0.Z, &x1.X, &x1.Z, (bit ^ prevBit))
+            //sProjectivePointConditionalSwap(&x0, &x1, (bit ^ prevBit))
+            tmp = DoubleProjFp2(&x0, &Aplus2C, &C4)
+            x1 = AddProjFp2(&x0, &x1, xP)
+            x0 = tmp
+            prevBit = bit
+        }
+    }
+    // now prevBit is the lowest bit of the scalar
+    kCurveOps.Params.Op.CondSwap(&x0.X, &x0.Z, &x1.X, &x1.Z, prevBit)
+    return x0
+}
+
+// Returns true if lhs = rhs.  Takes variable time.
+func VartimeEqFp2(lhs, rhs *Fp2Element) bool {
+    a := *lhs
+    b := *rhs
+
+    fp751StrongReduce(&a.A)
+    fp751StrongReduce(&a.B)
+    fp751StrongReduce(&b.A)
+    fp751StrongReduce(&b.B)
+
+    eq := true
+    for i := 0; i < len(a.A)&&eq; i++ {
+        eq = eq && (a.A[i] == b.A[i])
+        eq = eq && (a.B[i] == b.B[i])
+    }
+    return eq
+}
+
+// Returns true if lhs = rhs.  Takes variable time.
+func VartimeEqProjFp2(lhs, rhs *ProjectivePoint) bool {
+    var t0,t1 Fp2Element
+    kFieldOps.Mul(&t0, &lhs.X, &rhs.Z)
+    kFieldOps.Mul(&t1, &lhs.Z, &rhs.X)
+    return VartimeEqFp2(&t0, &t1)
+}
+
+func (GeneratedTestParams) generateFp2p751(rand *rand.Rand) Fp2Element {
+    // Generation strategy: low limbs taken from [0,2^64); high limb
+    // taken from smaller range
+    //
+    // Size hint is ignored since all elements are fixed size.
+    //
+    // Field elements taken in range [0,2p).  Emulate this by capping
+    // the high limb by the top digit of 2*p-1:
+    //
+    // sage: (2*p-1).digits(2^64)[-1]
+    // 246065832128056
+    //
+    // This still allows generating values >= 2p, but hopefully that
+    // excess is OK (and if it's not, we'll find out, because it's for
+    // testing...)
+    //
+    highLimb := rand.Uint64() % 246065832128056
+    fpElementGen := func() FpElement {
+        return FpElement{
+            rand.Uint64(),
+            rand.Uint64(),
+            rand.Uint64(),
+            rand.Uint64(),
+            rand.Uint64(),
+            rand.Uint64(),
+            rand.Uint64(),
+            rand.Uint64(),
+            rand.Uint64(),
+            rand.Uint64(),
+            rand.Uint64(),
+            highLimb,
+        }
+    }
+    return Fp2Element{A: fpElementGen(), B: fpElementGen()}
+}
+
+func (c GeneratedTestParams) Generate(rand *rand.Rand, size int) reflect.Value {
+    return reflect.ValueOf(
+        GeneratedTestParams{
+            ProjectivePoint{
+                X: c.generateFp2p751(rand),
+                Z: c.generateFp2p751(rand),
+            },
+            ProjectiveCurveParameters{
+                A: c.generateFp2p751(rand),
+                C: c.generateFp2p751(rand),
+            },
+            c.generateFp2p751(rand),
+        })
+}
+
+func (x primeFieldElement) Generate(rand *rand.Rand, size int) reflect.Value {
+       return reflect.ValueOf(primeFieldElement{
+        A: new(GeneratedTestParams).generateFp2p751(rand).A})
+}
+
+// Convert an FpElement to a big.Int for testing.  Because this is only
+// for testing, no big.Int to FpElement conversion is provided.
+func radix64ToBigInt(x []uint64) *big.Int {
+    radix := new(big.Int)
+    // 2^64
+    radix.UnmarshalText(([]byte)("18446744073709551616"))
+
+    base := new(big.Int).SetUint64(1)
+    val := new(big.Int).SetUint64(0)
+    tmp := new(big.Int)
+
+    for _, xi := range x {
+        tmp.SetUint64(xi)
+        tmp.Mul(tmp, base)
+        val.Add(val, tmp)
+        base.Mul(base, radix)
+    }
+
+    return val
+}
+
+// Converts number from Montgomery domain and returns big.Int
+func toBigInt(x *FpElement) *big.Int {
+    var fp Fp2Element
+    var in = Fp2Element{A:*x}
+    kFieldOps.FromMontgomery(&in, &fp)
+    return radix64ToBigInt(fp.A[:])
+}

p751toolbox/curve.go

@@ -1,293 +0,0 @@
-package p751toolbox
-
-// A point on the projective line P^1(F_{p^2}).
-//
-// This is used to work projectively with the curve coefficients.
-type ProjectiveCurveParameters struct {
-	A ExtensionFieldElement
-	C ExtensionFieldElement
-}
-
-// Stores curve projective parameters equivalent to A/C. Meaning of the
-// values depends on the context. When working with isogenies over
-// subgroup that are powers of:
-// * three then  (A:C) ~ (A+2C:A-2C)
-// * four then   (A:C) ~ (A+2C:  4C)
-// See Appendix A of SIKE for more details
-type CurveCoefficientsEquiv struct {
-	A ExtensionFieldElement
-	C ExtensionFieldElement
-}
-
-// A point on the projective line P^1(F_{p^2}).
-//
-// This represents a point on the Kummer line of a Montgomery curve.  The
-// curve is specified by a ProjectiveCurveParameters struct.
-type ProjectivePoint struct {
-	X ExtensionFieldElement
-	Z ExtensionFieldElement
-}
-
-func (params *ProjectiveCurveParameters) FromAffine(a *ExtensionFieldElement) {
-	params.A = *a
-	params.C.One()
-}
-
-// Computes j-invariant for a curve y2=x3+A/Cx+x with A,C in F_(p^2). Result
-// is returned in jBytes buffer, encoded in little-endian format. Caller
-// provided jBytes buffer has to be big enough to j-invariant value. In case
-// of SIDH, buffer size must be at least size of shared secret.
-// Implementation corresponds to Algorithm 9 from SIKE.
-func (cparams *ProjectiveCurveParameters) Jinvariant(jBytes []byte) {
-	var j, t0, t1 ExtensionFieldElement
-
-	j.Square(&cparams.A)  // j  = A^2
-	t1.Square(&cparams.C) // t1 = C^2
-	t0.Add(&t1, &t1)      // t0 = t1 + t1
-	t0.Sub(&j, &t0)       // t0 = j - t0
-	t0.Sub(&t0, &t1)      // t0 = t0 - t1
-	j.Sub(&t0, &t1)       // t0 = t0 - t1
-	t1.Square(&t1)        // t1 = t1^2
-	j.Mul(&j, &t1)        // t0 = t0 * t1
-	t0.Add(&t0, &t0)      // t0 = t0 + t0
-	t0.Add(&t0, &t0)      // t0 = t0 + t0
-	t1.Square(&t0)        // t1 = t0^2
-	t0.Mul(&t0, &t1)      // t0 = t0 * t1
-	t0.Add(&t0, &t0)      // t0 = t0 + t0
-	t0.Add(&t0, &t0)      // t0 = t0 + t0
-	j.Inv(&j)             // j  = 1/j
-	j.Mul(&t0, &j)        // j  = t0 * j
-
-	j.ToBytes(jBytes)
-}
-
-// Given affine points x(P), x(Q) and x(Q-P) in a extension field F_{p^2}, function
-// recorvers projective coordinate A of a curve. This is Algorithm 10 from SIKE.
-func (curve *ProjectiveCurveParameters) RecoverCoordinateA(xp, xq, xr *ExtensionFieldElement) {
-	var t0, t1 ExtensionFieldElement
-
-	t1.Add(xp, xq)                            // t1 = Xp + Xq
-	t0.Mul(xp, xq)                            // t0 = Xp * Xq
-	curve.A.Mul(xr, &t1)                      // A  = X(q-p) * t1
-	curve.A.Add(&curve.A, &t0)                // A  = A + t0
-	t0.Mul(&t0, xr)                           // t0 = t0 * X(q-p)
-	curve.A.Sub(&curve.A, &oneExtensionField) // A  = A - 1
-	t0.Add(&t0, &t0)                          // t0 = t0 + t0
-	t1.Add(&t1, xr)                           // t1 = t1 + X(q-p)
-	t0.Add(&t0, &t0)                          // t0 = t0 + t0
-	curve.A.Square(&curve.A)                  // A  = A^2
-	t0.Inv(&t0)                               // t0 = 1/t0
-	curve.A.Mul(&curve.A, &t0)                // A  = A * t0
-	curve.A.Sub(&curve.A, &t1)                // A  = A - t1
-}
-
-// Computes equivalence (A:C) ~ (A+2C : A-2C)
-func (curve *ProjectiveCurveParameters) CalcCurveParamsEquiv3() CurveCoefficientsEquiv {
-	var coef CurveCoefficientsEquiv
-	var c2 ExtensionFieldElement
-
-	c2.Add(&curve.C, &curve.C)
-	// A24p = A+2*C
-	coef.A.Add(&curve.A, &c2)
-	// A24m = A-2*C
-	coef.C.Sub(&curve.A, &c2)
-	return coef
-}
-
-// Computes equivalence (A:C) ~ (A+2C : 4C)
-func (cparams *ProjectiveCurveParameters) CalcCurveParamsEquiv4() CurveCoefficientsEquiv {
-	var coefEq CurveCoefficientsEquiv
-
-	coefEq.C.Add(&cparams.C, &cparams.C)
-	// A24p = A+2C
-	coefEq.A.Add(&cparams.A, &coefEq.C)
-	// C24 = 4*C
-	coefEq.C.Add(&coefEq.C, &coefEq.C)
-	return coefEq
-}
-
-// Helper function for RightToLeftLadder(). Returns A+2C / 4.
-func (cparams *ProjectiveCurveParameters) calcAplus2Over4() (ret ExtensionFieldElement) {
-	var tmp ExtensionFieldElement
-	// 2C
-	tmp.Add(&cparams.C, &cparams.C)
-	// A+2C
-	ret.Add(&cparams.A, &tmp)
-	// 1/4C
-	tmp.Add(&tmp, &tmp).Inv(&tmp)
-	// A+2C/4C
-	ret.Mul(&ret, &tmp)
-	return
-}
-
-// Recovers (A:C) curve parameters from projectively equivalent (A+2C:A-2C).
-func (cparams *ProjectiveCurveParameters) RecoverCurveCoefficients3(coefEq *CurveCoefficientsEquiv) {
-	cparams.A.Add(&coefEq.A, &coefEq.C)
-	// cparams.A = 2*(A+2C+A-2C) = 4A
-	cparams.A.Add(&cparams.A, &cparams.A)
-	// cparams.C = (A+2C-A+2C) = 4C
-	cparams.C.Sub(&coefEq.A, &coefEq.C)
-	return
-}
-
-// Recovers (A:C) curve parameters from projectively equivalent (A+2C:4C).
-func (cparams *ProjectiveCurveParameters) RecoverCurveCoefficients4(coefEq *CurveCoefficientsEquiv) {
-	var half = ExtensionFieldElement{
-		A: Fp751Element{
-			0x00000000000124D6, 0x0000000000000000, 0x0000000000000000,
-			0x0000000000000000, 0x0000000000000000, 0xB8E0000000000000,
-			0x9C8A2434C0AA7287, 0xA206996CA9A378A3, 0x6876280D41A41B52,
-			0xE903B49F175CE04F, 0x0F8511860666D227, 0x00004EA07CFF6E7F},
-		B: Fp751Element{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	}
-	// cparams.C = (4C)*1/2=2C
-	cparams.C.Mul(&coefEq.C, &half)
-	// cparams.A = A+2C - 2C = A
-	cparams.A.Sub(&coefEq.A, &cparams.C)
-	// cparams.C = 2C * 1/2 = C
-	cparams.C.Mul(&cparams.C, &half)
-	return
-}
-
-func (point *ProjectivePoint) FromAffine(x *ExtensionFieldElement) {
-	point.X = *x
-	point.Z = oneExtensionField
-}
-
-func (point *ProjectivePoint) ToAffine() *ExtensionFieldElement {
-	affine_x := new(ExtensionFieldElement)
-	affine_x.Inv(&point.Z).Mul(affine_x, &point.X)
-	return affine_x
-}
-
-func (lhs *ProjectivePoint) VartimeEq(rhs *ProjectivePoint) bool {
-	var t0, t1 ExtensionFieldElement
-	t0.Mul(&lhs.X, &rhs.Z)
-	t1.Mul(&lhs.Z, &rhs.X)
-	return t0.VartimeEq(&t1)
-}
-
-func ProjectivePointConditionalSwap(xP, xQ *ProjectivePoint, choice uint8) {
-	ExtensionFieldConditionalSwap(&xP.X, &xQ.X, choice)
-	ExtensionFieldConditionalSwap(&xP.Z, &xQ.Z, choice)
-}
-
-// Combined coordinate doubling and differential addition. Takes projective points
-// P,Q,Q-P and (A+2C)/4C curve E coefficient. Returns 2*P and P+Q calculated on E.
-// Function is used only by RightToLeftLadder. Corresponds to Algorithm 5 of SIKE
-func xDblAdd(P, Q, QmP *ProjectivePoint, a24 *ExtensionFieldElement) (dblP, PaQ ProjectivePoint) {
-	var t0, t1, t2 ExtensionFieldElement
-	xQmP, zQmP := &QmP.X, &QmP.Z
-	xPaQ, zPaQ := &PaQ.X, &PaQ.Z
-	x2P, z2P := &dblP.X, &dblP.Z
-	xP, zP := &P.X, &P.Z
-	xQ, zQ := &Q.X, &Q.Z
-
-	t0.Add(xP, zP)       // t0   = Xp+Zp
-	t1.Sub(xP, zP)       // t1   = Xp-Zp
-	x2P.Square(&t0)      // 2P.X = t0^2
-	t2.Sub(xQ, zQ)       // t2   = Xq-Zq
-	xPaQ.Add(xQ, zQ)     // Xp+q = Xq+Zq
-	t0.Mul(&t0, &t2)     // t0   = t0 * t2
-	z2P.Mul(&t1, &t1)    // 2P.Z = t1 * t1
-	t1.Mul(&t1, xPaQ)    // t1   = t1 * Xp+q
-	t2.Sub(x2P, z2P)     // t2   = 2P.X - 2P.Z
-	x2P.Mul(x2P, z2P)    // 2P.X = 2P.X * 2P.Z
-	xPaQ.Mul(a24, &t2)   // Xp+q = A24 * t2
-	zPaQ.Sub(&t0, &t1)   // Zp+q = t0 - t1
-	z2P.Add(xPaQ, z2P)   // 2P.Z = Xp+q + 2P.Z
-	xPaQ.Add(&t0, &t1)   // Xp+q = t0 + t1
-	z2P.Mul(z2P, &t2)    // 2P.Z = 2P.Z * t2
-	zPaQ.Square(zPaQ)    // Zp+q = Zp+q ^ 2
-	xPaQ.Square(xPaQ)    // Xp+q = Xp+q ^ 2
-	zPaQ.Mul(xQmP, zPaQ) // Zp+q = Xq-p * Zp+q
-	xPaQ.Mul(zQmP, xPaQ) // Xp+q = Zq-p * Xp+q
-	return
-}
-
-// Given the curve parameters, xP = x(P), and k >= 0, compute x2P = x([2^k]P).
-//
-// Returns x2P to allow chaining.  Safe to overlap xP, x2P.
-func (x2P *ProjectivePoint) Pow2k(params *CurveCoefficientsEquiv, xP *ProjectivePoint, k uint32) *ProjectivePoint {
-	var t0, t1 ExtensionFieldElement
-
-	*x2P = *xP
-	x, z := &x2P.X, &x2P.Z
-
-	for i := uint32(0); i < k; i++ {
-		t0.Sub(x, z)           // t0  = Xp - Zp
-		t1.Add(x, z)           // t1  = Xp + Zp
-		t0.Square(&t0)         // t0  = t0 ^ 2
-		t1.Square(&t1)         // t1  = t1 ^ 2
-		z.Mul(&params.C, &t0)  // Z2p = C24 * t0
-		x.Mul(z, &t1)          // X2p = Z2p * t1
-		t1.Sub(&t1, &t0)       // t1  = t1 - t0
-		t0.Mul(&params.A, &t1) // t0  = A24+ * t1
-		z.Add(z, &t0)          // Z2p = Z2p + t0
-		z.Mul(z, &t1)          // Zp  = Z2p * t1
-	}
-
-	return x2P
-}
-
-// Given the curve parameters, xP = x(P), and k >= 0, compute x3P = x([3^k]P).
-//
-// Returns x3P to allow chaining.  Safe to overlap xP, xR.
-func (x3P *ProjectivePoint) Pow3k(params *CurveCoefficientsEquiv, xP *ProjectivePoint, k uint32) *ProjectivePoint {
-	var t0, t1, t2, t3, t4, t5, t6 ExtensionFieldElement
-
-	*x3P = *xP
-	x, z := &x3P.X, &x3P.Z
-
-	for i := uint32(0); i < k; i++ {
-		t0.Sub(x, z)           // t0  = Xp - Zp
-		t2.Square(&t0)         // t2  = t0^2
-		t1.Add(x, z)           // t1  = Xp + Zp
-		t3.Square(&t1)         // t3  = t1^2
-		t4.Add(&t1, &t0)       // t4  = t1 + t0
-		t0.Sub(&t1, &t0)       // t0  = t1 - t0
-		t1.Square(&t4)         // t1  = t4^2
-		t1.Sub(&t1, &t3)       // t1  = t1 - t3
-		t1.Sub(&t1, &t2)       // t1  = t1 - t2
-		t5.Mul(&t3, &params.A) // t5  = t3 * A24+
-		t3.Mul(&t3, &t5)       // t3  = t5 * t3
-		t6.Mul(&t2, &params.C) // t6  = t2 * A24-
-		t2.Mul(&t2, &t6)       // t2  = t2 * t6
-		t3.Sub(&t2, &t3)       // t3  = t2 - t3
-		t2.Sub(&t5, &t6)       // t2  = t5 - t6
-		t1.Mul(&t2, &t1)       // t1  = t2 * t1
-		t2.Add(&t3, &t1)       // t2  = t3 + t1
-		t2.Square(&t2)         // t2  = t2^2
-		x.Mul(&t2, &t4)        // X3p = t2 * t4
-		t1.Sub(&t3, &t1)       // t1  = t3 - t1
-		t1.Square(&t1)         // t1  = t1^2
-		z.Mul(&t1, &t0)        // Z3p = t1 * t0
-	}
-	return x3P
-}
-
-// RightToLeftLadder is a right-to-left point multiplication that given the
-// x-coordinate of P, Q and P-Q calculates the x-coordinate of R=Q+[scalar]P.
-// nbits must be smaller or equal to len(scalar).
-func RightToLeftLadder(c *ProjectiveCurveParameters, P, Q, PmQ *ProjectivePoint,
-	nbits uint, scalar []uint8) ProjectivePoint {
-	var R0, R2, R1 ProjectivePoint
-
-	aPlus2Over4 := c.calcAplus2Over4()
-	R1 = *P
-	R2 = *PmQ
-	R0 = *Q
-
-	// Iterate over the bits of the scalar, bottom to top
-	prevBit := uint8(0)
-	for i := uint(0); i < nbits; i++ {
-		bit := (scalar[i>>3] >> (i & 7) & 1)
-		swap := prevBit ^ bit
-		prevBit = bit
-		ProjectivePointConditionalSwap(&R1, &R2, swap)
-		R0, R2 = xDblAdd(&R0, &R2, &R1, &aPlus2Over4)
-	}
-
-	ProjectivePointConditionalSwap(&R1, &R2, prevBit)
-	return R1
-}

p751toolbox/curve_test.go

@@ -1,372 +0,0 @@
-package p751toolbox
-
-import (
-	"bytes"
-	"math/rand"
-	"reflect"
-	"testing"
-	"testing/quick"
-)
-
-// Sage script for generating test vectors:
-// sage: p = 2^372 * 3^239 - 1; Fp = GF(p)
-// sage: R.<x> = Fp[]
-// sage: Fp2 = Fp.extension(x^2 + 1, 'i')
-// sage: i = Fp2.gen()
-// sage: A = 4385300808024233870220415655826946795549183378139271271040522089756750951667981765872679172832050962894122367066234419550072004266298327417513857609747116903999863022476533671840646615759860564818837299058134292387429068536219*i + 1408083354499944307008104531475821995920666351413327060806684084512082259107262519686546161682384352696826343970108773343853651664489352092568012759783386151707999371397181344707721407830640876552312524779901115054295865393760
-// sage: C = 933177602672972392833143808100058748100491911694554386487433154761658932801917030685312352302083870852688835968069519091048283111836766101703759957146191882367397129269726925521881467635358356591977198680477382414690421049768*i + 9088894745865170214288643088620446862479558967886622582768682946704447519087179261631044546285104919696820250567182021319063155067584445633834024992188567423889559216759336548208016316396859149888322907914724065641454773776307
-// sage: E = EllipticCurve(Fp2, [0,A/C,0,1,0])
-// sage: XP, YP, ZP = (8172151271761071554796221948801462094972242987811852753144865524899433583596839357223411088919388342364651632180452081960511516040935428737829624206426287774255114241789158000915683252363913079335550843837650671094705509470594*i + 9326574858039944121604015439381720195556183422719505497448541073272720545047742235526963773359004021838961919129020087515274115525812121436661025030481584576474033630899768377131534320053412545346268645085054880212827284581557, 2381174772709336084066332457520782192315178511983342038392622832616744048226360647551642232950959910067260611740876401494529727990031260499974773548012283808741733925525689114517493995359390158666069816204787133942283380884077*i + 5378956232034228335189697969144556552783858755832284194802470922976054645696324118966333158267442767138528227968841257817537239745277092206433048875637709652271370008564179304718555812947398374153513738054572355903547642836171, 1)
-// sage: XQ, YQ, ZQ = (58415083458086949460774631288254059520198050925769753726965257584725433404854146588020043440874302516770401547098747946831855304325931998303167820332782040123488940125615738529117852871565495005635582483848203360379709783913*i + 5253691516070829381103946367549411712646323371150903710970418364086793242952116417060222617507143384179331507713214357463046698181663711723374230655213762988214127964685312538735038552802853885461137159568726585741937776693865 : 5158855522131677856751020091923296498106687563952623634063620100943451774747711264382913073742659403103540173656392945587315798626256921403472228775673617431559609635074280718249654789362069112718012186738875122436282913402060*i + 6207982879267706771450280080677554401823496760317291100742934673980817977221212113336809046153999230028113719924529054308271353271269929392876016936836430864080079684360816968551808465413552611322477628209203066093559492022329 : 1)
-// sage: P = E((XP,YP,ZP))
-// sage: X2, Y2, Z2 = 2*P
-// sage: X3, Y3, Z3 = 3*P
-// sage: m = 96550223052359874398280314003345143371473380422728857598463622014420884224892
-
-// A = 4385300808024233870220415655826946795549183378139271271040522089756750951667981765872679172832050962894122367066234419550072004266298327417513857609747116903999863022476533671840646615759860564818837299058134292387429068536219*i + 1408083354499944307008104531475821995920666351413327060806684084512082259107262519686546161682384352696826343970108773343853651664489352092568012759783386151707999371397181344707721407830640876552312524779901115054295865393760
-var curve_A = ExtensionFieldElement{
-	A: Fp751Element{0x8319eb18ca2c435e, 0x3a93beae72cd0267, 0x5e465e1f72fd5a84, 0x8617fa4150aa7272, 0x887da24799d62a13, 0xb079b31b3c7667fe, 0xc4661b150fa14f2e, 0xd4d2b2967bc6efd6, 0x854215a8b7239003, 0x61c5302ccba656c2, 0xf93194a27d6f97a2, 0x1ed9532bca75},
-	B: Fp751Element{0xb6f541040e8c7db6, 0x99403e7365342e15, 0x457e9cee7c29cced, 0x8ece72dc073b1d67, 0x6e73cef17ad28d28, 0x7aed836ca317472, 0x89e1de9454263b54, 0x745329277aa0071b, 0xf623dfc73bc86b9b, 0xb8e3c1d8a9245882, 0x6ad0b3d317770bec, 0x5b406e8d502b}}
-
-// C = 933177602672972392833143808100058748100491911694554386487433154761658932801917030685312352302083870852688835968069519091048283111836766101703759957146191882367397129269726925521881467635358356591977198680477382414690421049768*i + 9088894745865170214288643088620446862479558967886622582768682946704447519087179261631044546285104919696820250567182021319063155067584445633834024992188567423889559216759336548208016316396859149888322907914724065641454773776307
-var curve_C = ExtensionFieldElement{
-	A: Fp751Element{0x4fb2358bbf723107, 0x3a791521ac79e240, 0x283e24ef7c4c922f, 0xc89baa1205e33cc, 0x3031be81cff6fee1, 0xaf7a494a2f6a95c4, 0x248d251eaac83a1d, 0xc122fca1e2550c88, 0xbc0451b11b6cfd3d, 0x9c0a114ab046222c, 0x43b957b32f21f6ea, 0x5b9c87fa61de},
-	B: Fp751Element{0xacf142afaac15ec6, 0xfd1322a504a071d5, 0x56bb205e10f6c5c6, 0xe204d2849a97b9bd, 0x40b0122202fe7f2e, 0xecf72c6fafacf2cb, 0x45dfc681f869f60a, 0x11814c9aff4af66c, 0x9278b0c4eea54fe7, 0x9a633d5baf7f2e2e, 0x69a329e6f1a05112, 0x1d874ace23e4}}
-
-var curve = ProjectiveCurveParameters{A: curve_A, C: curve_C}
-
-// x(P) = 8172151271761071554796221948801462094972242987811852753144865524899433583596839357223411088919388342364651632180452081960511516040935428737829624206426287774255114241789158000915683252363913079335550843837650671094705509470594*i + 9326574858039944121604015439381720195556183422719505497448541073272720545047742235526963773359004021838961919129020087515274115525812121436661025030481584576474033630899768377131534320053412545346268645085054880212827284581557
-var affine_xP = ExtensionFieldElement{
-	A: Fp751Element{0xe8d05f30aac47247, 0x576ec00c55441de7, 0xbf1a8ec5fe558518, 0xd77cb17f77515881, 0x8e9852837ee73ec4, 0x8159634ad4f44a6b, 0x2e4eb5533a798c5, 0x9be8c4354d5bc849, 0xf47dc61806496b84, 0x25d0e130295120e0, 0xdbef54095f8139e3, 0x5a724f20862c},
-	B: Fp751Element{0x3ca30d7623602e30, 0xfb281eddf45f07b7, 0xd2bf62d5901a45bc, 0xc67c9baf86306dd2, 0x4e2bd93093f538ca, 0xcfd92075c25b9cbe, 0xceafe9a3095bcbab, 0x7d928ad380c85414, 0x37c5f38b2afdc095, 0x75325899a7b779f4, 0xf130568249f20fdd, 0x178f264767d1}}
-
-// x([2]P) = 1476586462090705633631615225226507185986710728845281579274759750260315746890216330325246185232948298241128541272709769576682305216876843626191069809810990267291824247158062860010264352034514805065784938198193493333201179504845*i + 3623708673253635214546781153561465284135688791018117615357700171724097420944592557655719832228709144190233454198555848137097153934561706150196041331832421059972652530564323645509890008896574678228045006354394485640545367112224
-var affine_xP2 = ExtensionFieldElement{
-	A: Fp751Element{0x2a77afa8576ce979, 0xab1360e69b0aeba0, 0xd79e3e3cbffad660, 0x5fd0175aa10f106b, 0x1800ebafce9fbdbc, 0x228fc9142bdd6166, 0x867cf907314e34c3, 0xa58d18c94c13c31c, 0x699a5bc78b11499f, 0xa29fc29a01f7ccf1, 0x6c69c0c5347eebce, 0x38ecee0cc57},
-	B: Fp751Element{0x43607fd5f4837da0, 0x560bad4ce27f8f4a, 0x2164927f8495b4dd, 0x621103fdb831a997, 0xad740c4eea7db2db, 0x2cde0442205096cd, 0x2af51a70ede8324e, 0x41a4e680b9f3466, 0x5481f74660b8f476, 0xfcb2f3e656ff4d18, 0x42e3ce0837171acc, 0x44238c30530c}}
-
-// x([3]P) = 9351941061182433396254169746041546943662317734130813745868897924918150043217746763025923323891372857734564353401396667570940585840576256269386471444236630417779544535291208627646172485976486155620044292287052393847140181703665*i + 9010417309438761934687053906541862978676948345305618417255296028956221117900864204687119686555681136336037659036201780543527957809743092793196559099050594959988453765829339642265399496041485088089691808244290286521100323250273
-var affine_xP3 = ExtensionFieldElement{
-	A: Fp751Element{0x2096e3f23feca947, 0xf36f635aa4ad8634, 0xdae3b1c6983c5e9a, 0xe08df6c262cb74b4, 0xd2ca4edc37452d3d, 0xfb5f3fe42f500c79, 0x73740aa3abc2b21f, 0xd535fd869f914cca, 0x4a558466823fb67f, 0x3e50a7a0e3bfc715, 0xf43c6da9183a132f, 0x61aca1e1b8b9},
-	B: Fp751Element{0x1e54ec26ea5077bd, 0x61380572d8769f9a, 0xc615170684f59818, 0x6309c3b93e84ef6e, 0x33c74b1318c3fcd0, 0xfe8d7956835afb14, 0x2d5a7b55423c1ecc, 0x869db67edfafea68, 0x1292632394f0a628, 0x10bba48225bfd141, 0x6466c28b408daba, 0x63cacfdb7c43}}
-
-// x([2^2]P) = 441719501189485559222919502512761433931671682884872259563221427434901842337947564993718830905758163254463901652874331063768876314142359813382575876106725244985607032091781306919778265250690045578695338669105227100119314831452*i + 6961734028200975729170216310486458180126343885294922940439352055937945948015840788921225114530454649744697857047401608073256634790353321931728699534700109268264491160589480994022419317695690866764726967221310990488404411684053
-var affine_xP4 = ExtensionFieldElement{
-	A: Fp751Element{0x6f9dbe4c39175153, 0xf2fec757eb99e88, 0x43d7361a93733d91, 0x3abd10ed19c85a3d, 0xc4de9ab9c5ef7181, 0x53e375901684c900, 0x68ffc3e7d71c41ff, 0x47adab62c8d942fe, 0x226a33fd6fbb381d, 0x87ef4c8fdd83309a, 0xaca1cf44c5fa8799, 0x6cbae86c755f},
-	B: Fp751Element{0x4c80c37fe68282a7, 0xbd8b9d7248bf553a, 0x1fb0e8e74d5e1762, 0xb63fa0e4e5f91482, 0xc675ab8a45a1439, 0xdfa6772deace7820, 0xf0d813d71d9a9255, 0x53a1a58c634534bd, 0x4ebfc6485fdfd888, 0x6991fe4358bcf169, 0xc0547bdaca85b6fd, 0xf461548d632}}
-
-// x([3^2]P) = 3957171963425208493644602380039721164492341594850197356580248639045894821895524981729970650520936632013218950972842867220898274664982599375786979902471523505057611521217523103474682939638645404445093536997296151472632038973463*i + 1357869545269286021642168835877253886774707209614159162748874474269328421720121175566245719916322684751967981171882659798149072149161259103020057556362998810229937432814792024248155991141511691087135859252304684633946087474060
-var affine_xP9 = ExtensionFieldElement{
-	A: Fp751Element{0x7c0daa0f04ded4e0, 0x52dc4f883d85e065, 0x91afbdc2c1714d0b, 0xb7b3db8e658cfeba, 0x43d4e72a692882f3, 0x535c56d83753da30, 0xc8a58724433cbf5d, 0x351153c0a5e74219, 0x2c81827d19f93dd5, 0x26ef8aca3370ea1a, 0x1cf939a6dd225dec, 0x3403cb28ad41},
-	B: Fp751Element{0x93e7bc373a9ff7b, 0x57b8cc47635ebc0f, 0x92eab55689106cf3, 0x93643111d421f24c, 0x1c58b519506f6b7a, 0xebd409fb998faa13, 0x5c86ed799d09d80e, 0xd9a1d764d6363562, 0xf95e87f92fb0c4cc, 0x6b2bbaf5632a5609, 0x2d9b6a809dfaff7f, 0x29c0460348b}}
-
-// m = 96550223052359874398280314003345143371473380422728857598463622014420884224892
-var mScalarBytes = [...]uint8{0x7c, 0x7b, 0x95, 0xfa, 0xb4, 0x75, 0x6c, 0x48, 0x8c, 0x17, 0x55, 0xb4, 0x49, 0xf5, 0x1e, 0xa3, 0xb, 0x31, 0xf0, 0xa4, 0xa6, 0x81, 0xad, 0x94, 0x51, 0x11, 0xe7, 0xf5, 0x5b, 0x7d, 0x75, 0xd5}
-
-// x([m]P) = 7893578558852400052689739833699289348717964559651707250677393044951777272628231794999463214496545377542328262828965953246725804301238040891993859185944339366910592967840967752138115122568615081881937109746463885908097382992642*i + 8293895847098220389503562888233557012043261770526854885191188476280014204211818299871679993460086974249554528517413590157845430186202704783785316202196966198176323445986064452630594623103149383929503089342736311904030571524837
-var affine_xaP = ExtensionFieldElement{
-	A: Fp751Element{0x2112f3c7d7f938bb, 0x704a677f0a4df08f, 0x825370e31fb4ef00, 0xddbf79b7469f902, 0x27640c899ea739fd, 0xfb7b8b19f244108e, 0x546a6679dd3baebc, 0xe9f0ecf398d5265f, 0x223d2b350e75e461, 0x84b322a0b6aff016, 0xfabe426f539f8b39, 0x4507a0604f50},
-	B: Fp751Element{0xac77737e5618a5fe, 0xf91c0e08c436ca52, 0xd124037bc323533c, 0xc9a772bf52c58b63, 0x3b30c8f38ef6af4d, 0xb9eed160e134f36e, 0x24e3836393b25017, 0xc828be1b11baf1d9, 0x7b7dab585df50e93, 0x1ca3852c618bd8e0, 0x4efa73bcb359fa00, 0x50b6a923c2d4}}
-
-// Inputs for testing 3-point-ladder
-var threePointLadderInputs = []ProjectivePoint{
-	// x(P)
-	ProjectivePoint{
-		X: ExtensionFieldElement{
-			A: Fp751Element{0xe8d05f30aac47247, 0x576ec00c55441de7, 0xbf1a8ec5fe558518, 0xd77cb17f77515881, 0x8e9852837ee73ec4, 0x8159634ad4f44a6b, 0x2e4eb5533a798c5, 0x9be8c4354d5bc849, 0xf47dc61806496b84, 0x25d0e130295120e0, 0xdbef54095f8139e3, 0x5a724f20862c},
-			B: Fp751Element{0x3ca30d7623602e30, 0xfb281eddf45f07b7, 0xd2bf62d5901a45bc, 0xc67c9baf86306dd2, 0x4e2bd93093f538ca, 0xcfd92075c25b9cbe, 0xceafe9a3095bcbab, 0x7d928ad380c85414, 0x37c5f38b2afdc095, 0x75325899a7b779f4, 0xf130568249f20fdd, 0x178f264767d1}},
-		Z: oneExtensionField,
-	},
-	// x(Q)
-	ProjectivePoint{
-		X: ExtensionFieldElement{
-			A: Fp751Element{0x2b71a2a93ad1e10e, 0xf0b9842a92cfb333, 0xae17373615a27f5c, 0x3039239f428330c4, 0xa0c4b735ed7dcf98, 0x6e359771ddf6af6a, 0xe986e4cac4584651, 0x8233a2b622d5518, 0xbfd67bf5f06b818b, 0xdffe38d0f5b966a6, 0xa86b36a3272ee00a, 0x193e2ea4f68f},
-			B: Fp751Element{0x5a0f396459d9d998, 0x479f42250b1b7dda, 0x4016b57e2a15bf75, 0xc59f915203fa3749, 0xd5f90257399cf8da, 0x1fb2dadfd86dcef4, 0x600f20e6429021dc, 0x17e347d380c57581, 0xc1b0d5fa8fe3e440, 0xbcf035330ac20e8, 0x50c2eb5f6a4f03e6, 0x86b7c4571}},
-		Z: oneExtensionField,
-	},
-	// x(P-Q)
-	ProjectivePoint{
-		X: ExtensionFieldElement{
-			A: Fp751Element{0x4aafa9f378f7b5ff, 0x1172a683aa8eee0, 0xea518d8cbec2c1de, 0xe191bcbb63674557, 0x97bc19637b259011, 0xdbeae5c9f4a2e454, 0x78f64d1b72a42f95, 0xe71cb4ea7e181e54, 0xe4169d4c48543994, 0x6198c2286a98730f, 0xd21d675bbab1afa5, 0x2e7269fce391},
-			B: Fp751Element{0x23355783ce1d0450, 0x683164cf4ce3d93f, 0xae6d1c4d25970fd8, 0x7807007fb80b48cf, 0xa005a62ec2bbb8a2, 0x6b5649bd016004cb, 0xbb1a13fa1330176b, 0xbf38e51087660461, 0xe577fddc5dd7b930, 0x5f38116f56947cd3, 0x3124f30b98c36fde, 0x4ca9b6e6db37}},
-		Z: oneExtensionField,
-	},
-}
-
-func (P ProjectivePoint) Generate(rand *rand.Rand, size int) reflect.Value {
-	f := ExtensionFieldElement{}
-	x, _ := f.Generate(rand, size).Interface().(ExtensionFieldElement)
-	z, _ := f.Generate(rand, size).Interface().(ExtensionFieldElement)
-	return reflect.ValueOf(ProjectivePoint{
-		X: x,
-		Z: z,
-	})
-}
-
-func (curve ProjectiveCurveParameters) Generate(rand *rand.Rand, size int) reflect.Value {
-	f := ExtensionFieldElement{}
-	A, _ := f.Generate(rand, size).Interface().(ExtensionFieldElement)
-	C, _ := f.Generate(rand, size).Interface().(ExtensionFieldElement)
-	return reflect.ValueOf(ProjectiveCurveParameters{
-		A: A,
-		C: C,
-	})
-}
-
-// Helpers
-
-// Given xP = x(P), xQ = x(Q), and xPmQ = x(P-Q), compute xR = x(P+Q).
-//
-// Returns xR to allow chaining.  Safe to overlap xP, xQ, xR.
-func (xR *ProjectivePoint) Add(xP, xQ, xPmQ *ProjectivePoint) *ProjectivePoint {
-	// Algorithm 1 of Costello-Smith.
-	var v0, v1, v2, v3, v4 ExtensionFieldElement
-	v0.Add(&xP.X, &xP.Z)               // X_P + Z_P
-	v1.Sub(&xQ.X, &xQ.Z).Mul(&v1, &v0) // (X_Q - Z_Q)(X_P + Z_P)
-	v0.Sub(&xP.X, &xP.Z)               // X_P - Z_P
-	v2.Add(&xQ.X, &xQ.Z).Mul(&v2, &v0) // (X_Q + Z_Q)(X_P - Z_P)
-	v3.Add(&v1, &v2).Square(&v3)       // 4(X_Q X_P - Z_Q Z_P)^2
-	v4.Sub(&v1, &v2).Square(&v4)       // 4(X_Q Z_P - Z_Q X_P)^2
-	v0.Mul(&xPmQ.Z, &v3)               // 4X_{P-Q}(X_Q X_P - Z_Q Z_P)^2
-	xR.Z.Mul(&xPmQ.X, &v4)             // 4Z_{P-Q}(X_Q Z_P - Z_Q X_P)^2
-	xR.X = v0
-	return xR
-}
-
-// Given xP = x(P) and cached curve parameters Aplus2C = A + 2*C, C4 = 4*C,
-// compute xQ = x([2]P).
-//
-// Returns xQ to allow chaining.  Safe to overlap xP, xQ.
-func (xQ *ProjectivePoint) Double(xP *ProjectivePoint, Aplus2C, C4 *ExtensionFieldElement) *ProjectivePoint {
-	// Algorithm 2 of Costello-Smith, amended to work with projective curve coefficients.
-	var v1, v2, v3, xz4 ExtensionFieldElement
-	v1.Add(&xP.X, &xP.Z).Square(&v1) // (X+Z)^2
-	v2.Sub(&xP.X, &xP.Z).Square(&v2) // (X-Z)^2
-	xz4.Sub(&v1, &v2)                // 4XZ = (X+Z)^2 - (X-Z)^2
-	v2.Mul(&v2, C4)                  // 4C(X-Z)^2
-	xQ.X.Mul(&v1, &v2)               // 4C(X+Z)^2(X-Z)^2
-	v3.Mul(&xz4, Aplus2C)            // 4XZ(A + 2C)
-	v3.Add(&v3, &v2)                 // 4XZ(A + 2C) + 4C(X-Z)^2
-	xQ.Z.Mul(&v3, &xz4)              // (4XZ(A + 2C) + 4C(X-Z)^2)4XZ
-	// Now (xQ.x : xQ.z)
-	//   = (4C(X+Z)^2(X-Z)^2 : (4XZ(A + 2C) + 4C(X-Z)^2)4XZ )
-	//   = ((X+Z)^2(X-Z)^2 : (4XZ((A + 2C)/4C) + (X-Z)^2)4XZ )
-	//   = ((X+Z)^2(X-Z)^2 : (4XZ((a + 2)/4) + (X-Z)^2)4XZ )
-	return xQ
-}
-
-// Given x(P) and a scalar m in little-endian bytes, compute x([m]P) using the
-// Montgomery ladder.  This is described in Algorithm 8 of Costello-Smith.
-//
-// This function's execution time is dependent only on the byte-length of the
-// input scalar.  All scalars of the same input length execute in uniform time.
-// The scalar can be padded with zero bytes to ensure a uniform length.
-//
-// Safe to overlap the source with the destination.
-func (xQ *ProjectivePoint) ScalarMult(curve *ProjectiveCurveParameters, xP *ProjectivePoint, scalar []uint8) *ProjectivePoint {
-	var x0, x1, tmp ProjectivePoint
-	var Aplus2C, C4 ExtensionFieldElement
-
-	Aplus2C.Add(&curve.C, &curve.C) // = 2*C
-	C4.Add(&Aplus2C, &Aplus2C)      // = 4*C
-	Aplus2C.Add(&Aplus2C, &curve.A) // = 2*C + A
-
-	x0.X.One()
-	x0.Z.Zero()
-	x1 = *xP
-
-	// Iterate over the bits of the scalar, top to bottom
-	prevBit := uint8(0)
-	for i := len(scalar) - 1; i >= 0; i-- {
-		scalarByte := scalar[i]
-		for j := 7; j >= 0; j-- {
-			bit := (scalarByte >> uint(j)) & 0x1
-			ProjectivePointConditionalSwap(&x0, &x1, (bit ^ prevBit))
-			tmp.Double(&x0, &Aplus2C, &C4)
-			x1.Add(&x0, &x1, xP)
-			x0 = tmp
-			prevBit = bit
-		}
-	}
-	// now prevBit is the lowest bit of the scalar
-	ProjectivePointConditionalSwap(&x0, &x1, prevBit)
-	*xQ = x0
-	return xQ
-}
-
-// Tests
-
-func TestOne(t *testing.T) {
-	var tmp ExtensionFieldElement
-
-	tmp.Mul(&oneExtensionField, &affine_xP)
-	if !tmp.VartimeEq(&affine_xP) {
-		t.Error("Not equal 1")
-	}
-}
-
-func TestScalarMultVersusSage(t *testing.T) {
-	var xP ProjectivePoint
-
-	xP.FromAffine(&affine_xP)
-	affine_xQ := xP.ScalarMult(&curve, &xP, mScalarBytes[:]).ToAffine() // = x([m]P)
-	if !affine_xaP.VartimeEq(affine_xQ) {
-		t.Error("\nExpected\n", affine_xaP, "\nfound\n", affine_xQ)
-	}
-}
-
-func Test_jInvariant(t *testing.T) {
-	var curve = ProjectiveCurveParameters{A: curve_A, C: curve_C}
-	var jbufRes = make([]byte, P751_SharedSecretSize)
-	var jbufExp = make([]byte, P751_SharedSecretSize)
-	// Computed using Sage
-	// j = 3674553797500778604587777859668542828244523188705960771798425843588160903687122861541242595678107095655647237100722594066610650373491179241544334443939077738732728884873568393760629500307797547379838602108296735640313894560419*i + 3127495302417548295242630557836520229396092255080675419212556702820583041296798857582303163183558315662015469648040494128968509467224910895884358424271180055990446576645240058960358037224785786494172548090318531038910933793845
-	var known_j = ExtensionFieldElement{
-		A: Fp751Element{0xc7a8921c1fb23993, 0xa20aea321327620b, 0xf1caa17ed9676fa8, 0x61b780e6b1a04037, 0x47784af4c24acc7a, 0x83926e2e300b9adf, 0xcd891d56fae5b66, 0x49b66985beb733bc, 0xd4bcd2a473d518f, 0xe242239991abe224, 0xa8af5b20f98672f8, 0x139e4d4e4d98},
-		B: Fp751Element{0xb5b52a21f81f359, 0x715e3a865db6d920, 0x9bac2f9d8911978b, 0xef14acd8ac4c1e3d, 0xe81aacd90cfb09c8, 0xaf898288de4a09d9, 0xb85a7fb88c5c4601, 0x2c37c3f1dd303387, 0x7ad3277fe332367c, 0xd4cbee7f25a8e6f8, 0x36eacbe979eaeffa, 0x59eb5a13ac33},
-	}
-
-	curve.Jinvariant(jbufRes)
-	known_j.ToBytes(jbufExp)
-
-	if !bytes.Equal(jbufRes, jbufExp) {
-		t.Error("Computed incorrect j-invariant: found\n", jbufRes, "\nexpected\n", jbufExp)
-	}
-}
-
-func TestProjectivePointVartimeEq(t *testing.T) {
-	var xP ProjectivePoint
-
-	xP.FromAffine(&affine_xP)
-	xQ := xP
-	// Scale xQ, which results in the same projective point
-	xQ.X.Mul(&xQ.X, &curve_A)
-	xQ.Z.Mul(&xQ.Z, &curve_A)
-	if !xQ.VartimeEq(&xP) {
-		t.Error("Expected the scaled point to be equal to the original")
-	}
-}
-
-func TestPointDoubleVersusSage(t *testing.T) {
-	var curve = ProjectiveCurveParameters{A: curve_A, C: curve_C}
-	var params = curve.CalcCurveParamsEquiv4()
-	var xP, xQ ProjectivePoint
-
-	xP.FromAffine(&affine_xP)
-	affine_xQ := xQ.Pow2k(&params, &xP, 1).ToAffine()
-	if !affine_xQ.VartimeEq(&affine_xP2) {
-		t.Error("\nExpected\n", affine_xP2, "\nfound\n", affine_xQ)
-	}
-}
-
-func TestPointMul4VersusSage(t *testing.T) {
-	var params = curve.CalcCurveParamsEquiv4()
-	var xP, xQ ProjectivePoint
-
-	xP.FromAffine(&affine_xP)
-	affine_xQ := xQ.Pow2k(&params, &xP, 2).ToAffine()
-	if !affine_xQ.VartimeEq(&affine_xP4) {
-		t.Error("\nExpected\n", affine_xP4, "\nfound\n", affine_xQ)
-	}
-}
-
-func TestPointMul9VersusSage(t *testing.T) {
-	var params = curve.CalcCurveParamsEquiv3()
-	var xP, xQ ProjectivePoint
-
-	xP.FromAffine(&affine_xP)
-	affine_xQ := xQ.Pow3k(&params, &xP, 2).ToAffine()
-	if !affine_xQ.VartimeEq(&affine_xP9) {
-		t.Error("\nExpected\n", affine_xP9, "\nfound\n", affine_xQ)
-	}
-}
-
-func TestPointPow2kVersusScalarMult(t *testing.T) {
-	var xP, xQ, xR ProjectivePoint
-	var params = curve.CalcCurveParamsEquiv4()
-
-	xP.FromAffine(&affine_xP)
-	affine_xQ := xQ.Pow2k(&params, &xP, 5).ToAffine()              // = x([32]P)
-	affine_xR := xR.ScalarMult(&curve, &xP, []byte{32}).ToAffine() // = x([32]P)
-
-	if !affine_xQ.VartimeEq(affine_xR) {
-		t.Error("\nExpected\n", affine_xQ, "\nfound\n", affine_xR)
-	}
-}
-
-func TestRecoverCoordinateA(t *testing.T) {
-	var cparam ProjectiveCurveParameters
-	// Vectors generated with SIKE reference implementation
-	var a = ExtensionFieldElement{
-		A: Fp751Element{0x9331D9C5AAF59EA4, 0xB32B702BE4046931, 0xCEBB333912ED4D34, 0x5628CE37CD29C7A2, 0x0BEAC5ED48B7F58E, 0x1FB9D3E281D65B07, 0x9C0CFACC1E195662, 0xAE4BCE0F6B70F7D9, 0x59E4E63D43FE71A0, 0xEF7CE57560CC8615, 0xE44A8FB7901E74E8, 0x000069D13C8366D1},
-		B: Fp751Element{0xF6DA1070279AB966, 0xA78FB0CE7268C762, 0x19B40F044A57ABFA, 0x7AC8EE6160C0C233, 0x93D4993442947072, 0x757D2B3FA4E44860, 0x073A920F8C4D5257, 0x2031F1B054734037, 0xDEFAA1D2406555CD, 0x26F9C70E1496BE3D, 0x5B3F335A0A4D0976, 0x000013628B2E9C59}}
-	var affine_xP = ExtensionFieldElement{
-		A: Fp751Element{0xea6b2d1e2aebb250, 0x35d0b205dc4f6386, 0xb198e93cb1830b8d, 0x3b5b456b496ddcc6, 0x5be3f0d41132c260, 0xce5f188807516a00, 0x54f3e7469ea8866d, 0x33809ef47f36286, 0x6fa45f83eabe1edb, 0x1b3391ae5d19fd86, 0x1e66daf48584af3f, 0xb430c14aaa87},
-		B: Fp751Element{0x97b41ebc61dcb2ad, 0x80ead31cb932f641, 0x40a940099948b642, 0x2a22fd16cdc7fe84, 0xaabf35b17579667f, 0x76c1d0139feb4032, 0x71467e1e7b1949be, 0x678ca8dadd0d6d81, 0x14445daea9064c66, 0x92d161eab4fa4691, 0x8dfbb01b6b238d36, 0x2e3718434e4e}}
-	var affine_xQ = ExtensionFieldElement{
-		A: Fp751Element{0xb055cf0ca1943439, 0xa9ff5de2fa6c69ed, 0x4f2761f934e5730a, 0x61a1dcaa1f94aa4b, 0xce3c8fadfd058543, 0xeac432aaa6701b8e, 0x8491d523093aea8b, 0xba273f9bd92b9b7f, 0xd8f59fd34439bb5a, 0xdc0350261c1fe600, 0x99375ab1eb151311, 0x14d175bbdbc5},
-		B: Fp751Element{0xffb0ef8c2111a107, 0x55ceca3825991829, 0xdbf8a1ccc075d34b, 0xb8e9187bd85d8494, 0x670aa2d5c34a03b0, 0xef9fe2ed2b064953, 0xc911f5311d645aee, 0xf4411f409e410507, 0x934a0a852d03e1a8, 0xe6274e67ae1ad544, 0x9f4bc563c69a87bc, 0x6f316019681e}}
-	var affine_xQmP = ExtensionFieldElement{
-		A: Fp751Element{0x6ffb44306a153779, 0xc0ffef21f2f918f3, 0x196c46d35d77f778, 0x4a73f80452edcfe6, 0x9b00836bce61c67f, 0x387879418d84219e, 0x20700cf9fc1ec5d1, 0x1dfe2356ec64155e, 0xf8b9e33038256b1c, 0xd2aaf2e14bada0f0, 0xb33b226e79a4e313, 0x6be576fad4e5},
-		B: Fp751Element{0x7db5dbc88e00de34, 0x75cc8cb9f8b6e11e, 0x8c8001c04ebc52ac, 0x67ef6c981a0b5a94, 0xc3654fbe73230738, 0xc6a46ee82983ceca, 0xed1aa61a27ef49f0, 0x17fe5a13b0858fe0, 0x9ae0ca945a4c6b3c, 0x234104a218ad8878, 0xa619627166104394, 0x556a01ff2e7e}}
-
-	cparam.RecoverCoordinateA(&affine_xP, &affine_xQ, &affine_xQmP)
-	cparam.C.One()
-
-	// Check A is correct
-	if !cparam.A.VartimeEq(&a) {
-		t.Error("\nExpected\n", a, "\nfound\n", cparam.A)
-	}
-
-	// Check C is not changed
-	if !cparam.C.VartimeEq(&oneExtensionField) {
-		t.Error("\nExpected\n", cparam.C, "\nfound\n", oneExtensionField)
-	}
-}
-
-func TestR2LVersusSage(t *testing.T) {
-	var xR ProjectivePoint
-
-	sageAffine_xR := ExtensionFieldElement{
-		A: Fp751Element{0x729465ba800d4fd5, 0x9398015b59e514a1, 0x1a59dd6be76c748e, 0x1a7db94eb28dd55c, 0x444686e680b1b8ec, 0xcc3d4ace2a2454ff, 0x51d3dab4ec95a419, 0xc3b0f33594acac6a, 0x9598a74e7fd44f8a, 0x4fbf8c638f1c2e37, 0x844e347033052f51, 0x6cd6de3eafcf},
-		B: Fp751Element{0x85da145412d73430, 0xd83c0e3b66eb3232, 0xd08ff2d453ec1369, 0xa64aaacfdb395b13, 0xe9cba211a20e806e, 0xa4f80b175d937cfc, 0x556ce5c64b1f7937, 0xb59b39ea2b3fdf7a, 0xc2526b869a4196b3, 0x8dad90bca9371750, 0xdfb4a30c9d9147a2, 0x346d2130629b}}
-	xR = RightToLeftLadder(&curve, &threePointLadderInputs[0], &threePointLadderInputs[1], &threePointLadderInputs[2], uint(len(mScalarBytes)*8), mScalarBytes[:])
-	affine_xR := xR.ToAffine()
-
-	if !affine_xR.VartimeEq(&sageAffine_xR) {
-		t.Error("\nExpected\n", sageAffine_xR, "\nfound\n", affine_xR)
-	}
-}
-
-func TestPointTripleVersusAddDouble(t *testing.T) {
-	tripleEqualsAddDouble := func(curve ProjectiveCurveParameters, P ProjectivePoint) bool {
-		var P2, P3, P2plusP ProjectivePoint
-
-		eqivParams4 := curve.CalcCurveParamsEquiv4()
-		eqivParams3 := curve.CalcCurveParamsEquiv3()
-		P2.Pow2k(&eqivParams4, &P, 1) // = x([2]P)
-		P3.Pow3k(&eqivParams3, &P, 1) // = x([3]P)
-		P2plusP.Add(&P2, &P, &P)      // = x([2]P + P)
-		return P3.VartimeEq(&P2plusP)
-	}
-
-	if err := quick.Check(tripleEqualsAddDouble, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-func BenchmarkThreePointLadder379BitScalar(b *testing.B) {
-	var mScalarBytes = [...]uint8{84, 222, 146, 63, 85, 18, 173, 162, 167, 38, 10, 8, 143, 176, 93, 228, 247, 128, 50, 128, 205, 42, 15, 137, 119, 67, 43, 3, 61, 91, 237, 24, 235, 12, 53, 96, 186, 164, 232, 223, 197, 224, 64, 109, 137, 63, 246, 4}
-
-	for n := 0; n < b.N; n++ {
-		RightToLeftLadder(&curve, &threePointLadderInputs[0], &threePointLadderInputs[1], &threePointLadderInputs[2], uint(len(mScalarBytes)*8), mScalarBytes[:])
-	}
-}
-
-func BenchmarkR2L379BitScalar(b *testing.B) {
-	var mScalarBytes = [...]uint8{84, 222, 146, 63, 85, 18, 173, 162, 167, 38, 10, 8, 143, 176, 93, 228, 247, 128, 50, 128, 205, 42, 15, 137, 119, 67, 43, 3, 61, 91, 237, 24, 235, 12, 53, 96, 186, 164, 232, 223, 197, 224, 64, 109, 137, 63, 246, 4}
-
-	for n := 0; n < b.N; n++ {
-		RightToLeftLadder(&curve, &threePointLadderInputs[0], &threePointLadderInputs[1], &threePointLadderInputs[2], uint(len(mScalarBytes)*8), mScalarBytes[:])
-	}
-}

p751toolbox/field.go

@@ -1,427 +0,0 @@
-package p751toolbox
-
-//------------------------------------------------------------------------------
-// Extension Field
-//------------------------------------------------------------------------------
-
-// Represents an element of the extension field F_{p^2}.
-type ExtensionFieldElement struct {
-	// This field element is in Montgomery form, so that the value `A` is
-	// represented by `aR mod p`.
-	A Fp751Element
-	// This field element is in Montgomery form, so that the value `B` is
-	// represented by `bR mod p`.
-	B Fp751Element
-}
-
-var zeroExtensionField = ExtensionFieldElement{
-	A: Fp751Element{0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0},
-	B: Fp751Element{0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0},
-}
-
-var oneExtensionField = ExtensionFieldElement{
-	A: Fp751Element{0x249ad, 0x0, 0x0, 0x0, 0x0, 0x8310000000000000, 0x5527b1e4375c6c66, 0x697797bf3f4f24d0, 0xc89db7b2ac5c4e2e, 0x4ca4b439d2076956, 0x10f7926c7512c7e9, 0x2d5b24bce5e2},
-	B: Fp751Element{0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0},
-}
-
-// 2*p751
-var p751x2 = Fp751Element{
-	0xFFFFFFFFFFFFFFFE, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF,
-	0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xDD5FFFFFFFFFFFFF,
-	0xC7D92D0A93F0F151, 0xB52B363427EF98ED, 0x109D30CFADD7D0ED,
-	0x0AC56A08B964AE90, 0x1C25213F2F75B8CD, 0x0000DFCBAA83EE38}
-
-// p751
-var p751 = Fp751Element{
-	0xffffffffffffffff, 0xffffffffffffffff, 0xffffffffffffffff,
-	0xffffffffffffffff, 0xffffffffffffffff, 0xeeafffffffffffff,
-	0xe3ec968549f878a8, 0xda959b1a13f7cc76, 0x084e9867d6ebe876,
-	0x8562b5045cb25748, 0x0e12909f97badc66, 0x00006fe5d541f71c}
-
-// p751 + 1
-var p751p1 = Fp751Element{
-	0x0000000000000000, 0x0000000000000000, 0x0000000000000000,
-	0x0000000000000000, 0x0000000000000000, 0xeeb0000000000000,
-	0xe3ec968549f878a8, 0xda959b1a13f7cc76, 0x084e9867d6ebe876,
-	0x8562b5045cb25748, 0x0e12909f97badc66, 0x00006fe5d541f71c}
-
-// Set dest = 0.
-//
-// Returns dest to allow chaining operations.
-func (dest *ExtensionFieldElement) Zero() *ExtensionFieldElement {
-	*dest = zeroExtensionField
-	return dest
-}
-
-// Set dest = 1.
-//
-// Returns dest to allow chaining operations.
-func (dest *ExtensionFieldElement) One() *ExtensionFieldElement {
-	*dest = oneExtensionField
-	return dest
-}
-
-// Set dest = lhs * rhs.
-//
-// Allowed to overlap lhs or rhs with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *ExtensionFieldElement) Mul(lhs, rhs *ExtensionFieldElement) *ExtensionFieldElement {
-	// Let (a,b,c,d) = (lhs.a,lhs.b,rhs.a,rhs.b).
-	a := &lhs.A
-	b := &lhs.B
-	c := &rhs.A
-	d := &rhs.B
-
-	// We want to compute
-	//
-	// (a + bi)*(c + di) = (a*c - b*d) + (a*d + b*c)i
-	//
-	// Use Karatsuba's trick: note that
-	//
-	// (b - a)*(c - d) = (b*c + a*d) - a*c - b*d
-	//
-	// so (a*d + b*c) = (b-a)*(c-d) + a*c + b*d.
-
-	var ac, bd fp751X2
-	fp751Mul(&ac, a, c) // = a*c*R*R
-	fp751Mul(&bd, b, d) // = b*d*R*R
-
-	var b_minus_a, c_minus_d Fp751Element
-	fp751SubReduced(&b_minus_a, b, a) // = (b-a)*R
-	fp751SubReduced(&c_minus_d, c, d) // = (c-d)*R
-
-	var ad_plus_bc fp751X2
-	fp751Mul(&ad_plus_bc, &b_minus_a, &c_minus_d) // = (b-a)*(c-d)*R*R
-	fp751X2AddLazy(&ad_plus_bc, &ad_plus_bc, &ac) // = ((b-a)*(c-d) + a*c)*R*R
-	fp751X2AddLazy(&ad_plus_bc, &ad_plus_bc, &bd) // = ((b-a)*(c-d) + a*c + b*d)*R*R
-
-	fp751MontgomeryReduce(&dest.B, &ad_plus_bc) // = (a*d + b*c)*R mod p
-
-	var ac_minus_bd fp751X2
-	fp751X2SubLazy(&ac_minus_bd, &ac, &bd)       // = (a*c - b*d)*R*R
-	fp751MontgomeryReduce(&dest.A, &ac_minus_bd) // = (a*c - b*d)*R mod p
-
-	return dest
-}
-
-// Set dest = 1/x
-//
-// Allowed to overlap dest with x.
-//
-// Returns dest to allow chaining operations.
-func (dest *ExtensionFieldElement) Inv(x *ExtensionFieldElement) *ExtensionFieldElement {
-	a := &x.A
-	b := &x.B
-
-	// We want to compute
-	//
-	//    1          1     (a - bi)	    (a - bi)
-	// -------- = -------- -------- = -----------
-	// (a + bi)   (a + bi) (a - bi)   (a^2 + b^2)
-	//
-	// Letting c = 1/(a^2 + b^2), this is
-	//
-	// 1/(a+bi) = a*c - b*ci.
-
-	var asq_plus_bsq PrimeFieldElement
-	var asq, bsq fp751X2
-	fp751Mul(&asq, a, a)                         // = a*a*R*R
-	fp751Mul(&bsq, b, b)                         // = b*b*R*R
-	fp751X2AddLazy(&asq, &asq, &bsq)             // = (a^2 + b^2)*R*R
-	fp751MontgomeryReduce(&asq_plus_bsq.A, &asq) // = (a^2 + b^2)*R mod p
-	// Now asq_plus_bsq = a^2 + b^2
-
-	// Invert asq_plus_bsq
-	inv := asq_plus_bsq
-	inv.Mul(&asq_plus_bsq, &asq_plus_bsq)
-	inv.P34(&inv)
-	inv.Mul(&inv, &inv)
-	inv.Mul(&inv, &asq_plus_bsq)
-
-	var ac fp751X2
-	fp751Mul(&ac, a, &inv.A)
-	fp751MontgomeryReduce(&dest.A, &ac)
-
-	var minus_b Fp751Element
-	fp751SubReduced(&minus_b, &minus_b, b)
-	var minus_bc fp751X2
-	fp751Mul(&minus_bc, &minus_b, &inv.A)
-	fp751MontgomeryReduce(&dest.B, &minus_bc)
-
-	return dest
-}
-
-// Set (y1, y2, y3)  = (1/x1, 1/x2, 1/x3).
-//
-// All xi, yi must be distinct.
-func ExtensionFieldBatch3Inv(x1, x2, x3, y1, y2, y3 *ExtensionFieldElement) {
-	var x1x2, t ExtensionFieldElement
-	x1x2.Mul(x1, x2)           // x1*x2
-	t.Mul(&x1x2, x3).Inv(&t)   // 1/(x1*x2*x3)
-	y1.Mul(&t, x2).Mul(y1, x3) // 1/x1
-	y2.Mul(&t, x1).Mul(y2, x3) // 1/x2
-	y3.Mul(&t, &x1x2)          // 1/x3
-}
-
-// Set dest = x * x
-//
-// Allowed to overlap dest with x.
-//
-// Returns dest to allow chaining operations.
-func (dest *ExtensionFieldElement) Square(x *ExtensionFieldElement) *ExtensionFieldElement {
-	a := &x.A
-	b := &x.B
-
-	// We want to compute
-	//
-	// (a + bi)*(a + bi) = (a^2 - b^2) + 2abi.
-
-	var a2, a_plus_b, a_minus_b Fp751Element
-	fp751AddReduced(&a2, a, a)        // = a*R + a*R = 2*a*R
-	fp751AddReduced(&a_plus_b, a, b)  // = a*R + b*R = (a+b)*R
-	fp751SubReduced(&a_minus_b, a, b) // = a*R - b*R = (a-b)*R
-
-	var asq_minus_bsq, ab2 fp751X2
-	fp751Mul(&asq_minus_bsq, &a_plus_b, &a_minus_b) // = (a+b)*(a-b)*R*R = (a^2 - b^2)*R*R
-	fp751Mul(&ab2, &a2, b)                          // = 2*a*b*R*R
-
-	fp751MontgomeryReduce(&dest.A, &asq_minus_bsq) // = (a^2 - b^2)*R mod p
-	fp751MontgomeryReduce(&dest.B, &ab2)           // = 2*a*b*R mod p
-
-	return dest
-}
-
-// Set dest = lhs + rhs.
-//
-// Allowed to overlap lhs or rhs with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *ExtensionFieldElement) Add(lhs, rhs *ExtensionFieldElement) *ExtensionFieldElement {
-	fp751AddReduced(&dest.A, &lhs.A, &rhs.A)
-	fp751AddReduced(&dest.B, &lhs.B, &rhs.B)
-
-	return dest
-}
-
-// Set dest = lhs - rhs.
-//
-// Allowed to overlap lhs or rhs with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *ExtensionFieldElement) Sub(lhs, rhs *ExtensionFieldElement) *ExtensionFieldElement {
-	fp751SubReduced(&dest.A, &lhs.A, &rhs.A)
-	fp751SubReduced(&dest.B, &lhs.B, &rhs.B)
-
-	return dest
-}
-
-// If choice = 1u8, set (x,y) = (y,x). If choice = 0u8, set (x,y) = (x,y).
-//
-// Returns dest to allow chaining operations.
-func ExtensionFieldConditionalSwap(x, y *ExtensionFieldElement, choice uint8) {
-	fp751ConditionalSwap(&x.A, &y.A, choice)
-	fp751ConditionalSwap(&x.B, &y.B, choice)
-}
-
-// Returns true if lhs = rhs.  Takes variable time.
-func (lhs *ExtensionFieldElement) VartimeEq(rhs *ExtensionFieldElement) bool {
-	return lhs.A.vartimeEq(rhs.A) && lhs.B.vartimeEq(rhs.B)
-}
-
-// Convert the input to wire format.
-//
-// The output byte slice must be at least 188 bytes long.
-func (x *ExtensionFieldElement) ToBytes(output []byte) {
-	if len(output) < 188 {
-		panic("output byte slice too short, need 188 bytes")
-	}
-
-	var a,b Fp751Element
-	FromMontgomery(x, &a, &b)
-
-	// convert to bytes in little endian form. 8*12 = 96, but we drop the last two bytes
-	//  since p is 751 < 752=94*8 bits.
-	for i := 0; i < 94; i++ {
-		// set i = j*8 + k
-		j := i / 8
-		k := uint64(i % 8)
-
-		output[i] = byte(a[j] >> (8 * k))
-		output[i+94] = byte(b[j] >> (8 * k))
-	}
-}
-
-// Read 188 bytes into the given ExtensionFieldElement.
-//
-// It is an error to call this function if the input byte slice is less than 188 bytes long.
-func (x *ExtensionFieldElement) FromBytes(input []byte) {
-	if len(input) < 188 {
-		panic("input byte slice too short, need 188 bytes")
-	}
-
-	for i:=0; i<94; i++ {
-		j := i / 8
-		k := uint64(i % 8)
-		x.A[j] |= uint64(input[i]) << (8 * k)
-		x.B[j] |= uint64(input[i+94]) << (8 * k)
-	}
-
-	ToMontgomery(x)
-}
-
-// Converts values in x.A and x.B to Montgomery domain
-// x.A = x.A * R mod p
-// x.B = x.B * R mod p
-func ToMontgomery(x *ExtensionFieldElement) {
-	var aRR fp751X2
-
-	// convert to montgomery domain
-	fp751Mul(&aRR, &x.A, &montgomeryRsq)   // = a*R*R
-	fp751MontgomeryReduce(&x.A, &aRR)      // = a*R mod p
-	fp751Mul(&aRR, &x.B, &montgomeryRsq)
-	fp751MontgomeryReduce(&x.B, &aRR)
-}
-
-// Converts values in x.A and x.B from Montgomery domain
-// a = x.A mod p
-// b = x.B mod p
-//
-// After returning from the call x is not modified.
-func FromMontgomery(x *ExtensionFieldElement, a,b *Fp751Element) {
-	var aR fp751X2
-
-	// convert from montgomery domain
-	copy(aR[:], x.A[:])
-	fp751MontgomeryReduce(a, &aR) 	 // = a mod p in [0, 2p)
-	fp751StrongReduce(a)          	 // = a mod p in [0, p)
-
-	for i:=range(aR) {
-		aR[i] = 0
-	}
-	copy(aR[:], x.B[:])
-	fp751MontgomeryReduce(b, &aR)
-	fp751StrongReduce(b)
-}
-
-//------------------------------------------------------------------------------
-// Prime Field
-//------------------------------------------------------------------------------
-
-// Represents an element of the prime field F_p.
-type PrimeFieldElement struct {
-	// This field element is in Montgomery form, so that the value `A` is
-	// represented by `aR mod p`.
-	A Fp751Element
-}
-
-var zeroPrimeField = PrimeFieldElement{
-	A: Fp751Element{0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0},
-}
-
-var onePrimeField = PrimeFieldElement{
-	A: Fp751Element{0x249ad, 0x0, 0x0, 0x0, 0x0, 0x8310000000000000, 0x5527b1e4375c6c66, 0x697797bf3f4f24d0, 0xc89db7b2ac5c4e2e, 0x4ca4b439d2076956, 0x10f7926c7512c7e9, 0x2d5b24bce5e2},
-}
-
-// Set dest = lhs * rhs.
-//
-// Allowed to overlap lhs or rhs with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) Mul(lhs, rhs *PrimeFieldElement) *PrimeFieldElement {
-	a := &lhs.A // = a*R
-	b := &rhs.A // = b*R
-
-	var ab fp751X2
-	fp751Mul(&ab, a, b)                 // = a*b*R*R
-	fp751MontgomeryReduce(&dest.A, &ab) // = a*b*R mod p
-
-	return dest
-}
-
-// Set dest = x^(2^k), for k >= 1, by repeated squarings.
-//
-// Allowed to overlap x with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) Pow2k(x *PrimeFieldElement, k uint8) *PrimeFieldElement {
-	dest.Mul(x, x)
-	for i := uint8(1); i < k; i++ {
-		dest.Mul(dest, dest)
-	}
-
-	return dest
-}
-
-// Set dest = x^((p-3)/4).  If x is square, this is 1/sqrt(x).
-//
-// Allowed to overlap x with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) P34(x *PrimeFieldElement) *PrimeFieldElement {
-	// Sliding-window strategy computed with Sage, awk, sed, and tr.
-	//
-	// This performs sum(powStrategy) = 744 squarings and len(mulStrategy)
-	// = 137 multiplications, in addition to 1 squaring and 15
-	// multiplications to build a lookup table.
-	//
-	// In total this is 745 squarings, 152 multiplications.  Since squaring
-	// is not implemented for the prime field, this is 897 multiplications
-	// in total.
-	powStrategy := [137]uint8{5, 7, 6, 2, 10, 4, 6, 9, 8, 5, 9, 4, 7, 5, 5, 4, 8, 3, 9, 5, 5, 4, 10, 4, 6, 6, 6, 5, 8, 9, 3, 4, 9, 4, 5, 6, 6, 2, 9, 4, 5, 5, 5, 7, 7, 9, 4, 6, 4, 8, 5, 8, 6, 6, 2, 9, 7, 4, 8, 8, 8, 4, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 2}
-	mulStrategy := [137]uint8{31, 23, 21, 1, 31, 7, 7, 7, 9, 9, 19, 15, 23, 23, 11, 7, 25, 5, 21, 17, 11, 5, 17, 7, 11, 9, 23, 9, 1, 19, 5, 3, 25, 15, 11, 29, 31, 1, 29, 11, 13, 9, 11, 27, 13, 19, 15, 31, 3, 29, 23, 31, 25, 11, 1, 21, 19, 15, 15, 21, 29, 13, 23, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 3}
-	initialMul := uint8(27)
-
-	// Build a lookup table of odd multiples of x.
-	lookup := [16]PrimeFieldElement{}
-	xx := &PrimeFieldElement{}
-	xx.Mul(x, x) // Set xx = x^2
-	lookup[0] = *x
-	for i := 1; i < 16; i++ {
-		lookup[i].Mul(&lookup[i-1], xx)
-	}
-	// Now lookup = {x, x^3, x^5, ... }
-	// so that lookup[i] = x^{2*i + 1}
-	// so that lookup[k/2] = x^k, for odd k
-
-	*dest = lookup[initialMul/2]
-	for i := uint8(0); i < 137; i++ {
-		dest.Pow2k(dest, powStrategy[i])
-		dest.Mul(dest, &lookup[mulStrategy[i]/2])
-	}
-
-	return dest
-}
-
-//------------------------------------------------------------------------------
-// Internals
-//------------------------------------------------------------------------------
-
-const fp751NumWords = 12
-
-// (2^768)^2 mod p
-// This can't be a constant because Go doesn't allow array constants, so try
-// not to modify it.
-var montgomeryRsq = Fp751Element{2535603850726686808, 15780896088201250090, 6788776303855402382, 17585428585582356230, 5274503137951975249, 2266259624764636289, 11695651972693921304, 13072885652150159301, 4908312795585420432, 6229583484603254826, 488927695601805643, 72213483953973}
-
-// Internal representation of an element of the base field F_p.
-//
-// This type is distinct from PrimeFieldElement in that no particular meaning
-// is assigned to the representation -- it could represent an element in
-// Montgomery form, or not.  Tracking the meaning of the field element is left
-// to higher types.
-type Fp751Element [fp751NumWords]uint64
-
-// Represents an intermediate product of two elements of the base field F_p.
-type fp751X2 [2 * fp751NumWords]uint64
-
-func (x Fp751Element) vartimeEq(y Fp751Element) bool {
-	fp751StrongReduce(&x)
-	fp751StrongReduce(&y)
-	eq := true
-	for i := 0; i < fp751NumWords; i++ {
-		eq = (x[i] == y[i]) && eq
-	}
-
-	return eq
-}

p751toolbox/field_test.go

@@ -1,419 +0,0 @@
-package p751toolbox
-
-import (
-	"math/big"
-	"math/rand"
-	"reflect"
-	"testing"
-	"testing/quick"
-)
-
-var quickCheckScaleFactor = uint8(3)
-var quickCheckConfig = &quick.Config{MaxCount: (1 << (12 + quickCheckScaleFactor))}
-
-var cln16prime, _ = new(big.Int).SetString("10354717741769305252977768237866805321427389645549071170116189679054678940682478846502882896561066713624553211618840202385203911976522554393044160468771151816976706840078913334358399730952774926980235086850991501872665651576831", 10)
-
-// Convert an Fp751Element to a big.Int for testing.  Because this is only
-// for testing, no big.Int to Fp751Element conversion is provided.
-
-func radix64ToBigInt(x []uint64) *big.Int {
-	radix := new(big.Int)
-	// 2^64
-	radix.UnmarshalText(([]byte)("18446744073709551616"))
-
-	base := new(big.Int).SetUint64(1)
-	val := new(big.Int).SetUint64(0)
-	tmp := new(big.Int)
-
-	for _, xi := range x {
-		tmp.SetUint64(xi)
-		tmp.Mul(tmp, base)
-		val.Add(val, tmp)
-		base.Mul(base, radix)
-	}
-
-	return val
-}
-
-func VartimeEq(x,y *PrimeFieldElement) bool {
-	return x.A.vartimeEq(y.A)
-}
-
-func (x *Fp751Element) toBigInt() *big.Int {
-	// Convert from Montgomery form
-	return x.toBigIntFromMontgomeryForm()
-}
-
-func (x *Fp751Element) toBigIntFromMontgomeryForm() *big.Int {
-	// Convert from Montgomery form
-	a := Fp751Element{}
-	aR := fp751X2{}
-	copy(aR[:], x[:])              // = a*R
-	fp751MontgomeryReduce(&a, &aR) // = a mod p  in [0,2p)
-	fp751StrongReduce(&a)          // = a mod p  in [0,p)
-	return radix64ToBigInt(a[:])
-}
-
-func TestPrimeFieldElementToBigInt(t *testing.T) {
-	// Chosen so that p < xR < 2p
-	x := PrimeFieldElement{A: Fp751Element{
-		1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 140737488355328,
-	}}
-	// Computed using Sage:
-	// sage: p = 2^372 * 3^239 - 1
-	// sage: R = 2^768
-	// sage: from_radix_64 = lambda xs: sum((xi * (2**64)**i for i,xi in enumerate(xs)))
-	// sage: xR = from_radix_64([1]*11 + [2^47])
-	// sage: assert(p < xR)
-	// sage: assert(xR < 2*p)
-	// sage: (xR / R) % p
-	xBig, _ := new(big.Int).SetString("4469946751055876387821312289373600189787971305258234719850789711074696941114031433609871105823930699680637820852699269802003300352597419024286385747737509380032982821081644521634652750355306547718505685107272222083450567982240", 10)
-	if xBig.Cmp(x.A.toBigInt()) != 0 {
-		t.Error("Expected", xBig, "found", x.A.toBigInt())
-	}
-}
-
-func generateFp751(rand *rand.Rand) Fp751Element {
-	// Generation strategy: low limbs taken from [0,2^64); high limb
-	// taken from smaller range
-	//
-	// Size hint is ignored since all elements are fixed size.
-	//
-	// Field elements taken in range [0,2p).  Emulate this by capping
-	// the high limb by the top digit of 2*p-1:
-	//
-	// sage: (2*p-1).digits(2^64)[-1]
-	// 246065832128056
-	//
-	// This still allows generating values >= 2p, but hopefully that
-	// excess is OK (and if it's not, we'll find out, because it's for
-	// testing...)
-	//
-	highLimb := rand.Uint64() % 246065832128056
-
-	return Fp751Element{
-		rand.Uint64(),
-		rand.Uint64(),
-		rand.Uint64(),
-		rand.Uint64(),
-		rand.Uint64(),
-		rand.Uint64(),
-		rand.Uint64(),
-		rand.Uint64(),
-		rand.Uint64(),
-		rand.Uint64(),
-		rand.Uint64(),
-		highLimb,
-	}
-}
-
-func (x PrimeFieldElement) Generate(rand *rand.Rand, size int) reflect.Value {
-	return reflect.ValueOf(PrimeFieldElement{A: generateFp751(rand)})
-}
-
-func (x ExtensionFieldElement) Generate(rand *rand.Rand, size int) reflect.Value {
-	return reflect.ValueOf(ExtensionFieldElement{A: generateFp751(rand), B: generateFp751(rand)})
-}
-
-//------------------------------------------------------------------------------
-// Extension Field
-//------------------------------------------------------------------------------
-
-func TestOneExtensionFieldToBytes(t *testing.T) {
-	var x ExtensionFieldElement
-	var xBytes [188]byte
-
-	x.One()
-	x.ToBytes(xBytes[:])
-	if xBytes[0] != 1 {
-		t.Error("Expected 1, got", xBytes[0])
-	}
-	for i := 1; i < 188; i++ {
-		if xBytes[i] != 0 {
-			t.Error("Expected 0, got", xBytes[0])
-		}
-	}
-}
-
-func TestExtensionFieldElementToBytesRoundTrip(t *testing.T) {
-	roundTrips := func(x ExtensionFieldElement) bool {
-		var xBytes [188]byte
-		var xPrime ExtensionFieldElement
-		x.ToBytes(xBytes[:])
-		xPrime.FromBytes(xBytes[:])
-
-		return x.VartimeEq(&xPrime)
-	}
-
-	if err := quick.Check(roundTrips, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-func TestExtensionFieldElementMulDistributesOverAdd(t *testing.T) {
-	mulDistributesOverAdd := func(x, y, z ExtensionFieldElement) bool {
-		// Compute t1 = (x+y)*z
-		t1 := new(ExtensionFieldElement)
-		t1.Add(&x, &y)
-		t1.Mul(t1, &z)
-
-		// Compute t2 = x*z + y*z
-		t2 := new(ExtensionFieldElement)
-		t3 := new(ExtensionFieldElement)
-		t2.Mul(&x, &z)
-		t3.Mul(&y, &z)
-		t2.Add(t2, t3)
-
-		return t1.VartimeEq(t2)
-	}
-
-	if err := quick.Check(mulDistributesOverAdd, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-func TestExtensionFieldElementMulIsAssociative(t *testing.T) {
-	isAssociative := func(x, y, z ExtensionFieldElement) bool {
-		// Compute t1 = (x*y)*z
-		t1 := new(ExtensionFieldElement)
-		t1.Mul(&x, &y)
-		t1.Mul(t1, &z)
-
-		// Compute t2 = (y*z)*x
-		t2 := new(ExtensionFieldElement)
-		t2.Mul(&y, &z)
-		t2.Mul(t2, &x)
-
-		return t1.VartimeEq(t2)
-	}
-
-	if err := quick.Check(isAssociative, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-func TestExtensionFieldElementSquareMatchesMul(t *testing.T) {
-	sqrMatchesMul := func(x ExtensionFieldElement) bool {
-		// Compute t1 = (x*x)
-		t1 := new(ExtensionFieldElement)
-		t1.Mul(&x, &x)
-
-		// Compute t2 = x^2
-		t2 := new(ExtensionFieldElement)
-		t2.Square(&x)
-
-		return t1.VartimeEq(t2)
-	}
-
-	if err := quick.Check(sqrMatchesMul, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-func TestExtensionFieldElementInv(t *testing.T) {
-	inverseIsCorrect := func(x ExtensionFieldElement) bool {
-		z := new(ExtensionFieldElement)
-		z.Inv(&x)
-
-		// Now z = (1/x), so (z * x) * x == x
-		z.Mul(z, &x)
-		z.Mul(z, &x)
-
-		return z.VartimeEq(&x)
-	}
-
-	// This is more expensive; run fewer tests
-	var quickCheckConfig = &quick.Config{MaxCount: (1 << (8 + quickCheckScaleFactor))}
-	if err := quick.Check(inverseIsCorrect, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-func TestExtensionFieldElementBatch3Inv(t *testing.T) {
-	batchInverseIsCorrect := func(x1, x2, x3 ExtensionFieldElement) bool {
-		var x1Inv, x2Inv, x3Inv ExtensionFieldElement
-		x1Inv.Inv(&x1)
-		x2Inv.Inv(&x2)
-		x3Inv.Inv(&x3)
-
-		var y1, y2, y3 ExtensionFieldElement
-		ExtensionFieldBatch3Inv(&x1, &x2, &x3, &y1, &y2, &y3)
-
-		return (y1.VartimeEq(&x1Inv) && y2.VartimeEq(&x2Inv) && y3.VartimeEq(&x3Inv))
-	}
-
-	// This is more expensive; run fewer tests
-	var quickCheckConfig = &quick.Config{MaxCount: (1 << (5 + quickCheckScaleFactor))}
-	if err := quick.Check(batchInverseIsCorrect, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-//------------------------------------------------------------------------------
-// Prime Field
-//------------------------------------------------------------------------------
-func TestPrimeFieldElementMulVersusBigInt(t *testing.T) {
-	mulMatchesBigInt := func(x, y PrimeFieldElement) bool {
-		z := new(PrimeFieldElement)
-		z.Mul(&x, &y)
-
-		check := new(big.Int)
-		check.Mul(x.A.toBigInt(), y.A.toBigInt())
-		check.Mod(check, cln16prime)
-
-		return check.Cmp(z.A.toBigInt()) == 0
-	}
-
-	if err := quick.Check(mulMatchesBigInt, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-func TestPrimeFieldElementP34VersusBigInt(t *testing.T) {
-	var p34, _ = new(big.Int).SetString("2588679435442326313244442059466701330356847411387267792529047419763669735170619711625720724140266678406138302904710050596300977994130638598261040117192787954244176710019728333589599932738193731745058771712747875468166412894207", 10)
-	p34MatchesBigInt := func(x PrimeFieldElement) bool {
-		z := new(PrimeFieldElement)
-		z.P34(&x)
-
-		check := x.A.toBigInt()
-		check.Exp(check, p34, cln16prime)
-
-		return check.Cmp(z.A.toBigInt()) == 0
-	}
-
-	// This is more expensive; run fewer tests
-	var quickCheckConfig = &quick.Config{MaxCount: (1 << (8 + quickCheckScaleFactor))}
-	if err := quick.Check(p34MatchesBigInt, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-// Package-level storage for this field element is intended to deter
-// compiler optimizations.
-var benchmarkFp751Element Fp751Element
-var benchmarkFp751X2 fp751X2
-var bench_x = Fp751Element{17026702066521327207, 5108203422050077993, 10225396685796065916, 11153620995215874678, 6531160855165088358, 15302925148404145445, 1248821577836769963, 9789766903037985294, 7493111552032041328, 10838999828319306046, 18103257655515297935, 27403304611634}
-var bench_y = Fp751Element{4227467157325093378, 10699492810770426363, 13500940151395637365, 12966403950118934952, 16517692605450415877, 13647111148905630666, 14223628886152717087, 7167843152346903316, 15855377759596736571, 4300673881383687338, 6635288001920617779, 30486099554235}
-var bench_z = fp751X2{1595347748594595712, 10854920567160033970, 16877102267020034574, 12435724995376660096, 3757940912203224231, 8251999420280413600, 3648859773438820227, 17622716832674727914, 11029567000887241528, 11216190007549447055, 17606662790980286987, 4720707159513626555, 12887743598335030915, 14954645239176589309, 14178817688915225254, 1191346797768989683, 12629157932334713723, 6348851952904485603, 16444232588597434895, 7809979927681678066, 14642637672942531613, 3092657597757640067, 10160361564485285723, 240071237}
-
-func BenchmarkExtensionFieldElementMul(b *testing.B) {
-	z := &ExtensionFieldElement{A: bench_x, B: bench_y}
-	w := new(ExtensionFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Mul(z, z)
-	}
-}
-
-func BenchmarkExtensionFieldElementInv(b *testing.B) {
-	z := &ExtensionFieldElement{A: bench_x, B: bench_y}
-	w := new(ExtensionFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Inv(z)
-	}
-}
-
-func BenchmarkExtensionFieldElementSquare(b *testing.B) {
-	z := &ExtensionFieldElement{A: bench_x, B: bench_y}
-	w := new(ExtensionFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Square(z)
-	}
-}
-
-func BenchmarkExtensionFieldElementAdd(b *testing.B) {
-	z := &ExtensionFieldElement{A: bench_x, B: bench_y}
-	w := new(ExtensionFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Add(z, z)
-	}
-}
-
-func BenchmarkExtensionFieldElementSub(b *testing.B) {
-	z := &ExtensionFieldElement{A: bench_x, B: bench_y}
-	w := new(ExtensionFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Sub(z, z)
-	}
-}
-
-func BenchmarkPrimeFieldElementMul(b *testing.B) {
-	z := &PrimeFieldElement{A: bench_x}
-	w := new(PrimeFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Mul(z, z)
-	}
-}
-
-// --- field operation functions
-
-func BenchmarkFp751Multiply(b *testing.B) {
-	for n := 0; n < b.N; n++ {
-		fp751Mul(&benchmarkFp751X2, &bench_x, &bench_y)
-	}
-}
-
-func BenchmarkFp751MontgomeryReduce(b *testing.B) {
-	z := bench_z
-
-	// This benchmark actually computes garbage, because
-	// fp751MontgomeryReduce mangles its input, but since it's
-	// constant-time that shouldn't matter for the benchmarks.
-	for n := 0; n < b.N; n++ {
-		fp751MontgomeryReduce(&benchmarkFp751Element, &z)
-	}
-}
-
-func BenchmarkFp751AddReduced(b *testing.B) {
-	for n := 0; n < b.N; n++ {
-		fp751AddReduced(&benchmarkFp751Element, &bench_x, &bench_y)
-	}
-}
-
-func BenchmarkFp751SubReduced(b *testing.B) {
-	for n := 0; n < b.N; n++ {
-		fp751SubReduced(&benchmarkFp751Element, &bench_x, &bench_y)
-	}
-}
-
-func BenchmarkFp751ConditionalSwap(b *testing.B) {
-	x, y := bench_x, bench_y
-	for n := 0; n < b.N; n++ {
-		fp751ConditionalSwap(&x, &y, 1)
-		fp751ConditionalSwap(&x, &y, 0)
-	}
-}
-
-func BenchmarkFp751StrongReduce(b *testing.B) {
-	x := bench_x
-	for n := 0; n < b.N; n++ {
-		fp751StrongReduce(&x)
-	}
-}
-
-func BenchmarkFp751AddLazy(b *testing.B) {
-	var z Fp751Element
-	x, y := bench_x, bench_y
-	for n := 0; n < b.N; n++ {
-		fp751AddLazy(&z, &x, &y)
-	}
-}
-
-func BenchmarkFp751X2AddLazy(b *testing.B) {
-	x, y, z := bench_z, bench_z, bench_z
-	for n := 0; n < b.N; n++ {
-		fp751X2AddLazy(&x, &y, &z)
-	}
-}
-
-func BenchmarkFp751X2SubLazy(b *testing.B) {
-	x, y, z := bench_z, bench_z, bench_z
-	for n := 0; n < b.N; n++ {
-		fp751X2SubLazy(&x, &y, &z)
-	}
-}

p751toolbox/isogeny.go

@@ -1,144 +0,0 @@
-package p751toolbox
-
-// Interface for working with isogenies.
-type Isogeny interface {
-	// Given a torsion point on a curve computes isogenous curve.
-	// Returns curve coefficients (A:C), so that E_(A/C) = E_(A/C)/<P>,
-	// where P is a provided projective point. Sets also isogeny constants
-	// that are needed for isogeny evaluation.
-	GenerateCurve(*ProjectivePoint) CurveCoefficientsEquiv
-	// Evaluates isogeny at caller provided point. Requires isogeny curve constants
-	// to be earlier computed by GenerateCurve.
-	EvaluatePoint(*ProjectivePoint) ProjectivePoint
-}
-
-// Stores Isogeny 4 curve constants
-type isogeny4 struct {
-	isogeny3
-	K3 ExtensionFieldElement
-}
-
-// Stores Isogeny 3 curve constants
-type isogeny3 struct {
-	K1 ExtensionFieldElement
-	K2 ExtensionFieldElement
-}
-
-// Constructs isogeny4 objects
-func NewIsogeny4() Isogeny {
-	return new(isogeny4)
-}
-
-// Constructs isogeny3 objects
-func NewIsogeny3() Isogeny {
-	return new(isogeny3)
-}
-
-// Given a three-torsion point p = x(PB) on the curve E_(A:C), construct the
-// three-isogeny phi : E_(A:C) -> E_(A:C)/<P_3> = E_(A':C').
-//
-// Input: (XP_3: ZP_3), where P_3 has exact order 3 on E_A/C
-// Output: * Curve coordinates (A' + 2C', A' - 2C') corresponding to E_A'/C' = A_E/C/<P3>
-//		   * Isogeny phi with constants in F_p^2
-func (phi *isogeny3) GenerateCurve(p *ProjectivePoint) CurveCoefficientsEquiv {
-	var t0, t1, t2, t3, t4 ExtensionFieldElement
-	var coefEq CurveCoefficientsEquiv
-	var K1, K2 = &phi.K1, &phi.K2
-
-	K1.Sub(&p.X, &p.Z)           // K1 = XP3 - ZP3
-	t0.Square(K1)                // t0 = K1^2
-	K2.Add(&p.X, &p.Z)           // K2 = XP3 + ZP3
-	t1.Square(K2)                // t1 = K2^2
-	t2.Add(&t0, &t1)             // t2 = t0 + t1
-	t3.Add(K1, K2)               // t3 = K1 + K2
-	t3.Square(&t3)               // t3 = t3^2
-	t3.Sub(&t3, &t2)             // t3 = t3 - t2
-	t2.Add(&t1, &t3)             // t2 = t1 + t3
-	t3.Add(&t3, &t0)             // t3 = t3 + t0
-	t4.Add(&t3, &t0)             // t4 = t3 + t0
-	t4.Add(&t4, &t4)             // t4 = t4 + t4
-	t4.Add(&t1, &t4)             // t4 = t1 + t4
-	coefEq.C.Mul(&t2, &t4)       // A24m = t2 * t4
-	t4.Add(&t1, &t2)             // t4 = t1 + t2
-	t4.Add(&t4, &t4)             // t4 = t4 + t4
-	t4.Add(&t0, &t4)             // t4 = t0 + t4
-	t4.Mul(&t3, &t4)             // t4 = t3 * t4
-	t0.Sub(&t4, &coefEq.C)       // t0 = t4 - A24m
-	coefEq.A.Add(&coefEq.C, &t0) // A24p = A24m + t0
-	return coefEq
-}
-
-// Given a 3-isogeny phi and a point pB = x(PB), compute x(QB), the x-coordinate
-// of the image QB = phi(PB) of PB under phi : E_(A:C) -> E_(A':C').
-//
-// The output xQ = x(Q) is then a point on the curve E_(A':C'); the curve
-// parameters are returned by the GenerateCurve function used to construct phi.
-func (phi *isogeny3) EvaluatePoint(p *ProjectivePoint) ProjectivePoint {
-	var t0, t1, t2 ExtensionFieldElement
-	var q ProjectivePoint
-	var K1, K2 = &phi.K1, &phi.K2
-	var px, pz = &p.X, &p.Z
-
-	t0.Add(px, pz)   // t0 = XQ + ZQ
-	t1.Sub(px, pz)   // t1 = XQ - ZQ
-	t0.Mul(K1, &t0)  // t2 = K1 * t0
-	t1.Mul(K2, &t1)  // t1 = K2 * t1
-	t2.Add(&t0, &t1) // t2 = t0 + t1
-	t0.Sub(&t1, &t0) // t0 = t1 - t0
-	t2.Square(&t2)   // t2 = t2 ^ 2
-	t0.Square(&t0)   // t0 = t0 ^ 2
-	q.X.Mul(px, &t2) // XQ'= XQ * t2
-	q.Z.Mul(pz, &t0) // ZQ'= ZQ * t0
-	return q
-}
-
-// Given a four-torsion point p = x(PB) on the curve E_(A:C), construct the
-// four-isogeny phi : E_(A:C) -> E_(A:C)/<P_4> = E_(A':C').
-//
-// Input: (XP_4: ZP_4), where P_4 has exact order 4 on E_A/C
-// Output: * Curve coordinates (A' + 2C', 4C') corresponding to E_A'/C' = A_E/C/<P4>
-//		   * Isogeny phi with constants in F_p^2
-func (phi *isogeny4) GenerateCurve(p *ProjectivePoint) CurveCoefficientsEquiv {
-	var coefEq CurveCoefficientsEquiv
-	var xp4, zp4 = &p.X, &p.Z
-	var K1, K2, K3 = &phi.K1, &phi.K2, &phi.K3
-
-	K2.Sub(xp4, zp4)
-	K3.Add(xp4, zp4)
-	K1.Square(zp4)
-	K1.Add(K1, K1)
-	coefEq.C.Square(K1)
-	K1.Add(K1, K1)
-	coefEq.A.Square(xp4)
-	coefEq.A.Add(&coefEq.A, &coefEq.A)
-	coefEq.A.Square(&coefEq.A)
-	return coefEq
-}
-
-// Given a 4-isogeny phi and a point xP = x(P), compute x(Q), the x-coordinate
-// of the image Q = phi(P) of P under phi : E_(A:C) -> E_(A':C').
-//
-// Input: Isogeny returned by GenerateCurve and point q=(Qx,Qz) from E0_A/C
-// Output: Corresponding point q from E1_A'/C', where E1 is 4-isogenous to E0
-func (phi *isogeny4) EvaluatePoint(p *ProjectivePoint) ProjectivePoint {
-	var t0, t1 ExtensionFieldElement
-	var q = *p
-	var xq, zq = &q.X, &q.Z
-	var K1, K2, K3 = &phi.K1, &phi.K2, &phi.K3
-
-	t0.Add(xq, zq)
-	t1.Sub(xq, zq)
-	xq.Mul(&t0, K2)
-	zq.Mul(&t1, K3)
-	t0.Mul(&t0, &t1)
-	t0.Mul(&t0, K1)
-	t1.Add(xq, zq)
-	zq.Sub(xq, zq)
-	t1.Square(&t1)
-	zq.Square(zq)
-	xq.Add(&t0, &t1)
-	t0.Sub(zq, &t0)
-	xq.Mul(xq, &t1)
-	zq.Mul(zq, &t0)
-	return q
-}

p751toolbox/isogeny_test.go

@@ -1,114 +0,0 @@
-package p751toolbox
-
-import (
-	"testing"
-)
-
-func TestFourIsogenyVersusSage(t *testing.T) {
-	var xR, xP4, resPhiXr, expPhiXr ProjectivePoint
-	var phi = NewIsogeny4()
-
-	// sage: p = 2^372 * 3^239 - 1; Fp = GF(p)
-	// sage: R.<x> = Fp[]
-	// sage: Fp2 = Fp.extension(x^2 + 1, 'i')
-	// sage: i = Fp2.gen()
-	// sage: E0Fp = EllipticCurve(Fp, [0,0,0,1,0])
-	// sage: E0Fp2 = EllipticCurve(Fp2, [0,0,0,1,0])
-	// sage: x_PA = 11
-	// sage: y_PA = -Fp(11^3 + 11).sqrt()
-	// sage: x_PB = 6
-	// sage: y_PB = -Fp(6^3 + 6).sqrt()
-	// sage: P_A = 3^239 * E0Fp((x_PA,y_PA))
-	// sage: P_B = 2^372 * E0Fp((x_PB,y_PB))
-	// sage: def tau(P):
-	// ....:     return E0Fp2( (-P.xy()[0], i*P.xy()[1]))
-	// ....:
-	// sage: m_B = 3*randint(0,3^238)
-	// sage: m_A = 2*randint(0,2^371)
-	// sage: R_A = E0Fp2(P_A) + m_A*tau(P_A)
-	// sage: def y_recover(x, a):
-	// ....:     return (x**3 + a*x**2 + x).sqrt()
-	// ....:
-	// sage: first_4_torsion_point = E0Fp2(1, y_recover(Fp2(1),0))
-	// sage: sage_first_4_isogeny = E0Fp2.isogeny(first_4_torsion_point)
-	// sage: a = Fp2(0)
-	// sage: E1A = EllipticCurve(Fp2, [0,(2*(a+6))/(a-2),0,1,0])
-	// sage: sage_isomorphism = sage_first_4_isogeny.codomain().isomorphism_to(E1A)
-	// sage: isogenized_R_A = sage_isomorphism(sage_first_4_isogeny(R_A))
-	// sage: P_4 = (2**(372-4))*isogenized_R_A
-	// sage: P_4._order = 4 #otherwise falls back to generic group methods for order
-	// sage: X4, Z4 = P_4.xy()[0], 1
-	// sage: phi4 = EllipticCurveIsogeny(E1A, P_4, None, 4)
-	// sage: E2A_sage = phi4.codomain() # not in monty form
-	// sage: Aprime, Cprime = 2*(2*X4^4 - Z4^4), Z4^4
-	// sage: E2A = EllipticCurve(Fp2, [0,Aprime/Cprime,0,1,0])
-	// sage: sage_iso = E2A_sage.isomorphism_to(E2A)
-	// sage: isogenized2_R_A = sage_iso(phi4(isogenized_R_A))
-
-	xP4.FromAffine(&ExtensionFieldElement{
-		A: Fp751Element{0x2afd75a913f3d5e7, 0x2918fba06f88c9ab, 0xa4ac4dc7cb526f05, 0x2d19e9391a607300, 0x7a79e2b34091b54, 0x3ad809dcb42f1792, 0xd46179328bd6402a, 0x1afa73541e2c4f3f, 0xf602d73ace9bdbd8, 0xd77ac58f6bab7004, 0x4689d97f6793b3b3, 0x4f26b00e42b7},
-		B: Fp751Element{0x6cdf918dafdcb890, 0x666f273cc29cfae2, 0xad00fcd31ba618e2, 0x5fbcf62bef2f6a33, 0xf408bb88318e5098, 0x84ab97849453d175, 0x501bbfcdcfb8e1ac, 0xf2370098e6b5542c, 0xc7dc73f5f0f6bd32, 0xdd76dcd86729d1cf, 0xca22c905029996e4, 0x5cf4a9373de3}})
-	xR.FromAffine(&ExtensionFieldElement{
-		A: Fp751Element{0xff99e76f78da1e05, 0xdaa36bd2bb8d97c4, 0xb4328cee0a409daf, 0xc28b099980c5da3f, 0xf2d7cd15cfebb852, 0x1935103dded6cdef, 0xade81528de1429c3, 0x6775b0fa90a64319, 0x25f89817ee52485d, 0x706e2d00848e697, 0xc4958ec4216d65c0, 0xc519681417f},
-		B: Fp751Element{0x742fe7dde60e1fb9, 0x801a3c78466a456b, 0xa9f945b786f48c35, 0x20ce89e1b144348f, 0xf633970b7776217e, 0x4c6077a9b38976e5, 0x34a513fc766c7825, 0xacccba359b9cd65, 0xd0ca8383f0fd0125, 0x77350437196287a, 0x9fe1ad7706d4ea21, 0x4d26129ee42d}})
-	expPhiXr.FromAffine(&ExtensionFieldElement{
-		A: Fp751Element{0x111efd8bd0b7a01e, 0x6ab75a4f3789ca9b, 0x939dbe518564cac4, 0xf9eeaba1601d0434, 0x8d41f8ba6edac998, 0xfcd2557efe9aa170, 0xb3c3549c098b7844, 0x52874fef6f81127c, 0xb2b9ac82aa518bb3, 0xee70820230520a86, 0xd4012b7f5efb184a, 0x573e4536329b},
-		B: Fp751Element{0xa99952281e932902, 0x569a89a571f2c7b1, 0x6150143846ba3f6b, 0x11fd204441e91430, 0x7f469bd55c9b07b, 0xb72db8b9de35b161, 0x455a9a37a940512a, 0xb0cff7670abaf906, 0x18c785b7583375fe, 0x603ab9ca403c9148, 0xab54ba3a6e6c62c1, 0x2726d7d57c4f}})
-
-	phi.GenerateCurve(&xP4)
-	resPhiXr = phi.EvaluatePoint(&xR)
-	if !expPhiXr.VartimeEq(&resPhiXr) {
-		t.Error("\nExpected\n", expPhiXr.ToAffine(), "\nfound\n", resPhiXr.ToAffine())
-	}
-}
-
-func TestThreeIsogenyVersusSage(t *testing.T) {
-	var xR, xP3, resPhiXr, expPhiXr ProjectivePoint
-	var phi = NewIsogeny3()
-
-	// sage: %colors Linux
-	// sage: p = 2^372 * 3^239 - 1; Fp = GF(p)
-	// sage: R.<x> = Fp[]
-	// sage: Fp2 = Fp.extension(x^2 + 1, 'i')
-	// sage: i = Fp2.gen()
-	// sage: E0Fp = EllipticCurve(Fp, [0,0,0,1,0])
-	// sage: E0Fp2 = EllipticCurve(Fp2, [0,0,0,1,0])
-	// sage: x_PA = 11
-	// sage: y_PA = -Fp(11^3 + 11).sqrt()
-	// sage: x_PB = 6
-	// sage: y_PB = -Fp(6^3 + 6).sqrt()
-	// sage: P_A = 3^239 * E0Fp((x_PA,y_PA))
-	// sage: P_B = 2^372 * E0Fp((x_PB,y_PB))
-	// sage: def tau(P):
-	// ....:     return E0Fp2( (-P.xy()[0], i*P.xy()[1]))
-	// ....:
-	// sage: m_B = 3*randint(0,3^238)
-	// sage: R_B = E0Fp2(P_B) + m_B*tau(P_B)
-	// sage: P_3 = (3^238)*R_B
-	// sage: def three_isog(P_3, P):
-	// ....:     X3, Z3 = P_3.xy()[0], 1
-	// ....:     XP, ZP = P.xy()[0], 1
-	// ....:     x = (XP*(X3*XP - Z3*ZP)^2)/(ZP*(Z3*XP - X3*ZP)^2)
-	// ....:     A3, C3 = (Z3^4 + 9*X3^2*(2*Z3^2 - 3*X3^2)), 4*X3*Z3^3
-	// ....:     cod = EllipticCurve(Fp2, [0,A3/C3,0,1,0])
-	// ....:     return cod.lift_x(x)
-	// ....:
-	// sage: isogenized_R_B = three_isog(P_3, R_B)
-
-	xR.FromAffine(&ExtensionFieldElement{
-		A: Fp751Element{0xbd0737ed5cc9a3d7, 0x45ae6d476517c101, 0x6f228e9e7364fdb2, 0xbba4871225b3dbd, 0x6299ccd2e5da1a07, 0x38488fe4af5f2d0e, 0xec23cae5a86e980c, 0x26c804ba3f1edffa, 0xfbbed81932df60e5, 0x7e00e9d182ae9187, 0xc7654abb66d05f4b, 0x262d0567237b},
-		B: Fp751Element{0x3a3b5b6ad0b2ac33, 0x246602b5179127d3, 0x502ae0e9ad65077d, 0x10a3a37237e1bf70, 0x4a1ab9294dd05610, 0xb0f3adac30fe1fa6, 0x341995267faf70cb, 0xa14dd94d39cf4ec1, 0xce4b7527d1bf5568, 0xe0410423ed45c7e4, 0x38011809b6425686, 0x28f52472ebed}})
-	xP3.FromAffine(&ExtensionFieldElement{
-		A: Fp751Element{0x7bb7a4a07b0788dc, 0xdc36a3f6607b21b0, 0x4750e18ee74cf2f0, 0x464e319d0b7ab806, 0xc25aa44c04f758ff, 0x392e8521a46e0a68, 0xfc4e76b63eff37df, 0x1f3566d892e67dd8, 0xf8d2eb0f73295e65, 0x457b13ebc470bccb, 0xfda1cc9efef5be33, 0x5dbf3d92cc02},
-		B: Fp751Element{0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0}})
-	expPhiXr.FromAffine(&ExtensionFieldElement{
-		A: Fp751Element{0x286db7d75913c5b1, 0xcb2049ad50189220, 0xccee90ef765fa9f4, 0x65e52ce2730e7d88, 0xa6b6b553bd0d06e7, 0xb561ecec14591590, 0x17b7a66d8c64d959, 0x77778cecbe1461e, 0x9405c9c0c41a57ce, 0x8f6b4847e8ca7d3d, 0xf625eb987b366937, 0x421b3590e345},
-		B: Fp751Element{0x566b893803e7d8d6, 0xe8c71a04d527e696, 0x5a1d8f87bf5eb51, 0x42ae08ae098724f, 0x4ee3d7c7af40ca2e, 0xd9f9ab9067bb10a7, 0xecd53d69edd6328c, 0xa581e9202dea107d, 0x8bcdfb6c8ecf9257, 0xe7cbbc2e5cbcf2af, 0x5f031a8701f0e53e, 0x18312d93e3cb}})
-
-	phi.GenerateCurve(&xP3)
-	resPhiXr = phi.EvaluatePoint(&xR)
-
-	if !expPhiXr.VartimeEq(&resPhiXr) {
-		t.Error("\nExpected\n", expPhiXr.ToAffine(), "\nfound\n", resPhiXr.ToAffine())
-	}
-}

p751toolbox/print_test.go

@@ -1,19 +0,0 @@
-package p751toolbox
-
-// Tools used for testing and debugging
-
-import (
-	"fmt"
-)
-
-func (primeElement PrimeFieldElement) String() string {
-   return fmt.Sprintf("%X", primeElement.A.toBigInt().String())
-}
-
-func (extElement ExtensionFieldElement) String() string {
-	return fmt.Sprintf("\nA: %X\nB: %X", extElement.A.toBigInt().String(), extElement.B.toBigInt().String())
-}
-
-func (point ProjectivePoint) String() string {
-	return fmt.Sprintf("X:\n%sZ:\n%s", point.X.String(), point.Z.String())
-}

sidh/api.go

@@ -2,14 +2,10 @@ package sidh
 
 import (
 	"errors"
-	. "github.com/cloudflare/p751sidh/p751toolbox"
+	. "github.com/cloudflare/p751sidh/internal/isogeny"
 	"io"
 )
 
-// I keep it bool in order to be able to apply logical NOT
-type KeyVariant uint
-type PrimeFieldId uint
-
 // Id's correspond to bitlength of the prime field characteristic
 // Currently FP_751 is the only one supported by this implementation
 const (
@@ -43,9 +39,9 @@ type key struct {
 // Defines operations on public key
 type PublicKey struct {
 	key
-	affine_xP   ExtensionFieldElement
-	affine_xQ   ExtensionFieldElement
-	affine_xQmP ExtensionFieldElement
+	affine_xP   Fp2Element
+	affine_xQ   Fp2Element
+	affine_xQmP Fp2Element
 }
 
 // Defines operations on private key
@@ -91,9 +87,11 @@ func (pub *PublicKey) Import(input []byte) error {
 	if len(input) != pub.Size() {
 		return errors.New("sidh: input to short")
 	}
-	pub.affine_xP.FromBytes(input[0:pub.params.SharedSecretSize])
-	pub.affine_xQ.FromBytes(input[pub.params.SharedSecretSize : 2*pub.params.SharedSecretSize])
-	pub.affine_xQmP.FromBytes(input[2*pub.params.SharedSecretSize : 3*pub.params.SharedSecretSize])
+	op := CurveOperations{Params: pub.params}
+	ssSz := pub.params.SharedSecretSize
+	op.Fp2FromBytes(&pub.affine_xP, input[0:ssSz])
+	op.Fp2FromBytes(&pub.affine_xQ, input[ssSz:2*ssSz])
+	op.Fp2FromBytes(&pub.affine_xQmP, input[2*ssSz:3*ssSz])
 	return nil
 }
 
@@ -101,9 +99,11 @@ func (pub *PublicKey) Import(input []byte) error {
 // returned byte string is filled with zeros.
 func (pub *PublicKey) Export() []byte {
 	output := make([]byte, pub.params.PublicKeySize)
-	pub.affine_xP.ToBytes(output[0:pub.params.SharedSecretSize])
-	pub.affine_xQ.ToBytes(output[pub.params.SharedSecretSize : 2*pub.params.SharedSecretSize])
-	pub.affine_xQmP.ToBytes(output[2*pub.params.SharedSecretSize : 3*pub.params.SharedSecretSize])
+	op := CurveOperations{Params: pub.params}
+	ssSz := pub.params.SharedSecretSize
+	op.Fp2ToBytes(output[0:ssSz], &pub.affine_xP)
+	op.Fp2ToBytes(output[ssSz:2*ssSz], &pub.affine_xQ)
+	op.Fp2ToBytes(output[2*ssSz:3*ssSz], &pub.affine_xQmP)
 	return output
 }
 

sidh/params.go

@@ -1,36 +1,9 @@
 package sidh
 
-import . "github.com/cloudflare/p751sidh/p751toolbox"
-
-type DomainParams struct {
-	// P, Q and R=P-Q base points
-	Affine_P, Affine_Q, Affine_R ExtensionFieldElement
-	// Max size of secret key for x-torsion group
-	SecretBitLen uint
-	// MaskBytes
-	MaskBytes []byte
-	// Size of a compuatation strategy for x-torsion group
-	IsogenyStrategy []uint32
-}
-
-type SidhParams struct {
-	Id PrimeFieldId
-	// The secret key size, in bytes.
-	SecretKeySize int
-	// The public key size, in bytes.
-	PublicKeySize int
-	// The shared secret size, in bytes.
-	SharedSecretSize uint
-	// 2- and 3-torsion group parameter definitions
-	A, B DomainParams
-	// Sample rate to obtain a value in [0,3^238]
-	SampleRate uint
-	// Length of SIKE secret message. Must be one of {24,32,40},
-	// depending on size of prime field used (see [SIKE], 1.4 and 5.1)
-	MsgLen uint
-	// Length of SIKE ephemeral KEM key (see [SIKE], 1.4 and 5.1)
-	KemSize uint
-}
+import (
+	p751 "github.com/cloudflare/p751sidh/p751"
+	. "github.com/cloudflare/p751sidh/internal/isogeny"
+)
 
 // Keeps mapping: SIDH prime field ID to domain parameters
 var sidhParams = make(map[PrimeFieldId]SidhParams)
@@ -47,29 +20,33 @@ func Params(id PrimeFieldId) *SidhParams {
 func init() {
 	p751 := SidhParams{
 		Id:               FP_751,
-		SecretKeySize:    P751_SecretKeySize,
-		PublicKeySize:    P751_PublicKeySize,
-		SharedSecretSize: P751_SharedSecretSize,
+		SecretKeySize:    p751.P751_SecretKeySize,
+		PublicKeySize:    p751.P751_PublicKeySize,
+		SharedSecretSize: p751.P751_SharedSecretSize,
 		A: DomainParams{
-			Affine_P:        P751_affine_PA,
-			Affine_Q:        P751_affine_QA,
-			Affine_R:        P751_affine_RA,
-			SecretBitLen:    P751_SecretBitLenA,
-			MaskBytes:       []byte{P751_MaskAliceByte1, P751_MaskAliceByte2, P751_MaskAliceByte3},
-			IsogenyStrategy: P751_AliceIsogenyStrategy[:],
+			Affine_P:        p751.P751_affine_PA,
+			Affine_Q:        p751.P751_affine_QA,
+			Affine_R:        p751.P751_affine_RA,
+			SecretBitLen:    p751.P751_SecretBitLenA,
+			MaskBytes:       []byte{p751.P751_MaskAliceByte1, p751.P751_MaskAliceByte2, p751.P751_MaskAliceByte3},
+			IsogenyStrategy: p751.P751_AliceIsogenyStrategy[:],
 		},
 		B: DomainParams{
-			Affine_P:        P751_affine_PB,
-			Affine_Q:        P751_affine_QB,
-			Affine_R:        P751_affine_RB,
-			SecretBitLen:    P751_SecretBitLenB,
-			MaskBytes:       []byte{P751_MaskBobByte},
-			IsogenyStrategy: P751_BobIsogenyStrategy[:],
+			Affine_P:        p751.P751_affine_PB,
+			Affine_Q:        p751.P751_affine_QB,
+			Affine_R:        p751.P751_affine_RB,
+			SecretBitLen:    p751.P751_SecretBitLenB,
+			MaskBytes:       []byte{p751.P751_MaskBobByte},
+			IsogenyStrategy: p751.P751_BobIsogenyStrategy[:],
 		},
+		OneFp2:  p751.P751_OneFp2,
+		HalfFp2: p751.P751_HalfFp2,
 		MsgLen: 32,
 		// SIKEp751 provides 192 bit of classical security ([SIKE], 5.1)
 		KemSize:    24,
-		SampleRate: P751_SampleRate,
+		SampleRate: p751.P751_SampleRate,
+		Bytelen: p751.P751_Bytelen,
+		Op: p751.FieldOperations(),
 	}
 
 	sidhParams[FP_751] = p751

sidh/sidh.go

@@ -4,9 +4,7 @@ import (
 	"errors"
 	"io"
 
-	// TODO: This is needed by ExtensionFieldElement struct, which itself
-	// 		 depends on implementation of p751.
-	. "github.com/cloudflare/p751sidh/p751toolbox"
+	. "github.com/cloudflare/p751sidh/internal/isogeny"
 )
 
 // -----------------------------------------------------------------------------
@@ -19,9 +17,10 @@ func traverseTreePublicKeyA(curve *ProjectiveCurveParameters, xR, phiP, phiQ, ph
 	var points = make([]ProjectivePoint, 0, 8)
 	var indices = make([]int, 0, 8)
 	var i, sidx int
+	var op = CurveOperations{Params: pub.params}
 
-	cparam := curve.CalcCurveParamsEquiv4()
-	phi := NewIsogeny4()
+	cparam := op.CalcCurveParamsEquiv4(curve)
+	phi := Newisogeny4(op.Params.Op)
 	strat := pub.params.A.IsogenyStrategy
 	stratSz := len(strat)
 
@@ -32,7 +31,7 @@ func traverseTreePublicKeyA(curve *ProjectiveCurveParameters, xR, phiP, phiQ, ph
 
 			k := strat[sidx]
 			sidx++
-			xR.Pow2k(&cparam, xR, 2*k)
+			op.Pow2k(xR, &cparam, 2*k)
 			i += int(k)
 		}
 
@@ -57,9 +56,10 @@ func traverseTreeSharedKeyA(curve *ProjectiveCurveParameters, xR *ProjectivePoin
 	var points = make([]ProjectivePoint, 0, 8)
 	var indices = make([]int, 0, 8)
 	var i, sidx int
+	var op = CurveOperations{Params: pub.params}
 
-	cparam := curve.CalcCurveParamsEquiv4()
-	phi := NewIsogeny4()
+	cparam := op.CalcCurveParamsEquiv4(curve)
+	phi := Newisogeny4(op.Params.Op)
 	strat := pub.params.A.IsogenyStrategy
 	stratSz := len(strat)
 
@@ -70,7 +70,7 @@ func traverseTreeSharedKeyA(curve *ProjectiveCurveParameters, xR *ProjectivePoin
 
 			k := strat[sidx]
 			sidx++
-			xR.Pow2k(&cparam, xR, 2*k)
+			op.Pow2k(xR, &cparam, 2*k)
 			i += int(k)
 		}
 
@@ -91,9 +91,10 @@ func traverseTreePublicKeyB(curve *ProjectiveCurveParameters, xR, phiP, phiQ, ph
 	var points = make([]ProjectivePoint, 0, 8)
 	var indices = make([]int, 0, 8)
 	var i, sidx int
+	var op = CurveOperations{Params: pub.params}
 
-	cparam := curve.CalcCurveParamsEquiv3()
-	phi := NewIsogeny3()
+	cparam := op.CalcCurveParamsEquiv3(curve)
+	phi := Newisogeny3(op.Params.Op)
 	strat := pub.params.B.IsogenyStrategy
 	stratSz := len(strat)
 
@@ -104,7 +105,7 @@ func traverseTreePublicKeyB(curve *ProjectiveCurveParameters, xR, phiP, phiQ, ph
 
 			k := strat[sidx]
 			sidx++
-			xR.Pow3k(&cparam, xR, k)
+			op.Pow3k(xR, &cparam, k)
 			i += int(k)
 		}
 
@@ -129,9 +130,10 @@ func traverseTreeSharedKeyB(curve *ProjectiveCurveParameters, xR *ProjectivePoin
 	var points = make([]ProjectivePoint, 0, 8)
 	var indices = make([]int, 0, 8)
 	var i, sidx int
+	var op = CurveOperations{Params: pub.params}
 
-	cparam := curve.CalcCurveParamsEquiv3()
-	phi := NewIsogeny3()
+	cparam := op.CalcCurveParamsEquiv3(curve)
+	phi := Newisogeny3(op.Params.Op)
 	strat := pub.params.B.IsogenyStrategy
 	stratSz := len(strat)
 
@@ -142,7 +144,7 @@ func traverseTreeSharedKeyB(curve *ProjectiveCurveParameters, xR *ProjectivePoin
 
 			k := strat[sidx]
 			sidx++
-			xR.Pow3k(&cparam, xR, k)
+			op.Pow3k(xR, &cparam, k)
 			i += int(k)
 		}
 
@@ -224,29 +226,30 @@ func (prv *PrivateKey) generatePrivateKeyB(rand io.Reader) error {
 func publicKeyGenA(prv *PrivateKey) (pub *PublicKey) {
 	var xPA, xQA, xRA ProjectivePoint
 	var xPB, xQB, xRB, xR ProjectivePoint
-	var invZP, invZQ, invZR ExtensionFieldElement
+	var invZP, invZQ, invZR Fp2Element
 	var tmp ProjectiveCurveParameters
-	var phi = NewIsogeny4()
+
 	pub = NewPublicKey(prv.params.Id, KeyVariant_SIDH_A)
+	var op = CurveOperations{Params: pub.params}
+	var phi = Newisogeny4(op.Params.Op)
 
 	// Load points for A
-	xPA.FromAffine(&prv.params.A.Affine_P)
-	xQA.FromAffine(&prv.params.A.Affine_Q)
-	xRA.FromAffine(&prv.params.A.Affine_R)
+	xPA = ProjectivePoint{X: prv.params.A.Affine_P, Z: prv.params.OneFp2}
+	xQA = ProjectivePoint{X: prv.params.A.Affine_Q, Z: prv.params.OneFp2}
+	xRA = ProjectivePoint{X: prv.params.A.Affine_R, Z: prv.params.OneFp2}
 
 	// Load points for B
-	xRB.FromAffine(&prv.params.B.Affine_R)
-	xQB.FromAffine(&prv.params.B.Affine_Q)
-	xPB.FromAffine(&prv.params.B.Affine_P)
+	xRB = ProjectivePoint{X: prv.params.B.Affine_R, Z: prv.params.OneFp2}
+	xQB = ProjectivePoint{X: prv.params.B.Affine_Q, Z: prv.params.OneFp2}
+	xPB = ProjectivePoint{X: prv.params.B.Affine_P, Z: prv.params.OneFp2}
 
 	// Find isogeny kernel
-	tmp.A.Zero()
-	tmp.C.One()
-	xR = RightToLeftLadder(&tmp, &xPA, &xQA, &xRA, prv.params.A.SecretBitLen, prv.Scalar)
+	tmp.C = pub.params.OneFp2
+	xR = op.ScalarMul3Pt(&tmp, &xPA, &xQA, &xRA, prv.params.A.SecretBitLen, prv.Scalar)
 
 	// Reset params object and travers isogeny tree
-	tmp.A.Zero()
-	tmp.C.One()
+	tmp.C = pub.params.OneFp2
+	tmp.A.Zeroize()
 	traverseTreePublicKeyA(&tmp, &xR, &xPB, &xQB, &xRB, pub)
 
 	// Secret isogeny
@@ -254,11 +257,11 @@ func publicKeyGenA(prv *PrivateKey) (pub *PublicKey) {
 	xPA = phi.EvaluatePoint(&xPB)
 	xQA = phi.EvaluatePoint(&xQB)
 	xRA = phi.EvaluatePoint(&xRB)
-	ExtensionFieldBatch3Inv(&xPA.Z, &xQA.Z, &xRA.Z, &invZP, &invZQ, &invZR)
+	op.Fp2Batch3Inv(&xPA.Z, &xQA.Z, &xRA.Z, &invZP, &invZQ, &invZR)
 
-	pub.affine_xP.Mul(&xPA.X, &invZP)
-	pub.affine_xQ.Mul(&xQA.X, &invZQ)
-	pub.affine_xQmP.Mul(&xRA.X, &invZR)
+	op.Params.Op.Mul(&pub.affine_xP, &xPA.X, &invZP)
+	op.Params.Op.Mul(&pub.affine_xQ, &xQA.X, &invZQ)
+	op.Params.Op.Mul(&pub.affine_xQmP, &xRA.X, &invZR)
 	return
 }
 
@@ -266,38 +269,39 @@ func publicKeyGenA(prv *PrivateKey) (pub *PublicKey) {
 func publicKeyGenB(prv *PrivateKey) (pub *PublicKey) {
 	var xPB, xQB, xRB, xR ProjectivePoint
 	var xPA, xQA, xRA ProjectivePoint
-	var invZP, invZQ, invZR ExtensionFieldElement
+	var invZP, invZQ, invZR Fp2Element
 	var tmp ProjectiveCurveParameters
-	var phi = NewIsogeny3()
+
 	pub = NewPublicKey(prv.params.Id, prv.keyVariant)
+	var op = CurveOperations{Params: pub.params}
+	var phi = Newisogeny3(op.Params.Op)
 
 	// Load points for B
-	xRB.FromAffine(&prv.params.B.Affine_R)
-	xQB.FromAffine(&prv.params.B.Affine_Q)
-	xPB.FromAffine(&prv.params.B.Affine_P)
+	xRB = ProjectivePoint{X: prv.params.B.Affine_R, Z: prv.params.OneFp2}
+	xQB = ProjectivePoint{X: prv.params.B.Affine_Q, Z: prv.params.OneFp2}
+	xPB = ProjectivePoint{X: prv.params.B.Affine_P, Z: prv.params.OneFp2}
 
 	// Load points for A
-	xPA.FromAffine(&prv.params.A.Affine_P)
-	xQA.FromAffine(&prv.params.A.Affine_Q)
-	xRA.FromAffine(&prv.params.A.Affine_R)
+	xPA = ProjectivePoint{X: prv.params.A.Affine_P, Z: prv.params.OneFp2}
+	xQA = ProjectivePoint{X: prv.params.A.Affine_Q, Z: prv.params.OneFp2}
+	xRA = ProjectivePoint{X: prv.params.A.Affine_R, Z: prv.params.OneFp2}
 
-	tmp.A.Zero()
-	tmp.C.One()
-	xR = RightToLeftLadder(&tmp, &xPB, &xQB, &xRB, prv.params.B.SecretBitLen, prv.Scalar)
+	tmp.C = pub.params.OneFp2
+	xR = op.ScalarMul3Pt(&tmp, &xPB, &xQB, &xRB, prv.params.B.SecretBitLen, prv.Scalar)
 
-	tmp.A.Zero()
-	tmp.C.One()
+	tmp.C = pub.params.OneFp2
+	tmp.A.Zeroize()
 	traverseTreePublicKeyB(&tmp, &xR, &xPA, &xQA, &xRA, pub)
 
 	phi.GenerateCurve(&xR)
 	xPB = phi.EvaluatePoint(&xPA)
 	xQB = phi.EvaluatePoint(&xQA)
 	xRB = phi.EvaluatePoint(&xRA)
-	ExtensionFieldBatch3Inv(&xPB.Z, &xQB.Z, &xRB.Z, &invZP, &invZQ, &invZR)
+	op.Fp2Batch3Inv(&xPB.Z, &xQB.Z, &xRB.Z, &invZP, &invZQ, &invZR)
 
-	pub.affine_xP.Mul(&xPB.X, &invZP)
-	pub.affine_xQ.Mul(&xQB.X, &invZQ)
-	pub.affine_xQmP.Mul(&xRB.X, &invZR)
+	op.Params.Op.Mul(&pub.affine_xP, &xPB.X, &invZP)
+	op.Params.Op.Mul(&pub.affine_xQ, &xQB.X, &invZQ)
+	op.Params.Op.Mul(&pub.affine_xQmP, &xRB.X, &invZR)
 	return
 }
 
@@ -311,25 +315,26 @@ func deriveSecretA(prv *PrivateKey, pub *PublicKey) []byte {
 	var cparam ProjectiveCurveParameters
 	var xP, xQ, xQmP ProjectivePoint
 	var xR ProjectivePoint
-	var phi = NewIsogeny4()
+	var op = CurveOperations{Params: prv.params}
+	var phi = Newisogeny4(op.Params.Op)
 
 	// Recover curve coefficients
-	cparam.RecoverCoordinateA(&pub.affine_xP, &pub.affine_xQ, &pub.affine_xQmP)
-	cparam.C.One()
+	cparam.C = pub.params.OneFp2
+	op.RecoverCoordinateA(&cparam, &pub.affine_xP, &pub.affine_xQ, &pub.affine_xQmP)
 
 	// Find kernel of the morphism
-	xP.FromAffine(&pub.affine_xP)
-	xQ.FromAffine(&pub.affine_xQ)
-	xQmP.FromAffine(&pub.affine_xQmP)
-	xR = RightToLeftLadder(&cparam, &xP, &xQ, &xQmP, pub.params.A.SecretBitLen, prv.Scalar)
+	xP = ProjectivePoint{X: pub.affine_xP, Z: pub.params.OneFp2}
+	xQ = ProjectivePoint{X: pub.affine_xQ, Z: pub.params.OneFp2}
+	xQmP = ProjectivePoint{X: pub.affine_xQmP, Z: pub.params.OneFp2}
+	xR = op.ScalarMul3Pt(&cparam, &xP, &xQ, &xQmP, pub.params.A.SecretBitLen, prv.Scalar)
 
 	// Traverse isogeny tree
 	traverseTreeSharedKeyA(&cparam, &xR, pub)
 
 	// Calculate j-invariant on isogeneus curve
 	c := phi.GenerateCurve(&xR)
-	cparam.RecoverCurveCoefficients4(&c)
-	cparam.Jinvariant(sharedSecret)
+	op.RecoverCurveCoefficients4(&cparam, &c)
+	op.Jinvariant(&cparam, sharedSecret)
 	return sharedSecret
 }
 
@@ -339,24 +344,25 @@ func deriveSecretB(prv *PrivateKey, pub *PublicKey) []byte {
 	var xP, xQ, xQmP ProjectivePoint
 	var xR ProjectivePoint
 	var cparam ProjectiveCurveParameters
-	var phi = NewIsogeny3()
+	var op = CurveOperations{Params: prv.params}
+	var phi = Newisogeny3(op.Params.Op)
 
 	// Recover curve coefficients
-	cparam.RecoverCoordinateA(&pub.affine_xP, &pub.affine_xQ, &pub.affine_xQmP)
-	cparam.C.One()
+	cparam.C = pub.params.OneFp2
+	op.RecoverCoordinateA(&cparam, &pub.affine_xP, &pub.affine_xQ, &pub.affine_xQmP)
 
 	// Find kernel of the morphism
-	xP.FromAffine(&pub.affine_xP)
-	xQ.FromAffine(&pub.affine_xQ)
-	xQmP.FromAffine(&pub.affine_xQmP)
-	xR = RightToLeftLadder(&cparam, &xP, &xQ, &xQmP, pub.params.B.SecretBitLen, prv.Scalar)
+	xP = ProjectivePoint{X: pub.affine_xP, Z: pub.params.OneFp2}
+	xQ = ProjectivePoint{X: pub.affine_xQ, Z: pub.params.OneFp2}
+	xQmP = ProjectivePoint{X: pub.affine_xQmP, Z: pub.params.OneFp2}
+	xR = op.ScalarMul3Pt(&cparam, &xP, &xQ, &xQmP, pub.params.B.SecretBitLen, prv.Scalar)
 
 	// Traverse isogeny tree
 	traverseTreeSharedKeyB(&cparam, &xR, pub)
 
 	// Calculate j-invariant on isogeneus curve
 	c := phi.GenerateCurve(&xR)
-	cparam.RecoverCurveCoefficients3(&c)
-	cparam.Jinvariant(sharedSecret)
+	op.RecoverCurveCoefficients3(&cparam, &c)
+	op.Jinvariant(&cparam, sharedSecret)
 	return sharedSecret
 }

sidh/sidh_test.go

@@ -5,6 +5,8 @@ import (
 	"crypto/rand"
 	"encoding/hex"
 	"testing"
+
+	. "github.com/cloudflare/p751sidh/internal/isogeny"
 )
 
 const (