BN254-Snarks: bad performance on Fp4 squaring

[mratsim]

For some reason (multiplication by its super large 9+i non-residue?)
BN254-Snarks is very slow on Fp4

Suspect 1:

• Multiplication by Non-Residue here:

  constantine/constantine/tower_field_extensions/tower_instantiation.nim

  Lines 118 to 129 in c4a2dee

  	else:
  	# BN254_Snarks, u = 9
  	# Full 𝔽p2 multiplication is cheaper than addition chains
  	# for u*c0 and u*c1
  	static:
  	doAssert u >= 0 and uint64(u) <= uint64(high(BaseType))
  	doAssert v >= 0 and uint64(v) <= uint64(high(BaseType))
  	# TODO: compile-time
  	var NR {.noInit.}: Fp2
  	NR.c0.fromUint(uint u)
  	NR.c1.fromUint(uint v)
  	r.prod(a, NR)

• Rewritten from:

  constantine/constantine/tower_field_extensions/tower_instantiation.nim

  Lines 88 to 94 in 2c5e12d

  	let a0 = a.c0
  	let a1 = a.c1
  	when a.fromComplexExtension():
  	a.c0.diff(u * a0, v * a1)
  	else:
  	a.c0.sum(u * a0, (Beta * v) * a1)
  	a.c1.sum(v * a0, u * a1)

Suspect 2:

The general squaring in quadratic field:

constantine/constantine/tower_field_extensions/extension_fields.nim

Lines 477 to 524 in c4a2dee

	func square_generic(r: var QuadraticExt, a: QuadraticExt) =
	## Return a² in ``r``
	## ``r`` is initialized/overwritten
	# Algorithm (with β the non-residue in the base field)
	#
	# (c0, c1)² => (c0 + c1 w)²
	# => c0² + 2 c0 c1 w + c1²w²
	# => c0² + β c1² + 2 c0 c1 w
	# => (c0² + β c1², 2 c0 c1)
	# We have 2 squarings and 1 multiplication in the base field
	# which are significantly more costly than additions.
	# For example when construction 𝔽p12 from 𝔽p6:
	# - 4 limbs like BN254: multiplication is 20x slower than addition/substraction
	# - 6 limbs like BLS12-381: multiplication is 28x slower than addition/substraction
	#
	# We can save operations with one of the following expressions
	# of c0² + β c1² and noticing that c0c1 is already computed for the "y" coordinate
	#
	# Alternative 1:
	# c0² + β c1² <=> (c0 - c1)(c0 - β c1) + β c0c1 + c0c1
	#
	# Alternative 2:
	# c0² + β c1² <=> (c0 + c1)(c0 + β c1) - β c0c1 - c0c1
	mixin prod
	var v0 {.noInit.}, v1 {.noInit.}: typeof(r.c0)
	# v1 <- (c0 + β c1)
	v1.prod(a.c1, NonResidue)
	v1 += a.c0
	# v0 <- (c0 + c1)(c0 + β c1)
	v0.sum(a.c0, a.c1)
	v0 *= v1
	# v1 <- c0 c1
	v1.prod(a.c0, a.c1)
	# aliasing: a unneeded now
	# r0 = (c0 + c1)(c0 + β c1) - c0c1
	v0 -= v1
	# r1 = 2 c0c1
	r.c1.double(v1)
	# r0 = (c0 + c1)(c0 + β c1) - c0c1 - β c0c1
	v1 *= NonResidue
	r.c0.diff(v0, v1)

The basic expression is (c0² + β c1², 2 c0 c1)

We can either minimize Mul/Squarings, requiring only 2 multiplications by rewriting to:

• r0 = (c0 + c1)(c0 + β c1) - c0c1 - β c0c1
• r1 = 2 c0c1

or only use squarings by rewritting to

• r0 = a0² + β a1²
• r1 = (a0 + a1)² - a0² - a1²

The first expression requires 2 multiplication by non-residue.