go-sike/p503/arith.go

package sike

// Set dest = x^((p-3)/4).  If x is square, this is 1/sqrt(x).
// Uses variation of sliding-window algorithm from with window size
// of 5 and least to most significant bit sliding (left-to-right)
// See HAC 14.85 for general description.
//
// Allowed to overlap x with dest.
// All values in Montgomery domains
// Set dest = x^(2^k), for k >= 1, by repeated squarings.
func p34(dest, x *Fp) {
	var lookup [16]Fp

	// Sliding-window strategy computed with etc/scripts/sliding_window_strat_calc.py
	//
	// This performs sum(powStrategy) + 1 squarings and len(lookup) + len(mulStrategy)
	// multiplications.
	powStrategy := []uint8{1, 12, 5, 5, 2, 7, 11, 3, 8, 4, 11, 4, 7, 5, 6, 3, 7, 5, 7, 2, 12, 5, 6, 4, 6, 8, 6, 4, 7, 5, 5, 8, 5, 8, 5, 5, 8, 9, 3, 6, 2, 10, 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, 3}
	mulStrategy := []uint8{0, 12, 11, 10, 0, 1, 8, 3, 7, 1, 8, 3, 6, 7, 14, 2, 14, 14, 9, 0, 13, 9, 15, 5, 12, 7, 13, 7, 15, 6, 7, 9, 0, 5, 7, 6, 8, 8, 3, 7, 0, 10, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 3}

	// Precompute lookup table of odd multiples of x for window
	// size k=5.
	var xx Fp
	fpMulRdc(&xx, x, x)
	lookup[0] = *x
	for i := 1; i < 16; i++ {
		fpMulRdc(&lookup[i], &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[mulStrategy[0]]
	for i := uint8(1); i < uint8(len(powStrategy)); i++ {
		fpMulRdc(dest, dest, dest)
		for j := uint8(1); j < powStrategy[i]; j++ {
			fpMulRdc(dest, dest, dest)
		}
		fpMulRdc(dest, dest, &lookup[mulStrategy[i]])
	}
}

func add(dest, lhs, rhs *Fp2) {
	fpAddRdc(&dest.A, &lhs.A, &rhs.A)
	fpAddRdc(&dest.B, &lhs.B, &rhs.B)
}

func sub(dest, lhs, rhs *Fp2) {
	fpSubRdc(&dest.A, &lhs.A, &rhs.A)
	fpSubRdc(&dest.B, &lhs.B, &rhs.B)
}

func mul(dest, lhs, rhs *Fp2) {
	// 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 FpX2
	fpMul(&ac, a, c) // = a*c*R*R
	fpMul(&bd, b, d) // = b*d*R*R

	var b_minus_a, c_minus_d Fp
	fpSubRdc(&b_minus_a, b, a) // = (b-a)*R
	fpSubRdc(&c_minus_d, c, d) // = (c-d)*R

	var ad_plus_bc FpX2
	fpMul(&ad_plus_bc, &b_minus_a, &c_minus_d) // = (b-a)*(c-d)*R*R
	fp2Add(&ad_plus_bc, &ad_plus_bc, &ac)      // = ((b-a)*(c-d) + a*c)*R*R
	fp2Add(&ad_plus_bc, &ad_plus_bc, &bd)      // = ((b-a)*(c-d) + a*c + b*d)*R*R

	fpMontRdc(&dest.B, &ad_plus_bc) // = (a*d + b*c)*R mod p

	var ac_minus_bd FpX2
	fp2Sub(&ac_minus_bd, &ac, &bd)   // = (a*c - b*d)*R*R
	fpMontRdc(&dest.A, &ac_minus_bd) // = (a*c - b*d)*R mod p
}

func inv(dest, x *Fp2) {
	var e1, e2 FpX2
	var f1, f2 Fp

	fpMul(&e1, &x.A, &x.A) // = a*a*R*R
	fpMul(&e2, &x.B, &x.B) // = b*b*R*R
	fp2Add(&e1, &e1, &e2)  // = (a^2 + b^2)*R*R
	fpMontRdc(&f1, &e1)    // = (a^2 + b^2)*R mod p
	// Now f1 = a^2 + b^2

	fpMulRdc(&f2, &f1, &f1)
	p34(&f2, &f2)
	fpMulRdc(&f2, &f2, &f2)
	fpMulRdc(&f2, &f2, &f1)

	fpMul(&e1, &x.A, &f2)
	fpMontRdc(&dest.A, &e1)

	fpSubRdc(&f1, &Fp{}, &x.B)
	fpMul(&e1, &f1, &f2)
	fpMontRdc(&dest.B, &e1)
}

func sqr(dest, x *Fp2) {
	var a2, aPlusB, aMinusB Fp
	var a2MinB2, ab2 FpX2

	a := &x.A
	b := &x.B

	// (a + bi)*(a + bi) = (a^2 - b^2) + 2abi.
	fpAddRdc(&a2, a, a)                // = a*R + a*R = 2*a*R
	fpAddRdc(&aPlusB, a, b)            // = a*R + b*R = (a+b)*R
	fpSubRdc(&aMinusB, a, b)           // = a*R - b*R = (a-b)*R
	fpMul(&a2MinB2, &aPlusB, &aMinusB) // = (a+b)*(a-b)*R*R = (a^2 - b^2)*R*R
	fpMul(&ab2, &a2, b)                // = 2*a*b*R*R
	fpMontRdc(&dest.A, &a2MinB2)       // = (a^2 - b^2)*R mod p
	fpMontRdc(&dest.B, &ab2)           // = 2*a*b*R mod p
}

// In case choice == 1, performs following swap in constant time:
// 	xPx <-> xQx
//	xPz <-> xQz
// Otherwise returns xPx, xPz, xQx, xQz unchanged
func condSwap(xPx, xPz, xQx, xQz *Fp2, choice uint8) {
	fpSwapCond(&xPx.A, &xQx.A, choice)
	fpSwapCond(&xPx.B, &xQx.B, choice)
	fpSwapCond(&xPz.A, &xQz.A, choice)
	fpSwapCond(&xPz.B, &xQz.B, choice)
}