1package prover
  2
  3import (
  4	"bytes"
  5	"math/big"
  6)
  7
  8func arrayOfZeroes(n int) []*big.Int {
  9	var r []*big.Int
 10	for i := 0; i < n; i++ {
 11		r = append(r, new(big.Int).SetInt64(0))
 12	}
 13	return r
 14}
 15
 16func fAdd(a, b *big.Int) *big.Int {
 17	ab := new(big.Int).Add(a, b)
 18	return new(big.Int).Mod(ab, R)
 19}
 20
 21func fSub(a, b *big.Int) *big.Int {
 22	ab := new(big.Int).Sub(a, b)
 23	return new(big.Int).Mod(ab, R)
 24}
 25
 26func fMul(a, b *big.Int) *big.Int {
 27	ab := new(big.Int).Mul(a, b)
 28	return new(big.Int).Mod(ab, R)
 29}
 30
 31func fDiv(a, b *big.Int) *big.Int {
 32	ab := new(big.Int).Mul(a, new(big.Int).ModInverse(b, R))
 33	return new(big.Int).Mod(ab, R)
 34}
 35
 36func fNeg(a *big.Int) *big.Int {
 37	return new(big.Int).Mod(new(big.Int).Neg(a), R)
 38}
 39
 40func fInv(a *big.Int) *big.Int {
 41	return new(big.Int).ModInverse(a, R)
 42}
 43
 44func fExp(base *big.Int, e *big.Int) *big.Int {
 45	res := big.NewInt(1)
 46	rem := new(big.Int).Set(e)
 47	exp := base
 48
 49	for !bytes.Equal(rem.Bytes(), big.NewInt(int64(0)).Bytes()) {
 50		// if BigIsOdd(rem) {
 51		if rem.Bit(0) == 1 { // .Bit(0) returns 1 when is odd
 52			res = fMul(res, exp)
 53		}
 54		exp = fMul(exp, exp)
 55		rem = new(big.Int).Rsh(rem, 1)
 56	}
 57	return res
 58}
 59
 60func max(a, b int) int {
 61	if a > b {
 62		return a
 63	}
 64	return b
 65}
 66
 67func polynomialSub(a, b []*big.Int) []*big.Int {
 68	r := arrayOfZeroes(max(len(a), len(b)))
 69	for i := 0; i < len(a); i++ {
 70		r[i] = fAdd(r[i], a[i])
 71	}
 72	for i := 0; i < len(b); i++ {
 73		r[i] = fSub(r[i], b[i])
 74	}
 75	return r
 76}
 77
 78func polynomialMul(a, b []*big.Int) []*big.Int {
 79	r := arrayOfZeroes(len(a) + len(b) - 1)
 80	for i := 0; i < len(a); i++ {
 81		for j := 0; j < len(b); j++ {
 82			r[i+j] = fAdd(r[i+j], fMul(a[i], b[j]))
 83		}
 84	}
 85	return r
 86}
 87
 88func polynomialDiv(a, b []*big.Int) ([]*big.Int, []*big.Int) {
 89	// https://en.wikipedia.org/wiki/Division_algorithm
 90	r := arrayOfZeroes(len(a) - len(b) + 1)
 91	rem := a
 92	for len(rem) >= len(b) {
 93		l := fDiv(rem[len(rem)-1], b[len(b)-1])
 94		pos := len(rem) - len(b)
 95		r[pos] = l
 96		aux := arrayOfZeroes(pos)
 97		aux1 := append(aux, l)
 98		aux2 := polynomialSub(rem, polynomialMul(b, aux1))
 99		rem = aux2[:len(aux2)-1]
100	}
101	return r, rem
102}