1package r1csqapFloat
  2
  3import "math/big"
  4
  5func Transpose(matrix [][]*big.Float) [][]*big.Float {
  6	var r [][]*big.Float
  7	for i := 0; i < len(matrix[0]); i++ {
  8		var row []*big.Float
  9		for j := 0; j < len(matrix); j++ {
 10			row = append(row, matrix[j][i])
 11		}
 12		r = append(r, row)
 13	}
 14	return r
 15}
 16
 17func ArrayOfBigZeros(num int) []*big.Float {
 18	bigZero := big.NewFloat(float64(0))
 19	var r []*big.Float
 20	for i := 0; i < num; i++ {
 21		r = append(r, bigZero)
 22	}
 23	return r
 24}
 25
 26func PolMul(a, b []*big.Float) []*big.Float {
 27	r := ArrayOfBigZeros(len(a) + len(b) - 1)
 28	for i := 0; i < len(a); i++ {
 29		for j := 0; j < len(b); j++ {
 30			r[i+j] = new(big.Float).Add(
 31				r[i+j],
 32				new(big.Float).Mul(a[i], b[j]))
 33		}
 34	}
 35	return r
 36}
 37
 38func PolDiv(a, b []*big.Float) ([]*big.Float, []*big.Float) {
 39	// https://en.wikipedia.org/wiki/Division_algorithm
 40	r := ArrayOfBigZeros(len(a) - len(b) + 1)
 41	rem := a
 42	for len(rem) >= len(b) {
 43		l := new(big.Float).Quo(rem[len(rem)-1], b[len(b)-1])
 44		pos := len(rem) - len(b)
 45		r[pos] = l
 46		aux := ArrayOfBigZeros(pos)
 47		aux1 := append(aux, l)
 48		aux2 := PolSub(rem, PolMul(b, aux1))
 49		rem = aux2[:len(aux2)-1]
 50	}
 51	return r, rem
 52}
 53
 54func max(a, b int) int {
 55	if a > b {
 56		return a
 57	}
 58	return b
 59}
 60
 61func PolAdd(a, b []*big.Float) []*big.Float {
 62	r := ArrayOfBigZeros(max(len(a), len(b)))
 63	for i := 0; i < len(a); i++ {
 64		r[i] = new(big.Float).Add(r[i], a[i])
 65	}
 66	for i := 0; i < len(b); i++ {
 67		r[i] = new(big.Float).Add(r[i], b[i])
 68	}
 69	return r
 70}
 71
 72func PolSub(a, b []*big.Float) []*big.Float {
 73	r := ArrayOfBigZeros(max(len(a), len(b)))
 74	for i := 0; i < len(a); i++ {
 75		r[i] = new(big.Float).Add(r[i], a[i])
 76	}
 77	for i := 0; i < len(b); i++ {
 78		bneg := new(big.Float).Mul(b[i], big.NewFloat(float64(-1)))
 79		r[i] = new(big.Float).Add(r[i], bneg)
 80	}
 81	return r
 82
 83}
 84
 85func FloatPow(a *big.Float, e int) *big.Float {
 86	if e == 0 {
 87		return big.NewFloat(float64(1))
 88	}
 89	result := new(big.Float).Copy(a)
 90	for i := 0; i < e-1; i++ {
 91		result = new(big.Float).Mul(result, a)
 92	}
 93	return result
 94}
 95
 96func PolEval(v []*big.Float, x *big.Float) *big.Float {
 97	r := big.NewFloat(float64(0))
 98	for i := 0; i < len(v); i++ {
 99		xi := FloatPow(x, i)
100		elem := new(big.Float).Mul(v[i], xi)
101		r = new(big.Float).Add(r, elem)
102	}
103	return r
104}
105
106func NewPolZeroAt(pointPos, totalPoints int, height *big.Float) []*big.Float {
107	fac := 1
108	for i := 1; i < totalPoints+1; i++ {
109		if i != pointPos {
110			fac = fac * (pointPos - i)
111		}
112	}
113	facBig := big.NewFloat(float64(fac))
114	hf := new(big.Float).Quo(height, facBig)
115	r := []*big.Float{hf}
116	for i := 1; i < totalPoints+1; i++ {
117		if i != pointPos {
118			ineg := big.NewFloat(float64(-i))
119			b1 := big.NewFloat(float64(1))
120			r = PolMul(r, []*big.Float{ineg, b1})
121		}
122	}
123	return r
124}
125
126func LagrangeInterpolation(v []*big.Float) []*big.Float {
127	// https://en.wikipedia.org/wiki/Lagrange_polynomial
128	var r []*big.Float
129	for i := 0; i < len(v); i++ {
130		r = PolAdd(r, NewPolZeroAt(i+1, len(v), v[i]))
131	}
132	//
133	return r
134}
135
136func R1CSToQAP(a, b, c [][]*big.Float) ([][]*big.Float, [][]*big.Float, [][]*big.Float, []*big.Float) {
137	aT := Transpose(a)
138	bT := Transpose(b)
139	cT := Transpose(c)
140	var alphas [][]*big.Float
141	for i := 0; i < len(aT); i++ {
142		alphas = append(alphas, LagrangeInterpolation(aT[i]))
143	}
144	var betas [][]*big.Float
145	for i := 0; i < len(bT); i++ {
146		betas = append(betas, LagrangeInterpolation(bT[i]))
147	}
148	var gammas [][]*big.Float
149	for i := 0; i < len(cT); i++ {
150		gammas = append(gammas, LagrangeInterpolation(cT[i]))
151	}
152	z := []*big.Float{big.NewFloat(float64(1))}
153	for i := 1; i < len(aT[0])+1; i++ {
154		ineg := big.NewFloat(float64(-i))
155		b1 := big.NewFloat(float64(1))
156		z = PolMul(z, []*big.Float{ineg, b1})
157	}
158	return alphas, betas, gammas, z
159}
160
161func CombinePolynomials(r []*big.Float, ap, bp, cp [][]*big.Float) ([]*big.Float, []*big.Float, []*big.Float, []*big.Float) {
162	var alpha []*big.Float
163	for i := 0; i < len(r); i++ {
164		m := PolMul([]*big.Float{r[i]}, ap[i])
165		alpha = PolAdd(alpha, m)
166	}
167	var beta []*big.Float
168	for i := 0; i < len(r); i++ {
169		m := PolMul([]*big.Float{r[i]}, bp[i])
170		beta = PolAdd(beta, m)
171	}
172	var gamma []*big.Float
173	for i := 0; i < len(r); i++ {
174		m := PolMul([]*big.Float{r[i]}, cp[i])
175		gamma = PolAdd(gamma, m)
176	}
177
178	px := PolSub(PolMul(alpha, beta), gamma)
179	return alpha, beta, gamma, px
180}
181
182func DivisorPolinomial(px, z []*big.Float) []*big.Float {
183	quo, _ := PolDiv(px, z)
184	return quo
185}