zkhack2021_puzzle2_writeup.py

#!/usr/bin/env python3
"""Solve zkHack puzzle #2 (https://www.zkhack.dev/puzzle2.html)

On https://github.com/kobigurk/zkhack-trusted-setup there is a Rust program
which displays:

    Alice has computed a trusted setup for a Groth16 proof scheme.
    She decided to use a 128-bit long secret, and she swears that she does not
    know the secret s needed to get this setup.
    The trusted setup is constructed as follows using two additional scalars α
    and β:
    * [s^i] G1 for 0 ⩽ i ⩽ 62,
    * [α s^i] G1 for 0 ⩽ i ⩽ 31,
    * [β s^i] G1 for 0 ⩽ i ⩽ 31,
    * [s^i] G2 for 0 ⩽ i ⩽ 31.

    Can you recover the secret anyway?

The file src/data.rs
(https://github.com/kobigurk/zkhack-trusted-setup/blob/76d366b6304c8eb34c9f2244bb96a63fd61f60ba/src/data.rs)
contains 63+32+32 = 127 points on the first curve of BLS12 381, and 32 points on
the second curve.

When trying to see whether the points are on the prime subgroups of order
52435875175126190479447740508185965837690552500527637822603658699938581184513,
it appears that working with ark_* crates in Rust is very hard. Indeed there
does not seem to be a simple way to perform a scalar multiplication with this
number. The first issue is that importing a number in decimal form seems to only
be possible in ark_bls12_381::Fr. But this subgroup is actually numbers modulo
52435875175126190479447740508185965837690552500527637822603658699938581184513,
so in fact the following code is incorrect:

    use ark_ff::biginteger::BigInteger256;

    fn main() {
        let (_ts1, _ts2) = puzzle_data();
        let BLS12_381_ORDER: BigInteger256 = Fr::from_str(
            "52435875175126190479447740508185965837690552500527637822603658699938581184513",
        ).unwrap().into();
        println!("ts1[0] * order = {}", _ts1[0].mul(BLS12_381_ORDER));
    }

This code computes "_ts1[0] * 0", which displays GroupAffine(Infinity), even
though the result of _ts1[0] * 52435875175126190479447740508185965837690552500527637822603658699938581184513
is not the infinity.

When trying to replace Fr::from_str with another method, it quickly appeared
that using BigInteger256 was not suitable when performing multiplications, and
that ark_ec really expects the scalar to be given modulo
52435875175126190479447740508185965837690552500527637822603658699938581184513.

To work around these issues, I decided to work on the puzzle using Python, and
a Python-only implementation of BLS12 381.
I used https://github.com/Chia-Network/bls-signatures/tree/b3acfe5ecd23660f8ce21a2aafb9d14c3263fef7/python-impl
as reference.

Here is my Python code, with some comments.
First, let's import some standard modules.
"""
import re
from pathlib import Path
import sys
import math

# This scripts expects that https://github.com/kobigurk/zkhack-trusted-setup and
# https://github.com/Chia-Network/bls-signatures are cloned next to it.
SCRIPT_DIR = Path(__file__).parent
BLS_SIGNATURES_PYTHON_DIR = SCRIPT_DIR / "bls-signatures" / "python-impl"

sys.path.insert(0, str(BLS_SIGNATURES_PYTHON_DIR))
import bls12381
import ec

# Load the puzzle data
PUZZLE2_DATA_DIR = SCRIPT_DIR / "zkhack-trusted-setup" / "src" / "data.rs"
ts1_points = []
ts2_points = []
with PUZZLE2_DATA_DIR.open("r") as rust_file:
    for line in rust_file:
        matches = re.match(r'^ *let ts1_([0-9]+)_bytes = hex::decode\("([0-9a-f]+)"', line)
        if matches:
            idx = int(matches.group(1))
            assert len(ts1_points) == idx
            pt_bytes = bytes.fromhex(matches.group(2))
            assert len(pt_bytes) == 96
            x = bls12381.Fq(bls12381.q, int.from_bytes(pt_bytes[:48], "little"))
            y = bls12381.Fq(bls12381.q, int.from_bytes(pt_bytes[48:], "little"))
            pt = ec.AffinePoint(x, y, False, ec.default_ec)
            assert pt.is_on_curve()
            ts1_points.append(pt)
        matches = re.match(r'^ *let ts2_([0-9]+)_bytes = hex::decode\("([0-9a-f]+)"', line)
        if matches:
            idx = int(matches.group(1))
            assert len(ts2_points) == idx
            pt_bytes = bytes.fromhex(matches.group(2))
            assert len(pt_bytes) == 192
            x = bls12381.Fq2(
                bls12381.q, int.from_bytes(pt_bytes[:48], "little"), int.from_bytes(pt_bytes[48:96], "little")
            )
            y = bls12381.Fq2(
                bls12381.q, int.from_bytes(pt_bytes[96:144], "little"), int.from_bytes(pt_bytes[144:], "little")
            )
            pt = ec.AffinePoint(x, y, False, ec.default_ec_twist)
            assert pt.is_on_curve()
            ts2_points.append(pt)

assert len(ts1_points) == 127
assert len(ts2_points) == 32


# The first generator is not in the prime subgroup!
BLS12_381_ORDER = 52435875175126190479447740508185965837690552500527637822603658699938581184513
print(f"G1 * order = {ts1_points[0] * BLS12_381_ORDER}")

"""
This displays:

    G1 * order = AffinePoint(x=Fq(0x110a5..7af0e), y=Fq(0x68dfa..c0ad3), i=False)

If the point was in the prime subgroup, it would display:

    AffinePoint(x=Fq(0x0), y=Fq(0x0), i=True)

(This is the result of ec.G1Generator() * BLS12_381_ORDER)

Now, what is the order of G1? To answer, I found the cofactors of the BLS12 381
curve on https://hackmd.io/@benjaminion/bls12-381
"""
G1_COFACTOR = 0x396C8C005555E1568C00AAAB0000AAAB
G2_COFACTOR = 0x5D543A95414E7F1091D50792876A202CD91DE4547085ABAA68A205B2E5A7DDFA628F1CB4D9E82EF21537E293A6691AE1616EC6E786F0C70CF1C38E31C7238E5

print(f"G1 * order * G1_cofactor = {ts1_points[0] * BLS12_381_ORDER * G1_COFACTOR}")

"""
This displays:

    G1 * order * G1_cofactor = AffinePoint(x=Fq(0x0), y=Fq(0x0), i=True)

So the result is the point at infinity. This means that the order of G1 is a
factor of order * G1_cofactor.

To compute ord(G1), we first need to factorize G1_COFACTOR
"""

g1_cofactor_factors = []
remaining = G1_COFACTOR
factor = 3
while factor * factor <= remaining:
    while (remaining % factor) == 0:
        remaining //= factor
        print(f"Found G1_cofactor factor {factor}, remaining {remaining}")
        g1_cofactor_factors.append(factor)
    factor += 2
if remaining != 1:
    g1_cofactor_factors.append(remaining)

print(f"g1_cofactor_factors = {g1_cofactor_factors}")

"""
This displays:

    g1_cofactor_factors = [3, 11, 11, 10177, 10177, 859267, 859267]

Using these prime factors, it is possible to compute ord(G1) by reduction
"""

puzzle_g1_order = BLS12_381_ORDER * G1_COFACTOR
assert (ts1_points[0] * puzzle_g1_order).infinity
for f in g1_cofactor_factors:
    while puzzle_g1_order % f == 0 and (ts1_points[0] * (puzzle_g1_order // f)).infinity:
        print(f"Reducing temporary G1 order by {f}")
        puzzle_g1_order //= f

print(f"ord(G1) = {puzzle_g1_order}")

"""
This displays:

    ord(G1) = 41608392149640139394573879782083739887453721694564462058295541683768164056223783970589006205076495868811

Now, how does this help recovering the secret s? In the data:

- the first point, ts1_points[0], is defined as "[s^0] G1", which is G1,
- the second point, ts1_points[1], is defined as "[s^1] G1", which is s * G1

As ord(G1) * G1 is the point at infinity (by definition of the order),
ord(G1) * (s * G1) is also this point, so a theorem (known as "Lagrange's theorem")
states that:

    ord(s * G1) divides ord(G1)

So there exists an integer k >= 1 such that ord(G1) = k * ord(s * G1). Then it
is possible to prove that k divides s in this situation.

In practice, computing the order of s * G1 leads to: ord(s * G1) = ord(G1)
"""
for f in g1_cofactor_factors:
    if puzzle_g1_order % f == 0 and (ts1_points[1] * (puzzle_g1_order // f)).infinity:
        print(f"Found reduction of ord(s * G1) by {f}")  # This message is never displayed

"""
So this way does not enable guessing a diviser of s.

Nevertheless, we can compute (s - v) * G1 for any integer v, as:

    (s - v) * G1 = s * G1 - v * G1 = ts1_points[1] - v * ts1_points[0]

If ord((s - v) * G1) < ord(G1), then using a similar reasoning, there exists
an integer k such that:

    ord(G1) = k * ord((s - v) * G1)

    k divides (s - v)

    s mod k = v mod k

Moreover by testing all prime factors of ord(G1), it is possible to construct
these k.

This leads to a first algorithm to compute s mod k with k a factor of ord(G1):
"""


def compute_s_mod_k_using_g1(k):
    assert puzzle_g1_order % k == 0
    complement_factor = puzzle_g1_order // k
    assert not (ts1_points[0] * k).infinity

    for value in range(0, k):
        test_pt = ts1_points[1] - ts1_points[0] * value
        if (test_pt * complement_factor).infinity:
            print(f"Found s % {k} == {value}")
            return value


assert puzzle_g1_order == 3 * 11 * 10177 * 859267 * 52437899 * BLS12_381_ORDER
compute_s_mod_k_using_g1(3)
compute_s_mod_k_using_g1(11)

"""
This algorithms leads to:

    Found s % 3 == 2
    Found s % 11 == 4

Nevertheless it is quite slow for larger number and it can be optimized. Indeed
instead of performing a scalar multiplication by ord(G1)/k at every iteration,
it is possible to precompute the result and the steps between the successive
points.
"""


def compute_s_mod_k_using_g1_optimized(k):
    assert puzzle_g1_order % k == 0
    complement_factor = puzzle_g1_order // k
    assert not (ts1_points[0] * k).infinity

    test_pt = ts1_points[1] * complement_factor
    step_pt = ts1_points[0] * complement_factor
    for value in range(0, k):
        if test_pt.infinity:
            print(f"Found (with optim) s % {k} == {value}")
            return value
        test_pt = test_pt - step_pt


compute_s_mod_k_using_g1_optimized(3)
compute_s_mod_k_using_g1_optimized(11)
compute_s_mod_k_using_g1_optimized(10177)
compute_s_mod_k_using_g1_optimized(859267)
compute_s_mod_k_using_g1_optimized(52437899)

"""
Actually this code raises an exception because of a bug in the used library.
With the patch from https://github.com/Chia-Network/bls-signatures/pull/288,
this code works fine and displays:

    Found (with optim) s % 3 == 2
    Found (with optim) s % 11 == 4
    Found (with optim) s % 10177 == 5487
    Found (with optim) s % 859267 == 809074
    Found (with optim) s % 52437899 == 3648959

It is possible to then use the Chinese Remainder Theorem (CRT) to combine all
these constraints into one
"""


def modinv(a, m):
    """Compute the inverse of a modulo m"""
    return pow(a, -1, m)


def combine_crt(mod1, val1, mod2, val2):
    assert math.gcd(mod1, mod2) == 1
    new_mod = mod1 * mod2
    base_1 = modinv(mod2, mod1) * mod2
    base_2 = modinv(mod1, mod2) * mod1
    new_val = (val1 * base_1 + val2 * base_2) % new_mod
    assert new_val % mod1 == val1
    assert new_val % mod2 == val2
    return new_mod, new_val


combined_mod, combined_val = 1, 0
for mod, val in ((3, 2), (11, 4), (10177, 5487), (859267, 809074), (52437899, 3648959)):
    combined_mod, combined_val = combine_crt(combined_mod, combined_val, mod, val)

print(f"Combined: s % {combined_mod} == {combined_val}")

"""
This displays:

    Combined: s % 15132376222941642753 == 2335387132884273659

Repeating the same process on G2 leads to more constraints on s.
"""
g2_cofactor_factors = []
remaining = G2_COFACTOR
factor = 3
while factor * factor <= remaining and factor < 300000:
    while (remaining % factor) == 0:
        remaining //= factor
        print(f"Found G2_cofactor factor {factor}, remaining {remaining}")
        g2_cofactor_factors.append(factor)
    factor += 2
if remaining != 2:
    g2_cofactor_factors.append(remaining)

print(f"g2_cofactor_factors = {g2_cofactor_factors}")

puzzle_g2_order = BLS12_381_ORDER * G2_COFACTOR
assert (ts2_points[0] * puzzle_g2_order).infinity
for f in g2_cofactor_factors:
    while puzzle_g2_order % f == 0 and (ts2_points[0] * (puzzle_g2_order // f)).infinity:
        print(f"Reducing temporary G2 order by {f}")
        puzzle_g2_order //= f

print(f"ord(G2) = {puzzle_g2_order}")

for f in g2_cofactor_factors:
    if puzzle_g2_order % f == 0 and (ts2_points[1] * (puzzle_g2_order // f)).infinity:
        print(f"Found reduction of ord(s * G2) by {f}")  # This message is never displayed


def compute_s_mod_k_using_g2_optimized(k):
    assert puzzle_g2_order % k == 0
    complement_factor = puzzle_g2_order // k
    assert not (ts2_points[0] * k).infinity

    test_pt = ts2_points[1] * complement_factor
    step_pt = ts2_points[0] * complement_factor
    for value in range(0, k):
        if test_pt.infinity:
            print(f"Found (with optim) s % {k} == {value}")
            return value
        test_pt = test_pt - step_pt


LARGE_G2_ORDER_PRIME = 402096035359507321594726366720466575392706800671181159425656785868777272553337714697862511267018014931937703598282857976535744623203249
assert puzzle_g2_order == 13 * 23 * 2713 * 11953 * 262069 * LARGE_G2_ORDER_PRIME * BLS12_381_ORDER
compute_s_mod_k_using_g2_optimized(13)
compute_s_mod_k_using_g2_optimized(23)
compute_s_mod_k_using_g2_optimized(2713)
compute_s_mod_k_using_g2_optimized(11953)
compute_s_mod_k_using_g2_optimized(262069)

"""
This displays:

    g2_cofactor_factors = [13, 13, 23, 23, 2713, 11953, 262069, 402096035359507321594726366720466575392706800671181159425656785868777272553337714697862511267018014931937703598282857976535744623203249]
    Reducing temporary G2 order by 13
    Reducing temporary G2 order by 23
    ord(G2) = 53576194808460553217203172723474824808344883275353018929305633902963132938098184426089150766896967993227272192442055963012237813401673323498414509448806171002635876817216134661849632732818260738617639084599158761113588926055983
    Found (with optim) s % 13 == 2
    Found (with optim) s % 23 == 15
    Found (with optim) s % 2713 == 1184
    Found (with optim) s % 11953 == 10669
    Found (with optim) s % 262069 == 204769
"""

for mod, val in ((13, 2), (23, 15), (2713, 1184), (11953, 10669), (262069, 204769)):
    combined_mod, combined_val = combine_crt(combined_mod, combined_val, mod, val)

print(f"Combined: s % {combined_mod} == {combined_val}")

"""
The new combined constraint is:

    Combined: s % 38452154918091875653578148163112927 == 5592216610550884993006174526481245

38452154918091875653578148163112927 is 115-bit wide. Apparently Alice used a
128-bit secret, so there remains approximately 128 - 115 = 13 bits to find.
This can be brute-forced by finding an integer b such that:

    s = b * 38452154918091875653578148163112927 + 5592216610550884993006174526481245
"""
test_pt = ts1_points[0] * combined_val
step_pt = ts1_points[0] * combined_mod
for b in range(2 ** 14):
    if test_pt == ts1_points[1]:
        s = b * combined_mod + combined_val
        print(f"Found b = {b}, s = {s}")
        break
    test_pt = test_pt + step_pt


# Check that s is right
assert ts1_points[0] * s == ts1_points[1]
assert ts1_points[1] * s == ts1_points[2]

assert ts2_points[0] * s == ts2_points[1]
assert ts2_points[1] * s == ts2_points[2]

for i in range(62):
    assert ts1_points[i] * s == ts1_points[i + 1]  # Check s^(i+1) G1
for i in range(63, 63 + 31):
    assert ts1_points[i] * s == ts1_points[i + 1]  # Check alpha * s^(i+1) G1
for i in range(63 + 32, 63 + 32 + 31):
    assert ts1_points[i] * s == ts1_points[i + 1]  # Check beta * s^(i+1) G1
for i in range(31):
    assert ts2_points[i] * s == ts2_points[i + 1]  # Check s^(i+1) G2

"""
This leads to:

    Found b = 2989, s = 114939083266787167213538091034071020048

Back to the Rust project, injecting s leads to a program which does not fail:

    fn main() {
        let (_ts1, _ts2) = puzzle_data();

        /* Your solution here! (s in decimal)*/
        let s = Fr::from_str("114939083266787167213538091034071020048").unwrap();

        assert_eq!(_ts1[0].mul(s), _ts1[1]);
        assert_eq!(_ts2[0].mul(s), _ts2[1]);
    }

In short, it was possible to recover s thanks to the use of points which orders
contained small prime factors. These small factors enabled performing a discrete
logarithm attack to recover the secret.
"""