polynomial Division, SolPolynomials, DivisorPolynomial

r1csqap/r1csqap.go

@@ -47,6 +47,21 @@ func (pf PolynomialField) Mul(a, b []*big.Int) []*big.Int {
 	}
 	return r
 }
+func (pf PolynomialField) Div(a, b []*big.Int) ([]*big.Int, []*big.Int) {
+	// https://en.wikipedia.org/wiki/Division_algorithm
+	r := ArrayOfBigZeros(len(a) - len(b) + 1)
+	rem := a
+	for len(rem) >= len(b) {
+		l := pf.F.Div(rem[len(rem)-1], b[len(b)-1])
+		pos := len(rem) - len(b)
+		r[pos] = l
+		aux := ArrayOfBigZeros(pos)
+		aux1 := append(aux, l)
+		aux2 := pf.Sub(rem, pf.Mul(b, aux1))
+		rem = aux2[:len(aux2)-1]
+	}
+	return r, rem
+}
 
 func max(a, b int) int {
 	if a > b {
@@ -72,24 +87,11 @@ func (pf PolynomialField) Sub(a, b []*big.Int) []*big.Int {
 		r[i] = pf.F.Add(r[i], a[i])
 	}
 	for i := 0; i < len(b); i++ {
-		// bneg := pf.F.Mul(b[i], big.NewInt(int64(-1)))
-		// r[i] = pf.F.Add(r[i], bneg)
 		r[i] = pf.F.Sub(r[i], b[i])
 	}
 	return r
 }
 
-// func FloatPow(a *big.Int, e int) *big.Int {
-//         if e == 0 {
-//                 return big.NewInt(int64(1))
-//         }
-//         result := new(big.Int).Copy(a)
-//         for i := 0; i < e-1; i++ {
-//                 result = new(big.Int).Mul(result, a)
-//         }
-//         return result
-// }
-
 func (pf PolynomialField) Eval(v []*big.Int, x *big.Int) *big.Int {
 	r := big.NewInt(int64(0))
 	for i := 0; i < len(v); i++ {
@@ -155,3 +157,29 @@ func (pf PolynomialField) R1CSToQAP(a, b, c [][]*big.Int) ([][]*big.Int, [][]*bi
 	}
 	return alpha, beta, gamma, z
 }
+
+func (pf PolynomialField) SolPolynomials(r []*big.Int, ap, bp, cp [][]*big.Int) ([]*big.Int, []*big.Int, []*big.Int, []*big.Int) {
+	var alpha []*big.Int
+	for i := 0; i < len(r); i++ {
+		m := pf.Mul([]*big.Int{r[i]}, ap[i])
+		alpha = pf.Add(alpha, m)
+	}
+	var beta []*big.Int
+	for i := 0; i < len(r); i++ {
+		m := pf.Mul([]*big.Int{r[i]}, bp[i])
+		beta = pf.Add(beta, m)
+	}
+	var gamma []*big.Int
+	for i := 0; i < len(r); i++ {
+		m := pf.Mul([]*big.Int{r[i]}, cp[i])
+		gamma = pf.Add(gamma, m)
+	}
+
+	px := pf.Sub(pf.Mul(alpha, beta), gamma)
+	return alpha, beta, gamma, px
+}
+
+func (pf PolynomialField) DivisorPolinomial(px, z []*big.Int) []*big.Int {
+	quo, _ := pf.Div(px, z)
+	return quo
+}

r1csqap/r1csqap_test.go

@@ -108,7 +108,12 @@ func TestR1CSToQAP(t *testing.T) {
 
 	b0 := big.NewInt(int64(0))
 	b1 := big.NewInt(int64(1))
+	b3 := big.NewInt(int64(3))
 	b5 := big.NewInt(int64(5))
+	b9 := big.NewInt(int64(9))
+	b27 := big.NewInt(int64(27))
+	b30 := big.NewInt(int64(30))
+	b35 := big.NewInt(int64(35))
 	a := [][]*big.Int{
 		[]*big.Int{b0, b1, b0, b0, b0, b0},
 		[]*big.Int{b0, b0, b0, b1, b0, b0},
@@ -127,10 +132,27 @@ func TestR1CSToQAP(t *testing.T) {
 		[]*big.Int{b0, b0, b0, b0, b0, b1},
 		[]*big.Int{b0, b0, b1, b0, b0, b0},
 	}
-	alpha, beta, gamma, z := pf.R1CSToQAP(a, b, c)
+	ap, bp, cp, z := pf.R1CSToQAP(a, b, c)
+	fmt.Println(ap)
+	fmt.Println(bp)
+	fmt.Println(cp)
+	fmt.Println(z)
+
+	w := []*big.Int{b1, b3, b35, b9, b27, b30}
+	alpha, beta, gamma, px := pf.SolPolynomials(w, ap, bp, cp)
 	fmt.Println(alpha)
 	fmt.Println(beta)
 	fmt.Println(gamma)
-	fmt.Println(z)
+	fmt.Println(px)
+
+	h := pf.DivisorPolinomial(px, z)
+	fmt.Println(h)
+
+	// h==px/z so px==h*z
+	assert.Equal(t, px, pf.Mul(h, z))
 
+	// a(x) * b(x) - c(x) == h * z(x)
+	abc := pf.Sub(pf.Mul(alpha, beta), gamma)
+	hz := pf.Mul(h, z)
+	assert.Equal(t, abc, hz)
 }

r1csqapFloat/r1csqapFloat.go

@@ -1,8 +1,6 @@
 package r1csqapFloat
 
-import (
-	"math/big"
-)
+import "math/big"
 
 func Transpose(matrix [][]*big.Float) [][]*big.Float {
 	var r [][]*big.Float
@@ -37,6 +35,22 @@ func PolMul(a, b []*big.Float) []*big.Float {
 	return r
 }
 
+func PolDiv(a, b []*big.Float) ([]*big.Float, []*big.Float) {
+	// https://en.wikipedia.org/wiki/Division_algorithm
+	r := ArrayOfBigZeros(len(a) - len(b) + 1)
+	rem := a
+	for len(rem) >= len(b) {
+		l := new(big.Float).Quo(rem[len(rem)-1], b[len(b)-1])
+		pos := len(rem) - len(b)
+		r[pos] = l
+		aux := ArrayOfBigZeros(pos)
+		aux1 := append(aux, l)
+		aux2 := PolSub(rem, PolMul(b, aux1))
+		rem = aux2[:len(aux2)-1]
+	}
+	return r, rem
+}
+
 func max(a, b int) int {
 	if a > b {
 		return a
@@ -143,3 +157,29 @@ func R1CSToQAP(a, b, c [][]*big.Float) ([][]*big.Float, [][]*big.Float, [][]*big
 	}
 	return alpha, beta, gamma, z
 }
+
+func SolPolynomials(r []*big.Float, ap, bp, cp [][]*big.Float) ([]*big.Float, []*big.Float, []*big.Float, []*big.Float) {
+	var alpha []*big.Float
+	for i := 0; i < len(r); i++ {
+		m := PolMul([]*big.Float{r[i]}, ap[i])
+		alpha = PolAdd(alpha, m)
+	}
+	var beta []*big.Float
+	for i := 0; i < len(r); i++ {
+		m := PolMul([]*big.Float{r[i]}, bp[i])
+		beta = PolAdd(beta, m)
+	}
+	var gamma []*big.Float
+	for i := 0; i < len(r); i++ {
+		m := PolMul([]*big.Float{r[i]}, cp[i])
+		gamma = PolAdd(gamma, m)
+	}
+
+	px := PolSub(PolMul(alpha, beta), gamma)
+	return alpha, beta, gamma, px
+}
+
+func DivisorPolinomial(px, z []*big.Float) []*big.Float {
+	quo, _ := PolDiv(px, z)
+	return quo
+}

r1csqapFloat/r1csqapFloat_test.go

@@ -84,7 +84,12 @@ func TestLagrangeInterpolation(t *testing.T) {
 func TestR1CSToQAP(t *testing.T) {
 	b0 := big.NewFloat(float64(0))
 	b1 := big.NewFloat(float64(1))
+	b3 := big.NewFloat(float64(3))
 	b5 := big.NewFloat(float64(5))
+	b9 := big.NewFloat(float64(9))
+	b27 := big.NewFloat(float64(27))
+	b30 := big.NewFloat(float64(30))
+	b35 := big.NewFloat(float64(35))
 	a := [][]*big.Float{
 		[]*big.Float{b0, b1, b0, b0, b0, b0},
 		[]*big.Float{b0, b0, b0, b1, b0, b0},
@@ -103,10 +108,21 @@ func TestR1CSToQAP(t *testing.T) {
 		[]*big.Float{b0, b0, b0, b0, b0, b1},
 		[]*big.Float{b0, b0, b1, b0, b0, b0},
 	}
-	alpha, beta, gamma, z := R1CSToQAP(a, b, c)
+	ap, bp, cp, z := R1CSToQAP(a, b, c)
+	fmt.Println(ap)
+	fmt.Println(bp)
+	fmt.Println(cp)
+	fmt.Println(z)
+	zexpected := []*big.Float{big.NewFloat(float64(24)), big.NewFloat(float64(-50)), big.NewFloat(float64(35)), big.NewFloat(float64(-10)), big.NewFloat(float64(1))}
+	assert.Equal(t, z, zexpected)
+
+	w := []*big.Float{b1, b3, b35, b9, b27, b30}
+	alpha, beta, gamma, px := SolPolynomials(w, ap, bp, cp)
 	fmt.Println(alpha)
 	fmt.Println(beta)
 	fmt.Println(gamma)
-	fmt.Println(z)
+	fmt.Println(px)
 
+	h := DivisorPolinomial(px, z)
+	fmt.Println(h)
 }