1package sn
  2
  3import (
  4	"math/big"
  5)
  6
  7func Transpose(matrix [][]*big.Float) [][]*big.Float {
  8	var r [][]*big.Float
  9	for i := 0; i < len(matrix[0]); i++ {
 10		var row []*big.Float
 11		for j := 0; j < len(matrix); j++ {
 12			row = append(row, matrix[j][i])
 13		}
 14		r = append(r, row)
 15	}
 16	return r
 17}
 18
 19func ArrayOfBigZeros(num int) []*big.Float {
 20	bigZero := big.NewFloat(float64(0))
 21	var r []*big.Float
 22	for i := 0; i < num; i++ {
 23		r = append(r, bigZero)
 24	}
 25	return r
 26}
 27
 28func PolMul(a, b []*big.Float) []*big.Float {
 29	r := ArrayOfBigZeros(len(a) + len(b) - 1)
 30	for i := 0; i < len(a); i++ {
 31		for j := 0; j < len(b); j++ {
 32			r[i+j] = new(big.Float).Add(
 33				r[i+j],
 34				new(big.Float).Mul(a[i], b[j]))
 35		}
 36	}
 37	return r
 38}
 39
 40func max(a, b int) int {
 41	if a > b {
 42		return a
 43	}
 44	return b
 45}
 46
 47func PolAdd(a, b []*big.Float) []*big.Float {
 48	r := ArrayOfBigZeros(max(len(a), len(b)))
 49	for i := 0; i < len(a); i++ {
 50		r[i] = new(big.Float).Add(r[i], a[i])
 51	}
 52	for i := 0; i < len(b); i++ {
 53		r[i] = new(big.Float).Add(r[i], b[i])
 54	}
 55	return r
 56}
 57
 58func PolSub(a, b []*big.Float) []*big.Float {
 59	r := ArrayOfBigZeros(max(len(a), len(b)))
 60	for i := 0; i < len(a); i++ {
 61		r[i] = new(big.Float).Add(r[i], a[i])
 62	}
 63	for i := 0; i < len(b); i++ {
 64		bneg := new(big.Float).Mul(b[i], big.NewFloat(float64(-1)))
 65		r[i] = new(big.Float).Add(r[i], bneg)
 66	}
 67	return r
 68
 69}
 70
 71func FloatPow(a *big.Float, e int) *big.Float {
 72	if e == 0 {
 73		return big.NewFloat(float64(1))
 74	}
 75	result := new(big.Float).Copy(a)
 76	for i := 0; i < e-1; i++ {
 77		result = new(big.Float).Mul(result, a)
 78	}
 79	return result
 80}
 81
 82func PolEval(v []*big.Float, x *big.Float) *big.Float {
 83	r := big.NewFloat(float64(0))
 84	for i := 0; i < len(v); i++ {
 85		xi := FloatPow(x, i)
 86		elem := new(big.Float).Mul(v[i], xi)
 87		r = new(big.Float).Add(r, elem)
 88	}
 89	return r
 90}
 91
 92func NewPolZeroAt(pointPos, totalPoints int, height *big.Float) []*big.Float {
 93	fac := 1
 94	for i := 1; i < totalPoints+1; i++ {
 95		if i != pointPos {
 96			fac = fac * (pointPos - i)
 97		}
 98	}
 99	facBig := big.NewFloat(float64(fac))
100	hf := new(big.Float).Quo(height, facBig)
101	r := []*big.Float{hf}
102	for i := 1; i < totalPoints+1; i++ {
103		if i != pointPos {
104			ineg := big.NewFloat(float64(-i))
105			b1 := big.NewFloat(float64(1))
106			r = PolMul(r, []*big.Float{ineg, b1})
107		}
108	}
109	return r
110}
111
112func LagrangeInterpolation(v []*big.Float) []*big.Float {
113	// https://en.wikipedia.org/wiki/Lagrange_polynomial
114	var r []*big.Float
115	for i := 0; i < len(v); i++ {
116		r = PolAdd(r, NewPolZeroAt(i+1, len(v), v[i]))
117	}
118	//
119	return r
120}
121
122func R1CSToQAP(a, b, c [][]*big.Float) ([][]*big.Float, [][]*big.Float, [][]*big.Float, []*big.Float) {
123	aT := Transpose(a)
124	bT := Transpose(b)
125	cT := Transpose(c)
126	var alpha [][]*big.Float
127	for i := 0; i < len(aT); i++ {
128		alpha = append(alpha, LagrangeInterpolation(aT[i]))
129	}
130	var beta [][]*big.Float
131	for i := 0; i < len(bT); i++ {
132		beta = append(beta, LagrangeInterpolation(bT[i]))
133	}
134	var gamma [][]*big.Float
135	for i := 0; i < len(cT); i++ {
136		gamma = append(gamma, LagrangeInterpolation(cT[i]))
137	}
138	z := []*big.Float{big.NewFloat(float64(1))}
139	for i := 1; i < len(aT[0])+1; i++ {
140		ineg := big.NewFloat(float64(-i))
141		b1 := big.NewFloat(float64(1))
142		z = PolMul(z, []*big.Float{ineg, b1})
143	}
144	return alpha, beta, gamma, z
145}