[research] add sage implementation of R1CS->QAP (#190)

Add sage code that runs through the Groth16 R1CS --> QAP conversion
contained in the article
"R1CS and QAP - From Zero to Hero with Finite Fields & sagemath"
https://risencrypto.github.io/R1CSQAP/

Co-authored-by: y <y>

script/research/zk/groth16/qap.sage

@@ -0,0 +1,117 @@
+#!/usr/bin/env sage
+
+"""
+The following code is taken from the examples outlined in the article:
+
+"R1CS and QAP - From Zero to Hero with Finite Fields & sagemath".
+(https://risencrypto.github.io/R1CSQAP/). The Sage examples begin at the 
+end of the section called "Conversion to Arithemetic Circuits & then R1CS".
+
+Note that this example uses the same polynomial as the file qap.py
+in this directory. 
+"""
+
+header = """
+This script demonstrates the conversion of a set of R1CS matrices to QAP form.
+For more information, read the code comments in the file for this script.
+"""
+print(header)
+
+# Define the R1CS matrices over the finite field GF(41)
+F41 = GF(41)
+L = Matrix(F41, [
+[0,1,0,0,0,0],
+[0,0,0,1,0,0],
+[0,1,0,0,1,0],
+[5,0,0,0,0,1],]
+)
+R = Matrix(F41, [
+[0,1,0,0,0,0],
+[0,1,0,0,0,0],
+[1,0,0,0,0,0],
+[1,0,0,0,0,0,],])
+O = Matrix(F41, [
+[0,0,0,1,0,0],
+[0,0,0,0,1,0],
+[0,0,0,0,0,1],
+[0,0,1,0,0,0,],])
+
+# Print R1CS matrices
+print(L)
+print()
+print(R)
+print()
+print(O)
+print()
+
+# Debug: this section should print 35*x^3 + 36*x^2 + 16*x + 36
+# It represents the Lagrange interpolation for the first column of matrix L.
+R41.<x> = PolynomialRing(F41)
+points = [(1,0),(2,0),(3,0),(4,5)]
+R41.lagrange_polynomial(points)
+
+
+# Convert matrices L, R, and O to QAP form using Lagrange interpolation, i.e.
+# repeat the above lines of code for every column of every matrix.
+M = [L, R, O]
+PolyM = []
+for m in M:
+	PolyList = []
+	for i in range(m.ncols()):
+		points = []
+		for j in range(m.nrows()):
+			points.append([j+1,m[j,i]])
+		Poly = R41.lagrange_polynomial(points).coefficients(sparse=False)
+
+		if(len(Poly) < m.nrows()):
+			# if degree of the polynomial is less than 4
+			# we add zeroes to represent the missed out terms
+			dif = m.nrows() - len(Poly)
+			for c in range(dif):
+				Poly.append(0);
+		PolyList.append(Poly)
+	PolyM.append(Matrix(F41, PolyList))
+
+# Print the new matrices containing QAP form polynomials.
+# Each matrix should have 6 polynomials. This corresponds to the number of elements
+# in the solution vector (which is 6 using the example from the article).
+print(PolyM[0])
+print()
+print(PolyM[1])
+print()
+print(PolyM[2])
+print()
+
+# Solution vector (defined over the finite field). This is known only to the prover.
+S = vector(F41,[1, 3, 35, 9, 27, 30]) 
+# Create the L, Rx & Ox polynomials. Perform the dot product of the solution vector
+# with each of the matrices obtained in the previous step.
+Lx = R41(list(S*PolyM[0]))
+Rx = R41(list(S*PolyM[1]))
+Ox = R41(list(S*PolyM[2]))
+print("Lx = " + str(Lx))
+print("Rx = " + str(Rx))
+print("Ox = " + str(Ox))
+
+# Multiply polynomials Lx and Rx and subtract Ox. This is how we check the
+# constraints of the circuit and corresponds to checking the vectors in
+# the R1CS step.
+# Note: the verifier does not know the polynomial T.
+T = Lx*Rx - Ox
+print("T(x) = ", end="")
+print(T)
+
+# Check Tx at x = 1, 2, 3, 4. The output of each should be 0.
+print("T(1) = " + str(T(1)))
+print("T(2) = " + str(T(2)))
+print("T(3) = " + str(T(3)))
+print("T(4) = " + str(T(4)))
+
+# This polynomial is known to both the prover and verifier
+Z = R41((x-1)*(x-2)*(x-3)*(x-4))
+H = T.quo_rem(Z)
+print("Quotient of Z/T = ", end="")
+print(H[0])
+# The remainder here must be 0, indicating that Z exactly divides T.
+print("Remainder of Z/T = ", end="")
+print(H[1])