Write-up for ZK Hack III puzzle #2: Power corrupts

Author: Nicolas IOOSS
Date: 2022-12-05
Puzzle: https://zkhack.dev/zkhackIII/puzzleT2.html, https://github.com/ZK-Hack/puzzle-power-corrupts

1. Subject

The GitHub repository contains this introduction.

Bob has invented a new pairing-friendly elliptic curve, which he wanted to use with Groth16. For that purpose, Bob has performed a trusted setup, which resulted in an SRS containing a secret $\tau$ raised to high powers multiplied by a specific generator in both source groups. The exact parameters of the curve and part of the output of the setup are described in the document linked below.

Alice wants to recover $\tau$ and she noticed a few interesting details about the curve and the setup. Specifically, she noticed that the sum $d$ of the highest power $d_1$ of $\tau$ in $\mathbb{G}_1$ portion of the SRS, meaning the SRS contains an element of the form $\tau^{d_1} G_1$ where $G_1$ is a generator of $\mathbb{G}_1$, and the highest power $d_2$ of $\tau$ in $\mathbb{G}_2$ divides $q-1$, where $q$ is the order of the groups.

Additionally, she managed to perform a social engineering attack on Bob and extract the following information: if you express $\tau$ as $\tau = 2^{k_0 + k_1((q-1/d))} \mod r$, where $r$ is the order of the scalar field, $k_0$ is 51 bits and its fifteen most significant bits are 10111101110 (15854 in decimal). That is A < k0 < B where A = 1089478584172543 and B = 1089547303649280.

Alice then remembered the Cheon attack...

NOTE: for exponentiating $\mathbb{F}_r$ elements, use the pow_sp and pow_sp2 functions in utils.rs.

The parameters of the curve and the setup are available at https://gist.github.com/kobigurk/352036cee6cb8e44ddf0e231ee9c3f9b

The last link defines some parameters which will be used throughout this write-up:

$p = 5739212180072054691886037078074598258234828766015465682035977006377650233269727206694786492328072119776769214299497$. It is a 382-bit prime number used as the finite field over which the elliptic curve is defined
$q = 1114157594638178892192613$. It is a 80-bit prime number and is the size of the subgroup of the curve.
$E: y^2 = x^3 + 4$ is the equation of the curve which is considered. It is a pairing-friendly curve.
$k = 12$ is the embedding degree of $E(\mathbb{F}_p)$.
$\mathbb{G}_1$ is the subgroup of $E(\mathbb{F}_p)$ generated by $P$ = (5666793736468521094298738103921693621508309431638529348089160865885867964805742107934338916926827890667749984768215, 3326251098281288448602352180414320622119338868949804322483594574847275370379159993188118130786632123868776051955196). This subgroup has order $q$.
$\mathbb{G}_2$ is the subgroup of $E(\mathbb{F}_{p^2})$ generated by $Q$ = (2847190178490156899798643792842723617787968359868175140038826869144776012793105029391523604954249120667126821536281 + 1513797577242500304678874752065526230408447782356629533374984043360635354098197307045487457331199798062484459984831 * u, 2398127858646538650279262747029238501121661957103909673770298065006753715123740323569605568913154172079135187452386 + 5444946257649901533268220138726124417824817651440748374257708320447300055543369665159277001725118567443194417165086 * u)

$\tau$ is a secret integer used in the Structured Reference String (SRS) of a trusted setup. This SRS contains several points:

$\tau P \in \mathbb{G}_1$ = (4098523714512767373909357151251054769972251070043521443432671846694698445310918888674563384159054413429407635337545, 5104199250567875649831070757916425527562742910121942289068977475282886267374972984853923971702383949149250750103817)
$\tau^{d_1} P \in \mathbb{G}_1$ = (5666793736468521094298738103921693621508309431638529348089160865885867964805742107934338916926827890667749984768215, 3326251098281288448602352180414320622119338868949804322483594574847275370379159993188118130786632123868776051955196) with $d_1 = 11726539$
$\tau^{d_2} Q \in \mathbb{G}_2$ = (3536383419772898871062064633012296862124372086039789534814905834326827967479599778599887194095392165550880125330266 + 2315075704417849395497347082310199859284883937672695000597201154920791698799875018503579607990866594389648640170976 * u, 58522461936731088989461032245338237080030815519467180488197672431529427745827070450637290003234818635632376291077 + 200313434320582884299950030908390796161004965251373896142196467499133624968316891420231529223691679020778320981956 * u) with $d_2 = 690320833$

The GitHub repository contains a Rust project which uses arkworks to implement the fields and curves. This project uses different notations:

Fr, the scalar field, represents the field of integers modulo q.
Fq represents the field of integers modulo p.
Fq2 represents the field of polynomials $\mathbb{F}_p[U]$ modulo $U^2 - 23$. This field has $p^2$ elements.
Fq6 represents the field of polynomials $\mathbb{F}_{p^2}[V]$ modulo $V^3 - U$. This field has $p^6$ elements.
Fq12 represents the field of polynomials $\mathbb{F}_{p^6}[W]$ modulo $W^2 - V$. This is similar to a field of polynomials $\mathbb{F}_p[W]$ modulo $W^{12} - 23$ and has $p^{12}$ elements.
G1Affine and G1Projective represents the curve $E(\mathbb{F}_p): y^2 = x^3 + 4$. The generator defined in this file is not $P$ but another point. When using this type, we need to be careful not to rely on this generator point.
G2Affine and G2Projective represents the curve $E(\mathbb{F}_{p^2}): y^2 = x^3 + 4/U$. This $4/U$ comes from a transformation called "sextic twist" (because it reduces the degree of the of the extension field by a factor of six). Understanding why this curve is slightly different is not needed to solve the puzzle. Anyway, Appendix B of this write-up contains some helpful information on this topic.
Bls12Cheon is a type which defines a pairing with some parameters. This pairing function, Bls12Cheon::pairing will be written $e$ in this write-up. This function takes two points, $A \in \mathbb{G}_1$ and $B \in \mathbb{G}_2$, and returns a number $e(A, B) \in \mathbb{F}_{p^{12}}$. This function produces elements in a subgroup $\mathbb{G}_T \subset \mathbb{F}_{p^{12}}^{\times}$ of order $q$ generated by $e(P, Q)$.

Writing the definitions of all these notations is quite boring but is necessary to understand how this puzzle can be solved.

By the way, if you were lost, the puzzle consists in finding $\tau$ from the provided points. Now, let's solve it!

2. Cheon's Attack

In the subject, Alice noticed that $d = d_1 + d_2$ divides $q - 1$. Indeed, we can compute:

$$d = 11726539 + 690320833 = 702047372$$

$$q - 1 = 1114157594638178892192612 = d \times 1587011986761171$$

The subject hints that Cheon's attack is possible. This attack is described in a document called "Discrete Logarithm Problems with Auxiliary Inputs" (https://www.math.snu.ac.kr/~jhcheon/publications/2010/StrongDH_JoC_Final2.pdf) from Jung Hee Cheon. Some websites present how this attack can be used against Groth16 trusted setups:

https://hackmd.io/2oUhPtzWSRulLQ83Ctoy_g "Cheon's discrete log attack, and its relevance to zk-SNARKs" from Ariel Gabizon.
https://ethresear.ch/t/cheons-attack-and-its-effect-on-the-security-of-big-trusted-setups/6692 "Cheon's attack and its effect on the security of big trusted setups" from Kobi Gurkan.

This attack is about finding $\alpha$ from three elements $g$, $g^\alpha$ and $g^{\alpha^d}$ of a group of order $p$ with $d | p-1$. In the puzzle, we do not have something similar to this $g$. Instead, we have some points on elliptic curves... which are pairing-friendly. This means we can use a pairing function $e$ to transform a point $A \in \mathbb{G}_1$ to an element of group $\mathbb{G}_T$, by computing $e(A, Q)$.

Moreover, this pairing is bilinear:

$$e(\tau P, Q) = e(P, Q)^\tau$$

$$e(\tau^{d_1} P, Q) = e(P, Q)^{\tau^{d_1}}$$

$$e(\tau^{d_1} P, \tau^{d_2} Q) = e(P, Q)^{\tau^{d_1} \tau^{d_2}} = e(P, Q)^{\tau^{d_1 + d_2}} = e(P, Q)^{\tau^d}$$

By defining $g = e(P, Q) \in \mathbb{G}_T$, we can compute $g$, $g^\tau$ and $g^{\tau^d}$ which are all in the group $\mathbb{G}_T$ of order $q$.

And $d$ divides $q - 1$. Nice! Cheon's attack could work :) How?

When using 2 as a generator of the multiplicative group of $\mathbb{F}_q$, $\tau$ can be written using two numbers $0 \le k_0 < \frac{q - 1}{d}$ and $0 \le k_1 < d$:

$$\tau = 2^{k_0 + k_1 \frac{q - 1}{d}} \mod q$$

As $\tau^{q - 1} = 1 \mod q$ (according to Fermat's little theorem),

$$\tau^d = 2^{d k_0 + k_1 (q - 1)} = 2^{d k_0} = (2^d)^{k_0} \mod q$$

Enumerating all possible values of $k_0$ is quite expensive: $\frac{q - 1}{d} = 1587011986761171$ is 51-bit long. Thankfully, the subject gives the fifteen most significant bits of $k_0$. With $A' = 1089478584172544$ (which is $A + 1$ in the subject), $k_0 - A$ is at most 36-bit long. Let's split these remaining bits in two by defining $i, j$ such that:

$$0 \le i < 2^{18}$$

$$0 \le j < 2^{18}$$

$$k_0 = A' + i 2^{18} + j$$

We want to find $i$ and $j$ such that:

$$\tau^d = (2^d)^{k_0} = (2^d)^{A' + i 2^{18} + j} = (2^d)^{A' + i 2^{18}}(2^d)^j$$

This can be rewritten as:

$$\tau^d (2^{-d})^{A' + i 2^{18}} = (2^d)^j$$

Adding an exponentiation from $g$:

$$(g^{\tau^d})^{(2^{-d})^{A' + i 2^{18}}} = g^{(2^d)^j}$$

Instead of searching a 36-bit number by brute-force, we can compute all the possible values of $g^{(2^d)^j}$, keep them in memory and iterate through all the possible $i$. Here comes Shanks's Baby-Step Giant-Step algorithm! The trick here is that the overall complexity of this algorithm is in $2^{18}$ in time and space, instead of $2^{36}$ in time.

3. Arkworks' Exponentiation Bug

To implement Cheon's attack, we first need to compute $2^d$. In Rust, we can do:

let d1 = 11726539;
let d2 = 690320833;
let d: u64 = d1 + d2;
let two_pow_d = Fr::from(2).pow([d]);
println!("2^d = {}", two_pow_d);

This displays 749844686810961635743610, which is wrong! For example by computing pow(2, 11726539 + 690320833, 1114157594638178892192613) in Python, we obtain 981957778667961401372462. The subject provides a very helpful note:

NOTE: for exponentiating $\mathbb{F}_r$ elements, use the pow_sp and pow_sp2 functions in utils.rs.

pow_sp is indeed a regular modular exponentiation algorithm. Let's give it a try!

let two_pow_d = pow_sp(Fr::from(2u64), d.into(), 31);
println!("2^d = {}", two_pow_d);

This displays 981957778667961401372462, which is the same a with Python. So Fr::pow is buggy and gives wrong result. This is due to a bug in the way 128-bit field elements are squared in arkworks 0.4.0-alpha.3. For example this code should compute $2^2 = 4$

println!("x^2 = {}", Fr::from(2).square_in_place());

It displays 120063879707056511493220. At the time of writing, a Pull Request has been opened to fix this issue: arkworks-rs/algebra#526. Without this fix, we need to avoid using Fr::pow when solving the puzzle.

By the way, this bug does not seem to affect Fq12::pow, which relies on Fp384 instead of Fp128. In practice, when using Fq12::pow to solve the puzzle, the solution is correct.

Now that we know we should use pow_sp instead of Fr::pow, how can be implement Cheon's attack?

4. Implementation of Baby-Step Giant-Step

A previous section started explaining Cheon's attack and its first step: using Baby-Step Giant-Step to recover $k_0$.

We first need to compute $g^{(2^d)^j}$ for many $j$ into a table which gives $j$ for each value which is computed. This is exactly what function pow_sp2 does!

There is comment which states that this function "assumes exp is 80bits". As $q$ has 80 bits, this is always true.

So the first part of Baby-Step Giant-Step algorithm is:

let g: Fq12 = Bls12Cheon::pairing(P, Q).0;
let g_tau_d: Fq12 = Bls12Cheon::pairing(tau_d1_P, tau_d2_Q).0; // = g^(tau^d)

let baby_steps = pow_sp2(g, two_pow_d, 1 << 18);

Then, the given steps consists in computing $(g^{\tau^d})^{(2^{-d})^{A' + i 2^{18}}}$ for every $i$ until one is found in baby_steps.

let mut giant = g_tau_d.pow(pow_sp(two_pow_d.inverse().unwrap(), k0_base, 51).into_bigint());
let giant_inc = pow_sp(two_pow_d.inverse().unwrap(), 1 << 18, 64).into_bigint();
for i in 0..q1_d / (1 << 18) {
    if let Some(j) = baby_steps.get(&giant) {
        k0 = k0_base + (i << 18) + u128::from(*j);
        println!("Found i = {}, j = {}, k0 = {}", i, j, k0);
        break;
    }
    giant = giant.pow(giant_inc);
}

On my computer (using an Intel i7-8565U CPU), Baby-Step takes 11 seconds and Giant-Step 156 seconds before displaying:

Found i = 233602, j = 226194, k0 = 1089539821761426

Now that we know $k_0$, we can compute $\tau^d = (2^d)^{k_0} \mod q$. This was the first part of Cheon's attack.

By the way, when using $2^{20}$ instead of $2^{18}$ to split $k_0$, the same algorithm is faster. It finds $i = 58400, j = 750482$ in 39 + 39 seconds. These parameters may be tweaked depending on the computer.

5. Finishing Cheon's Attack

The second part of Cheon's attack consists in finding $k_1$.

The idea is to use the same Baby-Step Giant-Step algorithm but with a different generator:

$$\eta = 2^\frac{q - 1}{d} \mod q$$

Indeed, $\tau$ can be written:

$$\tau = 2^{k_0 + k_1 \frac{q - 1}{d}} = 2^{k_0} \eta^{k_1} \mod q$$

$$\tau 2^{-k_0} = \eta^{k_1} \mod q$$

$$(g^\tau)^{2^{-k_0}} = g^{\eta^{k_1}}$$

As $k_1 < d$ and $d$ has 30 bits, we can split $k_1 = i 2^{15} + j$ with two 15-bit integers $(i, j)$ and apply the Baby-Step Giant-Step algorithm. It finds $i = 21075$ and $j = 14120$ in 3 + 14 seconds. So $k_1 = 690599720$.

Finally,

$$\tau = 2^{k_0 + k_1 \frac{q - 1}{d}} = 284865198031253921498207 \mod q$$

This breaks the trusted setup! And this concludes this puzzle with: $\tau$ = 284865198031253921498207.

Appendix A: Full Rust Code of the Attack

use crate::utils::{bigInt_to_u128, pow_sp, pow_sp2};
use ark_bls12_cheon::{Bls12Cheon, Fq12, Fr, G1Projective as G1, G2Projective as G2};
use ark_ec::pairing::Pairing;
use ark_ff::{Field, PrimeField};
use std::collections::HashMap;
use std::ops::Mul;
use std::time::Instant;

pub fn attack(P: G1, tau_P: G1, tau_d1_P: G1, Q: G2, tau_d2_Q: G2) -> i128 {
    let d1 = 11726539;
    let d2 = 690320833;
    let d: u64 = d1 + d2;
    let q = 1114157594638178892192613u128;
    let q1_d = (q - 1) / u128::from(d);

    // Translate the problem into the target group
    // (".0" needs to be used here to extract the Fq12 from the PairingOutput type)
    let g: Fq12 = Bls12Cheon::pairing(P, Q).0;
    let g_tau: Fq12 = Bls12Cheon::pairing(tau_P, Q).0; // = g^tau
    let g_tau_d: Fq12 = Bls12Cheon::pairing(tau_d1_P, tau_d2_Q).0; // = g^(tau^d)

    // Find k0 such that tau^d = (2^d)^k0 mod q
    let two_pow_d = pow_sp(Fr::from(2u64), d.into(), 31);
    let g_2_d = g_tau_d.pow(two_pow_d.into_bigint());
    println!("2^d = {}", two_pow_d);

    // We know the high bits of k0 and the subject provides A and B such as A < k0 < B
    let k0_base: u128 = 0x3dee000000000;
    assert_eq!(1089478584172543, k0_base - 1); // A
    assert_eq!(1089547303649280, k0_base + (1 << 36)); // B

    // Shanks's Baby-Step Giant-Step algorithm: with 36 bits to find, split in two
    // To find i, j such that g^((2^d)^k0) = g^(tau^d) and k0 = k0_base + i * 2^20 + j,
    // Compute g^((2^d)^j) and for each i, search (g^(tau^d))^(2^(-d*k0_A))^((2^(-d*2^20))^i)
    let mut k0 = 0u128;
    let bsgs_split = 1u64 << 20;
    println!("Baby step for k0...");
    let start = Instant::now();
    let baby_steps = pow_sp2(g, two_pow_d, bsgs_split);
    println!(
        "Baby step for k0 done in {} seconds",
        start.elapsed().as_secs()
    );

    let mut giant = g_tau_d.pow(pow_sp(two_pow_d.inverse().unwrap(), k0_base, 51).into_bigint());
    let giant_inc = pow_sp(two_pow_d.inverse().unwrap(), bsgs_split.into(), 64).into_bigint();
    println!("Giant step for k0...");
    let start = Instant::now();
    for i in 0..q1_d / u128::from(bsgs_split) {
        if let Some(j) = baby_steps.get(&giant) {
            let duration = start.elapsed().as_secs();
            k0 = k0_base + i * u128::from(bsgs_split) + u128::from(*j);
            println!(
                "Found i = {}, j = {}, k0 = {} in {} seconds",
                i, j, k0, duration
            );
            break;
        }
        giant = giant.pow(giant_inc);
    }
    assert_eq!(k0, 1089539821761426);

    // Ensure that tau^d == 2^(d*k0) modulo q by checking the pairing
    assert_eq!(g_tau_d, g.pow(pow_sp(two_pow_d, k0, 51).into_bigint()));

    // Compute 2^((q - 1) / d)
    let eta = pow_sp(Fr::from(2u64), q1_d, 51);

    // Shanks's Baby-Step Giant-Step algorithm: with 30 bits to find
    // Find i, j such that (g^tau)^(2^(-k0)) = g^(eta^k1) with k1 = i * 2^17 + j
    // Compute g^(eta^j) and for each i, search (g^tau)^(2^(-k0))((eta^(-2^17))^i)
    let mut k1 = 0u128;
    let bsgs_split = 1u64 << 16;
    println!("Baby step for k1...");
    let start = Instant::now();
    let baby_steps = pow_sp2(g, eta, bsgs_split);
    println!(
        "Baby step for k1 done in {} seconds",
        start.elapsed().as_secs()
    );

    let mut giant = g_tau.pow(pow_sp(Fr::from(2u64).inverse().unwrap(), k0, 51).into_bigint());
    let giant_inc = pow_sp(eta.inverse().unwrap(), bsgs_split.into(), 64).into_bigint();
    println!("Giant step for k1...");
    let start = Instant::now();
    for i in 0..d / bsgs_split {
        if let Some(j) = baby_steps.get(&giant) {
            let duration = start.elapsed().as_secs();
            k1 = u128::from(i * bsgs_split + *j);
            println!(
                "Found i = {}, j = {}, k1 = {} in {} seconds",
                i, j, k1, duration
            );
            break;
        }
        giant = giant.pow(giant_inc);
    }
    assert_eq!(k1, 690599720);

    let tau = pow_sp(Fr::from(2u64), k0 + k1 * q1_d, 80);
    println!("Found tau = {}", tau);
    assert_eq!(tau, Fr::from(284865198031253921498207u128));
    assert_eq!(P.mul(tau), tau_P);

    return bigInt_to_u128(tau.into_bigint()) as i128;
}

This outputs:

2^d = 981957778667961401372462
Baby step for k0...
Baby step for k0 done in 39 seconds
Giant step for k0...
Found i = 58400, j = 750482, k0 = 1089539821761426 in 41 seconds
Baby step for k1...
Baby step for k1 done in 4 seconds
Giant step for k1...
Found i = 10537, j = 46888, k1 = 690599720 in 7 seconds
Found tau = 284865198031253921498207

Appendix B: Using SageMath

SageMath is a powerful software. It greatly helps working on mathematics problems and implementing cryptographic attacks. It can be easily installed on Ubuntu using apt install sagemath and running sage in a terminal. However, working with curves like BLS12-381 with SageMath can be quite challenging due to the difficulty of defining a pairing function.

This Appendix aims at understanding the structure of the curve used in the puzzle. Indeed, it used a specific prime number and a curve which is not BLS12-381. But the curve is inspired by BLS12-381 and websites such as https://hackmd.io/@benjaminion/bls12-381 ("BLS12-381 For The Rest Of Us"), https://gist.github.com/HarryR/eb5ad0e5de51633678e015a6b06969a1 ("Sage script to derive all necessary parameters for BLS12-381 curve"), https://eprint.iacr.org/2020/875.pdf ("Efficient Final Exponentiation via Cyclotomic Structure for Pairings over Families of Elliptic Curves") and https://research.nccgroup.com/2020/07/13/pairing-over-bls12-381-part-2-curves/ ("Pairing over BLS12-381" series) can be much helpful.

The curve in the puzzle uses a parameter x = 0xdf017120670c6e6a, defined in bls12_cheon/src/curves/mod.rs. This parameter actually defines the prime number:

param_x = 0xdf017120670c6e6a
p = (param_x - 1)^2 * (param_x^4 - param_x^2 + 1) // 3 + param_x
assert p == 5739212180072054691886037078074598258234828766015465682035977006377650233269727206694786492328072119776769214299497

With a curve of embedding degree 12, we are interested in a curve defined in $\mathbb{F}_{p^{12}}$. In the puzzle, the coefficients are defined in bls12_cheon/src/curves/g1.rs as a = 0, b = 4. So the curve equation is $E: y^2 = x^3 + 4$. Here is how we can create such a curve in SageMath:

Fp = GF(p)
G1 = EllipticCurve(Fp, (0, 4))
P = G1.point((
    4366845981406663127346140105392043296067620632748305894915559567990751463871846461571751242076416842353760718219463,
    2322936086289479068066490612801140015408721761579169972917817367391045591350600034547617432914931369293155123935728,
))

# Check some properties
assert G1.order() == p - param_x

Two groups are defined on $E(\mathbb{F}_{p^{12}})$: $\mathbb{G}_1$ which uses coordinates in $\mathbb{F}_p$ and $\mathbb{G}'_2$ which is twisted into a group $\mathbb{G}_2$ with coordinates in $\mathbb{F}_{p^2}$. Before presenting the sextic twist, how can the fields be defined in SageMath?

To define a field $\mathbb{F}_{p^k}$, we need to specify which polynomial to use to extend the field appropriately. The puzzle defines these polynomials in each file of bls12_cheon/src/fields/. This leads to a tower of field extensions:

$$\mathbb{F}_{p^2} = \mathbb{F}_p[U] / (U^2 - 23)$$

$$\mathbb{F}{p^6} = \mathbb{F}{p^2}[V] / (V^3 - U)$$

$$\mathbb{F}{p^{12}} = \mathbb{F}{p^6}[W] / (W^2 - V)$$

In SageMath, this can be written:

K2.<U> = PolynomialRing(Fp)
Fp2.<u> = Fp.extension(U^2 - 23, 'u')

K6.<V> = PolynomialRing(Fp2)
Fp6.<v> = Fp2.extension(V2^3 - u)

K12.<W> = PolynomialRing(Fp6)
Fp12.<w> = Fp6.extension(W^2 - v)

However this explodes: SageMath tries to perform some computations which take a very long time to complete. This issue was reported two years ago on https://ask.sagemath.org/question/49663/efficiently-computing-tower-fields-for-pairings/ ("Efficiently computing tower fields for pairings"). It is a known problem and the way to work around it consists in skipping Fp6 when defining Fp12. Instead of using a polynomial which makes $w^2 = v$, we can use a polynomial which makes $w^6 = v^3 = u$.

K2.<U> = PolynomialRing(Fp)
Fp2.<u> = Fp.extension(U^2 - 23)

K12.<W> = PolynomialRing(Fp2)
Fp12.<w> = Fp2.extension(W^6 - u)
v = w^2

This takes few seconds instead and it can be used in a similar way: it all boils down to $\mathbb{F}_{p^{12}}$ being a field of polynomials with 12 coefficients in $\mathbb{F}_p$ taken modulo $W^{12} - 23$.

Now, what is the sextic twist of $\mathbb{G}'_2$? The parameters in bls12_cheon/src/curves/mod.rs defines a divisive twist type (according to the documentation of ark_ec::models::bls12::TwistType ; this is the opposite of BLS12-381 which uses a multiplicative twist type). This means that points $(x', y') \in \mathbb{G}'_2$ are mapped to points $(x, y) \in \mathbb{G}_2$ with:

$$x = \frac{x'}{w^2}$$

$$y = \frac{y'}{w^3}$$

This new point is on a curve with equation:

$$y^2 = \frac{y'^2}{w^6} = \frac{x'^3 + 4}{w^6} = \frac{x'^3}{w^6} + \frac{4}{u} = x^3 + \frac{4}{u}$$

And voila! This is the equation of points in $\mathbb{G}_2$. The coordinates of the points in this group are in $\mathbb{F}_{p^2}$ instead of $\mathbb{F}_{p^{12}}$ but this is trivial projection of the polynomials, as the group is generated by a point $Q \in \mathbb{F}_{p^2}^2$ and the curve operations make it stay in this set.

G2 = EllipticCurve(Fp2, (0, 4 / u))
Q = G2.point((
    2847190178490156899798643792842723617787968359868175140038826869144776012793105029391523604954249120667126821536281 +
    1513797577242500304678874752065526230408447782356629533374984043360635354098197307045487457331199798062484459984831 * u,
    2398127858646538650279262747029238501121661957103909673770298065006753715123740323569605568913154172079135187452386 +
    5444946257649901533268220138726124417824817651440748374257708320447300055543369665159277001725118567443194417165086 * u,
))

assert G2.order() == (p - param_x) * (p + param_x + 1) + (param_x^4 - param_x^2 + 1)

By the way, bls12_cheon/src/curves/g2.rs uses a large number to define this curve. This verifies that it is $4 / u$.

assert 4 / u == 499061928701917799294438006789095500716072066610040494090084957076317411588671931016937955854614966937110366460826 * u

The generators $P$ and $Q$ were chosen so that they generate groups of $q = 1114157594638178892192613$ elements. We can check that q divides $p^{12} - 1$, and in fact divides the 12th cyclotomic polynomial $\Phi_{12}(p) = p^4 - p^2 + 1$

q = 1114157594638178892192613
assert q.is_prime()
assert (p^12 - 1) % q == 0
assert (p^4 - p^2 + 1) % q == 0

In SageMath, checking that $P$ and $Q$ are of order $q$ is simple.

assert q * P == 0
assert q * Q == 0
# These are slow
assert P.order() == q
assert Q.order() == q

With this, it is possible to define the pairing function used by arkworks, in SageMath.

# Define the pairing as similar as possible to arkworks:
# - https://github.com/arkworks-rs/algebra/blob/v0.4.0-alpha.3/ec/src/models/bls12/mod.rs
# - https://gist.github.com/cygnusv/f41ee4026dc56cec693a306b1926b472
# - https://research.nccgroup.com/2020/08/13/pairing-over-bls12-381-part-3-pairing/
# - https://crypto.stanford.edu/pbc/notes/ep/miller.html
def double_eval_for_miller_loop(A, P):
    # Compute the tangent of A in G2, untwisted to the curve on Fp12, then evaluate at P in G1
    # This function is usually written g_{A,A}(P)
    (x_A, y_A) = A.xy()  # in Fp2.<u>
    (x_P, y_P) = P.xy()  # in Fp.<u>
    return (3 * x_A^2) * x_P * w + - (2 * y_A) * y_P + (2 * y_A^2 - 3 * x_A^3) * w^3


def add_eval_for_miller_loop(A, B, P):
    # Compute the line between A and B in G2, untwisted to the curve on Fp12, then evaluate at P in G1
    # This function is usually written g_{A,B}(P)
    (x_A, y_A) = A.xy()  # in Fp2.<u>
    (x_B, y_B) = B.xy()  # in Fp2.<u>
    (x_P, y_P) = P.xy()  # in Fp.<u>
    return (y_B - y_A) * x_P * w + (x_A - x_B) * y_P + (x_B * y_A - y_B * x_A) * w^3


def miller_loop(P, Q):
    bits_x = "{:b}".format(param_x)[1:]
    V = Q
    f = 1
    for bit_char in bits_x:
        f = f^2 * double_eval_for_miller_loop(V, P)
        V = V + V
        if bit_char == "1":
            f = f * add_eval_for_miller_loop(V, Q, P)
            V = V + Q
    return f


def final_exponentiation(f):
    # Compute f^((p^6 - 1)(p^2 + 1)(3 + (x - 1)(x - 1)(p + x)(x^2 + p^2 - 1))
    f_p6 = (((((f ^ p) ^ p) ^ p) ^ p) ^ p) ^ p
    r = f_p6 / f
    r_p2 = (r ^ p) ^ p
    r = r_p2 * r
    r3 = r^3
    tmp = ((r ^ (param_x - 1)) ^ (param_x - 1)) ^ (p + param_x)
    tmp = ((tmp ^ param_x) ^ param_x) * ((tmp ^ p) ^ p) / tmp
    return r3 * tmp


def pairing(P, Q):
    return final_exponentiation(miller_loop(P, Q))

For example computing pairing(P, Q) gives:

(2596687636486138420394379509962158672058256403889805443963896318232572236714713485252941701728591712624927032780492*u +
 4303501139115967476262904356705480817784306099580583090821857746172182684903701393121516196830521897114032826974427)*w^5 +
(5008010838243441695640882112202762672290378615201925731231185094071542277263809397079549568681504922053250863170771*u +
 1596997608270217820284409577980476028723719785708631423309006109395942036464559736380502982206926355075767746142974)*w^4 +
(73361568819745165664407692475320932461064678353803870409113048257281260434507678817330471775366990061990861123025*u +
 1339223782441700725691244882967696375519190880575999438334442105140927207899445186184202053352487251883706158858093)*w^3 +
(3084949435723987792361351554140461672700881784628158409122838513843817713388798825566183591987564920907356542558160*u +
 504870203720958150703083229714841289668311437599225812754632447890282487631780993024470871008366357453753673068376)*w^2 +
(232865636422510956473097645204868117990681890931583054092657258701339079744586374465102200728629622348904206102206*u +
 1465499049673074406040935342411290095942936938383984208191127200925432090417683448390466041787571632691901526482424)*w +
 4931839319440249149390561440646934939741776044920779405513298363034651764899270310506725401216615002132343166816510*u +
 2539716661100746027014011411866504595638032217324382986916999470045068691715591410927522528893435259478098135288250

In Rust, Bls12Cheon::pairing(P, Q) gives the same result, formatted differently:

QuadExtField(
  CubicExtField(
    QuadExtField(
      2539716661100746027014011411866504595638032217324382986916999470045068691715591410927522528893435259478098135288250 +
      4931839319440249149390561440646934939741776044920779405513298363034651764899270310506725401216615002132343166816510 * u
    ), QuadExtField(
      504870203720958150703083229714841289668311437599225812754632447890282487631780993024470871008366357453753673068376 +
      3084949435723987792361351554140461672700881784628158409122838513843817713388798825566183591987564920907356542558160 * u
    ), QuadExtField(
      1596997608270217820284409577980476028723719785708631423309006109395942036464559736380502982206926355075767746142974 +
      5008010838243441695640882112202762672290378615201925731231185094071542277263809397079549568681504922053250863170771 * u
  )) + CubicExtField(
    QuadExtField(
      1465499049673074406040935342411290095942936938383984208191127200925432090417683448390466041787571632691901526482424 +
      232865636422510956473097645204868117990681890931583054092657258701339079744586374465102200728629622348904206102206 * u
    ), QuadExtField(
      1339223782441700725691244882967696375519190880575999438334442105140927207899445186184202053352487251883706158858093 +
      73361568819745165664407692475320932461064678353803870409113048257281260434507678817330471775366990061990861123025 * u
    ), QuadExtField(
      4303501139115967476262904356705480817784306099580583090821857746172182684903701393121516196830521897114032826974427 +
      2596687636486138420394379509962158672058256403889805443963896318232572236714713485252941701728591712624927032780492 * u
  )) * u)

With all this, it is possible to check that the solution of the puzzle was correct, despite possible bugs in arkworks.

tau_P = G1.point((
    4098523714512767373909357151251054769972251070043521443432671846694698445310918888674563384159054413429407635337545,
    5104199250567875649831070757916425527562742910121942289068977475282886267374972984853923971702383949149250750103817,
))
tau_d1_P = G1.point((
    5666793736468521094298738103921693621508309431638529348089160865885867964805742107934338916926827890667749984768215,
    3326251098281288448602352180414320622119338868949804322483594574847275370379159993188118130786632123868776051955196,
))
tau_d2_Q = G2.point((
    3536383419772898871062064633012296862124372086039789534814905834326827967479599778599887194095392165550880125330266 +
    2315075704417849395497347082310199859284883937672695000597201154920791698799875018503579607990866594389648640170976 * u,
    58522461936731088989461032245338237080030815519467180488197672431529427745827070450637290003234818635632376291077 +
    200313434320582884299950030908390796161004965251373896142196467499133624968316891420231529223691679020778320981956 * u,
))

tau = 284865198031253921498207  # Puzzle solution

Fr = GF(q)
assert tau_P == tau * P
assert tau_d1_P == (Fr(tau)^11726539).lift() * P
assert tau_d2_Q == (Fr(tau)^690320833).lift() * Q

g = pairing(P, Q)
g_tau = pairing(tau_P, Q)
g_tau_d = pairing(tau_d1_P, tau_d2_Q)

# Check that the points are not the one and are q-th roots of unity in Fp12
assert g != 1 and g^q == 1
assert g_tau != 1 and g_tau^q == 1
assert g_tau_d != 1 and g_tau_d^q == 1

# Verify the solution of the puzzle in Fp12
assert g_tau == g^tau
assert g_tau_d == g^((Fr(tau)^702047372).lift())

niooss-ledger/zkhack2022_III_puzzle2_writeup.md