Ethereum KZG / EIP-4844 / Proto-danksharding followup (#270)

* Pass all verify_kzg_proof test cases

* pass blob_to_commitment tests

* move tests

* KZG: WIP on compute_proof

* eip4844: Pass all compute_kzg_proof tests

* pass compute_blob_kzg_proof tests

* pass all verify_blob_kzg_proof tests

* CI needs yaml

* fix memory leaks and add effect tags

* CI: lock yaml version too pre Nim 2.0

.github/workflows/ci.yml

@@ -211,21 +211,21 @@ jobs:
         run: |
           pacman -S --needed --noconfirm mingw-w64-x86_64-gmp
           nimble refresh --verbose -y
-          nimble install --verbose -y gmp jsony asynctools
+          nimble install --verbose -y gmp jsony asynctools yaml@1.1.0
 
       - name: Install test dependencies
         if: runner.os != 'Windows'
         shell: bash
         run: |
           nimble refresh --verbose -y
-          nimble install --verbose -y gmp jsony asynctools
+          nimble install --verbose -y gmp jsony asynctools yaml@1.1.0
 
       - name: Print Nim & compiler versions
         shell: bash
         # gcc is an alias to Clang on MacOS
         run: |
           nim -v
-          gcc -v 
+          gcc -v
 
       - name: Run Constantine tests (UNIX with Assembly)
         if: runner.os != 'Windows' && matrix.target.BACKEND == 'ASM'

constantine.nimble

@@ -496,6 +496,7 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   ("tests/t_ethereum_evm_precompiles.nim", false),
   ("tests/t_ethereum_bls_signatures.nim", false),
   ("tests/t_ethereum_eip2333_bls12381_key_derivation.nim", false),
+  ("tests/t_ethereum_eip4844_deneb_kzg.nim", false),
 ]
 
 const testDescNvidia: seq[string] = @[

constantine/commitments/kzg_polynomial_commitments.nim

@@ -104,7 +104,7 @@ import
 ## We have 2 parties, a Prover and a Verifier.
 ##
 ## They share a public Structured Reference String (SRS), also called trusted setup:
-##   srs_g1: [[1]₁, [τ]₁, [τ²]₁, ... [τⁿ]₁] also called powers of tau, with a bounded degree n
+##   srs_g1: [[1]₁, [τ]₁, [τ²]₁, ... [τⁿ⁻¹]₁] also called powers of tau, with a bounded degree n-1
 ##   srs_g2: [[1]₂, [τ]₂]
 ##
 ## τ and its powers are secrets that no one know, we only work with [τⁱ]₁ and [τ]₂
@@ -120,7 +120,7 @@ import
 ##
 ## 0. A data blob is interpreted as up to n 𝔽r elements
 ##    corresponding to a polynomial p(x) = blob₀ + blob₁ x + blob₂ x² + ... + blobₙ₋₁ xⁿ⁻¹
-##                                  p(x) = ∑ blobᵢ xⁱ
+##                                  p(x) = ∑₀ⁿ⁻¹ blobᵢ xⁱ
 ##
 ##    So we can commit/prove up to 4096*log₂(r) bits of data
 ##    For Ethereum, n = 4096 and log₂(r) = 255 bits
@@ -167,42 +167,7 @@ import
 ##   [(ω⁰, p(ω⁰)), (ω¹, p(ω¹)), (ω², p(ω²)), ..., (ωⁿ⁻¹, p(ωⁿ⁻¹))]
 ##   with ω ∈ 𝔽r a root of unity of order n, i.e. ωⁿ = 1
 
-type
+type G1aff[C: static Curve] = ECP_ShortW_Aff[Fp[C], G1]
-  PowersOfTauCoef[D: static int, F; G: static Subgroup] = object
-    coefs: array[D, ECP_ShortW_Aff[F, G]]
-
-  PowersOfTauEval[D: static int, F; G: static Subgroup] = object
-    evals: array[D, ECP_ShortW_Aff[F, G]]
-
-  G1aff[C: static Curve] = ECP_ShortW_Aff[Fp[C], G1]
-  G1jac[C: static Curve] = ECP_ShortW_Jac[Fp[C], G1]
-
-# Helper functions
-# ------------------------------------------------------------
-
-func g1_lincomb[C: static Curve](r: var G1jac[C],
-                points: ptr UncheckedArray[G1aff[C]],
-                scalars: ptr UncheckedArray[matchingOrderBigInt(C)],
-                len: int) =
-  ## Multi-scalar-multiplication / linear combination
-  r.raw.multiScalarMul_vartime(
-    scalars,
-    cast[ptr UncheckedArray[typeof points[0].raw]](points),
-    len)
-
-func g1_lincomb[C: static Curve](r: var G1jac[C],
-                points: ptr UncheckedArray[G1aff[C]],
-                scalars: ptr UncheckedArray[Fr[C]],
-                len: int) =
-  ## Multi-scalar-multiplication / linear combination
-  let scalars2 = allocHeapArray(matchingOrderBigInt(C), len)
-
-  for i in 0 ..< len:
-    scalars2[i].fromField(scalars[i])
-
-  r.g1_lincomb(points, scalars2, len)
-
-  scalars2.freeHeap()
 
 # KZG - Prover - Lagrange basis
 # ------------------------------------------------------------
@@ -211,54 +176,74 @@ func g1_lincomb[C: static Curve](r: var G1jac[C],
 # as the powers of τ
 
 func kzg_commit*[N: static int, C: static Curve](
-       commitment: var ECP_ShortW_Jac[Fp[C], G1],
+       commitment: var ECP_ShortW_Aff[Fp[C], G1],
-       poly_evals: array[N, matchingOrderBigInt(C)],
+       poly_evals: array[N, BigInt],
-       powers_of_tau: PowersOfTauEval[N, Fp[C], G1]) =
+       powers_of_tau: PolynomialEval[N, G1aff[C]]) {.tags:[Alloca, HeapAlloc, Vartime].} =
-  commitment.g1_lincomb(powers_of_tau.evals.asUnchecked(), poly_evals.asUnchecked(), N)
+
+  var commitmentJac {.noInit.}: ECP_ShortW_Jac[Fp[C], G1]
+  commitmentJac.multiScalarMul_vartime(poly_evals, powers_of_tau.evals)
+  commitment.affine(commitmentJac)
 
 func kzg_prove*[N: static int, C: static Curve](
-       proof: var ECP_ShortW_Jac[Fp[C], G1],
+       proof: var ECP_ShortW_Aff[Fp[C], G1],
        eval_at_challenge: var Fr[C],
        poly: PolynomialEval[N, Fr[C]],
        domain: PolyDomainEval[N, Fr[C]],
        challenge: Fr[C],
-       powers_of_tau: PowersOfTauEval[N, Fp[C], G1]) =
+       powers_of_tau: PolynomialEval[N, G1aff[C]],
+       isBitReversedDomain: static bool) {.tags:[Alloca, HeapAlloc, Vartime].} =
 
   # Note:
   #   The order of inputs in
-  #  `kzg_prove`, `evalPolyAt_vartime`, `differenceQuotientEvalOffDomain`, `differenceQuotientEvalInDomain`
+  #  `kzg_prove`, `evalPolyAt`, `differenceQuotientEvalOffDomain`, `differenceQuotientEvalInDomain`
   #  minimizes register changes when parameter passing.
+  #
+  # z = challenge in the following code
 
-  # z = challenge
+  let diffQuotientPolyFr = allocHeapAligned(PolynomialEval[N, Fr[C]], alignment = 64)
+  let invRootsMinusZ = allocHeapAligned(array[N, Fr[C]], alignment = 64)
 
-  let invRootsMinusZ = allocHeap(array[N, Fr[C]])
+  # Compute 1/(ωⁱ - z) with ω a root of unity, i in [0, N).
-  let diffQuotientPoly = allocHeap(PolynomialEval[N, Fr[C]])
+  # zIndex = i if ωⁱ - z == 0 (it is the i-th root of unity) and -1 otherwise.
-
+  let zIndex = invRootsMinusZ[].inverseRootsMinusZ_vartime(
-  let zIndex = invRootsMinusZ.inverseRootsMinusZ_vartime(domain, challenge)
+                                  domain, challenge,
+                                  earlyReturnOnZero = false)
 
   if zIndex == -1:
     # p(z)
-    eval_at_challenge.evalPolyAt_vartime(
+    eval_at_challenge.evalPolyAt(
-      invRootsMinusZ,
+      poly, challenge,
-      poly, domain,
+      invRootsMinusZ[],
-      challenge)
+      domain)
 
     # q(x) = (p(x) - p(z)) / (x - z)
-    diffQuotientPoly.differenceQuotientEvalOffDomain(
+    diffQuotientPolyFr[].differenceQuotientEvalOffDomain(
-      invRootsMinusZ, poly, eval_at_challenge)
+      poly, eval_at_challenge, invRootsMinusZ[])
   else:
     # p(z)
     # But the challenge z is equal to one of the roots of unity (how likely is that?)
-    eval_at_challenge = poly[zIndex]
+    eval_at_challenge = poly.evals[zIndex]
 
     # q(x) = (p(x) - p(z)) / (x - z)
-    diffQuotientPoly.differenceQuotientEvalInDomain(
+    diffQuotientPolyFr[].differenceQuotientEvalInDomain(
-      invRootsMinusZ, poly, domain, zIndex)
+      poly, uint32 zIndex, invRootsMinusZ[], domain, isBitReversedDomain)
 
-  proof.g1_lincomb(powers_of_tau.evals.asUnchecked(), diffQuotientPoly.asUnchecked(), N)
+  freeHeapAligned(invRootsMinusZ)
+
+  const orderBits = C.getCurveOrderBitwidth()
+  let diffQuotientPolyBigInt = allocHeapAligned(array[N, BigInt[orderBits]], alignment = 64)
+
+  for i in 0 ..< N:
+    diffQuotientPolyBigInt[i].fromField(diffQuotientPolyFr.evals[i])
+
+  freeHeapAligned(diffQuotientPolyFr)
+
+  var proofJac {.noInit.}: ECP_ShortW_Jac[Fp[C], G1]
+  proofJac.multiScalarMul_vartime(diffQuotientPolyBigInt[], powers_of_tau.evals)
+  proof.affine(proofJac)
+
+  freeHeapAligned(diffQuotientPolyBigInt)
 
-  freeHeap(diffQuotientPoly)
-  freeHeap(invRootsMinusZ)
 
 # KZG - Verifier
 # ------------------------------------------------------------
@@ -268,7 +253,7 @@ func kzg_verify*[F2; C: static Curve](
        challenge: BigInt, # matchingOrderBigInt(C),
        eval_at_challenge: BigInt, # matchingOrderBigInt(C),
        proof: ECP_ShortW_Aff[Fp[C], G1],
-       tauG2: ECP_ShortW_Aff[F2, G2]): bool =
+       tauG2: ECP_ShortW_Aff[F2, G2]): bool {.tags:[Alloca, Vartime].} =
   ## Verify a short KZG proof that ``p(challenge) = eval_at_challenge``
   ## without doing the whole p(challenge) computation
   #
@@ -306,10 +291,10 @@ func kzg_verify*[F2; C: static Curve](
   tauG2Jac.fromAffine(tauG2)
   commitmentJac.fromAffine(commitment)
 
-  tau_minus_challenge_G2.scalarMul(challenge)
+  tau_minus_challenge_G2.scalarMul_vartime(challenge)
   tau_minus_challenge_G2.diff(tauG2Jac, tau_minus_challenge_G2)
 
-  commitment_minus_eval_at_challenge_G1.scalarMul(eval_at_challenge)
+  commitment_minus_eval_at_challenge_G1.scalarMul_vartime(eval_at_challenge)
   commitment_minus_eval_at_challenge_G1.diff(commitmentJac, commitment_minus_eval_at_challenge_G1)
 
   var tmzG2 {.noInit.}: ECP_ShortW_Aff[F2, G2]

constantine/ethereum_eip4844_kzg_polynomial_commitments.nim

@@ -7,17 +7,17 @@
 # at your option. This file may not be copied, modified, or distributed except according to those terms.
 
 import
-  math/config/curves,
+  ./math/config/curves,
-  math/io/io_bigints,
+  ./math/io/[io_bigints, io_fields],
-  math/[ec_shortweierstrass, arithmetic, extension_fields],
+  ./math/[ec_shortweierstrass, arithmetic, extension_fields],
-  math/arithmetic/limbs_montgomery,
+  ./math/arithmetic/limbs_montgomery,
-  math/elliptic/ec_multi_scalar_mul,
+  ./math/elliptic/ec_multi_scalar_mul,
-  math/polynomials/polynomials,
+  ./math/polynomials/polynomials,
-  commitments/kzg_polynomial_commitments,
+  ./commitments/kzg_polynomial_commitments,
-  hashes,
+  ./hashes,
-  platforms/[abstractions, views, allocs],
+  ./platforms/[abstractions, views, allocs],
-  serialization/[codecs_bls12_381, endians],
+  ./serialization/[codecs_bls12_381, endians],
-  trusted_setups/ethereum_kzg_srs
+  ./trusted_setups/ethereum_kzg_srs
 
 export loadTrustedSetup, TrustedSetupStatus, EthereumKZGContext
 
@@ -110,9 +110,9 @@ func fiatShamirChallenge(dst: var Fr[BLS12_381], blob: Blob, commitmentBytes: ar
 
   transcript.update(FIAT_SHAMIR_PROTOCOL_DOMAIN)
 
-  # Append the degree of polynomial as a domain separator
+  # Append the degree of polynomial as 16-byte big-endian integer as a domain separator
-  transcript.update(FIELD_ELEMENTS_PER_BLOB.uint64.toBytes(bigEndian))
   transcript.update(default(array[16-sizeof(uint64), byte]))
+  transcript.update(FIELD_ELEMENTS_PER_BLOB.uint64.toBytes(bigEndian))
 
   transcript.update(blob)
   transcript.update(commitmentBytes)
@@ -183,7 +183,7 @@ func blob_to_bigint_polynomial(
     doAssert sizeof(dst[]) == sizeof(Blob)
     doAssert sizeof(array[FIELD_ELEMENTS_PER_BLOB, array[32, byte]]) == sizeof(Blob)
 
-  let view = cast[ptr array[FIELD_ELEMENTS_PER_BLOB, array[32, byte]]](blob.unsafeAddr())
+  let view = cast[ptr array[FIELD_ELEMENTS_PER_BLOB, array[32, byte]]](blob)
 
   for i in 0 ..< FIELD_ELEMENTS_PER_BLOB:
     let status = dst.evals[i].bytes_to_bls_bigint(view[i])
@@ -201,7 +201,7 @@ func blob_to_field_polynomial(
     doAssert sizeof(dst[]) == sizeof(Blob)
     doAssert sizeof(array[FIELD_ELEMENTS_PER_BLOB, array[32, byte]]) == sizeof(Blob)
 
-  let view = cast[ptr array[FIELD_ELEMENTS_PER_BLOB, array[32, byte]]](blob.unsafeAddr())
+  let view = cast[ptr array[FIELD_ELEMENTS_PER_BLOB, array[32, byte]]](blob)
 
   for i in 0 ..< FIELD_ELEMENTS_PER_BLOB:
     let status = dst.evals[i].bytes_to_bls_field(view[i])
@@ -214,6 +214,8 @@ func blob_to_field_polynomial(
 # ------------------------------------------------------------
 
 template check(evalExpr: CttCodecScalarStatus): untyped =
+  # Translate codec status code to KZG status code
+  # Beware of resource cleanup like heap allocation, this can early exit the caller.
   block:
     let status = evalExpr # Ensure single evaluation
     case status
@@ -222,6 +224,8 @@ template check(evalExpr: CttCodecScalarStatus): untyped =
     of cttCodecScalar_ScalarLargerThanCurveOrder:       return cttEthKZG_ScalarLargerThanCurveOrder
 
 template check(evalExpr: CttCodecEccStatus): untyped =
+  # Translate codec status code to KZG status code
+  # Beware of resource cleanup like heap allocation, this can early exit the caller.
   block:
     let status = evalExpr # Ensure single evaluation
     case status
@@ -236,21 +240,90 @@ func blob_to_kzg_commitment*(
        ctx: ptr EthereumKZGContext,
        dst: var array[48, byte],
        blob: ptr Blob): CttEthKzgStatus =
+  ## Compute a commitment to the `blob`.
+  ## The commitment can be verified without needing the full `blob`
+  ##
+  ## Mathematical description
+  ##   commitment = [p(τ)]₁
+  ##
+  ##   The blob data is used as a polynomial,
+  ##   the polynomial is evaluated at powers of tau τ, a trusted setup.
+  ##
+  ##   Verification can be done by verifying the relation:
+  ##     proof.(τ - z) = p(τ)-p(z)
+  ##   which doesn't require the full blob but only evaluations of it
+  ##   - at τ, p(τ) is the commitment
+  ##   - and at the verification challenge z.
+  ##
+  ##   with proof = [(p(τ) - p(z)) / (τ-z)]₁
   let poly = allocHeapAligned(PolynomialEval[FIELD_ELEMENTS_PER_BLOB, matchingOrderBigInt(BLS12_381)], 64)
   let status = poly.blob_to_bigint_polynomial(blob)
-  if status == cttCodecScalar_Zero:
+  if status == cttCodecScalar_ScalarLargerThanCurveOrder:
-    return cttEthKZG_ScalarZero
+    freeHeapAligned(poly)
-  elif status == cttCodecScalar_ScalarLargerThanCurveOrder:
+    return cttEthKZG_ScalarLargerThanCurveOrder
+  elif status != cttCodecScalar_Success:
+    debugEcho "Unreachable status in blob_to_kzg_commitment: ", status
+    debugEcho "Panicking ..."
+    quit 1
+
+  var r {.noinit.}: ECP_ShortW_Aff[Fp[BLS12_381], G1]
+  kzg_commit(r, poly.evals, ctx.srs_lagrange_g1) # symbol resolution need explicit generics
+  discard dst.serialize_g1_compressed(r)
+
+  freeHeapAligned(poly)
+  return cttEthKZG_Success
+
+func compute_kzg_proof*(
+       ctx: ptr EthereumKZGContext,
+       proof_bytes: var array[48, byte],
+       y_bytes: var array[32, byte],
+       blob: ptr Blob,
+       z_bytes: array[32, byte]): CttEthKzgStatus {.tags:[Alloca, HeapAlloc, Vartime].} =
+  ## Generate:
+  ## - y = p(z), the evaluation of p at the challenge z, with p being the Blob interpreted as a polynomial.
+  ## - A zero-knowledge proof of correct evaluation.
+  ##
+  ## Mathematical description
+  ##   [proof]₁ = [(p(τ) - p(z)) / (τ-z)]₁, with p(τ) being the commitment, i.e. the evaluation of p at the powers of τ
+  ##   The notation [a]₁ corresponds to the scalar multiplication of a by the generator of 𝔾1
+  ##
+  ##   Verification can be done by verifying the relation:
+  ##     proof.(τ - z) = p(τ)-p(z)
+  ##   which doesn't require the full blob but only evaluations of it
+  ##   - at τ, p(τ) is the commitment
+  ##   - and at the verification challenge z.
+
+  # Random or Fiat-Shamir challenge
+  var z {.noInit.}: Fr[BLS12_381]
+  var status = bytes_to_bls_field(z, z_bytes)
+  if status != cttCodecScalar_Success:
+    # cttCodecScalar_Zero is not possible
     return cttEthKZG_ScalarLargerThanCurveOrder
 
-  var r {.noInit.}: ECP_ShortW_Jac[Fp[BLS12_381], G1]
+  # Blob -> Polynomial
-  r.multiScalarMul_vartime(poly.evals, ctx.srs_lagrange_g1)
+  let poly = allocHeapAligned(PolynomialEval[FIELD_ELEMENTS_PER_BLOB, Fr[BLS12_381]], 64)
+  status = poly.blob_to_field_polynomial(blob)
+  if status == cttCodecScalar_ScalarLargerThanCurveOrder:
+    freeHeapAligned(poly)
+    return cttEthKZG_ScalarLargerThanCurveOrder
+  elif status != cttCodecScalar_Success:
+    debugEcho "Unreachable status in compute_kzg_proof: ", status
+    debugEcho "Panicking ..."
+    quit 1
 
-  var r_aff {.noinit.}: ECP_ShortW_Aff[Fp[BLS12_381], G1]
+  var y {.noInit.}: Fr[BLS12_381]                         # y = p(z), eval at challenge z
-  r_aff.affine(r)
+  var proof {.noInit.}: ECP_ShortW_Aff[Fp[BLS12_381], G1] # [proof]₁ = [(p(τ) - p(z)) / (τ-z)]₁
-  discard dst.serialize_g1_compressed(r_aff)
 
-  freeHeap(poly)
+  kzg_prove(
+    proof, y,
+    poly[], ctx.domain,
+    z, ctx.srs_lagrange_g1,
+    isBitReversedDomain = true)
+
+  discard proof_bytes.serialize_g1_compressed(proof) # cannot fail
+  y_bytes.marshal(y, bigEndian) # cannot fail
+
+  freeHeapAligned(poly)
   return cttEthKZG_Success
 
 func verify_kzg_proof*(
@@ -258,7 +331,7 @@ func verify_kzg_proof*(
        commitment_bytes: array[48, byte],
        z_bytes: array[32, byte],
        y_bytes: array[32, byte],
-       proof_bytes: array[48, byte]): CttEthKzgStatus =
+       proof_bytes: array[48, byte]): CttEthKzgStatus {.tags:[Alloca, Vartime].} =
   ## Verify KZG proof that p(z) == y where p(z) is the polynomial represented by "polynomial_kzg"
 
   var commitment {.noInit.}: KZGCommitment
@@ -273,7 +346,103 @@ func verify_kzg_proof*(
   var proof {.noInit.}: KZGProof
   check proof.bytes_to_kzg_proof(proof_bytes)
 
-  let verif = kzg_verify(commitment.raw, challenge, eval_at_challenge, proof.raw, ctx.srs_monomial_g2[1])
+  let verif = kzg_verify(commitment.raw, challenge, eval_at_challenge, proof.raw, ctx.srs_monomial_g2.coefs[1])
+  if verif:
+    return cttEthKZG_Success
+  else:
+    return cttEthKZG_VerificationFailure
+
+func compute_blob_kzg_proof*(
+       ctx: ptr EthereumKZGContext,
+       proof_bytes: var array[48, byte],
+       blob: ptr Blob,
+       commitment_bytes: array[48, byte]): CttEthKzgStatus {.tags:[Alloca, HeapAlloc, Vartime].} =
+  ## Given a blob, return the KZG proof that is used to verify it against the commitment.
+  ## This method does not verify that the commitment is correct with respect to `blob`.
+
+  var commitment {.noInit.}: KZGCommitment
+  check commitment.bytes_to_kzg_commitment(commitment_bytes)
+
+  # Blob -> Polynomial
+  let poly = allocHeapAligned(PolynomialEval[FIELD_ELEMENTS_PER_BLOB, Fr[BLS12_381]], 64)
+  var status = poly.blob_to_field_polynomial(blob)
+  if status == cttCodecScalar_ScalarLargerThanCurveOrder:
+    freeHeapAligned(poly)
+    return cttEthKZG_ScalarLargerThanCurveOrder
+  elif status != cttCodecScalar_Success:
+    debugEcho "Unreachable status in compute_kzg_proof: ", status
+    debugEcho "Panicking ..."
+    quit 1
+
+  var challenge {.noInit.}: Fr[BLS12_381]
+  challenge.fiatShamirChallenge(blob[], commitment_bytes)
+
+  var y {.noInit.}: Fr[BLS12_381]                         # y = p(z), eval at challenge z
+  var proof {.noInit.}: ECP_ShortW_Aff[Fp[BLS12_381], G1] # [proof]₁ = [(p(τ) - p(z)) / (τ-z)]₁
+
+  kzg_prove(
+    proof, y,
+    poly[], ctx.domain,
+    challenge, ctx.srs_lagrange_g1,
+    isBitReversedDomain = true)
+
+  discard proof_bytes.serialize_g1_compressed(proof) # cannot fail
+
+  freeHeapAligned(poly)
+  return cttEthKZG_Success
+
+func verify_blob_kzg_proof*(
+       ctx: ptr EthereumKZGContext,
+       blob: ptr Blob,
+       commitment_bytes: array[48, byte],
+       proof_bytes: array[48, byte]): CttEthKzgStatus {.tags:[Alloca, HeapAlloc, Vartime].} =
+  ## Given a blob and a KZG proof, verify that the blob data corresponds to the provided commitment.
+
+  var commitment {.noInit.}: KZGCommitment
+  check commitment.bytes_to_kzg_commitment(commitment_bytes)
+
+  var proof {.noInit.}: KZGProof
+  check proof.bytes_to_kzg_proof(proof_bytes)
+
+  # Blob -> Polynomial
+  let poly = allocHeapAligned(PolynomialEval[FIELD_ELEMENTS_PER_BLOB, Fr[BLS12_381]], 64)
+  var status = poly.blob_to_field_polynomial(blob)
+  if status == cttCodecScalar_ScalarLargerThanCurveOrder:
+    freeHeapAligned(poly)
+    return cttEthKZG_ScalarLargerThanCurveOrder
+  elif status != cttCodecScalar_Success:
+    debugEcho "Unreachable status in compute_kzg_proof: ", status
+    debugEcho "Panicking ..."
+    quit 1
+
+  var challengeFr {.noInit.}: Fr[BLS12_381]
+  challengeFr.fiatShamirChallenge(blob[], commitment_bytes)
+
+  var challenge, eval_at_challenge {.noInit.}: matchingOrderBigInt(BLS12_381)
+  challenge.fromField(challengeFr)
+
+  let invRootsMinusZ = allocHeapAligned(array[FIELD_ELEMENTS_PER_BLOB, Fr[BLS12_381]], alignment = 64)
+
+  # Compute 1/(ωⁱ - z) with ω a root of unity, i in [0, N).
+  # zIndex = i if ωⁱ - z == 0 (it is the i-th root of unity) and -1 otherwise.
+  let zIndex = invRootsMinusZ[].inverseRootsMinusZ_vartime(
+                                  ctx.domain, challengeFr,
+                                  earlyReturnOnZero = true)
+
+  if zIndex == -1:
+    var eval_at_challenge_fr{.noInit.}: Fr[BLS12_381]
+    eval_at_challenge_fr.evalPolyAt(
+      poly[], challengeFr,
+      invRootsMinusZ[],
+      ctx.domain)
+    eval_at_challenge.fromField(eval_at_challenge_fr)
+  else:
+    eval_at_challenge.fromField(poly.evals[zIndex])
+
+  freeHeapAligned(invRootsMinusZ)
+  freeHeapAligned(poly)
+
+  let verif = kzg_verify(commitment.raw, challenge, eval_at_challenge, proof.raw, ctx.srs_monomial_g2.coefs[1])
   if verif:
     return cttEthKZG_Success
   else:

constantine/math/polynomials/polynomials.nim

@@ -18,15 +18,15 @@ import
 ## ############################################################
 
 type
-  PolynomialCoef*[N: static int, Field] = object
+  PolynomialCoef*[N: static int, Group] = object
     ## A polynomial in monomial basis
     ## [a₀, a₁, a₂, ..., aₙ]
     ##
     ## mapping to the canonical formula
-    ## p(x) = a₀ + a₁ x + a₂ x² + ... + aₙ xⁿ
+    ## p(x) = a₀ + a₁ x + a₂ x² + ... + aₙ₋₁ xⁿ⁻¹
-    coefs*{.align: 64.}: array[N, Field]
+    coefs*{.align: 64.}: array[N, Group]
 
-  PolynomialEval*[N: static int, Field] = object
+  PolynomialEval*[N: static int, Group] = object
     ## A polynomial in Lagrange basis (evaluation form)
     ## [f(0), f(ω), ..., f(ωⁿ⁻¹)]
     ## with n < 2³² and ω a root of unity
@@ -36,62 +36,82 @@ type
     ##
     ## https://ethresear.ch/t/kate-commitments-from-the-lagrange-basis-without-ffts/6950
     ## https://en.wikipedia.org/wiki/Lagrange_polynomial#Barycentric_form
-    evals*{.align: 64.}: array[N, Field]
+    evals*{.align: 64.}: array[N, Group]
 
   PolyDomainEval*[N: static int, Field] = object
     ## Metadata for polynomial in Lagrange basis (evaluation form)
+    ##
+    ## Note on inverses
+    ##   1/ωⁱ (mod N) = ωⁿ⁻ⁱ (mod N)
+    ## Hence in canonical representation rootsOfUnity[(N-i) and (N-1)] contains the inverse of rootsOfUnity[i]
+    ## This translates into rootsOfUnity[brp((N-brp(i)) and (N-1))] when bit-reversal permuted
     rootsOfUnity*{.align: 64.}: array[N, Field]
     invMaxDegree*: Field
 
 func inverseRootsMinusZ_vartime*[N: static int, Field](
        invRootsMinusZ: var array[N, Field],
        domain: PolyDomainEval[N, Field],
-       z: Field): int =
+       z: Field,
+       earlyReturnOnZero: static bool): int =
   ## Compute 1/(ωⁱ-z) for i in [0, N)
   ##
   ## Returns -1 if z ∉ {1, ω, ω², ... , ωⁿ⁻¹}
   ## Returns the index of ωⁱ==z otherwise
   ##
-  ## If ωⁱ-z == 0, the other inverses are still computed
+  ## If ωⁱ-z == 0 AND earlyReturnOnZero is false
-  ## and 0 is returned at that index.
+  ##   the other inverses are still computed
+  ##   and 0 is returned at that index
+  ## If ωⁱ-z == 0 AND earlyReturnOnZero is true
+  ##   the index of ωⁱ==z is returned
+  ##   the content of invRootsMinusZ is undefined
 
   # Mongomery's batch inversion
   # ω is a root of unity of order N,
   # so if ωⁱ-z == 0, it can only happen in one place
   var accInv{.noInit.}: Field
+  var rootsMinusZ{.noInit.}: array[N, Field]
+
+  accInv.setOne()
   var index0 = -1
 
+  when earlyReturnOnZero: # Split computation in 2 phases
+    for i in 0 ..< N:
+      rootsMinusZ[i].diff(domain.rootsOfUnity[i], z)
+      if rootsMinusZ[i].isZero().bool():
+        return i
+
   for i in 0 ..< N:
-    invRootsMinusZ[i].diff(domain.rootsOfUnity[i], z)
+    when not earlyReturnOnZero: # Fused substraction and batch inversion
+      rootsMinusZ[i].diff(domain.rootsOfUnity[i], z)
+      if rootsMinusZ[i].isZero().bool():
+        index0 = i
+        invRootsMinusZ[i].setZero()
+        continue
 
-    if invRootsMinusZ[i].isZero().bool():
+    invRootsMinusZ[i] = accInv
-      index0 = i
+    accInv *= rootsMinusZ[i]
-      continue
-
-    if i == 0:
-      accInv = invRootsMinusZ[i]
-    else:
-      accInv *= invRootsMinusZ[i]
 
   accInv.inv_vartime()
 
   for i in countdown(N-1, 1):
     if i == index0:
-      invRootsMinusZ[i].setZero()
       continue
 
     invRootsMinusZ[i] *= accInv
-    accInv *= domain.rootsOfUnity[i]
+    accInv *= rootsMinusZ[i]
 
-  invRootsMinusZ[0] *= accInv
+  if index0 == 0:
+    invRootsMinusZ[0].setZero()
+  else: # invRootsMinusZ[0] was init to accInv=1
+    invRootsMinusZ[0] = accInv
   return index0
 
-func evalPolyAt_vartime*[N: static int, Field](
+func evalPolyAt*[N: static int, Field](
        r: var Field,
        poly: PolynomialEval[N, Field],
-       domain: PolyDomainEval[N, Field],
+       z: Field,
        invRootsMinusZ: array[N, Field],
-       z: Field) =
+       domain: PolyDomainEval[N, Field]) =
   ## Evaluate a polynomial in evaluation form
   ## at the point z
   ## z MUST NOT be one of the roots of unity
@@ -103,22 +123,22 @@ func evalPolyAt_vartime*[N: static int, Field](
   for i in 0 ..< N:
     var summand {.noInit.}: Field
     summand.prod(domain.rootsOfUnity[i], invRootsMinusZ[i])
-    summand *= poly[i]
+    summand *= poly.evals[i]
     r += summand
 
   var t {.noInit.}: Field
   t = z
-  const numDoublings = log2_vartime(N) # N is a power of 2
+  const numDoublings = log2_vartime(uint32 N) # N is a power of 2
-  t.square_repeated(numDoublings)      # exponentiation by a power of 2
+  t.square_repeated(int numDoublings)         # exponentiation by a power of 2
   t.diff(Field(mres: Field.getMontyOne()), t) # TODO: refactor getMontyOne to getOne and return a field element.
   r *= t
   r *= domain.invMaxDegree
 
 func differenceQuotientEvalOffDomain*[N: static int, Field](
        r: var PolynomialEval[N, Field],
-       invRootsMinusZ: array[N, Field],
        poly: PolynomialEval[N, Field],
-       pZ: Field) =
+       pZ: Field,
+       invRootsMinusZ: array[N, Field]) =
   ## Compute r(x) = (p(x) - p(z)) / (x - z)
   ##
   ## for z != ωⁱ a power of a root of unity
@@ -131,15 +151,16 @@ func differenceQuotientEvalOffDomain*[N: static int, Field](
   for i in 0 ..< N:
     # qᵢ = (p(ωⁱ) - p(z))/(ωⁱ-z)
     var qi {.noinit.}: Field
-    qi.diff(poly[i], pZ)
+    qi.diff(poly.evals[i], pZ)
-    r[i].prod(qi, invRootsMinusZ[i])
+    r.evals[i].prod(qi, invRootsMinusZ[i])
 
 func differenceQuotientEvalInDomain*[N: static int, Field](
        r: var PolynomialEval[N, Field],
-       invRootsMinusZ: array[N, Field],
        poly: PolynomialEval[N, Field],
+       zIndex: uint32,
+       invRootsMinusZ: array[N, Field],
        domain: PolyDomainEval[N, Field],
-       zIndex: int) =
+       isBitReversedDomain: static bool) =
   ## Compute r(x) = (p(x) - p(z)) / (x - z)
   ##
   ## for z = ωⁱ a power of a root of unity
@@ -149,18 +170,14 @@ func differenceQuotientEvalInDomain*[N: static int, Field](
   ##   - rootsOfUnity:    ωⁱ
   ##   - invRootsMinusZ:  1/(ωⁱ-z)
   ##   - zIndex:          the index of the root of unity power that matches z = ωⁱᵈˣ
-  r[zIndex].setZero()
-  template invZ(): untyped =
-    # 1/z
-    #  from ωⁿ = 1 and z = ωⁱᵈˣ
-    #  hence ωⁿ⁻ⁱᵈˣ = 1/z
-    #  Note if using bit-reversal permutation (BRP):
-    #    BRP maintains the relationship
-    #    that the inverse of ωⁱ is at position n-i (mod n) in the array of roots of unity
-    static: doAssert N.isPowerOf2_vartime()
-    domain.rootsOfUnity[(N-zIndex) and (N-1)]
 
-  for i in 0 ..< N:
+  static:
+    # For powers of 2: x mod N == x and (N-1)
+    doAssert N.isPowerOf2_vartime()
+
+  r.evals[zIndex].setZero()
+
+  for i in 0'u32 ..< N:
     if i == zIndex:
       # https://dankradfeist.de/ethereum/2021/06/18/pcs-multiproofs.html
       # section "Dividing when one of the points is zero".
@@ -168,25 +185,37 @@ func differenceQuotientEvalInDomain*[N: static int, Field](
 
     # qᵢ = (p(ωⁱ) - p(z))/(ωⁱ-z)
     var qi {.noinit.}: Field
-    qi.diff(poly[i], poly[zIndex])
+    qi.diff(poly.evals[i], poly.evals[zIndex])
-    r[i].prod(qi, invRootsMinusZ[i])
+    r.evals[i].prod(qi, invRootsMinusZ[i])
 
     # q'ᵢ = -qᵢ * ωⁱ/z
     # q'idx = ∑ q'ᵢ
-    # since z is a power of ω, ωⁱ/z = ωⁱ⁻ⁱᵈˣ
-    # However some protocols use bit-reversal permutation (brp) to store the ωⁱ
-    # Hence retrieving the data would require roots[brp((brp(i)-brp(index)) mod n)] for those
-    # But is this fast? There is no single instruction for reversing bits of an integer.
-    # and the reversal depends on N.
-    # - https://stackoverflow.com/questions/746171/efficient-algorithm-for-bit-reversal-from-msb-lsb-to-lsb-msb-in-c
-    # - https://stackoverflow.com/questions/52226858/bit-reversal-algorithm-by-rutkowska
-    # - https://www.hpl.hp.com/techreports/93/HPL-93-89.pdf
-    # - https://graphics.stanford.edu/~seander/bithacks.html#BitReverseObvious
-    # The C version from Stanford's bithacks need log₂(n) loop iterations
-    # A 254~255-bit multiplication takes 38 cycles, we need 3 brp so at most ~13 cycles per brp
-    # For small Ethereum KZG, n = 2¹² = 4096, we're already at the breaking point
-    # even if an iteration takes a single cycle with instruction-level parallelism
     var ri {.noinit.}: Field
-    ri.neg(domain.rootsOfUnity[i])
+    ri.neg(r.evals[i])                                  # -qᵢ
-    ri *= invZ
+    when isBitReversedDomain:
-    r[zIndex].prod(ri, qi)
+      const logN = log2_vartime(uint32 N)
+      let invZidx = N - reverseBits(uint32 zIndex, logN)
+      let canonI = reverseBits(i, logN)
+      let idx = reverseBits((canonI + invZidx) and (N-1), logN)
+      ri *= domain.rootsOfUnity[idx]                    # -qᵢ * ωⁱ/z  (explanation at the bottom)
+    else:
+      ri *= domain.rootsOfUnity[(i+N-zIndex) and (N-1)] # -qᵢ * ωⁱ/z  (explanation at the bottom)
+    r.evals[zIndex] += ri                               # r[zIndex] = ∑ -qᵢ * ωⁱ/z
+
+    # * 1/z computation detail
+    #    from ωⁿ = 1 and z = ωⁱᵈˣ
+    #    hence ωⁿ⁻ⁱᵈˣ = 1/z
+    #    However our z may be in bit-reversal permuted
+    #
+    # * We want ωⁱ/z which translate to ωⁱ*ωⁿ⁻ⁱᵈˣ hence ωⁱ⁺ⁿ⁻ⁱᵈˣ
+    #   with the roots of unity being a cyclic group of order N so we compute i+N-zIndex (mod N)
+    #
+    #   However some protocols use bit-reversal permutation (brp) to store the ωⁱ
+    #   Hence retrieving the data requires roots[brp((brp(i)-n-brp(idx)) mod n)] for those (note: n = brp(n))
+    #
+    #   For Ethereum:
+    #     A 254~255-bit multiplication takes 11ns / 38 cycles (Fr[BLS12-381]),
+    #     A brp with n = 2¹² = 4096 (for EIP4844) takes about 6ns
+    #   We could also cache either ωⁿ⁻ⁱ or a map i' = brp(n - brp(i))
+    #   in non-brp order but cache misses are expensive
+    #   and brp can benefits from instruction-level parallelism

constantine/serialization/codecs.nim

@@ -97,6 +97,9 @@ func toHex*(bytes: openarray[byte]): string =
     result[2 + 2*i] = hexChars.secretLookup(SecretWord bi shr 4 and 0xF)
     result[2 + 2*i+1] = hexChars.secretLookup(SecretWord bi and 0xF)
 
+func fromHex*(dst: var openArray[byte], hex: string) =
+  dst.paddedFromHex(hex, bigEndian)
+
 func fromHex*[N: static int](T: type array[N, byte], hex: string): T =
   result.paddedFromHex(hex, bigEndian)
 

constantine/serialization/codecs_bls12_381.nim

@@ -61,7 +61,6 @@ type
     cttCodecEcc_PointNotInSubgroup
     cttCodecEcc_PointAtInfinity
 
-  Scalar* = matchingOrderBigInt(BLS12_381)
   G1P* = ECP_ShortW_Aff[Fp[BLS12_381], G1]
   G2P* = ECP_ShortW_Aff[Fp2[BLS12_381], G2]
 
@@ -69,7 +68,7 @@ type
 # Input validation
 # ------------------------------------------------------------------------------------------------
 
-func validate_scalar*(scalar: Scalar): CttCodecScalarStatus =
+func validate_scalar*(scalar: matchingOrderBigInt(BLS12_381)): CttCodecScalarStatus =
   ## Validate a scalar
   ## Regarding timing attacks, this will leak information
   ## if the scalar is 0 or larger than the curve order.
@@ -104,13 +103,13 @@ func validate_g2*(g2Point: G2P): CttCodecEccStatus =
 # Codecs
 # ------------------------------------------------------------------------------------------------
 
-func serialize_scalar*(dst: var array[32, byte], scalar: Scalar): CttCodecScalarStatus =
+func serialize_scalar*(dst: var array[32, byte], scalar: matchingOrderBigInt(BLS12_381)): CttCodecScalarStatus =
   ## Serialize a scalar
   ## Returns cttCodecScalar_Success if successful
   dst.marshal(scalar, bigEndian)
   return cttCodecScalar_Success
 
-func deserialize_scalar*(dst: var Scalar, src: array[32, byte]): CttCodecScalarStatus =
+func deserialize_scalar*(dst: var matchingOrderBigInt(BLS12_381), src: array[32, byte]): CttCodecScalarStatus =
   ## Deserialize a scalar
   ## Also validates the scalar range
   ##

constantine/trusted_setups/ethereum_kzg_srs.nim

@@ -23,10 +23,6 @@ import
 # Aliases
 # ------------------------------------------------------------
 
-type
-  G1Point = ECP_ShortW_Aff[Fp[BLS12_381], G1]
-  G2Point = ECP_ShortW_Aff[Fp2[BLS12_381], G2]
-
 # Presets
 # ------------------------------------------------------------
 const FIELD_ELEMENTS_PER_BLOB* {.intdefine.} = 4096
@@ -62,7 +58,7 @@ type
 
     # Trusted setup, see https://vitalik.ca/general/2022/03/14/trustedsetup.html
 
-    srs_lagrange_g1*{.align: 64.}: array[FIELD_ELEMENTS_PER_BLOB, G1Point]
+    srs_lagrange_g1*{.align: 64.}: PolynomialEval[FIELD_ELEMENTS_PER_BLOB, ECP_ShortW_Aff[Fp[BLS12_381], G1]]
     # Part of the Structured Reference String (SRS) holding the 𝔾1 points
     # This is used for committing to polynomials and producing an opening proof at
     # a random value (chosen via Fiat-Shamir heuristic)
@@ -80,7 +76,7 @@ type
     #
     # Conversion can be done with a discrete Fourier transform.
 
-    srs_monomial_g2*{.align: 64.}: array[KZG_SETUP_G2_LENGTH, G2Point]
+    srs_monomial_g2*{.align: 64.}: PolynomialCoef[KZG_SETUP_G2_LENGTH, ECP_ShortW_Aff[Fp2[BLS12_381], G2]]
     # Part of the SRS holding the 𝔾2 points
     #
     # Referring to the 𝔾2 generator as H, we store
@@ -257,7 +253,7 @@ proc loadTrustedSetup*(ctx: ptr EthereumKZGContext, filePath: string): TrustedSe
     block: # Last sanity check
       # When the srs is in monomial form we can check that
       # the first point is the generator
-      if bool(ctx.srs_monomial_g2[0] != BLS12_381.getGenerator"G2"):
+      if bool(ctx.srs_monomial_g2.coefs[0] != BLS12_381.getGenerator"G2"):
         return tsWrongPreset
 
     return tsSuccess

constantine/trusted_setups/ethereum_kzg_testing_setups_generator.nim

@@ -125,7 +125,7 @@ func getLagrange[EC](fftDesc: ECFFT_Descriptor[EC], monomial: seq[EC]): seq[EC]
   ## The polynomial is also bit-reversal permuted.
 
   result.setLen(monomial.len)
-  let status = fftDesc.ifft(result, monomial)
+  let status = fftDesc.ifft_vartime(result, monomial)
   doAssert status == FFTS_Success, "Ethereum testing trusted setup failure during Lagrange form: " & $status
 
   result.bit_reversal_permutation()
@@ -208,6 +208,7 @@ proc genEthereumKzgTestingTrustedSetup(filepath: string, secret: auto, length: i
   # Projective coordinates are slightly faster than jacobian on 𝔾1
   var fftDesc = ECFFTDescriptor[ECP_ShortW_Prj[Fp[BLS12_381], G1]].new(
     order = length, ctt_eth_kzg_fr_pow2_roots_of_unity[log2_vartime(length.uint)])
+  defer: fftDesc.delete()
 
   block: # Metadata 3 - roots of unity - bit-reversal permuted
     var meta: array[32, byte]
@@ -229,7 +230,7 @@ proc genEthereumKzgTestingTrustedSetup(filepath: string, secret: auto, length: i
 
   f.padNUL64()
 
-  block: # Data 2 - srs 𝔾2 points - bit-reversal permuted
+  block: # Data 2 - srs 𝔾2 points
     const g2Length = 65
     let ts2 = ECP_ShortW_Aff[Fp2[BLS12_381], G2].newTrustedSetupMonomial(secret, g2Length)
     # Raw dump requires little-endian
@@ -238,7 +239,7 @@ proc genEthereumKzgTestingTrustedSetup(filepath: string, secret: auto, length: i
   f.padNUL64()
 
   bit_reversal_permutation(fftDesc.rootsOfUnity.toOpenArray(0, fftDesc.order-1))
-  block: # Data 2 - roots of unity - bit-reversal permuted
+  block: # Data 3 - roots of unity - bit-reversal permuted
     # Raw dump requires little-endian
     # and we convert them all to Montgomery form
     for i in 0 ..< fftDesc.order:

constantine/trusted_setups/gen_eth_kzg_testing_setups.nim

@@ -1,249 +0,0 @@
-# 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
-  ../math/polynomials/fft,
-  ../math/[arithmetic, extension_fields, ec_shortweierstrass],
-  ../math/config/curves,
-  ../math/elliptic/[ec_scalar_mul_vartime, ec_shortweierstrass_batch_ops],
-  ../math/io/io_fields,
-  ../math/constants/zoo_generators,
-  ../platforms/abstractions,
-  ../serialization/endians,
-  std/streams
-
-# This tool generates the same testing setups that are used in Ethereum consensus-spec
-# in a Constantine-specific format specified in README.md
-
-# Trusted setup source:
-#
-# - Minimal preset: https://github.com/ethereum/consensus-specs/blob/v1.3.0/presets/minimal/trusted_setups/testing_trusted_setups.json
-# - Mainnet preset: https://github.com/ethereum/consensus-specs/blob/v1.3.0/presets/mainnet/trusted_setups/testing_trusted_setups.json
-#
-# The upstream trusted setups are stored in `./tests/protocol_ethereum_deneb_kzg`
-#
-# The upstream testing setup generator is:
-# - dump_kzg_trusted_setup_files
-#   https://github.com/ethereum/consensus-specs/blob/v1.3.0/tests/core/pyspec/eth2spec/utils/kzg.py#L96-L123
-# Called with
-# - python3 ./gen_kzg_trusted_setups.py --secret=1337 --g1-length=4 --g2-length=65
-#   python3 ./gen_kzg_trusted_setups.py --secret=1337 --g1-length=4096 --g2-length=65
-#   https://github.com/ethereum/consensus-specs/blob/v1.3.0/Makefile#L209-L210
-
-# Roots of unity
-# ------------------------------------------------------------
-#
-# Computation:
-#   Reference: https://crypto.stanford.edu/pbc/notes/numbertheory/gen.html
-#
-#   1. Find a primitive root of the finite field of modulus q
-#      i.e. root^k != 1 for all k < q-1 so powers of root generate the field.
-#
-#   sagemath: GF(r).multiplicative_generator()
-#
-#   2. primitive_root⁽ᵐᵒᵈᵘˡᵘˢ⁻¹⁾/⁽²^ⁱ⁾ for i in [0, 32)
-#
-#   sagemath: [primitive_root^((r-1)//(1 << i)) for i in range(32)]
-#
-# Usage:
-#   The roots of unity ω allow usage of polynomials in evaluation form (Lagrange basis)
-#   see ω https://dankradfeist.de/ethereum/2021/06/18/pcs-multiproofs.html
-#
-# Where does the 32 come from?
-#   Recall the definition of the BLS12-381 curve:
-#   sagemath:
-#     x = -(2^63 + 2^62 + 2^60 + 2^57 + 2^48 + 2^16)
-#     order = x^4 - x^2 + 1
-#
-#   and check the 2-adicity
-#     factor(order-1)
-#     => 2^32 * 3 * 11 * 19 * 10177 * 125527 * 859267 * 906349^2 * 2508409 * 2529403 * 52437899 * 254760293^2
-#
-#   BLS12-381 was chosen for its high 2-adicity, as 2^32 is a factor of its order-1
-
-const ctt_eth_kzg_fr_pow2_roots_of_unity = [
-  # primitive_root⁽ᵐᵒᵈᵘˡᵘˢ⁻¹⁾/⁽²^ⁱ⁾ for i in [0, 32)
-  # The primitive root chosen is 7
-  Fr[BLS12_381].fromHex"0x1",
-  Fr[BLS12_381].fromHex"0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000000",
-  Fr[BLS12_381].fromHex"0x8d51ccce760304d0ec030002760300000001000000000000",
-  Fr[BLS12_381].fromHex"0x345766f603fa66e78c0625cd70d77ce2b38b21c28713b7007228fd3397743f7a",
-  Fr[BLS12_381].fromHex"0x20b1ce9140267af9dd1c0af834cec32c17beb312f20b6f7653ea61d87742bcce",
-  Fr[BLS12_381].fromHex"0x50e0903a157988bab4bcd40e22f55448bf6e88fb4c38fb8a360c60997369df4e",
-  Fr[BLS12_381].fromHex"0x45af6345ec055e4d14a1e27164d8fdbd2d967f4be2f951558140d032f0a9ee53",
-  Fr[BLS12_381].fromHex"0x6898111413588742b7c68b4d7fdd60d098d0caac87f5713c5130c2c1660125be",
-  Fr[BLS12_381].fromHex"0x4f9b4098e2e9f12e6b368121ac0cf4ad0a0865a899e8deff4935bd2f817f694b",
-  Fr[BLS12_381].fromHex"0x95166525526a65439feec240d80689fd697168a3a6000fe4541b8ff2ee0434e",
-  Fr[BLS12_381].fromHex"0x325db5c3debf77a18f4de02c0f776af3ea437f9626fc085e3c28d666a5c2d854",
-  Fr[BLS12_381].fromHex"0x6d031f1b5c49c83409f1ca610a08f16655ea6811be9c622d4a838b5d59cd79e5",
-  Fr[BLS12_381].fromHex"0x564c0a11a0f704f4fc3e8acfe0f8245f0ad1347b378fbf96e206da11a5d36306",
-  Fr[BLS12_381].fromHex"0x485d512737b1da3d2ccddea2972e89ed146b58bc434906ac6fdd00bfc78c8967",
-  Fr[BLS12_381].fromHex"0x56624634b500a166dc86b01c0d477fa6ae4622f6a9152435034d2ff22a5ad9e1",
-  Fr[BLS12_381].fromHex"0x3291357ee558b50d483405417a0cbe39c8d5f51db3f32699fbd047e11279bb6e",
-  Fr[BLS12_381].fromHex"0x2155379d12180caa88f39a78f1aeb57867a665ae1fcadc91d7118f85cd96b8ad",
-  Fr[BLS12_381].fromHex"0x224262332d8acbf4473a2eef772c33d6cd7f2bd6d0711b7d08692405f3b70f10",
-  Fr[BLS12_381].fromHex"0x2d3056a530794f01652f717ae1c34bb0bb97a3bf30ce40fd6f421a7d8ef674fb",
-  Fr[BLS12_381].fromHex"0x520e587a724a6955df625e80d0adef90ad8e16e84419c750194e8c62ecb38d9d",
-  Fr[BLS12_381].fromHex"0x3e1c54bcb947035a57a6e07cb98de4a2f69e02d265e09d9fece7e0e39898d4b",
-  Fr[BLS12_381].fromHex"0x47c8b5817018af4fc70d0874b0691d4e46b3105f04db5844cd3979122d3ea03a",
-  Fr[BLS12_381].fromHex"0xabe6a5e5abcaa32f2d38f10fbb8d1bbe08fec7c86389beec6e7a6ffb08e3363",
-  Fr[BLS12_381].fromHex"0x73560252aa0655b25121af06a3b51e3cc631ffb2585a72db5616c57de0ec9eae",
-  Fr[BLS12_381].fromHex"0x291cf6d68823e6876e0bcd91ee76273072cf6a8029b7d7bc92cf4deb77bd779c",
-  Fr[BLS12_381].fromHex"0x19fe632fd3287390454dc1edc61a1a3c0ba12bb3da64ca5ce32ef844e11a51e",
-  Fr[BLS12_381].fromHex"0xa0a77a3b1980c0d116168bffbedc11d02c8118402867ddc531a11a0d2d75182",
-  Fr[BLS12_381].fromHex"0x23397a9300f8f98bece8ea224f31d25db94f1101b1d7a628e2d0a7869f0319ed",
-  Fr[BLS12_381].fromHex"0x52dd465e2f09425699e276b571905a7d6558e9e3f6ac7b41d7b688830a4f2089",
-  Fr[BLS12_381].fromHex"0xc83ea7744bf1bee8da40c1ef2bb459884d37b826214abc6474650359d8e211b",
-  Fr[BLS12_381].fromHex"0x2c6d4e4511657e1e1339a815da8b398fed3a181fabb30adc694341f608c9dd56",
-  Fr[BLS12_381].fromHex"0x4b5371495990693fad1715b02e5713b5f070bb00e28a193d63e7cb4906ffc93f"
-]
-
-func newTrustedSetupImpl(
-       EC: typedesc, secret: auto, length: int): seq[EC] =
-
-  result.setLen(length)
-
-  var P {.noInit.}: EC
-  P.fromAffine(EC.F.C.getGenerator($EC.G))
-  result[0] = P
-  for i in 1 ..< length:
-    P.scalarMul_minHammingWeight_windowed_vartime(secret, window = 5)
-    result[i] = P
-
-func newTrustedSetupMonomial(EC: typedesc, secret: auto, length: int): seq[EC] =
-  let ts = newTrustedSetupImpl(projective(EC), secret, length)
-  result.setLen(length)
-  batchAffine(result.asUnchecked(), ts.asUnchecked(), length)
-
-func getLagrange[EC](fftDesc: ECFFT_Descriptor[EC], monomial: seq[EC]): seq[EC] =
-  ## Get a polynomial in lagrange basis from a polynomial in monomial form.
-  ## The polynomial is also bit-reversal permuted.
-
-  result.setLen(monomial.len)
-  let status = fftDesc.ifft(result, monomial)
-  doAssert status == FFTS_Success, "Ethereum testing trusted setup failure during Lagrange form: " & $status
-
-  result.bit_reversal_permutation()
-
-func newTrustedSetupLagrange[EC](fftDesc: ECFFT_Descriptor[EC], secret: auto, length: int): auto =
-  let ts = newTrustedSetupImpl(EC, secret, length)
-  let ts2 = fftDesc.getLagrange(ts)
-
-  let tsAffine = newSeq[affine(EC)](length)
-  batchAffine(tsAffine.asUnchecked(), ts2.asUnchecked(), length)
-  return tsAffine
-
-proc padNUL64(f: FileStream) =
-  ## Pad NUL bytes until we reach a 64-byte boundary
-  let pos = f.getPosition()
-  let posMod64 = pos and 63
-
-  let pad = default(array[63, byte])
-  if posMod64 != 0:
-    f.writeData(pad[0].unsafeAddr, 64-posMod64)
-
-proc genEthereumKzgTestingTrustedSetup(filepath: string, secret: auto, length: int) =
-  ## Generate an Ethereum KZG testing trusted setup
-  ## in the Trusted Setup Interchange Format
-  ## `length` is the length of the SRS 𝔾1
-  ## the SRS 𝔾2 is fixed at 65.
-  ## SRS 𝔾1 and roots of unity are bit-reversal permuted
-
-  static: doAssert cpuEndian == littleEndian, "Trusted setup creation is only supported on little-endian CPUs at the moment."
-  doAssert length.uint.isPowerOf2_vartime(), "Expected power of 2 but found length " & $length
-
-  let f = openFileStream(filepath, fmWrite)
-  f.write"∃⋃∈∎" # ∃⋃∈∎ in UTF-8. (magic bytes)
-
-  # v1.0
-  f.write 'v'
-  f.write uint8 1
-  f.write '.'
-  f.write uint8 0
-
-  # Protocol
-  const protocol = "ethereum_deneb_kzg"
-  f.write protocol
-  let padProtocol = default(array[32 - protocol.len, byte]) # zero-init padding
-  f.writeData(padProtocol[0].unsafeAddr, padProtocol.len)
-
-  # Curve
-  const curve = "bls12_381"
-  f.write curve
-  let padCurve = default(array[15 - curve.len, byte]) # zero-init padding
-  f.writeData(padCurve[0].unsafeAddr, padCurve.len)
-
-  # Number of fields
-  f.write uint8 3
-
-  block: # Metadata 1 - srs 𝔾1 points - bit-reversal permuted
-    var meta: array[32, byte]
-    meta[0..<12] = asBytes"srs_lagrange"
-
-    meta[15..<17] = asBytes"g1"
-    meta[17..<20] = asBytes"brp"
-    meta[20..<24] = toBytes(uint32 sizeof(ECP_ShortW_Aff[Fp[BLS12_381], G1]), littleEndian)
-    meta[24..<32] = toBytes(uint64 length, littleEndian)
-
-    f.write meta
-
-  block: # Metadata 2 - srs 𝔾2 points (hardcoded to 65)
-    var meta: array[32, byte]
-    meta[0..<12] = asBytes"srs_monomial"
-
-    meta[15..<17] = asBytes"g2"
-    meta[17..<20] = asBytes"asc"
-    meta[20..<24] = toBytes(uint32 sizeof(ECP_ShortW_Aff[Fp2[BLS12_381], G2]), littleEndian)
-    meta[24..<32] = toBytes(65'u64, littleEndian)
-
-    f.write meta
-
-  # Projective coordinates are slightly faster than jacobian on 𝔾1
-  var fftDesc = ECFFTDescriptor[ECP_ShortW_Prj[Fp[BLS12_381], G1]].new(
-    order = length, ctt_eth_kzg_fr_pow2_roots_of_unity[log2_vartime(length.uint)])
-  defer: fftDesc.delete()
-
-  block: # Metadata 3 - roots of unity - bit-reversal permuted
-    var meta: array[32, byte]
-    meta[0..<11] = asBytes"roots_unity"
-
-    meta[15..<17] = asBytes"fr"
-    meta[17..<20] = asBytes"brp"
-    meta[20..<24] = toBytes(uint32 sizeof(fftDesc.rootsOfUnity[0]), littleEndian)
-    meta[24..<32] = toBytes(fftDesc.order.uint64, littleEndian)
-
-    f.write meta
-
-  f.padNUL64()
-
-  block: # Data 1 - srs 𝔾1 points - bit-reversal permuted
-    let ts1 = fftDesc.newTrustedSetupLagrange(secret, length)
-    # Raw dump requires little-endian
-    f.writeData(ts1[0].unsafeAddr, sizeof(ts1[0]) * length)
-
-  f.padNUL64()
-
-  block: # Data 2 - srs 𝔾2 points
-    const g2Length = 65
-    let ts2 = ECP_ShortW_Aff[Fp2[BLS12_381], G2].newTrustedSetupMonomial(secret, g2Length)
-    # Raw dump requires little-endian
-    f.writeData(ts2[0].unsafeAddr, sizeof(ts2[0]) * g2Length)
-
-  f.padNUL64()
-
-  bit_reversal_permutation(fftDesc.rootsOfUnity.toOpenArray(0, fftDesc.order-1))
-  block: # Data 3 - roots of unity - bit-reversal permuted
-    # Raw dump requires little-endian
-    f.writeData(fftDesc.rootsOfUnity, sizeof(fftDesc.rootsOfUnity[0]) * fftDesc.order)
-
-when isMainModule:
-  import ../math/io/io_bigints
-
-  let testSecret = BigInt[11].fromUint(1337'u64)
-  genEthereumKzgTestingTrustedSetup("trusted_setup_ethereum_kzg_test_minimal.tsif", testSecret, 4)
-  genEthereumKzgTestingTrustedSetup("trusted_setup_ethereum_kzg_test_mainnet.tsif", testSecret, 4096)