Implementation of PolyCommit_{DL} from "Constant-Size Commitments to Polynomials and Their Applications" https://www.cypherpunks.ca/~iang/pubs/PolyCommit-AsiaCrypt.pdf

batch_KZG10.py

from functools import reduce
import operator

from KZG10 import TrustedSetup, CommitDivision, curve, GF, polynomial, CommitSum


"""
The batch commitment scheme,
using a `gamma` variable randomly chosen by the verifier

We then multiply the commitments by the powers of gamma to 
commit to multiple polynomials.

Say, for example we have two polynomials with 3 coefficients each

sage:

phi = lambda k, poly: poly[0] + ((k^1)*poly[1]) + ((k^2)*poly[2])
psi = lambda a, b, poly: ((phi(a, poly) - phi(b, poly)) / (a - b))

poly0 = var(['p0c%d' % _ for _ in range(3)])
poly1 = var(['p1c%d' % _ for _ in range(3)])
poly2 = var(['p2c%d' % _ for _ in range(3)])

x, z, gamma = var('x z gamma')

# Then make witnesses for the two polynomials evaluated at the same point 

h_0 = (gamma^1) * psi(x, z, poly0)
h_1 = (gamma^2) * psi(x, z, poly1)
h_2 = (gamma^3) * psi(x, z, poly2)
W = H = h_0 + h_1 + h_2
# Where `H` is equivalent to the batch witness `W`

# We can see the form the batch witness polynomial takes
sage: factor(h_0)
(p0c2*x + p0c2*z + p0c1)*gamma
sage: factor(h_0 + h_1)
(gamma*p1c2*x + gamma*p1c2*z + gamma*p1c1 + p0c2*x + p0c2*z + p0c1)*gamma
sage: factor(h_0 + h_1 + h_2)
(gamma^2*p2c2*x + gamma^2*p2c2*z + gamma^2*p2c1 + gamma*p1c2*x + gamma*p1c2*z + gamma*p1c1 + p0c2*x + p0c2*z + p0c1)*gamma

# Commit to the polynomials
cm_0 = phi(x, poly0)
cm_1 = phi(x, poly1)
cm_2 = phi(x, poly2)

# Then verifier computes the elements F and v
F_0 = (gamma^1) * cm_0
F_1 = (gamma^2) * cm_1
F_2 = (gamma^3) * cm_2
F = F_0 + F_1 + F_2

# For each of the commitments
(gamma^1) * (cm_0 - phi(z, poly0)) == h_0 * (x-z)
(gamma^2) * (cm_1 - phi(z, poly1)) == h_1 * (x-z)
(gamma^3) * (cm_2 - phi(z, poly2)) == h_2 * (x-z)

# Symbolic computations of the evaluations... 
s0 = gamma^1 * phi(z, poly0)
s1 = gamma^2 * phi(z, poly1)
s2 = gamma^3 * phi(z, poly2)
v = s0 + s1 + s2

# The pairing equation verifies that
# F-v == W * (x-z)
# Looking at the expansion, we can see they are equivalent
sage: factor(expand(W * (x-z)))
(gamma^2*p2c2*x + gamma^2*p2c2*z + gamma^2*p2c1 + gamma*p1c2*x + gamma*p1c2*z + gamma*p1c1 + p0c2*x + p0c2*z + p0c1)*gamma*(x - z)
sage: factor(expand(F-v))
(gamma^2*p2c2*x + gamma^2*p2c2*z + gamma^2*p2c1 + gamma*p1c2*x + gamma*p1c2*z + gamma*p1c1 + p0c2*x + p0c2*z + p0c1)*gamma*(x - z)

sage: factor((F-v) / (W * (x-z)))
1
"""

def main():
	F = GF(curve.curve_order)
	n_polynomials = 3
	n_coeffs = 3

	poly_coeffs = [ [F.random() for _ in range(n_coeffs)] for i in range(n_polynomials) ]

	PK = TrustedSetup.generate(F, n_coeffs, True)

	# Commit to polynomials, phi(f_i, X)
	cm = []
	for coeffs_i in poly_coeffs:
		cm.append( CommitSum(PK, coeffs_i) )

	# Prover creates a random scalar
	z = F.random()

	# Verifier sends random scalar
	gamma = F.random()

	# Prover computes the polynomial commitment h(X)
	s = []
	W = None
	for i, coeffs_i in enumerate(poly_coeffs):
		s.append(polynomial(z, coeffs_i))  # s_i = f_i(z)
		h_i = CommitDivision(PK, z, coeffs_i) # h_i = psi_z(f_i, X)
		g1_h_i = curve.multiply(h_i, int(gamma**(i+1))) # y^i * h_i
		W = g1_h_i if W is None else curve.add(W, g1_h_i)
	# Prover sends W, z and s to verifier

	# Verifier computes the element F
	g1_F = None
	for i, cm_i in enumerate(cm):
		F_i = curve.multiply(cm_i, int(gamma**(i+1)))
		g1_F = F_i if g1_F is None else curve.add(g1_F, F_i)

	# And verifier computes the element `v`
	v = reduce(operator.add, [(gamma**(i+1)) * s_i for i, s_i in enumerate(s)])
	g1_v = curve.multiply(curve.G1, int(v))

	# Then performs the pairing equation
	g1_F_minus_v = curve.add(g1_F, curve.neg(g1_v))
	g2_z = curve.multiply(curve.G2, int(z))
	g2_x_minus_z = curve.add(PK.g2_powers[1], curve.neg(g2_z))
	a = curve.pairing(curve.G2, g1_F_minus_v)
	b = curve.pairing(g2_x_minus_z, W)
	print('a', a)
	print('b', b)
	c = curve.pairing(curve.neg(g2_x_minus_z), W)
	print('a*c', a*c)


if __name__ == "__main__":
	main()

KZG10.py

from typing import List, NamedTuple, Tuple, Union
from math import ceil, log2
from random import randint
from functools import reduce
import operator
from py_ecc import bn128 as curve

"""
Implementation of PolyCommit_{DL} from:

"Constant-Size Commitments to Polynomials and Their Applications"
 - Aniket Kate, Gregory M. Zaverucha, and Ian Goldberg
 - Max Planck Institute for Software Systems (MPI-SWS)
 - https://www.cypherpunks.ca/~iang/pubs/PolyCommit-AsiaCrypt.pdf
 - http://cacr.uwaterloo.ca/techreports/2010/cacr2010-10.pdf (extended version)

Section 3.2
"""



def as_bits_bytes(n):
	n_bits = ceil(log2(n))
	n_bytes = (n_bits + (8 - (n_bits % 8))) // 8
	return n_bits, n_bytes

FIELD_MODULUS = curve.field_modulus
FIELD_MODULUS_bits, FIELD_MODULUS_bytes = as_bits_bytes(FIELD_MODULUS)

CURVE_ORDER = curve.curve_order
CURVE_ORDER_bits, CURVE_ORDER_bytes = as_bits_bytes(CURVE_ORDER)

assert FIELD_MODULUS_bytes == 32
assert CURVE_ORDER_bytes == 32

PointG1 = Tuple[curve.FQ, curve.FQ]
PointG2 = Tuple[curve.FQ2, curve.FQ2]


class Field(object):
	def __init__(self, value, modulus: int):
		if isinstance(value, Field):
			value = value.v
		else:
			value = value % modulus
		self.v = value
		self.m = modulus

	def __eq__(self, other):
		if isinstance(other, int):
			return self.v == other % self.m
		elif isinstance(other, Field):
			return self.v == other.v and self.m == other.m
		else:
			raise ValueError(f'Cannot compare {self} with {other}')

	def __int__(self):
		return self.v

	def __add__(self, other):
		if isinstance(other, Field):
			other = other.v
		return Field(self.v + other, self.m)

	def __neg__(self):
		return Field(-self.v, self.m)

	def __mul__(self, other):
		if isinstance(other, Field):
			other = other.v
		return Field(self.v * other, self.m)

	def __repr__(self):
		return f"Field<{self.v}>"

	def __sub__(self, other):
		if isinstance(other, Field):
			other = other.v
		return Field(self.v - other, self.m)

	def __truediv__(self, other):
		return self * other.inverse()

	def __pow__(self, other):
		assert isinstance(other, int)
		return Field(pow(self.v, other, self.m), self.m)

	def inverse(self):
		return Field(pow(self.v, self.m-2, self.m), self.m)


class GF(object):
	def __init__(self, modulus: int):
		self.m = modulus

	def primitive_root(self, n: int):
		"""
		Find a primitive n-th root of unity
		 - http://www.csd.uwo.ca/~moreno//AM583/Lectures/Newton2Hensel.html/node9.html
		 - https://crypto.stackexchange.com/questions/63614/finding-the-n-th-root-of-unity-in-a-finite-field
		"""
		assert n >= 2
		x = 2
		while True:
			# x != 0, g = x^(q-1/n)
			# if g^(n/2) != 1, then it is a primitive n-th root
			g = pow(int(x), (self.m-1)//n, self.m)
			if pow(g, n//2, self.m) != 1:
				return self(g)
			x += 1

	def random(self):
		return self(randint(0, self.m-1))

	def __call__(self, value) -> Field:
		return Field(value, self.m)


class TrustedSetup(NamedTuple):
	F: GF
	t: int
	g1_powers: List[PointG1]
	g2_powers: List[PointG2]
	alpha_powers: List[Field]

	@classmethod
	def generate(cls, F: GF, t: int, g1andg2: bool):
		"""
		Simulate the trusted setup
		Generate a random `alpha`
		Then return `t` powers of `alpha` in G1, and the representation of `alpha` in G2

		 g, g^a, ... g^{a^t}
		"""
		alpha = F.random()
		alpha_powers = [F(1)]
		g1_powers = [curve.G1]
		g2_powers = [curve.G2]
		for i in range(t+1):
			alpha_powers.append(alpha_powers[-1] * alpha)
			if g1andg2:
				g1_powers.append(curve.multiply(g1_powers[-1], int(alpha)))
				g2_powers.append(curve.multiply(g2_powers[-1], int(alpha)))
		
		return cls(F, t, g1_powers, g2_powers, alpha_powers)


def polynomial(x: Field, coeffs: List[Field]):
	result = coeffs[0]
	for c_i in coeffs[1:]:
		result += c_i * x
		x = x*x
	return result


def CommitProduct(PK: TrustedSetup, coeff: List[Field]):
	"""
	XXX: unsure if we need this, but it looks useful

	Computes commitment 'C' at `a`, where `a` is part of the trusted setup

		C = reduce(operator.mul, [(a**j) * c_j for j, c_j in enumerate(coeff)])

	For example:

		sage: factor((a^0 * phi0) * ((a^1) * phi1) * ((a^2) * phi2) * ((a^3)*phi3) * ((a^4)*phi4))
		a^10*(phi0*phi1*phi2*phi3*phi4)

	Where the element 'k' is the n-th triangle number
	"""
	n = len(coeff)
	k = (n*(n+1))//2  # equivalent to: reduce(operator.add, range(n))
	element = PK.g1_powers[k]
	product = PK.F(1)
	for c_i in coeffs:
		product *= c_i
	return curve.multiply(element, product)


def CommitSumTrusted(PK: TrustedSetup, coeff: List[Field]):
	return reduce(operator.add, [PK.alpha_powers[i] * c_i for i, c_i in enumerate(coeff)])


def CommitSum(PK: TrustedSetup, coeff: List[Field]):
	"""
	Copute commitment to the evaluation of a polynomial with coefficients
	At `x`, where `x` is part of the trusted setup
	"""
	return reduce(curve.add, [curve.multiply(PK.g1_powers[i], int(c_i))
							  for i, c_i in enumerate(coeff)])


def CommitRemainder(PK: TrustedSetup, y: Field, coeff: List[Field]):
	# TODO: implement
	"""
	f(x) = phi(x)
sage: f
x |--> c2*x^2 + c1*x + c0
sage: f.maxima_methods().divide((x-a) * (x-b))[1]
-a*b*c2 + (a*c2 + b*c2 + c1)*x + c0
sage: f.maxima_methods().divide((x-a))[1]
a^2*c2 + a*c1 + c0
sage: f.maxima_methods().divide((x-a) * (x-b) - (x-c))[1]
-a*b*c2 - c*c2 + (a*c2 + b*c2 + c1 + c2)*x + c0

	"""
	pass


def CommitDivisionTrusted(PK: TrustedSetup, y: Field, coeff: List[Field]):
	n = len(coeff)
	y_powers = [PK.F(1)]
	for i in range(n):
		y_powers.append(y_powers[-1] * y)

	result = PK.F(0)
	for i in range(0, n-1):
		for j in range(i, -1, -1):
			a = PK.alpha_powers[j]
			b = y_powers[i-j]
			c = coeff[i+1]
			result += a*b*c
	return result
	"""
	return reduce(operator.add, [PK.alpha_powers[j] * (y_powers[i-j] * coeff[i+1])
							     for j in range(i, -1, -1)
							     for i in range(0, n-1)])
	"""


def CommitDivision(PK: TrustedSetup, y: Field, coeff: List[Field]):
	"""
	Compute commitment to the division: `(phi(x) - phi(y)) / (x - y)`
	Using the trusted setup secret `x`

	Example:

		sage: poly4 = lambda k: c0 + ((k^1)*c1) + ((k^2)*c2) + ((k^3)*c3)
		sage: factor((poly4(x) - poly4(y)) / (x - y))
		c3*x^2 + c3*x*y + c3*y^2 + c2*x + c2*y + c1

	TODO: number of multiplications can be reduced significantly
	      number of additions can also be reduced significantly
	"""
	n = len(coeff)
	y_powers = [PK.F(1)]
	for i in range(n):
		y_powers.append(y_powers[-1] * y)

	result = None
	for i in range(0, n-1):
		for j in range(i, -1, -1):
			a = PK.g1_powers[j]
			b = y_powers[i-j]
			c = coeff[i+1]
			term = curve.multiply(a, int(b*c))
			result = term if result is None else curve.add(result, term)
	return result


def Prove():
	F = GF(curve.curve_order)
	coeff = [F.random() for _ in range(3)]
	PK = TrustedSetup.generate(F, len(coeff), True)

	"""
	The pairing equation works to verify that the witness and evaluation
	match the committed polynomial.

	sage:
		var('x y c0 c1 c2 c3')
		phi = lambda k: c0 + ((k^1)*c1) + ((k^2)*c2)
		psi = lambda a, b: ((phi(a) - phi(b)) / (a - b))
		psi(x, y) * (x-y) == phi(x) - phi(y)
		(psi3(x, y) * (x-y)) + phi(y) == phi(x)

	sage: psi(x, y) * (x-y)
	c2*x^2 - c2*y^2 + c1*x - c1*y

	sage: psi(x, y)
	(c2*x^2 - c2*y^2 + c1*x - c1*y)/(x - y)

	sage: factor(psi(x, y))
	c2*x + c2*y + c1

	sage: phi(x) - phi(y)
	c2*x^2 - c2*y^2 + c1*x - c1*y

	sage: factor(phi(x) - phi(y))
	(c2*x + c2*y + c1)*(x - y)
	"""

	# Verify with trusted information
	x = PK.alpha_powers[1]
	phi_at_x = polynomial(x, coeff)
	assert phi_at_x == CommitSumTrusted(PK, coeff)

	i = F(3)
	phi_at_i = polynomial(i, coeff)
	a = polynomial(x, coeff) - phi_at_i
	b = a / (x - i)

	psi_i_at_x = CommitDivisionTrusted(PK, i, coeff)
	assert psi_i_at_x == b
	assert psi_i_at_x * (x-i) == phi_at_x - phi_at_i


	# Then make commitment without access to trusted setup secrets
	g1_phi_at_x = CommitSum(PK, coeff)  # Commit to polynomial

	# Commit to an evaluation of the same polynomial at i
	i = F.random()  # randomly sampled i
	phi_at_i = polynomial(i, coeff)	
	g1_psi_i_at_x = CommitDivision(PK, i, coeff)

	# Compute `x - i` in G2
	g2_i = curve.multiply(curve.G2, int(i))
	g2_x_sub_i = curve.add(PK.g2_powers[1], curve.neg(g2_i)) # x-i

	# Verifier
	g1_phi_at_i = curve.multiply(curve.G1, int(phi_at_i))
	g1_phi_at_x_sub_i = curve.add(g1_phi_at_x, curve.neg(g1_phi_at_i))
	a = curve.pairing(g2_x_sub_i, g1_psi_i_at_x)
	b = curve.pairing(curve.G2, curve.neg(g1_phi_at_x_sub_i))
	ab = a*b
	print('ab', ab, ab == curve.FQ12.one())

if __name__ == "__main__":
	Prove()



def derp(n):
	"""
	Example:

	sage: poly5 = lambda k: c0 + ((k^1)*c1) + ((k^2)*c2) + ((k^3)*c3) + ((k^4)*c4) + ((k^5)*c5)
	sage: poly6 = lambda k: c0 + ((k^1)*c1) + ((k^2)*c2) + ((k^3)*c3) + ((k^4)*c4) + ((k^5)*c5)
	sage: poly5 = lambda k: c0 + ((k^1)*c1) + ((k^2)*c2) + ((k^3)*c3) + ((k^4)*c4)
	sage: poly4 = lambda k: c0 + ((k^1)*c1) + ((k^2)*c2) + ((k^3)*c3)
	sage: poly3 = lambda k: c0 + ((k^1)*c1) + ((k^2)*c2)
	sage: poly2 = lambda k: c0 + ((k^1)*c1)
	sage: poly2(x) - poly2(y)
	c1*x - c1*y
	sage: factor(poly2(x) - poly2(y))
	c1*(x - y)
	sage: factor((poly2(x) - poly2(y)) / (x - y))
	c1
	sage: factor((poly3(x) - poly3(y)) / (x - y))
	c2*x + c2*y + c1
	sage: factor((poly4(x) - poly4(y)) / (x - y))
	c3*x^2 + c3*x*y + c3*y^2 + c2*x + c2*y + c1
	sage: factor((poly5(x) - poly5(y)) / (x - y))
	c4*x^3 + c4*x^2*y + c4*x*y^2 + c4*y^3 + c3*x^2 + c3*x*y + c3*y^2 + c2*x + c2*y + c1
	sage: factor((poly6(x) - poly6(y)) / (x - y))
	c5*x^4 + c5*x^3*y + c5*x^2*y^2 + c5*x*y^3 + c5*y^4 + c4*x^3 + c4*x^2*y + c4*x*y^2 + c4*y^3 + c3*x^2 + c3*x*y + c3*y^2 + c2*x + c2*y + c1

	Emit this sequence:

	>>> derp(4)
	c[1] * x^0 * y^0
	c[2] * x^1 * y^0
	c[2] * x^0 * y^1
	c[3] * x^2 * y^0
	c[3] * x^1 * y^1
	c[3] * x^0 * y^2

	>>> derp(5)
	c[1] * x^0 * y^0
	c[2] * x^1 * y^0
	c[2] * x^0 * y^1
	c[3] * x^2 * y^0
	c[3] * x^1 * y^1
	c[3] * x^0 * y^2
	c[4] * x^3 * y^0
	c[4] * x^2 * y^1
	c[4] * x^1 * y^2
	c[4] * x^0 * y^3
	"""
	for i in range(0, n-1):
		for j in range(i, -1, -1):
			print(f'c[{i+1}] * x^{j} * y^{i-j}')

multibatch_KZG10.py

from typing import List, NamedTuple, Tuple
from functools import reduce
import operator

from KZG10 import TrustedSetup, CommitDivision, curve, GF, polynomial, CommitSum, Field, PointG1


"""
As per the batch version of KZG10
This opens many polynomials at many points

sage:

phi = lambda k, poly: poly[0] + ((k^1)*poly[1]) + ((k^2)*poly[2])
psi = lambda a, b, poly: ((phi(a, poly) - phi(b, poly)) / (a - b))

# Proving two polynomials, opened at different positions
polynomials = [var(['p%dc%d' % (_, j) for j in range(3)])
			   for _ in range(2)]

x = var('x')


z1, z2 = var('z1 z2')
s1 = phi(z1, polynomials[0])
s2 = phi(z2, polynomials[1])

cm1 = phi(x, polynomials[0])
cm2 = phi(x, polynomials[1])

gamma1, gamma2 = var('gamma1 gamma2')

W1 = gamma1 * psi(x, z1, polynomials[0])
W2 = gamma2 * psi(x, z2, polynomials[1])

r1, r2 = var('r1 r2')

v1 = (gamma1*s1)
v2 = (gamma2*s2)
F_1 = (gamma1 * cm1) - v1
F_2 = (gamma2 * cm2) - v2

F = (r1*F_1) + (r2*F_2)

lhs = F + ((r1*z1)*W1) + ((r2*z2)*W2)

rhs = ((r1 * W1) + (r2*W2)) * x

sage: factor((r1*W1+r2*W2)*x)
(gamma1*p0c2*r1*x + gamma2*p1c2*r2*x + gamma1*p0c2*r1*z1 + gamma2*p1c2*r2*z2 + gamma1*p0c1*r1 + gamma2*p1c1*r2)*x
sage: factor(expand(lhs))
(gamma1*p0c2*r1*x + gamma2*p1c2*r2*x + gamma1*p0c2*r1*z1 + gamma2*p1c2*r2*z2 + gamma1*p0c1*r1 + gamma2*p1c1*r2)*x

"""

Coefficients = List[Field]
PolynomialGroup = List[Coefficients]


class Commitment(NamedTuple):
	c: PolynomialGroup
	cm: List[PointG1]
	z: Field
	s: List[Field]

	@classmethod
	def create(cls, PK: TrustedSetup, polysys: PolynomialGroup):
		z = PK.F.random()
		s = [polynomial(z, _) for _ in polysys]
		cm = [CommitSum(PK, _) for _ in polysys]
		return cls(polysys, cm, z, s)

	def witness(self, PK: TrustedSetup, gamma: Field) -> PointG1:
		return reduce(curve.add, [curve.multiply(CommitDivision(PK, self.z, coeffs_i), int(gamma**(i+1)))
								  for i, coeffs_i in enumerate(self.c)])

	def prepare_F_and_v(self, gamma: Field) -> Tuple[PointG1, Field]:
		g1_F = reduce(curve.add, [curve.multiply(cm_i, int(gamma**(i+1))) for i, cm_i in enumerate(self.cm)])
		v = reduce(operator.add, [(gamma**(j+1)) * s_i for j, s_i in enumerate(self.s)])
		return g1_F, v

	def prepare_F_minus_v(self, gamma: Field) -> PointG1:
		g1_F, v = self.prepare_F_and_v(gamma)
		g1_v = curve.multiply(curve.G1, int(v))
		return curve.add(g1_F, curve.neg(g1_v))


def Step1_Prover(PK: TrustedSetup, many_polynomials: List[PolynomialGroup]) -> List[Commitment]:
	"""
	Commit to many groups of polynomials
	"""
	return [Commitment.create(PK, _) for _ in many_polynomials]


def Step2_Verifier(PK: TrustedSetup, step1: List[Commitment]) -> List[Field]:
	"""
	Verifier challenges prover, by providing random gamma values for each commitment
	"""
	return [PK.F.random() for _ in range(len(step1))]


def Step3_Prover(PK: TrustedSetup, step1: List[Commitment], step2: List[Field]) -> List[PointG1]:
	"""
	Prover commits to the groups of polynomials, using the gamma values from the verifier in step2
	"""
	assert len(step1) == len(step2)
	return [_.witness(PK, gamma) for gamma, _ in zip(step2, step1)]


def Step4_Verifier_Individually(PK: TrustedSetup, step1: List[Commitment], step2: List[Field], step3: List[PointG1]):
	results = []
	for group, gamma, W in zip(step1, step2, step3):
		g1_F_minus_v = group.prepare_F_minus_v(gamma)
		g2_z = curve.multiply(curve.G2, int(group.z))
		g2_x_minus_z = curve.add(PK.g2_powers[1], curve.neg(g2_z))
		a = curve.pairing(curve.G2, g1_F_minus_v)
		b = curve.pairing(curve.neg(g2_x_minus_z), W)
		results.append(a*b == curve.FQ12.one())
	return all(results)


def Step4_Verifier(PK: TrustedSetup, step1: List[Commitment], step2: List[Field], step3: List[PointG1]):
	lhs_args = list()
	rhs_args = list()
	effs = list()
	for group, gamma, W in zip(step1, step2, step3):
		r = PK.F.random()
		g1_F_minus_v = group.prepare_F_minus_v(gamma)
		effs.append(curve.multiply(g1_F_minus_v, int(r)))
		lhs_args.append(curve.multiply(W, int(r*group.z)))
		rhs_args.append(curve.multiply(W, int(r)))

	lhs_args.insert(0, reduce(curve.add, effs))

	lhs = reduce(curve.add, lhs_args)
	rhs = reduce(curve.add, rhs_args)

	a = curve.pairing(curve.G2, lhs)
	b = curve.pairing(PK.g2_powers[1], curve.neg(rhs))
	print('a', a)
	print('b', b)
	return a*b == curve.FQ12.one()


def main():
	F = GF(curve.curve_order)
	n_polynomials = 3
	n_coeffs = 3
	n_openings = 2

	groups = [[[F.random() for _ in range(n_coeffs)]
			    for i in range(n_polynomials)]
			    for j in range(n_openings)]

	PK = TrustedSetup.generate(F, n_coeffs, True)
	step1 = Step1_Prover(PK, groups)
	step2 = Step2_Verifier(PK, step1)
	step3 = Step3_Prover(PK, step1, step2)

	print('Verifying each opening individually')
	print(Step4_Verifier_Individually(PK, step1, step2, step3))

	print('Verifying all openings together in a batch')
	step4 = Step4_Verifier(PK, step1, step2, step3)
	print(step4)


if __name__ == "__main__":
	main()

polynomial.py

from typing import NamedTuple, Union, List
from functools import reduce
import operator

from KZG10 import PointG1, curve, Field, GF


class DivisionResult(NamedTuple):
    q: 'Polynomial'
    r: 'Polynomial'


class Polynomial(object):
    def __init__(self, F, coeffs=None):
        self.F = F
        coeffs = coeffs or []
        self.coeffs = [F(_) for _ in coeffs]

    def eval_g1(self, point_powers: List[PointG1]) -> PointG1:
        """
        Evaluate the polynomial at an abstract point

        Provide powers of the point, e.g.

            [g^0, g^1, g^2...]

        returns

            sum([(g^0) * c_0, g*c_1, (g^i)*c_i, ...])
        """
        assert len(point_powers) >= len(self.coeffs)
        result = None
        for i, c_i in enumerate(self.coeffs):
            x = point_powers[i]         
            y = curve.multiply(x, c_i)              
            result = y if result is None else curve.add(result, y)
        return result

    def eval_scalar(self, point: Field) -> Field:
        """
        Evaluate the polynomial at a known scalar point
        """
        return reduce(operator.add, [(point**i) * c_i for i, c_i in enumerate(self.coeffs)])

    def __repr__(self):
        x = [f'{int(c_i)}*x^{i}' for i, c_i in enumerate(self.coeffs) if i != 0]
        terms = ', '.join(x)
        return f'<Polynomial[{terms}]>'

    @classmethod
    def _compact(cls, coeffs) -> List[Field]:
        """Removes all trailing zero coefficients"""
        new_coeffs = list()
        skipping = True
        # TODO: use slicing instead of copying
        for c_i in coeffs[::-1]:
            if skipping:
                skipping = c_i == 0
            if not skipping:
                new_coeffs.append(c_i)
        return new_coeffs[::-1] 

    def compact(self) -> 'Polynomial':
        return Polynomial(self.F, self._compact(self.coeffs))

    def __call__(self, point: Union[List[PointG1], Field]) -> Union[PointG1, Field]:
        """
        Evaluate the polynomial at a specific point
        """
        if isinstance(point, list):
            assert isinstance(point[0], PointG1)
            return self.eval_g1(point)
        elif isinstance(point, Field):
            return self.eval_scalar(point)
        raise TypeError(f"Unknown how to evaluate at {type(point)}: {point}")

    def mul_poly(self, other: 'Polynomial') -> 'Polynomial':
        """
        From: https://jeremykun.com/tag/division-algorithm/
        """
        newCoeffs = [self.F(0) for _ in range(len(self) + len(other) - 1)] 
        for i,a in enumerate(self.coeffs):
            for j,b in enumerate(other.coeffs):
                newCoeffs[i+j] += a*b
        return Polynomial(self.F, newCoeffs).compact()

    def mul_scalar(self, other: Field) -> 'Polynomial':
        return Polynomial(self.F, [_ * other for _ in self.coeffs])

    def __mul__(self, other: Union[Field, 'Polynomial']) -> 'Polynomial':
        if isinstance(other, Polynomial):
            return self.mul_poly(other)
        elif isinstance(other, Field):
            return self.mul_scalar(other)
        raise TypeError(f"Unknown how to evaluate at {type(point)}: {point}")

    def __add__(self, other: 'Polynomial') -> 'Polynomial':
        new_length = max(len(self.coeffs), len(other.coeffs))
        result = [self.F(0)] * new_length
        for i in range(new_length):
            if len(self) > i:
                result[i] += self.coeffs[i]
            if len(other) > i:
                result[i] += other.coeffs[i]
        return Polynomial(self.F, result).compact()

    def __neg__(self) -> 'Polynomial':
        return Polynomial(self.F, [-c_i for c_i in self.coeffs])

    def __sub__(self, other: 'Polynomial') -> 'Polynomial':
        return self + -other

    def is_zero(self):
        return len(self.coeffs) == 0 or all([c_i == 0 for c_i in self.coeffs])

    @classmethod
    def random(cls, F: GF, n: int) -> 'Polynomial':
        return cls(F, [F.random() for _ in range(n)])

    def inverse(self):
        return Polynomial(self.F, [_.inverse() for _ in self.coeffs])

    def __iter__(self):
        return iter(self.coeffs)

    def __len__(self):
        return len(self.coeffs)

    def __getitem__(self, index: int) -> Field:
        return self.coeffs[index]

    def degree(self):
        return len(self) - 1

    def __eq__(self, other: 'Polynomial') -> 'Polynomial':
        return self.degree() == other.degree() and all([x==y for x, y in zip(self.coeffs, other.coeffs)])

    def leadingCoefficient(self) -> Field:
        return self.coeffs[-1]

    @classmethod
    def zero(cls, F):
        return cls(F, [])

    def __divmod__(self, other: 'Polynomial') -> 'Polynomial':
        # Synthetic division of polynomials
        # https://rosettacode.org/wiki/Polynomial_synthetic_division#Python
        """
        dividend = list(self.compact().coeffs[::-1])
        divisor = list(other.compact().coeffs[::-1])
        out = list(dividend)
        normalizer = divisor[0].inverse()
        for i in range(len(dividend) - (len(divisor)-1)):
            out[i] *= normalizer
            coef = out[i]
            if coef != 0:
                for j in range(1, len(divisor)):
                    out[i + j] += -divisor[j] * coef
        separator = -(len(divisor)-1)
        q = Polynomial(self.F, out[:separator][::-1]).compact()
        r = Polynomial(self.F, out[separator:][::-1]).compact()
        return DivisionResult(q, r)
        """
        # From: https://jeremykun.com/tag/division-algorithm/
        quotient, remainder = self.zero(self.F), self
        divisorDeg = other.degree()
        divisorLC = other.leadingCoefficient().inverse()

        while remainder.degree() >= divisorDeg:
            monomialExponent = remainder.degree() - divisorDeg
            monomialZeros = [self.F(0) for _ in range(monomialExponent)]
            monomialDivisor = Polynomial(self.F, monomialZeros + [remainder.leadingCoefficient() * divisorLC])

            quotient += monomialDivisor
            remainder -= monomialDivisor * other

        return DivisionResult(quotient, remainder)

    def __truediv__(self, other: 'Polynomial') -> 'Polynomial':
        q, r = divmod(self, other)
        return q


def test_division():
    Fp = GF(curve.curve_order)
    a = Polynomial(Fp, [14, 9, 1])
    b = Polynomial(Fp, [7, 1])
    q, r = divmod(a, b)
    assert q[0] == 2
    assert q[1] == 1
    assert r.is_zero() == True

    # (2x^2 - 5x - 1) / x-3
    a = Polynomial(Fp, [-1, -5, 2])
    b = Polynomial(Fp, [-3, 1])
    q, r = divmod(a, b)
    assert len(q) == 2
    assert len(r) == 1
    assert q[0] == 1
    assert q[1] == 2
    assert r[0] == 2

    # sage: f(x) = 18*x^3 - x^2 - 5*x - 8
    # sage: g(x) = x^2 - 3
    # sage: f.maxima_methods().divide(g)
    # [18*x - 1, 49*x - 11] ... ?
    # Instead, just check quotient remainder theorem
    a = Polynomial(Fp, [-8, -5, -1, 18])
    b = Polynomial(Fp, [3, 1])
    q, r = divmod(a, b)
    c = (b*q) + r
    assert a == c


def test_multiplication():
    Fp = GF(curve.curve_order)
    a = Polynomial(Fp, [Fp(_) for _ in [-5, 1]])
    b = Polynomial(Fp, [Fp(_) for _ in [2, 1]])
    c = a * b
    assert len(c) == 3
    assert c[0] == Fp(-10)
    assert c[1] == Fp(-3)
    assert c[2] == Fp(1)


def main():
    test_division()
    test_multiplication()


if __name__ == "__main__":
    main()