WIP

The method is implemented in exactly same way as multiplication.
Currently there is no explenation for it's existance. Benchmarking shows
no change

.travis.yml

@@ -8,6 +8,7 @@ matrix:
   fast_finish: true
 
 script:
+  - make test
   - make bench
   - make cover
   - NOASM=1 make cover-sidh

Makefile

@@ -3,37 +3,49 @@ MK_FILE_PATH = $(lastword $(MAKEFILE_LIST))
 PRJ_DIR      = $(abspath $(dir $(MK_FILE_PATH)))
 GOPATH_LOCAL = $(PRJ_DIR)/build
 GOPATH_DIR   = github.com/cloudflare/p751sidh
-CSHAKE_PKG   = github.com/henrydcase/nobs/hash/sha3
+CSHAKE_PKG   ?= github.com/henrydcase/nobs/hash/sha3
 TARGETS      = p751toolbox sidh sike
 GOARCH       ?=
 OPTS_GCCGO   ?= -compiler gccgo -O2 -g
-OPTS_TAGS    ?= -tags=noasm
 OPTS         ?=
+OPTS_TAGS    ?= -tags=noasm
 NOASM        ?=
+# -run="NonExistent" is set to make sure tests are not run before benchmarking
+BENCH_OPTS   ?= -bench=. -run="NonExistent"
+# whether to be verbose
+V            ?= 1
 
 ifeq ($(NOASM),1)
 	OPTS+=$(OPTS_TAGS)
 endif
 
+ifeq ($(V),1)
+	OPTS += -v              # Be verbose
+	OPTS += -gcflags=-m     # Show results from inlining
+endif
+
+
 clean:
 	rm -rf $(GOPATH_LOCAL)
 	rm -rf coverage*.txt
 
-prep:
+build_env:
 	GOPATH=$(GOPATH_LOCAL) go get $(CSHAKE_PKG)
 	mkdir -p $(GOPATH_LOCAL)/src/$(GOPATH_DIR)
-	cp -rf p751toolbox $(GOPATH_LOCAL)/src/$(GOPATH_DIR)
-	cp -rf sidh $(GOPATH_LOCAL)/src/$(GOPATH_DIR)
-	cp -rf sike $(GOPATH_LOCAL)/src/$(GOPATH_DIR)
 	cp -rf etc $(GOPATH_LOCAL)/src/$(GOPATH_DIR)
 
-test-%: prep
-	GOPATH=$(GOPATH_LOCAL) go test -v $(OPTS) $(GOPATH_DIR)/$*
+copy-target-%:
+	cp -rf $* $(GOPATH_LOCAL)/src/$(GOPATH_DIR)
+
+prep_targets: build_env $(addprefix copy-target-, $(TARGETS))
+
+test-%: prep_targets
+	GOPATH=$(GOPATH_LOCAL) go test $(OPTS) $(GOPATH_DIR)/$*
 
-bench-%: prep
-	cd $*; GOPATH=$(GOPATH_LOCAL) go test -v $(OPTS) -bench=.
+bench-%: prep_targets
+	cd $*; GOPATH=$(GOPATH_LOCAL) go test $(OPTS) $(BENCH_OPTS)
 
-cover-%: prep
+cover-%: prep_targets
 	GOPATH=$(GOPATH_LOCAL) go test \
 		-race -coverprofile=coverage_$*.txt -covermode=atomic $(OPTS) $(GOPATH_DIR)/$*
 	cat coverage_$*.txt >> coverage.txt

p751toolbox/curve.go

@@ -28,15 +28,6 @@ type ProjectivePoint struct {
 	Z ExtensionFieldElement
 }
 
-// A point on the projective line P^1(F_p).
-//
-// This represents a point on the (Kummer line) of the prime-field subgroup of
-// the base curve E_0(F_p), defined by E_0 : y^2 = x^3 + x.
-type ProjectivePrimeFieldPoint struct {
-	X PrimeFieldElement
-	Z PrimeFieldElement
-}
-
 func (params *ProjectiveCurveParameters) FromAffine(a *ExtensionFieldElement) {
 	params.A = *a
 	params.C.One()
@@ -158,34 +149,17 @@ func (cparams *ProjectiveCurveParameters) RecoverCurveCoefficients4(coefEq *Curv
 	return
 }
 
-func (point *ProjectivePoint) FromAffinePrimeField(x *PrimeFieldElement) {
-	point.X.A = x.A
-	point.X.B = zeroExtensionField.B
-	point.Z = oneExtensionField
-}
-
 func (point *ProjectivePoint) FromAffine(x *ExtensionFieldElement) {
 	point.X = *x
 	point.Z = oneExtensionField
 }
 
-func (point *ProjectivePrimeFieldPoint) FromAffine(x *PrimeFieldElement) {
-	point.X = *x
-	point.Z = onePrimeField
-}
-
 func (point *ProjectivePoint) ToAffine() *ExtensionFieldElement {
 	affine_x := new(ExtensionFieldElement)
 	affine_x.Inv(&point.Z).Mul(affine_x, &point.X)
 	return affine_x
 }
 
-func (point *ProjectivePrimeFieldPoint) ToAffine() *PrimeFieldElement {
-	affine_x := new(PrimeFieldElement)
-	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)
@@ -193,23 +167,11 @@ func (lhs *ProjectivePoint) VartimeEq(rhs *ProjectivePoint) bool {
 	return t0.VartimeEq(&t1)
 }
 
-func (lhs *ProjectivePrimeFieldPoint) VartimeEq(rhs *ProjectivePrimeFieldPoint) bool {
-	var t0, t1 PrimeFieldElement
-	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)
 }
 
-func ProjectivePrimeFieldPointConditionalSwap(xP, xQ *ProjectivePrimeFieldPoint, choice uint8) {
-	PrimeFieldConditionalSwap(&xP.X, &xQ.X, choice)
-	PrimeFieldConditionalSwap(&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
@@ -268,6 +230,89 @@ func (x2P *ProjectivePoint) Pow2k(params *CurveCoefficientsEquiv, xP *Projective
 	return x2P
 }
 
+//type Mul_t func(res, lhs, rhs *ExtensionFieldElement)
+//type Square_t func (res, x *ExtensionFieldElement)
+//type Sub_t func(res, lhs, rhs *ExtensionFieldElement)
+//type Add_t func(res, lhs, rhs *ExtensionFieldElement)
+type Ops interface {
+	MulFp2(res, lhs, rhs *ExtensionFieldElement)
+	SquareFp2(res, x *ExtensionFieldElement)
+	SubFp2(res, lhs, rhs *ExtensionFieldElement)
+	AddFp2(res, lhs, rhs *ExtensionFieldElement)
+}
+
+type OpsCtx struct {
+	op [2]Ops;
+}
+var op OpsCtx
+type Ops1 struct{};
+
+func (ctx Ops1) MulFp2(res, lhs, rhs *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(&res.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(&res.A, &ac_minus_bd) // = (a*c - b*d)*R mod p
+}
+func (ctx Ops1) SquareFp2(res, x *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(&res.A, &asq_minus_bsq) // = (a^2 - b^2)*R mod p
+       fp751MontgomeryReduce(&res.B, &ab2)           // = 2*a*b*R mod p
+}
+func (ctx Ops1) AddFp2(res, lhs, rhs *ExtensionFieldElement) {
+       fp751AddReduced(&res.A, &lhs.A, &rhs.A)
+       fp751AddReduced(&res.B, &lhs.B, &rhs.B)
+}
+func (ctx Ops1) SubFp2(res, lhs, rhs *ExtensionFieldElement) {
+       fp751SubReduced(&res.A, &lhs.A, &rhs.A)
+       fp751SubReduced(&res.B, &lhs.B, &rhs.B)
+}
+
+
 // 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.
@@ -277,29 +322,31 @@ func (x3P *ProjectivePoint) Pow3k(params *CurveCoefficientsEquiv, xP *Projective
 	*x3P = *xP
 	x, z := &x3P.X, &x3P.Z
 
+	pp := op.op[0]
+
 	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
+        pp.SubFp2(&t0, x, z)           // t0  = Xp - Zp
+        pp.SquareFp2(&t2, &t0)         // t2  = t0^2
+        pp.AddFp2(&t1,x, z)           // t1  = Xp + Zp
+        pp.SquareFp2(&t3, &t1)         // t3  = t1^2
+        pp.AddFp2(&t4,&t1, &t0)       // t4  = t1 + t0
+        pp.SubFp2(&t0, &t1, &t0)       // t0  = t1 - t0
+        pp.SquareFp2(&t1, &t4)         // t1  = t4^2
+        pp.SubFp2(&t1, &t1, &t3)       // t1  = t1 - t3
+        pp.SubFp2(&t1, &t1, &t2)       // t1  = t1 - t2
+        pp.MulFp2(&t5,&t3, &params.A) // t5  = t3 * A24+
+        pp.MulFp2(&t3,&t3, &t5)       // t3  = t5 * t3
+        pp.MulFp2(&t6,&t2, &params.C) // t6  = t2 * A24-
+        pp.MulFp2(&t2,&t2, &t6)       // t2  = t2 * t6
+        pp.SubFp2(&t3, &t2, &t3)       // t3  = t2 - t3
+        pp.SubFp2(&t2, &t5, &t6)       // t2  = t5 - t6
+        pp.MulFp2(&t1,&t2, &t1)       // t1  = t2 * t1
+        pp.AddFp2(&t2,&t3, &t1)       // t2  = t3 + t1
+        pp.SquareFp2(&t2, &t2)         // t2  = t2^2
+        pp.MulFp2(x,&t2, &t4)        // X3p = t2 * t4
+        pp.SubFp2(&t1, &t3, &t1)       // t1  = t3 - t1
+        pp.SquareFp2(&t1, &t1)         // t1  = t1^2
+        pp.MulFp2(z,&t1, &t0)        // Z3p = t1 * t0
 	}
 	return x3P
 }
@@ -329,3 +376,7 @@ func RightToLeftLadder(c *ProjectiveCurveParameters, P, Q, PmQ *ProjectivePoint,
 	ProjectivePointConditionalSwap(&R1, &R2, prevBit)
 	return R1
 }
+
+func init() {
+	op.op[0] = new(Ops1)
+}

p751toolbox/field.go

@@ -105,16 +105,6 @@ func (dest *ExtensionFieldElement) Mul(lhs, rhs *ExtensionFieldElement) *Extensi
 	return dest
 }
 
-// Set dest = -x
-//
-// Allowed to overlap dest with x.
-//
-// Returns dest to allow chaining operations.
-func (dest *ExtensionFieldElement) Neg(x *ExtensionFieldElement) *ExtensionFieldElement {
-	dest.Sub(&zeroExtensionField, x)
-	return dest
-}
-
 // Set dest = 1/x
 //
 // Allowed to overlap dest with x.
@@ -142,18 +132,21 @@ func (dest *ExtensionFieldElement) Inv(x *ExtensionFieldElement) *ExtensionField
 	fp751MontgomeryReduce(&asq_plus_bsq.A, &asq) // = (a^2 + b^2)*R mod p
 	// Now asq_plus_bsq = a^2 + b^2
 
-	var asq_plus_bsq_inv PrimeFieldElement
-	asq_plus_bsq_inv.Inv(&asq_plus_bsq)
-	c := &asq_plus_bsq_inv.A
+	// 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, c)
+	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, c)
+	fp751Mul(&minus_bc, &minus_b, &inv.A)
 	fp751MontgomeryReduce(&dest.B, &minus_bc)
 
 	return dest
@@ -277,35 +270,6 @@ var onePrimeField = PrimeFieldElement{
 	A: Fp751Element{0x249ad, 0x0, 0x0, 0x0, 0x0, 0x8310000000000000, 0x5527b1e4375c6c66, 0x697797bf3f4f24d0, 0xc89db7b2ac5c4e2e, 0x4ca4b439d2076956, 0x10f7926c7512c7e9, 0x2d5b24bce5e2},
 }
 
-// Set dest = 0.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) Zero() *PrimeFieldElement {
-	*dest = zeroPrimeField
-	return dest
-}
-
-// Set dest = 1.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) One() *PrimeFieldElement {
-	*dest = onePrimeField
-	return dest
-}
-
-// Set dest to x.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) SetUint64(x uint64) *PrimeFieldElement {
-	var xRR fp751X2
-	dest.A = Fp751Element{}                 // = 0
-	dest.A[0] = x                           // = x
-	fp751Mul(&xRR, &dest.A, &montgomeryRsq) // = x*R*R
-	fp751MontgomeryReduce(&dest.A, &xRR)    // = x*R mod p
-
-	return dest
-}
-
 // Set dest = lhs * rhs.
 //
 // Allowed to overlap lhs or rhs with dest.
@@ -328,104 +292,14 @@ func (dest *PrimeFieldElement) Mul(lhs, rhs *PrimeFieldElement) *PrimeFieldEleme
 //
 // Returns dest to allow chaining operations.
 func (dest *PrimeFieldElement) Pow2k(x *PrimeFieldElement, k uint8) *PrimeFieldElement {
-	dest.Square(x)
+	dest.Mul(x, x)
 	for i := uint8(1); i < k; i++ {
-		dest.Square(dest)
+		dest.Mul(dest, dest)
 	}
 
 	return dest
 }
 
-// Set dest = x^2
-//
-// Allowed to overlap x with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) Square(x *PrimeFieldElement) *PrimeFieldElement {
-	a := &x.A // = a*R
-	b := &x.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
-//
-// Allowed to overlap x with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) Neg(x *PrimeFieldElement) *PrimeFieldElement {
-	dest.Sub(&zeroPrimeField, x)
-	return dest
-}
-
-// Set dest = lhs + rhs.
-//
-// Allowed to overlap lhs or rhs with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) Add(lhs, rhs *PrimeFieldElement) *PrimeFieldElement {
-	fp751AddReduced(&dest.A, &lhs.A, &rhs.A)
-
-	return dest
-}
-
-// Set dest = lhs - rhs.
-//
-// Allowed to overlap lhs or rhs with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) Sub(lhs, rhs *PrimeFieldElement) *PrimeFieldElement {
-	fp751SubReduced(&dest.A, &lhs.A, &rhs.A)
-
-	return dest
-}
-
-// Returns true if lhs = rhs.  Takes variable time.
-func (lhs *PrimeFieldElement) VartimeEq(rhs *PrimeFieldElement) bool {
-	return lhs.A.vartimeEq(rhs.A)
-}
-
-// If choice = 1u8, set (x,y) = (y,x). If choice = 0u8, set (x,y) = (x,y).
-//
-// Returns dest to allow chaining operations.
-func PrimeFieldConditionalSwap(x, y *PrimeFieldElement, choice uint8) {
-	fp751ConditionalSwap(&x.A, &y.A, choice)
-}
-
-// Set dest = sqrt(x), if x is a square.  If x is nonsquare dest is undefined.
-//
-// Allowed to overlap x with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) Sqrt(x *PrimeFieldElement) *PrimeFieldElement {
-	tmp_x := *x // Copy x in case dest == x
-	// Since x is assumed to be square, x = y^2
-	dest.P34(x)            // dest = (y^2)^((p-3)/4) = y^((p-3)/2)
-	dest.Mul(dest, &tmp_x) // dest = y^2 * y^((p-3)/2) = y^((p+1)/2)
-	// Now dest^2 = y^(p+1) = y^2 = x, so dest = sqrt(x)
-
-	return dest
-}
-
-// Set dest = 1/x.
-//
-// Allowed to overlap x with dest.
-//
-// Returns dest to allow chaining operations.
-func (dest *PrimeFieldElement) Inv(x *PrimeFieldElement) *PrimeFieldElement {
-	tmp_x := *x            // Copy x in case dest == x
-	dest.Square(x)         // dest = x^2
-	dest.P34(dest)         // dest = (x^2)^((p-3)/4) = x^((p-3)/2)
-	dest.Square(dest)      // dest = x^(p-3)
-	dest.Mul(dest, &tmp_x) // dest = x^(p-2)
-
-	return dest
-}
-
 // Set dest = x^((p-3)/4).  If x is square, this is 1/sqrt(x).
 //
 // Allowed to overlap x with dest.
@@ -448,7 +322,7 @@ func (dest *PrimeFieldElement) P34(x *PrimeFieldElement) *PrimeFieldElement {
 	// Build a lookup table of odd multiples of x.
 	lookup := [16]PrimeFieldElement{}
 	xx := &PrimeFieldElement{}
-	xx.Square(x) // Set xx = x^2
+	xx.Mul(x, x) // Set xx = x^2
 	lookup[0] = *x
 	for i := 1; i < 16; i++ {
 		lookup[i].Mul(&lookup[i-1], xx)

p751toolbox/field_test.go

@@ -35,6 +35,10 @@ func radix64ToBigInt(x []uint64) *big.Int {
 	return val
 }
 
+func VartimeEq(x,y *PrimeFieldElement) bool {
+	return x.A.vartimeEq(y.A)
+}
+
 func (x *PrimeFieldElement) toBigInt() *big.Int {
 	// Convert from Montgomery form
 	return x.A.toBigIntFromMontgomeryForm()
@@ -249,91 +253,6 @@ func TestExtensionFieldElementBatch3Inv(t *testing.T) {
 //------------------------------------------------------------------------------
 // Prime Field
 //------------------------------------------------------------------------------
-
-func TestPrimeFieldElementSetUint64VersusBigInt(t *testing.T) {
-	setUint64RoundTrips := func(x uint64) bool {
-		z := new(PrimeFieldElement).SetUint64(x).toBigInt().Uint64()
-		return x == z
-	}
-
-	if err := quick.Check(setUint64RoundTrips, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-func TestPrimeFieldElementAddVersusBigInt(t *testing.T) {
-	addMatchesBigInt := func(x, y PrimeFieldElement) bool {
-		z := new(PrimeFieldElement)
-		z.Add(&x, &y)
-
-		check := new(big.Int)
-		check.Add(x.toBigInt(), y.toBigInt())
-		check.Mod(check, cln16prime)
-
-		return check.Cmp(z.toBigInt()) == 0
-	}
-
-	if err := quick.Check(addMatchesBigInt, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-func TestPrimeFieldElementSubVersusBigInt(t *testing.T) {
-	subMatchesBigInt := func(x, y PrimeFieldElement) bool {
-		z := new(PrimeFieldElement)
-		z.Sub(&x, &y)
-
-		check := new(big.Int)
-		check.Sub(x.toBigInt(), y.toBigInt())
-		check.Mod(check, cln16prime)
-
-		return check.Cmp(z.toBigInt()) == 0
-	}
-
-	if err := quick.Check(subMatchesBigInt, quickCheckConfig); err != nil {
-		t.Error(err)
-	}
-}
-
-func TestPrimeFieldElementInv(t *testing.T) {
-	inverseIsCorrect := func(x PrimeFieldElement) bool {
-		z := new(PrimeFieldElement)
-		z.Inv(&x)
-
-		// Now z = (1/x), so (z * x) * x == x
-		z.Mul(z, &x).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 TestPrimeFieldElementSqrt(t *testing.T) {
-	inverseIsCorrect := func(x PrimeFieldElement) bool {
-		// Construct y = x^2 so we're sure y is square.
-		y := new(PrimeFieldElement)
-		y.Square(&x)
-
-		z := new(PrimeFieldElement)
-		z.Sqrt(y)
-
-		// Now z = sqrt(y), so z^2 == y
-		z.Square(z)
-		return z.VartimeEq(y)
-	}
-
-	// 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 TestPrimeFieldElementMulVersusBigInt(t *testing.T) {
 	mulMatchesBigInt := func(x, y PrimeFieldElement) bool {
 		z := new(PrimeFieldElement)
@@ -370,26 +289,6 @@ func TestPrimeFieldElementP34VersusBigInt(t *testing.T) {
 	}
 }
 
-func TestFp751ElementConditionalSwap(t *testing.T) {
-	var one = Fp751Element{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
-	var two = Fp751Element{2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2}
-
-	var x = one
-	var y = two
-
-	fp751ConditionalSwap(&x, &y, 0)
-
-	if !(x == one && y == two) {
-		t.Error("Found", x, "expected", one)
-	}
-
-	fp751ConditionalSwap(&x, &y, 1)
-
-	if !(x == two && y == one) {
-		t.Error("Found", x, "expected", two)
-	}
-}
-
 // Package-level storage for this field element is intended to deter
 // compiler optimizations.
 var benchmarkFp751Element Fp751Element
@@ -452,51 +351,6 @@ func BenchmarkPrimeFieldElementMul(b *testing.B) {
 	}
 }
 
-func BenchmarkPrimeFieldElementInv(b *testing.B) {
-	z := &PrimeFieldElement{A: bench_x}
-	w := new(PrimeFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Inv(z)
-	}
-}
-
-func BenchmarkPrimeFieldElementSqrt(b *testing.B) {
-	z := &PrimeFieldElement{A: bench_x}
-	w := new(PrimeFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Sqrt(z)
-	}
-}
-
-func BenchmarkPrimeFieldElementSquare(b *testing.B) {
-	z := &PrimeFieldElement{A: bench_x}
-	w := new(PrimeFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Square(z)
-	}
-}
-
-func BenchmarkPrimeFieldElementAdd(b *testing.B) {
-	z := &PrimeFieldElement{A: bench_x}
-	w := new(PrimeFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Add(z, z)
-	}
-}
-
-func BenchmarkPrimeFieldElementSub(b *testing.B) {
-	z := &PrimeFieldElement{A: bench_x}
-	w := new(PrimeFieldElement)
-
-	for n := 0; n < b.N; n++ {
-		w.Sub(z, z)
-	}
-}
-
 // --- field operation functions
 
 func BenchmarkFp751Multiply(b *testing.B) {

p751toolbox/isogeny.go

@@ -79,16 +79,18 @@ func (phi *isogeny3) EvaluatePoint(p *ProjectivePoint) 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
+	pp := op.op[0]
+
+	pp.AddFp2(&t0, px, pz)   // t0 = XQ + ZQ
+	pp.SubFp2(&t1,px, pz)   // t1 = XQ - ZQ
+	pp.MulFp2(&t0,K1, &t0)  // t2 = K1 * t0
+	pp.MulFp2(&t1,K2, &t1)  // t1 = K2 * t1
+	pp.AddFp2(&t2, &t0, &t1) // t2 = t0 + t1
+	pp.SubFp2(&t0,&t1, &t0) // t0 = t1 - t0
+	pp.SquareFp2(&t2, &t2)   // t2 = t2 ^ 2
+	pp.SquareFp2(&t0,&t0)   // t0 = t0 ^ 2
+	pp.MulFp2(&q.X,px, &t2) // XQ'= XQ * t2
+	pp.MulFp2(&q.Z, pz, &t0) // ZQ'= ZQ * t0
 	return q
 }
 
@@ -126,19 +128,21 @@ func (phi *isogeny4) EvaluatePoint(p *ProjectivePoint) ProjectivePoint {
 	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)
+	pp := op.op[0]
+
+	pp.AddFp2(&t0, xq, zq)
+	pp.SubFp2(&t1,xq, zq)
+	pp.MulFp2(xq, &t0, K2)
+	pp.MulFp2(zq, &t1, K3)
+	pp.MulFp2(&t0, &t0, &t1)
+	pp.MulFp2(&t0, &t0, K1)
+	pp.AddFp2(&t1,xq, zq)
+	pp.SubFp2(zq,xq, zq)
+	pp.SquareFp2(&t1,&t1)
+	pp.SquareFp2(zq,zq)
+	pp.AddFp2(xq,&t0, &t1)
+	pp.SubFp2(&t0,zq, &t0)
+	pp.MulFp2(xq, xq, &t1)
+	pp.MulFp2(zq, zq, &t0)
 	return q
 }

p751toolbox/print_test.go

@@ -38,10 +38,6 @@ func (point ProjectivePoint) String() string {
 	return fmt.Sprintf("X:\n%sZ:\n%s", point.X.String(), point.Z.String())
 }
 
-func (point ProjectivePrimeFieldPoint) String() string {
-	return fmt.Sprintf("X:\n%sZ:\n%s", point.X.String(), point.Z.String())
-}
-
 func (point Fp751Element) PrintHex() {
 	fmt.Printf("\t")
 	for i := 0; i < fp751NumWords/2; i++ {