Rareskills zk Homework 8

groth16.py

from py_ecc import bn128
import numpy as np
import galois
import sys

n = bn128.curve_order
G1 = bn128.G1
G2 = bn128.G2

def mod_scalar(a):
    n = bn128.curve_order
    return int((n - (-a % n) if a < 0 else a) % n )

def mod_array(a):
    return np.array([mod_scalar(u) for u in a.flatten()]).reshape(a.shape)

def ecmul(point, scalar):
    return bn128.multiply(point, mod_scalar(scalar))

def ecadd(s, *args):
    for x in args:
        try:
            s = bn128.add(s, x)
        except Exception as e:
            try:
                s = bn128.add(s, ecmul(G1, mod_scalar(x)))   
            except Exception as e:
                s = bn128.add(s, ecmul(G2, mod_scalar(x)))
    return s

alpha_secret = 0xe3f19b7abac5dc3ea3dde7f841ae35c97383dc333895371b4242c16f1381db20
beta_secret = 0x96990440feea840e8ee46d5ded5b9e9c348fab00c4f3535fbd515e190db7dfa9
gamma_secret = 0x2fbc33181cb4c00ea71830b0f3acb01ef8ec67da0f474a57768fcab4ce1cd667
delta_secret = 0x137d011e68df5bec3073a97d45f0362762a5da287526773c60ed000c6401453a
tau = 0x31185f42d90d02b050cc19feb19c569062c326eeafd25c45986e2850e690cae7

alpha1 = ecmul(G1, alpha_secret) 
alpha2 = ecmul(G2, alpha_secret) 
beta1 = ecmul(G1, beta_secret) 
beta2 = ecmul(G2, beta_secret) 
delta1 = ecmul(G1, delta_secret) 
delta2 = ecmul(G2, delta_secret) 
gamma1 = ecmul(G1, gamma_secret)
gamma2 = ecmul(G2, gamma_secret)
delta_inv1 = ecmul(G1, pow(delta_secret, -1, n))
gamma_inv1 = ecmul(G1, pow(gamma_secret, -1, n))
alpha_times_delta_inv1 = ecmul(G1, alpha_secret * pow(delta_secret, -1, n))
beta_times_delta_inv1 = ecmul(G1, beta_secret * pow(delta_secret, -1, n))
alpha_times_gamma_inv1 = ecmul(G1, alpha_secret * pow(gamma_secret, -1, n))
beta_times_gamma_inv1 = ecmul(G1, beta_secret * pow(gamma_secret, -1, n))

L = np.array([
    [0, 0, 1, 0, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, -5, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 1],
])

R = np.array([
    [0, 0, 1, 0, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0],
])

O = np.array([
    [0, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, 0, 0, 1, 0],
    [0, 0, 0, 0, 0, 0, 1],
    [0, 1, 0, 0, 0, -1, 0],
])

GF = galois.GF(n)

L_galois = GF(mod_array(L))
R_galois = GF(mod_array(R))
O_galois = GF(mod_array(O))

x = GF(mod_scalar(4))
y = GF(mod_scalar(-2))
v1 = x * x
v2 = v1 * v1         # x^4
v3 = GF(mod_scalar(-5))*y * y
out = v3*v1 + v2    # -5y^2 * x^2

a = GF(np.array([1, out, x, y, v1, v2, v3]))

assert all(np.equal(np.matmul(L_galois, a) * np.matmul(R_galois, a), np.matmul(O_galois, a))), "not equal"

def interpolate_column(col):
    xs = GF(np.array([1, 2, 3, 4]))
    return galois.lagrange_poly(xs, col)

U_polys = np.apply_along_axis(interpolate_column, 0, L_galois)
V_polys = np.apply_along_axis(interpolate_column, 0, R_galois)
W_polys = np.apply_along_axis(interpolate_column, 0, O_galois)

def print_polys(name, polys):
    print(name, "polys:")
    for poly in polys:
        print("  ", poly)

# print_polys('U', U_polys)
# print_polys('V', V_polys)
# print_polys('W', W_polys)

def multiply_polynomials_with_witness(polys, a):
    return list(map(lambda x, y: x * y, polys, a))

def reduce_sum_polynomials(products, start, stop):
    import functools
    return functools.reduce(lambda x, y: x + y, products[start:stop])

r = 0x117160873cfc83474bfce2186d0f35f76da677fb894480ea9a549816f8b0c4e8
s = 0x3f6143c35ee2618a3a049fc7d398147c73c183fb85d2ee4429b43828d604f655

au_products = multiply_polynomials_with_witness(U_polys, a)
av_products = multiply_polynomials_with_witness(V_polys, a)
aw_products = multiply_polynomials_with_witness(W_polys, a)

au_private = reduce_sum_polynomials(au_products, 4, 7) 
av_private = reduce_sum_polynomials(av_products, 4, 7) 
aw_private = reduce_sum_polynomials(aw_products, 4, 7) 

au_public = reduce_sum_polynomials(au_products, 0, 4) 
av_public = reduce_sum_polynomials(av_products, 0, 4) 
aw_public = reduce_sum_polynomials(aw_products, 0, 4) 

au = reduce_sum_polynomials(au_products, 0, 7) 
av = reduce_sum_polynomials(av_products, 0, 7) 
aw = reduce_sum_polynomials(aw_products, 0, 7) 

def make_t(num_rows, n, GF):
    t = galois.Poly([0, 1], field=GF)
    for i in range(num_rows):
        t = t * galois.Poly([1, mod_scalar(-(i + 1))], field=GF)
    return t 

t = make_t(4, n, GF)
h = (au * av - aw) // t
ht = h * t

assert au * av == aw + ht, "division has a remainder"

au = int(au(mod_scalar(tau)))
av = int(av(mod_scalar(tau)))
aw = int(aw(mod_scalar(tau)))
ht = int(ht(mod_scalar(tau)))

au_private = int(au_private(mod_scalar(tau)))
av_private = int(av_private(mod_scalar(tau)))
aw_private = int(aw_private(mod_scalar(tau)))

au_public = int(au_public(mod_scalar(tau)))
av_public = int(av_public(mod_scalar(tau)))
aw_public = int(aw_public(mod_scalar(tau)))

A1 = ecadd(alpha1, au, ecmul(delta1, r))
B2 = ecadd(beta2, av, ecmul(delta2, s))
B1 = ecadd(beta1, av, ecmul(delta1, s))
C1 = ecadd(
    ecmul(beta_times_delta_inv1, au_private), 
    ecmul(alpha_times_delta_inv1, av_private),
    ecmul(delta_inv1, aw_private),
    ecmul(delta_inv1, ht), 
    ecmul(A1, s), 
    ecmul(B1, r), 
    ecmul(delta1, -r * s)
)

def print_G1Point(name, point):
    print(name + ".x = " + str(point[0]) + ";")
    print(name + ".y = " + str(point[1]) + ";")

def print_G2Point(name, point):
    print(name + ".x1 = " + str(eval(str(point[0]))[0]) + ";")
    print(name + ".x2 = " + str(eval(str(point[0]))[1]) + ";")
    print(name + ".y1 = " + str(eval(str(point[1]))[0]) + ";")
    print(name + ".y2 = " + str(eval(str(point[1]))[1]) + ";")

# print_G1Point('alpha1', alpha1)
# print_G1Point('beta2', beta2)
# print_G1Point('A1', A1)
# print_G2Point('B2', B2)
# print_G1Point('C1', C1)
# print_G2Point('G2', G2)

def verify(pairings):
    exponent = bn128.FQ12.one()
    for p in pairings:
        exponent = exponent * p
    return bn128.final_exponentiate(exponent) == bn128.FQ12.one()

print(verify([
    bn128.pairing(B2, bn128.neg(A1)),
    bn128.pairing(beta2, alpha1),
    bn128.pairing(
        gamma2,
        ecadd(
            ecmul(beta_times_gamma_inv1, au_public), 
            ecmul(alpha_times_gamma_inv1, av_public),
            ecmul(gamma_inv1, aw_public)
        )
    ),
    bn128.pairing(delta2, C1)
]))

impressive, very nice