<gh_stars>0
import numpy as np

# R1CS for a polynomial using the Horner algoirhtm
# QAP using the Newton interpolation

def R1CS(poly):
    d = len(poly) - 1
    A = np.zeros((d, d+2))
    B = np.zeros((d, d+2))
    C = np.diag(np.ones(d), k=2)
    A[0,:2] = poly[:-3:-1]
    A[1:,0] = poly[-3::-1]
    A[1:,2:-1] = np.diag(np.ones(d-1))
    B[:,0] = 1
    C = C[:-2,:]
    return A, B, C

def _newtonInterpolation(x, y):
    m = len(x)
    x = np.copy(x)
    a = np.copy(y)
    for k in range(1,m):
        a[k:m] = (a[k:m] - a[k - 1])/(x[k:m] - x[k - 1])
    return a

def QAP(A, B, C):
    n, m = A.shape
    x = list(range(1,n+1))
    PA = np.zeros((n,m))
    PB = np.zeros((n,m))
    PC = np.zeros((n,m))
    for M, PM in zip([A, B, C], [PA, PB, PC]):
        for i in range(m):
            PM[:,i] = _newtonInterpolation(x, M[:,i])
    return PA, PB, PC

def _polymul(p, q):
    m, n = len(p), len(q)
    r = np.zeros(m + n - 1)
    for i in range(m + n):
        for j in range(i+1):
            if j < m and i < n + j:
                r[i] += p[j] * q[i-j]
    return r

def getProof(PA, PB, PC, s):
    p = _polymul(PA @ s, PB @ s) + (PC @ s)
    return p, s
    
poly = [5,1,0,1]
arrays = R1CS(poly)
for arr in arrays:
    print(arr)

polys = QAP(*arrays)
for p in polys:
    print(p)

x = [0, 1]
print(_polymul(x, poly))
