constantine G2 scalar multiplication bug

bug.nim

import strutils

import constantine/math/arithmetic    
import constantine/math/io/io_fields  
import constantine/math/io/io_bigints
import constantine/math/config/curves  
import constantine/math/extension_fields/towers 
import constantine/math/elliptic/ec_shortweierstrass_affine     as aff 
import constantine/math/elliptic/ec_shortweierstrass_projective as prj 
import constantine/math/elliptic/ec_scalar_mul                

#-------------------------------------------------------------------------------

type FFr  = Fr[BN254Snarks]
type FFp  = Fp[BN254Snarks]
type FFp2 = QuadraticExt[FFp]

type GG1 = ECP_ShortW_Aff[FFp , aff.G1]
type GG2 = ECP_ShortW_Aff[FFp2, aff.G2]

type ProjG1  = ECP_ShortW_Prj[FFp , prj.G1]
type ProjG2  = ECP_ShortW_Prj[FFp2, prj.G2]

const oneFr : FFr = fromHex( FFr, "0x01" )
const twoFr : FFr = fromHex( FFr, "0x02" )

#-------------------------------------------------------------------------------

func isZeroFp (x: FFp ): bool = bool(isZero(x))
func isZeroFp2(x: FFp2): bool = bool(isZero(x))
func isZeroFr (x: FFr ): bool = bool(isZero(x))

func isEqualFp (x, y: FFp ): bool = bool(x == y)
func isEqualFp2(x, y: FFp2): bool = bool(x == y)
func isEqualFr (x, y: FFr ): bool = bool(x == y)

func `+`*(x, y: FFp ): FFp  =  ( var z : FFp  = x ; z += y ; return z )
func `+`*(x, y: FFp2): FFp2 =  ( var z : FFp2 = x ; z += y ; return z )
func `+`*(x, y: FFr ): FFr  =  ( var z : FFr  = x ; z += y ; return z )

func `-`*(x, y: FFp ): FFp  =  ( var z : FFp  = x ; z -= y ; return z )
func `-`*(x, y: FFp2): FFp2 =  ( var z : FFp2 = x ; z -= y ; return z )
func `-`*(x, y: FFr ): FFr  =  ( var z : FFr  = x ; z -= y ; return z )

func `*`*(x, y: FFp ): FFp  =  ( var z : FFp  = x ; z *= y ; return z )
func `*`*(x, y: FFp2): FFp2 =  ( var z : FFp2 = x ; z *= y ; return z )
func `*`*(x, y: FFr ): FFr  =  ( var z : FFr  = x ; z *= y ; return z )

#-------------------------------------------------------------------------------

func mkFp2 (i: FFp, u: FFp) : FFp2 =
  let c : array[2, FFp] = [i,u]
  return QuadraticExt[FFp]( coords: c )

func squareFp2(y: FFp2): FFp2 =  ( var z : FFp2 = y ; square(z) ; return z )

func unsafeMkG2 ( X, Y: FFp2 ) : GG2 =
  return aff.ECP_ShortW_Aff[FFp2, aff.G2](x: X, y: Y)

const zeroFp  : FFp = fromHex( FFp, "0x00" )
const zeroFr  : FFr = fromHex( FFr, "0x00" )
const zeroFp2 : FFp2 = mkFp2( zeroFp, zeroFp )
const infG2   : GG2  = unsafeMkG2( zeroFp2 , zeroFp2 )

#-------------------------------------------------------------------------------

# y^2 = x^3 + B
# B = b1 + bu*u
# b1 = 19485874751759354771024239261021720505790618469301721065564631296452457478373
# b2 = 266929791119991161246907387137283842545076965332900288569378510910307636690
const twistCoeffB_1 : FFp  = fromHex(FFp, "0x2b149d40ceb8aaae81be18991be06ac3b5b4c5e559dbefa33267e6dc24a138e5")
const twistCoeffB_u : FFp  = fromHex(FFp, "0x009713b03af0fed4cd2cafadeed8fdf4a74fa084e52d1852e4a2bd0685c315d2")
const twistCoeffB   : FFp2 = mkFp2( twistCoeffB_1 , twistCoeffB_u )

func checkCurveEqG2( x, y: FFp2 ) : bool =
  if isZeroFp2(x) and isZeroFp2(y):
    # the point at infinity is on the curve by definition
    return true
  else:
    var x2 : FFp2 = squareFp2(x)
    var y2 : FFp2 = squareFp2(y)
    var x3 : FFp2 = x2 * x;
    var eq : FFp2
    eq =  x3
    eq += twistCoeffB
    eq -= y2
    return isZeroFp2(eq)

func isOnCurveG2 ( p: GG2 ) : bool =
  return checkCurveEqG2( p.x, p.y )

#-------------------------------------------------------------------------------

func mkG2( x, y: FFp2 ) : GG2 =
  if isZeroFp2(x) and isZeroFp2(y):
    return infG2
  else:
    assert( checkCurveEqG2(x,y) , "mkG2: not a G2 curve point" )
    return unsafeMkG2(x,y)

func negG2*(p: GG2): GG2 =
  var r : GG2 = p
  neg(r)
  return r

func addG2(p,q: GG2): GG2 =
  var r, x, y : ProjG2
  prj.fromAffine(x, p)
  prj.fromAffine(y, q)
  prj.sum(r, x, y)
  var s : GG2
  prj.affine(s, r)
  return s

func `+`(p,q: GG2): GG2 = addG2(p,q)

func `-`(p,q: GG2): GG2 = addG2(p,negG2(q))

func `**`( coeff: FFr , point: GG2 ) : GG2 =
  var q : ProjG2
  prj.fromAffine( q , point )
  scalarMul(  q , coeff.toBig() )
  var r : GG2
  prj.affine( r, q )
  return r

#-------------------------------------------------------------------------------

func toDecimalFp*(a : FFp): string =
  var s : string = toDecimal(a)
  s = s.strip( leading=true, trailing=false, chars={'0'} )
  if s.len == 0: s="0"
  return s

func toDecimalFr*(a : FFr): string =
  var s : string = toDecimal(a)
  s = s.strip( leading=true, trailing=false, chars={'0'} )
  if s.len == 0: s="0"
  return s

proc debugPrintFp*(prefix: string, x: FFp) =
  echo(prefix & toDecimalFp(x))

proc debugPrintFp2*(prefix: string, z: FFp2) =
  echo(prefix & " 1 ~> " & toDecimalFp(z.coords[0]))
  echo(prefix & " u ~> " & toDecimalFp(z.coords[1]))

proc debugPrintFr*(prefix: string, x: FFr) =
  echo(prefix & toDecimalFr(x))

proc debugPrintG1*(msg: string, pt: GG1) =
  echo(msg & ":")
  debugPrintFp( " x = ", pt.x )
  debugPrintFp( " y = ", pt.y )

proc debugPrintG2*(msg: string, pt: GG2) =
  echo(msg & ":")
  debugPrintFp2( " x = ", pt.x )
  debugPrintFp2( " y = ", pt.y )

#-------------------------------------------------------------------------------

let coeff  : FFr = fromHex( FFr, "0x7b17fcc286b01af79176aa7da3a8615020eacda89a90e4ff5d0a085483f0448" )
#let coeff_good  : FFr = fromHex( FFr, "0x7b17fcc286b01af79176aa7da3a8615020eacda89a90e4ff5d0a085483f044" )     # most coefficients work

let coeffA : FFr = fromHex( FFr, "0x12345678901234567890" )
let coeffB : FFr = coeff - coeffA

let gen2_xi : FFp = fromHex(FFp, "0x1adcd0ed10df9cb87040f46655e3808f98aa68a570acf5b0bde23fab1f149701")
let gen2_xu : FFp = fromHex(FFp, "0x09e847e9f05a6082c3cd2a1d0a3a82e6fbfbe620f7f31269fa15d21c1c13b23b")
let gen2_yi : FFp = fromHex(FFp, "0x056c01168a5319461f7ca7aa19d4fcfd1c7cdf52dbfc4cbee6f915250b7f6fc8")
let gen2_yu : FFp = fromHex(FFp, "0x0efe500a2d02dd77f5f401329f30895df553b878fc3c0dadaaa86456a623235c")

let gen2_x : FFp2 = mkFp2(gen2_xi, gen2_xu)
let gen2_y : FFp2 = mkFp2(gen2_yi, gen2_yu)

let gen2 : GG2 = mkG2(gen2_x, gen2_y)

#-------------------------------------------------------------------------------

debugPrintFr("coeff", coeff)
debugPrintG2("gen2 ", gen2 )

let lhs = coeff ** gen2                              # WRONG
let rhs = (coeffA ** gen2) + (coeffB ** gen2)        # CORRECT
let xxx = ((coeff - oneFr)**gen2) + gen2             # WRONG
let yyy = ((coeff - twoFr)**gen2) + (twoFr**gen2 )   # WRONG

debugPrintG2("coeff*gen2", lhs)

debugPrintG2("A*gen2 + (coeff-A)*gen2 ", rhs )

debugPrintG2("(coeff-1)*gen2 + gen2  ", xxx)
debugPrintG2("(coeff-2)*gen2 + 2*gen2", yyy)


echo("isOnCurve(lhs) = ",isOnCurveG2(lhs))
echo("isOnCurve(rhs) = ",isOnCurveG2(rhs))

debugPrintG2("lhs - rhs = ", lhs - rhs )

echo("")
echo("and the result according to SageMath:")
echo(" x1 = 17216390949661727229956939928583223226083668728437958793715435751523027888005 ")
echo(" xu =  3082945034329785101034278215941854680789766318859358488904629243958221738137 ")
echo(" y1 = 20108673238932196920264801702661201943173809015346082727725783869161803474440 ")
echo(" yu = 10405477402946058176045590740070709500904395284580129777629727895349459816649 ")