diff --git a/bindings/c_kzg.swg b/bindings/c_kzg.swg
index 36dcf08..60be872 100644
--- a/bindings/c_kzg.swg
+++ b/bindings/c_kzg.swg
@@ -54,6 +54,7 @@
 }
 
 %array_class(fr_t, frArray)
+%array_class(byte, byteArray)
 %array_class(g1_t, g1Array)
 %array_class(g2_t, g2Array)
 %pointer_class(poly, polyp)
diff --git a/bindings/python/fft.py b/bindings/python/fft.py
new file mode 100644
index 0000000..7ff24ac
--- /dev/null
+++ b/bindings/python/fft.py
@@ -0,0 +1,120 @@
+
+from py_ecc import optimized_bls12_381 as b
+
+def _simple_ft(vals, modulus, roots_of_unity):
+    L = len(roots_of_unity)
+    o = []
+    for i in range(L):
+        last = b.Z1 if type(vals[0]) == tuple else 0
+        for j in range(L):
+            if type(vals[0]) == tuple:
+                last = b.add(last, b.multiply(vals[j], roots_of_unity[(i*j)%L]))
+            else:
+                last += vals[j] * roots_of_unity[(i*j)%L]
+        o.append(last if type(last) == tuple else last % modulus)
+    return o
+
+def _fft(vals, modulus, roots_of_unity):
+    if len(vals) <= 4 and type(vals[0]) != tuple:
+        #return vals
+        return _simple_ft(vals, modulus, roots_of_unity)
+    elif len(vals) == 1 and type(vals[0]) == tuple:
+        return vals
+    L = _fft(vals[::2], modulus, roots_of_unity[::2])
+    R = _fft(vals[1::2], modulus, roots_of_unity[::2])
+    o = [0 for i in vals]
+    for i, (x, y) in enumerate(zip(L, R)):
+        y_times_root = b.multiply(y, roots_of_unity[i]) if type(y) == tuple else y*roots_of_unity[i]
+        o[i] = b.add(x, y_times_root) if type(x) == tuple else (x+y_times_root) % modulus
+        o[i+len(L)] = b.add(x, b.neg(y_times_root)) if type(x) == tuple else (x-y_times_root) % modulus
+    return o
+
+def expand_root_of_unity(root_of_unity, modulus):
+    # Build up roots of unity
+    rootz = [1, root_of_unity]
+    while rootz[-1] != 1:
+        rootz.append((rootz[-1] * root_of_unity) % modulus)
+    return rootz
+
+def fft(vals, modulus, root_of_unity, inv=False):
+    rootz = expand_root_of_unity(root_of_unity, modulus)
+    # Fill in vals with zeroes if needed
+    if len(rootz) > len(vals) + 1:
+        vals = vals + [0] * (len(rootz) - len(vals) - 1)
+    if inv:
+        # Inverse FFT
+        invlen = pow(len(vals), modulus-2, modulus)
+        if type(vals[0]) == tuple:
+            return [b.multiply(x, invlen) for x in
+                    _fft(vals, modulus, rootz[:0:-1])]
+        else:
+            return [(x*invlen) % modulus for x in
+                    _fft(vals, modulus, rootz[:0:-1])]
+    else:
+        # Regular FFT
+        return _fft(vals, modulus, rootz[:-1])
+
+# Evaluates f(x) for f in evaluation form
+def inv_fft_at_point(vals, modulus, root_of_unity, x):
+    if len(vals) == 1:
+        return vals[0]
+    # 1/2 in the field
+    half = (modulus + 1)//2
+    # 1/w
+    inv_root = pow(root_of_unity, len(vals)-1, modulus)
+    # f(-x) in evaluation form
+    f_of_minus_x_vals = vals[len(vals)//2:] + vals[:len(vals)//2]
+    # e(x) = (f(x) + f(-x)) / 2 in evaluation form
+    evens = [(f+g) * half % modulus for f,g in zip(vals, f_of_minus_x_vals)]
+    # o(x) = (f(x) - f(-x)) / 2 in evaluation form
+    odds = [(f-g) * half % modulus for f,g in zip(vals, f_of_minus_x_vals)]
+    # e(x^2) + coordinate * x * o(x^2) in evaluation form
+    comb = [(o * x * inv_root**i + e) % modulus for i, (o, e) in enumerate(zip(odds, evens))]
+    return inv_fft_at_point(comb[:len(comb)//2], modulus, root_of_unity ** 2 % modulus, x**2 % modulus)
+
+def shift_domain(vals, modulus, root_of_unity, factor):
+    if len(vals) == 1:
+        return vals
+    # 1/2 in the field
+    half = (modulus + 1)//2
+    # 1/w
+    inv_factor = pow(factor, modulus - 2, modulus)
+    half_length = len(vals)//2
+    # f(-x) in evaluation form
+    f_of_minus_x_vals = vals[half_length:] + vals[:half_length]
+    # e(x) = (f(x) + f(-x)) / 2 in evaluation form
+    evens = [(f+g) * half % modulus for f,g in zip(vals, f_of_minus_x_vals)]
+    print('e', evens)
+    # o(x) = (f(x) - f(-x)) / 2 in evaluation form
+    odds = [(f-g) * half % modulus for f,g in zip(vals, f_of_minus_x_vals)]
+    print('o', odds)
+    shifted_evens = shift_domain(evens[:half_length], modulus, root_of_unity ** 2 % modulus, factor ** 2 % modulus)
+    print('se', shifted_evens)
+    shifted_odds = shift_domain(odds[:half_length], modulus, root_of_unity ** 2 % modulus, factor ** 2 % modulus)
+    print('so', shifted_odds)
+    return (
+        [(e + inv_factor * o) % modulus for e, o in zip(shifted_evens, shifted_odds)] + 
+        [(e - inv_factor * o) % modulus for e, o in zip(shifted_evens, shifted_odds)]
+    )
+
+def shift_poly(poly, modulus, factor):
+    factor_power = 1
+    inv_factor = pow(factor, modulus - 2, modulus)
+    o = []
+    for p in poly:
+        o.append(p * factor_power % modulus)
+        factor_power = factor_power * inv_factor % modulus
+    return o
+
+def mul_polys(a, b, modulus, root_of_unity):
+    rootz = [1, root_of_unity]
+    while rootz[-1] != 1:
+        rootz.append((rootz[-1] * root_of_unity) % modulus)
+    if len(rootz) > len(a) + 1:
+        a = a + [0] * (len(rootz) - len(a) - 1)
+    if len(rootz) > len(b) + 1:
+        b = b + [0] * (len(rootz) - len(b) - 1)
+    x1 = _fft(a, modulus, rootz[:-1])
+    x2 = _fft(b, modulus, rootz[:-1])
+    return _fft([(v1*v2)%modulus for v1,v2 in zip(x1,x2)],
+               modulus, rootz[:0:-1])
diff --git a/bindings/python/kzg_proofs.py b/bindings/python/kzg_proofs.py
new file mode 100644
index 0000000..9b3afb0
--- /dev/null
+++ b/bindings/python/kzg_proofs.py
@@ -0,0 +1,252 @@
+from py_ecc import optimized_bls12_381 as b
+from fft import fft
+from multicombs import lincomb
+
+# Generatore for the field
+PRIMITIVE_ROOT = 5
+MODULUS = b.curve_order
+
+assert pow(PRIMITIVE_ROOT, (MODULUS - 1) // 2, MODULUS) != 1
+assert pow(PRIMITIVE_ROOT, MODULUS - 1, MODULUS) == 1
+
+#########################################################################################
+#
+# Helpers
+#
+#########################################################################################
+
+def is_power_of_two(x):
+    return x > 0 and x & (x-1) == 0
+
+
+def generate_setup(s, size):
+    """
+    # Generate trusted setup, in coefficient form.
+    # For data availability we always need to compute the polynomials anyway, so it makes little sense to do things in Lagrange space
+    """
+    return (
+        [b.multiply(b.G1, pow(s, i, MODULUS)) for i in range(size + 1)],
+        [b.multiply(b.G2, pow(s, i, MODULUS)) for i in range(size + 1)],
+    )
+
+#########################################################################################
+#
+# Field operations
+#
+#########################################################################################
+
+
+def get_root_of_unity(order):
+    """
+    Returns a root of unity of order "order"
+    """
+    assert (MODULUS - 1) % order == 0
+    return pow(PRIMITIVE_ROOT, (MODULUS - 1) // order, MODULUS)
+
+def inv(a):
+    """
+    Modular inverse using eGCD algorithm
+    """
+    if a == 0:
+        return 0
+    lm, hm = 1, 0
+    low, high = a % MODULUS, MODULUS
+    while low > 1:
+        r = high // low
+        nm, new = hm - lm * r, high - low * r
+        lm, low, hm, high = nm, new, lm, low
+    return lm % MODULUS
+
+def div(x, y):
+    return x * inv(y) % MODULUS
+
+#########################################################################################
+#
+# Polynomial operations
+#
+#########################################################################################
+
+def eval_poly_at(p, x):
+    """
+    Evaluate polynomial p (coefficient form) at point x
+    """
+    y = 0
+    power_of_x = 1
+    for i, p_coeff in enumerate(p):
+        y += power_of_x * p_coeff
+        power_of_x = (power_of_x * x) % MODULUS
+    return y % MODULUS
+
+def div_polys(a, b):
+    """
+    Long polynomial difivion for two polynomials in coefficient form
+    """
+    a = [x for x in a]
+    o = []
+    apos = len(a) - 1
+    bpos = len(b) - 1
+    diff = apos - bpos
+    while diff >= 0:
+        quot = div(a[apos], b[bpos])
+        o.insert(0, quot)
+        for i in range(bpos, -1, -1):
+            a[diff + i] -= b[i] * quot
+        apos -= 1
+        diff -= 1
+    return [x % MODULUS for x in o]
+
+#########################################################################################
+#
+# Utils for reverse bit order
+#
+#########################################################################################
+
+
+def reverse_bit_order(n, order):
+    """
+    Reverse the bit order of an integer n
+    """
+    assert is_power_of_two(order)
+    # Convert n to binary with the same number of bits as "order" - 1, then reverse its bit order
+    return int(('{:0' + str(order.bit_length() - 1) + 'b}').format(n)[::-1], 2)
+    
+
+def list_to_reverse_bit_order(l):
+    """
+    Convert a list between normal and reverse bit order. This operation is idempotent.
+    """
+    return [l[reverse_bit_order(i, len(l))] for i in range(len(l))]
+
+
+#########################################################################################
+#
+# Converting between polynomials (in coefficient form) and data (in reverse bit order)
+# and extending data
+#
+#########################################################################################
+
+def get_polynomial(data):
+    """
+    Interpolate a polynomial (coefficients) from data in reverse bit order
+    """
+    assert is_power_of_two(len(data))
+    root_of_unity = get_root_of_unity(len(data))
+    return fft(list_to_reverse_bit_order(data), MODULUS, root_of_unity, True)
+
+def get_data(polynomial):
+    """
+    Get data (in reverse bit order) from polynomial in coefficient form
+    """
+    assert is_power_of_two(len(polynomial))
+    root_of_unity = get_root_of_unity(len(polynomial))
+    return list_to_reverse_bit_order(fft(polynomial, MODULUS, root_of_unity, False))
+
+def get_extended_data(polynomial):
+    """
+    Get extended data (expanded by 2x, reverse bit order) from polynomial in coefficient form
+    """
+    assert is_power_of_two(len(polynomial))
+    extended_polynomial = polynomial + [0] * len(polynomial)
+    root_of_unity = get_root_of_unity(len(extended_polynomial))
+    return list_to_reverse_bit_order(fft(extended_polynomial, MODULUS, root_of_unity, False))
+
+#########################################################################################
+#
+# Kate single proofs
+#
+#########################################################################################
+
+def commit_to_poly(polynomial, setup):
+    """
+    Kate commitment to polynomial in coefficient form
+    """
+    return lincomb(setup[0][:len(polynomial)], polynomial, b.add, b.Z1)
+
+def compute_proof_single(polynomial, x, setup):
+    """
+    Compute Kate proof for polynomial in coefficient form at position x
+    """
+    quotient_polynomial = div_polys(polynomial, [-x, 1])
+    return lincomb(setup[0][:len(quotient_polynomial)], quotient_polynomial, b.add, b.Z1)
+
+def check_proof_single(commitment, proof, x, y, setup):
+    """
+    Check a proof for a Kate commitment for an evaluation f(x) = y
+    """
+    # Verify the pairing equation
+    #
+    # e([commitment - y], [1]) = e([proof],  [s - x])
+    #    equivalent to
+    # e([commitment - y]^(-1), [1]) * e([proof],  [s - x]) = 1_T
+    #
+
+    s_minus_x = b.add(setup[1][1], b.multiply(b.neg(b.G2), x))
+    commitment_minus_y = b.add(commitment, b.multiply(b.neg(b.G1), y))
+
+    pairing_check = b.pairing(b.G2, b.neg(commitment_minus_y), False)
+    pairing_check *= b.pairing(s_minus_x, proof, False)
+    pairing = b.final_exponentiate(pairing_check)
+
+    return pairing == b.FQ12.one()
+
+#########################################################################################
+#
+# Kate multiproofs on a coset
+#
+#########################################################################################
+
+def compute_proof_multi(polynomial, x, n, setup):
+    """
+    Compute Kate proof for polynomial in coefficient form at positions x * w^y where w is
+    an n-th root of unity (this is the proof for one data availability sample, which consists
+    of several polynomial evaluations)
+    """
+    quotient_polynomial = div_polys(polynomial, [-pow(x, n, MODULUS)] + [0] * (n - 1) + [1])
+    return lincomb(setup[0][:len(quotient_polynomial)], quotient_polynomial, b.add, b.Z1)
+
+def check_proof_multi(commitment, proof, x, ys, setup):
+    """
+    Check a proof for a Kate commitment for an evaluation f(x w^i) = y_i
+    """
+    n = len(ys)
+    root_of_unity = get_root_of_unity(n)
+    
+    # Interpolate at a coset. Note because it is a coset, not the subgroup, we have to multiply the
+    # polynomial coefficients by x^i
+    interpolation_polynomial = fft(ys, MODULUS, root_of_unity, True)
+    interpolation_polynomial = [div(c, pow(x, i, MODULUS)) for i, c in enumerate(interpolation_polynomial)]
+
+    # Verify the pairing equation
+    #
+    # e([commitment - interpolation_polynomial(s)], [1]) = e([proof],  [s^n - x^n])
+    #    equivalent to
+    # e([commitment - interpolation_polynomial]^(-1), [1]) * e([proof],  [s^n - x^n]) = 1_T
+    #
+
+    xn_minus_yn = b.add(setup[1][n], b.multiply(b.neg(b.G2), pow(x, n, MODULUS)))
+    commitment_minus_interpolation = b.add(commitment, b.neg(lincomb(
+        setup[0][:len(interpolation_polynomial)], interpolation_polynomial, b.add, b.Z1)))
+    pairing_check = b.pairing(b.G2, b.neg(commitment_minus_interpolation), False)
+    pairing_check *= b.pairing(xn_minus_yn, proof, False)
+    pairing = b.final_exponentiate(pairing_check)
+    return pairing == b.FQ12.one()
+
+if __name__ == "__main__":
+    polynomial = [1, 2, 3, 4, 7, 7, 7, 7, 13, 13, 13, 13, 13, 13, 13, 13]
+    n = len(polynomial)
+
+    setup = generate_setup(1927409816240961209460912649124, n)
+
+    commitment = commit_to_poly(polynomial, setup)
+    proof = compute_proof_single(polynomial, 17, setup)
+    value = eval_poly_at(polynomial, 17)
+    assert check_proof_single(commitment, proof, 17, value, setup)
+    print("Single point check passed")
+
+    root_of_unity = get_root_of_unity(8)
+    x = 5431
+    coset = [x * pow(root_of_unity, i, MODULUS) for i in range(8)]
+    ys = [eval_poly_at(polynomial, z) for z in coset]
+    proof = compute_proof_multi(polynomial, x, 8, setup)
+    assert check_proof_multi(commitment, proof, x, ys, setup)
+    print("Coset check passed")
\ No newline at end of file
diff --git a/bindings/python/multicombs.py b/bindings/python/multicombs.py
new file mode 100644
index 0000000..972526c
--- /dev/null
+++ b/bindings/python/multicombs.py
@@ -0,0 +1,131 @@
+import random, time, sys, math
+
+# For each subset in `subsets` (provided as a list of indices into `numbers`),
+# compute the sum of that subset of `numbers`. More efficient than the naive method.
+def multisubset(numbers, subsets, adder=lambda x,y: x+y, zero=0):
+    numbers = numbers[::]
+    subsets = {i: {x for x in subset} for i, subset in enumerate(subsets)}
+    output = [zero for _ in range(len(subsets))]
+    
+    for roundcount in range(9999999):
+        # Compute counts of every pair of indices in the subset list
+        pair_count = {}
+        for index, subset in subsets.items():
+            for x in subset:
+                for y in subset:
+                    if y > x:
+                        pair_count[(x, y)] = pair_count.get((x, y), 0) + 1
+
+        # Determine pairs with highest count. The cutoff parameter [:len(numbers)]
+        # determines a tradeoff between group operation count and other forms of overhead
+        pairs_by_count = sorted([el for el in pair_count.keys()], key=lambda el: pair_count[el], reverse=True)[:len(numbers)*int(math.log(len(numbers)))]
+
+        # Exit condition: all subsets have size 1, no pairs
+        if len(pairs_by_count) == 0:
+            for key, subset in subsets.items():
+                for index in subset:
+                    output[key] = adder(output[key], numbers[index])
+            return output
+
+        # In each of the highest-count pairs, take the sum of the numbers at those indices,
+        # and add the result as a new value, and modify `subsets` to include the new value
+        # wherever possible
+        used = set()
+        for maxx, maxy in pairs_by_count:
+            if maxx in used or maxy in used:
+                continue
+            used.add(maxx)
+            used.add(maxy)
+            numbers.append(adder(numbers[maxx], numbers[maxy]))
+            for key, subset in list(subsets.items()):
+                if maxx in subset and maxy in subset:
+                    subset.remove(maxx)
+                    subset.remove(maxy)
+                    if not subset:
+                        output[key] = numbers[-1]
+                        del subsets[key]
+                    else:
+                        subset.add(len(numbers)-1)
+
+# Alternative algorithm. Less optimal than the above, but much lower bit twiddling
+# overhead and much simpler.
+def multisubset2(numbers, subsets, adder=lambda x,y: x+y, zero=0):
+    # Split up the numbers into partitions
+    partition_size = 1 + int(math.log(len(subsets) + 1))
+    # Align number count to partition size (for simplicity)
+    numbers = numbers[::]
+    while len(numbers) % partition_size != 0:
+        numbers.append(zero)
+    # Compute power set for each partition (eg. a, b, c -> {0, a, b, a+b, c, a+c, b+c, a+b+c})
+    power_sets = []
+    for i in range(0, len(numbers), partition_size):
+        new_power_set = [zero]
+        for dimension, value in enumerate(numbers[i:i+partition_size]):
+            new_power_set += [adder(n, value) for n in new_power_set]
+        power_sets.append(new_power_set)
+    # Compute subset sums, using elements from power set for each range of values
+    # ie. with a single power set lookup you can get the sum of _all_ elements in
+    # the range partition_size*k...partition_size*(k+1) that are in that subset
+    subset_sums = []
+    for subset in subsets:
+        o = zero
+        for i in range(len(power_sets)):
+            index_in_power_set = 0
+            for j in range(partition_size):
+                if i * partition_size + j in subset:
+                    index_in_power_set += 2 ** j
+            o = adder(o, power_sets[i][index_in_power_set])
+        subset_sums.append(o)
+    return subset_sums
+
+# Reduces a linear combination `numbers[0] * factors[0] + numbers[1] * factors[1] + ...`
+# into a multi-subset problem, and computes the result efficiently
+def lincomb(numbers, factors, adder=lambda x,y: x+y, zero=0):
+    # Maximum bit length of a number; how many subsets we need to make
+    maxbitlen = max((len(bin(f))-2 for f in factors), default=0)
+    # Compute the subsets: the ith subset contains the numbers whose corresponding factor
+    # has a 1 at the ith bit
+    subsets = [{i for i in range(len(numbers)) if factors[i] & (1 << j)} for j in range(maxbitlen+1)]
+    subset_sums = multisubset2(numbers, subsets, adder=adder, zero=zero)
+    # For example, suppose a value V has factor 6 (011 in increasing-order binary). Subset 0
+    # will not have V, subset 1 will, and subset 2 will. So if we multiply the output of adding
+    # subset 0 with twice the output of adding subset 1, with four times the output of adding
+    # subset 2, then V will be represented 0 + 2 + 4 = 6 times. This reasoning applies for every
+    # value. So `subset_0_sum + 2 * subset_1_sum + 4 * subset_2_sum` gives us the result we want.
+    # Here, we compute this as `((subset_2_sum * 2) + subset_1_sum) * 2 + subset_0_sum` for
+    # efficiency: an extra `maxbitlen * 2` group operations.
+    o = zero
+    for i in range(len(subsets)-1, -1, -1):
+        o = adder(adder(o, o), subset_sums[i])
+    return o
+
+# Tests go here
+def make_mock_adder():
+    counter = [0]
+    def adder(x, y):
+        if x and y:
+            counter[0] += 1
+        return x+y
+    return adder, counter
+
+def test_multisubset(numcount, setcount):
+    numbers = [random.randrange(10**20) for _ in range(numcount)]
+    subsets = [{i for i in range(numcount) if random.randrange(2)} for i in range(setcount)]
+    adder, counter = make_mock_adder()
+    o = multisubset(numbers, subsets, adder=adder)
+    for output, subset in zip(o, subsets):
+        assert output == sum([numbers[x] for x in subset])
+
+def test_lincomb(numcount, bitlength=256):
+    numbers = [random.randrange(10**20) for _ in range(numcount)]
+    factors = [random.randrange(2**bitlength) for _ in range(numcount)]
+    adder, counter = make_mock_adder()
+    o = lincomb(numbers, factors, adder=adder)
+    assert o == sum([n*f for n,f in zip(numbers, factors)])
+    total_ones = sum(bin(f).count('1') for f in factors)
+    print("Naive operation count: %d" % (bitlength * numcount + total_ones))
+    print("Optimized operation count: %d" % (bitlength * 2 + counter[0]))
+    print("Optimization factor: %.2f" % ((bitlength * numcount + total_ones) / (bitlength * 2 + counter[0])))
+
+if __name__ == '__main__':
+    test_lincomb(int(sys.argv[1]) if len(sys.argv) >= 2 else 80)
diff --git a/bindings/python/py_ecc_tests.py b/bindings/python/py_ecc_tests.py
new file mode 100644
index 0000000..a1a5660
--- /dev/null
+++ b/bindings/python/py_ecc_tests.py
@@ -0,0 +1,161 @@
+import atexit
+import ckzg
+import kzg_proofs
+
+def int_from_uint64s(digits):
+    """Convert a 4-tuple of base64 digits to the int it denotes"""
+    res, mult = 0, 1
+    for x in digits:
+        res += mult * x
+        mult *= 2 ** 64
+    return res
+
+def fr_from_int(x):
+    r = []
+    while x > 0:
+        r.append(x % 2**64)
+        x //= 2**64
+    assert len(r) <= 4
+    while len(r) < 4:
+        r.append(0)
+    return ckzg.fr_from_uint64s(tuple(r))
+
+def fr_to_int(fr):
+    return int_from_uint64s(ckzg.fr_to_uint64s(fr))
+
+import random
+
+polynomial = [random.randint(0, 0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001) for i in range(1024)]
+n = len(polynomial)
+
+secret_value = 1927409816240961209460912649124348975783457236447839525141982745761361378419
+
+cfa = ckzg.frArray(len(polynomial))
+for i, c in enumerate(polynomial):
+    cfa[i] = fr_from_int(c)
+
+# Build the polynomial
+
+ret, p = ckzg.new_poly_with_coeffs(cfa.cast(), len(polynomial))
+assert ret == ckzg.C_KZG_OK
+
+# Build a trusted setup with an arbitrary secret s
+# and max scale 4 (so 16 secret values)
+
+max_scale = n.bit_length() - 1
+
+ret, fs = ckzg.new_fft_settings(max_scale)
+assert ret == 0
+
+sa = ckzg.byteArray(32)
+for i, c in enumerate(secret_value.to_bytes(32, "little")):
+    sa[i] = c
+
+ret, secret_s = ckzg.blst_scalar_from_le_bytes(sa.cast(), 32)
+assert ret == 1
+
+test_secret = ckzg.fr_from_scalar(secret_s)
+assert fr_to_int(test_secret) == secret_value
+
+
+num_secrets = 2 ** max_scale
+
+g1s = ckzg.g1Array(num_secrets)
+g2s = ckzg.g2Array(num_secrets)
+
+ckzg.generate_trusted_setup(g1s.cast(), g2s.cast(), secret_s, num_secrets)
+
+ret, ks = ckzg.new_kzg_settings(g1s.cast(), g2s.cast(), num_secrets, fs)
+assert ret == 0
+
+# Compute the Lagrange form of our polynomial in this setup
+
+ret, p_l = ckzg.new_poly_l_from_poly(p, ks)
+assert ret == 0
+
+x = 832877253762587406983796272890571809375809175809315245
+fr_x = fr_from_int(x)
+assert fr_to_int(fr_x) == x
+
+ret, y_l = ckzg.eval_poly_l(p_l, fr_x, fs)
+assert ret == 0
+
+y = ckzg.eval_poly(p, fr_x)
+assert ckzg.fr_equal(y, y_l)
+
+y_p = kzg_proofs.eval_poly_at(polynomial, x)
+assert fr_to_int(y) == y_p
+
+# Commit to the polynomial, in both Lagrange and coefficient form
+# The commitment should be the same
+
+ret, commitment = ckzg.commit_to_poly(p, ks)
+assert ret == 0
+ret, commitment_l = ckzg.commit_to_poly_l(p_l, ks)
+assert ret == 0
+assert ckzg.g1_equal(commitment, commitment_l)
+
+# Compute proof at an arbitrary point (for both forms)
+
+ret, π = ckzg.compute_proof_single(p, fr_x, ks)
+assert ret == 0
+ret, v = ckzg.eval_poly_l(p_l, fr_x, fs)
+assert ret == 0
+ret, π_l = ckzg.compute_proof_single_l(p_l, fr_x, v, ks)
+assert ret == 0
+
+# Compute proof using py_ecc
+
+pyecc_setup = kzg_proofs.generate_setup(secret_value, n)
+
+pyecc_commitment = kzg_proofs.commit_to_poly(polynomial, pyecc_setup)
+pyecc_proof = kzg_proofs.compute_proof_single(polynomial, x, pyecc_setup)
+assert kzg_proofs.check_proof_single(pyecc_commitment, pyecc_proof, x, y_p, pyecc_setup)
+
+from py_ecc.bls.point_compression import compress_G1
+
+pyecc_commitment_compressed = compress_G1(pyecc_commitment)
+pyecc_proof_compressed = compress_G1(pyecc_proof)
+
+commitment_compressed = ckzg.byteArray(48)
+
+ckzg.blst_p1_compress(commitment_compressed.cast(), commitment)
+
+commitment_compressed = bytes([commitment_compressed[i] for i in range(48)])
+assert commitment_compressed == pyecc_commitment_compressed.to_bytes(48, "big")
+
+
+proof_compressed = ckzg.byteArray(48)
+
+ckzg.blst_p1_compress(proof_compressed.cast(), π)
+
+proof_compressed = bytes([proof_compressed[i] for i in range(48)])
+assert proof_compressed == pyecc_proof_compressed.to_bytes(48, "big")
+
+# Check the proofs using the commitments
+
+ret, res = ckzg.check_proof_single(commitment, π, fr_x, v, ks)
+assert ret == 0
+assert res
+
+ret, res = ckzg.check_proof_single(commitment_l, π_l, fr_x, v, ks)
+assert ret == 0
+assert res
+
+# Check the proof fails with the wrong value
+
+w = ckzg.fr_add(v, ckzg.fr_one)
+ret, res = ckzg.check_proof_single(commitment_l, π_l, fr_x, w, ks)
+assert ret == 0
+assert not res
+
+print("All tests passed.")
+
+# We need to manually free the C allocated arrays
+# Use atexit so this file can be loaded interactively before freeing
+def cleanup():
+    ckzg.free_poly(p)
+    ckzg.free_poly_l(p_l)
+    ckzg.free_fft_settings(fs)
+    ckzg.free_kzg_settings(ks)
+atexit.register(cleanup)
diff --git a/bindings/python/requirements.txt b/bindings/python/requirements.txt
new file mode 100644
index 0000000..9565b69
--- /dev/null
+++ b/bindings/python/requirements.txt
@@ -0,0 +1 @@
+py_ecc
\ No newline at end of file