Groth16 prover

constantine.nimble

@@ -614,6 +614,7 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   # Polynomials
   # ----------------------------------------------------------
   ("tests/math_polynomials/t_polynomials.nim", false),
+  ("tests/math_polynomials/t_finite_field_fft.nim", false),
 
   # Protocols
   # ----------------------------------------------------------
@@ -633,6 +634,9 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   # Proof systems
   # ----------------------------------------------------------
   ("tests/proof_systems/t_r1cs_parser.nim", false),
+  ("tests/proof_systems/t_wtns_parser.nim", false),
+  ("tests/proof_systems/t_zkey_parser.nim", false),
+  ("tests/proof_systems/t_groth16_prover.nim", false),
   ("tests/interactive_proofs/t_multilinear_extensions.nim", false),
 ]
 

constantine/math/elliptic/ec_shortweierstrass_jacobian.nim

@@ -337,6 +337,8 @@ template sumImpl[F; G: static Subgroup](
 
   # if P or R were infinity points they would have spread 0 with Z₁Z₂
   block: # Infinity points
+    bind isNeutral
+    bind ccopy
     o.ccopy(Q, P.isNeutral())
     o.ccopy(P, Q.isNeutral())
 

constantine/math/io/io_ec.nim

@@ -56,6 +56,31 @@ func toHex*[EC: EC_ShortW_Prj or EC_ShortW_Jac or EC_ShortW_Aff or EC_ShortW_Jac
   result.appendHex(aff.y)
   result &= "\n" & sp & ")"
 
+func toDecimal*[EC: EC_ShortW_Prj or EC_ShortW_Jac or EC_ShortW_Aff or EC_ShortW_JacExt](P: EC, indent: static int = 0): string =
+  ## Stringify an elliptic curve point to Hex
+  ## Note. Leading zeros are not removed.
+  ## Output as decimal.
+  ##
+  ## WARNING: NOT constant time!
+  ##
+  ## This proc output may change format in the future
+
+  var aff {.noInit.}: EC_ShortW_Aff[EC.F, EC.G]
+  when EC isnot EC_ShortW_Aff:
+    aff.affine(P)
+  else:
+    aff = P
+
+  const sp = spaces(indent)
+
+  result = sp & $EC & "(\n" & sp & "  x: "
+  result.add toDecimal(aff.x)
+  result &= ",\n" & sp & "  y: "
+  result.add toDecimal(aff.y)
+  result &= "\n" & sp & ")"
+
+
+
 func toHex*[EC: EC_TwEdw_Aff or EC_TwEdw_Prj](P: EC, indent: static int = 0): string =
   ## Stringify an elliptic curve point to Hex for Twisted Edwards Curve
   ## Note, leading zeros are not removed.

constantine/math/io/io_extfields.nim

@@ -52,6 +52,40 @@ func toHex*(f: ExtensionField, indent = 0, order: static Endianness = bigEndian)
   ##   - no leaks
   result.appendHex(f, indent, order)
 
+func appendDecimal*(accum: var string, f: Fp, indent = 0, order: static Endianness = bigEndian) =
+  accum.add toDecimal(f)
+
+func appendDecimal*(accum: var string, f: ExtensionField, indent = 0, order: static Endianness = bigEndian) =
+  ## Stringify a tower field element to hex.
+  ## Note. Leading zeros are not removed.
+  ## Result is prefixed with 0x
+  ##
+  ## Output will be padded with 0s to maintain constant-time.
+  ##
+  ## CT:
+  ##   - no leaks
+  accum.add static($f.typeof.genericHead() & '(')
+  staticFor i, 0, f.coords.len:
+    when i != 0:
+      accum.add ", "
+    accum.add "\n" & spaces(indent+2) & "c" & $i & ": "
+    when f is Fp2:
+      accum.appendDecimal(f.coords[i], order = order)
+    else:
+      accum.appendDecimal(f.coords[i], indent+2, order)
+  accum.add ")"
+
+func toDecimal*(f: ExtensionField, indent = 0, order: static Endianness = bigEndian): string =
+  ## Stringify a tower field element to hex.
+  ## Note. Leading zeros are not removed.
+  ## Result is prefixed with 0x
+  ##
+  ## Output will be padded with 0s to maintain constant-time.
+  ##
+  ## CT:
+  ##   - no leaks
+  result.appendDecimal(f, indent, order)
+
 func fromHex*(dst: var Fp2, c0, c1: string) =
   ## Convert 2 coordinates to an element of 𝔽p2
   ## with dst = c0 + β * c1

constantine/math/io/io_fields.nim

@@ -116,7 +116,9 @@ func toDecimal*(f: FF): string =
   ## This function is NOT constant-time at the moment.
   f.toBig().toDecimal()
 
-func fromDecimal*(dst: var FF, decimalString: string) =
+{.pop.} # `fromDecimal` can raise `ValueError`
+
+func fromDecimal*(dst: var FF, decimalString: string) {.raises: [ValueError].} =
   ## Convert a decimal string. The input must be packed
   ## with no spaces or underscores.
   ## This assumes that bits and decimal length are **public.**
@@ -136,7 +138,7 @@ func fromDecimal*(dst: var FF, decimalString: string) =
   let raw {.noinit.} = fromDecimal(dst.mres.typeof, decimalString)
   dst.fromBig(raw)
 
-func fromDecimal*(T: type FF, hexString: string): T {.noInit.}=
+func fromDecimal*(T: type FF, hexString: string): T {.raises: [ValueError], noInit.}=
   ## Convert a decimal string. The input must be packed
   ## with no spaces or underscores.
   ## This assumes that bits and decimal length are **public.**

constantine/math/polynomials/fft.nim

@@ -98,7 +98,7 @@ func fft_internal[EC; bits: static int](
 
   for i in 0 ..< half:
     # FFT Butterfly
-    y_times_root   .scalarMul_vartime(output[i+half], rootsOfUnity[i])
+    y_times_root   .scalarMul_vartime(rootsOfUnity[i], output[i+half])
     output[i+half] .diff_vartime(output[i], y_times_root)
     output[i]      .sum_vartime(output[i], y_times_root)
 

constantine/math/polynomials/fft_fields.nim

@@ -0,0 +1,225 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import
+  constantine/named/algebras,
+  constantine/math/arithmetic,
+  constantine/math/io/io_bigints,
+  constantine/platforms/[abstractions, allocs, views],
+  ./fft_lut,
+  constantine/platforms/bithacks # for nextPowerOf2
+
+# ############################################################
+#
+#               Fast Fourier Transform
+#
+# ############################################################
+
+# Fast Fourier Transform (Number Theoretic Transform - NTT) over finite fields
+# ----------------------------------------------------------------
+
+type
+  FFTStatus* = enum
+    FFTS_Success
+    FFTS_TooManyValues = "Input length greater than the field 2-adicity (number of roots of unity)"
+    FFTS_SizeNotPowerOfTwo = "Input must be of a power of 2 length"
+
+  FFT_Descriptor*[F] = object # `F` is either `Fp[Name]` or `Fr[Name]`
+    ## Metadata for FFT on Elliptic Curve
+    order*: int
+    rouGen*: F ## Roots of unity generator based on primitive root of `F`: `ω = g^( (p - 1) // n )`
+    rootsOfUnity*: ptr UncheckedArray[getBigInt(F)] # `getBigInt` gives us the right type depending on Fr/Fp
+      ## domain, starting and ending with 1, length is cardinality+1
+      ## This allows FFT and inverse FFT to use the same buffer for roots.
+
+func computeRootsOfUnity[F](ctx: var FFT_Descriptor[F], generatorRootOfUnity: auto) =
+  static:
+    doAssert typeof(generatorRootOfUnity) is Fr[F.Name] or typeof(generatorRootOfUnity) is Fp[F.Name]
+
+  ctx.rootsOfUnity[0].setOne()
+
+  var cur = generatorRootOfUnity
+  for i in 1 .. ctx.order:
+    ctx.rootsOfUnity[i].fromField(cur)
+    cur *= generatorRootOfUnity
+
+  doAssert ctx.rootsOfUnity[ctx.order].isOne().bool(), "The given generator does not seem to be a root of unity " &
+    "of " & $F & " for order: " & $ctx.order & "."
+
+func init*[Name: static Algebra](T: type FFT_Descriptor, order: int, generatorRootOfUnity: FF[Name]): T =
+  result.order = order
+  result.rouGen = generatorRootOfUnity
+  result.rootsOfUnity = allocHeapArrayAligned(T.F.getBigInt(), order+1, alignment = 64)
+
+  result.computeRootsOfUnity(generatorRootOfUnity)
+
+proc rootOfUnityGenerator*[F](_: typedesc[F], order: int): F =
+  ## Computes a root of unity generator for the order `n`, using
+  ## `ω = g^( (p - 1) // n )` where `g` is the primitive root
+  ## of the field `F`.
+  ##
+  ## `p` = prime of the field (or order of subgroup)
+  ## `n` = FFT order
+  ## Highlighted the part we compute in each comment.
+  # ω = g^( `(p - 1)` // n )
+  var exponentI {.noInit.}: BigInt[F.bits()]
+  exponentI = F.getModulus()
+  exponentI -= One
+  var exponent = F.fromBig(exponentI)
+
+  # ω = g^( (p - 1) // `n` )
+  var n = F.fromInt(order.uint64)
+  # ω = g^( (p - 1) `// n` )
+  n.inv()
+  # ω = g^( `(p - 1) // n` )
+  exponent *= n
+
+  var g: F = F.fromUint(primitiveRoot(F.Name).uint64)
+  # ω = `g^( (p - 1) // n )`
+  g.pow_vartime(toBig(exponent))
+  result = g
+
+proc init*(T: typedesc[FFT_Descriptor], order: int): T =
+  ## Initialize an `FFT_Descriptor` for the given `order`. The root of unity generator is
+  ## computed automatically. However, a primitive root is required for the field over which
+  ## the FFT is to be done. See `fft_lut.nim` for definitions and more information.
+  let g = rootOfUnityGenerator(T.F, order)
+  result = T.init(order, g)
+
+func delete*(ctx: FFT_Descriptor) =
+  ctx.rootsOfUnity.freeHeapAligned()
+
+proc toFr[S: static int, Name: static Algebra](x: BigInt[S]): Fr[Name] =
+  result.fromBig(x)
+
+proc toFp[S: static int, Name: static Algebra](x: BigInt[S]): Fp[Name] =
+  result.fromBig(x)
+
+proc toF[F; S: static int](T: typedesc[F], x: BigInt[S]): auto =
+  when T is Fr:
+    toFr[S, T.Name](x)
+  else:
+    toFp[S, T.Name](x)
+
+func simpleFT[F; bits: static int](
+       output: var StridedView[F],
+       vals: StridedView[F],
+       rootsOfUnity: StridedView[BigInt[bits]]) =
+  # FFT is a recursive algorithm
+  # This is the base-case using a O(n²) algorithm
+
+  let L = output.len
+  var last {.noInit.}, v {.noInit.}: F
+
+  var v0w0 {.noinit}: F
+  v0w0.prod(vals[0], F.toF(rootsOfUnity[0]))
+
+  for i in 0 ..< L:
+    last = v0w0
+    for j in 1 ..< L:
+      v.prod(F.toF(rootsOfUnity[(i*j) mod L]), vals[j])
+      last.sum(last, v)
+    output[i] = last
+
+func fft_internal[F; bits: static int](
+       output: var StridedView[F],
+       vals: StridedView[F],
+       rootsOfUnity: StridedView[BigInt[bits]]) =
+  if output.len <= 4:
+    simpleFT(output, vals, rootsOfUnity)
+    return
+
+  let (evenVals, oddVals) = vals.splitAlternate()
+  var (outLeft, outRight) = output.splitHalf()
+  let halfROI = rootsOfUnity.skipHalf()
+
+  fft_internal(outLeft, evenVals, halfROI)
+  fft_internal(outRight, oddVals, halfROI)
+
+  let half = outLeft.len
+  var y_times_root{.noinit.}: F
+
+  for i in 0 ..< half:
+    # FFT Butterfly
+    y_times_root.prod(F.toF(rootsOfUnity[i]), output[i+half])
+    output[i+half].diff(output[i], y_times_root)
+    output[i].sum(output[i], y_times_root)
+
+
+func fft_vartime*[F](
+       desc: FFT_Descriptor[F],
+       output: var openarray[F],
+       vals: openarray[F]): FFT_Status =
+  if vals.len > desc.order:
+    return FFTS_TooManyValues
+  if not vals.len.uint64.isPowerOf2_vartime():
+    return FFTS_SizeNotPowerOfTwo
+
+  let rootz = desc.rootsOfUnity
+                  .toStridedView(desc.order)
+                  .slice(0, desc.order-1, desc.order div vals.len)
+
+  var voutput = output.toStridedView()
+  fft_internal(voutput, vals.toStridedView(), rootz)
+  return FFTS_Success
+
+# Similar adjustments would be made for ifft_vartime
+
+func ifft_vartime*[F](
+       desc: FFT_Descriptor[F],
+       output: var openarray[F],
+       vals: openarray[F]): FFT_Status =
+  ## Inverse FFT
+  if vals.len > desc.order:
+    return FFTS_TooManyValues
+  if not vals.len.uint64.isPowerOf2_vartime():
+    return FFTS_SizeNotPowerOfTwo
+
+  let rootz = desc.rootsOfUnity
+                  .toStridedView(desc.order+1) # Extra 1 at the end so that when reversed the buffer starts with 1
+                  .reversed()
+                  .slice(0, desc.order-1, desc.order div vals.len)
+
+  var voutput = output.toStridedView()
+  fft_internal(voutput, vals.toStridedView(), rootz)
+
+  var invLen {.noInit.}: F.getBigInt()
+  invLen.fromUint(vals.len.uint64)
+  invLen.invmod_vartime(invLen, F.getModulus())
+
+  for i in 0 ..< output.len:
+    let inp = output[i]
+    output[i].prod(inp, F.toF(invLen))
+
+  return FFTS_Success
+
+func fft_vartime*[F](vals: openarray[F]): seq[F] =
+  ## Performs an FFT on the given values and returns a seq of the result.
+  ##
+  ## For convenience only!
+  let order = nextPowerOfTwo_vartime(vals.len.uint64)
+  var fftDesc = FFTDescriptor[F].init(order.int)
+  defer: fftDesc.delete()
+
+  result = newSeq[F](order)
+  let status = fftDesc.fft_vartime(result, vals)
+
+  doAssert (status == FFTS_Success).bool, "FFT failed with " & $status
+
+proc ifft_vartime*[F](vals: openarray[F]): seq[F] =
+  ## Performs an inverse FFT on the given values and returns a seq of the result.
+  ##
+  ## For convenience only!
+  let order = nextPowerOfTwo_vartime(vals.len.uint64)
+  var fftDesc = FFTDescriptor[F].init(order.int)
+  defer: fftDesc.delete()
+
+  result = newSeq[F](order)
+  let status = fftDesc.ifft_vartime(result, vals)
+
+  doAssert (status == FFTS_Success).bool, "FFT failed with " & $status