Succinct Non-Interactive Zero Knowledge for a von Neumann Architecture, Eli Ben-Sasson, Alessandro Chiesa, Eran Tromer, Madars Virza https://eprint.iacr.org/2013/879.pdf
Pinocchio: Nearly practical verifiable computation, Bryan Parno, Craig Gentry, Jon Howell, Mariana Raykova https://eprint.iacr.org/2013/279.pdf

Caution

Implementation from scratch in Go to understand the concepts. Do not use in production.

Not finished, implementing this in my free time to understand it better, so I don't have much time.

Current implementation status:

Finite Fields (1, 2, 6, 12) operations
G1 and G2 curve operations
BN128 Pairing
circuit code compiler
- code to flat code
- flat code compiler
circuit to R1CS
polynomial operations
R1CS to QAP
generate trusted setup
generate proofs
verify proofs with BN128 pairing

Usage

zkSnark
bn128 (more details: https://github.com/arnaucube/go-snark/tree/master/bn128)
Finite Fields operations
R1CS to QAP (more details: https://github.com/arnaucube/go-snark/tree/master/r1csqap)
Circuit Compiler

Example:

// compile circuit and get the R1CS
flatCode := `
func test(x):
	aux = x*x
	y = aux*x
	z = x + y
	out = z + 5
`

// parse the code
parser := circuitcompiler.NewParser(strings.NewReader(flatCode))
circuit, err := parser.Parse()
assert.Nil(t, err)
fmt.Println(circuit)

// witness
b3 := big.NewInt(int64(3))
inputs := []*big.Int{b3}
w := circuit.CalculateWitness(inputs)
fmt.Println("\nwitness", w)
/*
now we have the witness:
w = [1 3 35 9 27 30]
*/

// flat code to R1CS
fmt.Println("generating R1CS from flat code")
a, b, c := circuit.GenerateR1CS()

/*
now we have the R1CS from the circuit:
a == [[0 1 0 0 0 0] [0 0 0 1 0 0] [0 1 0 0 1 0] [5 0 0 0 0 1]]
b == [[0 1 0 0 0 0] [0 1 0 0 0 0] [1 0 0 0 0 0] [1 0 0 0 0 0]]
c == [[0 0 0 1 0 0] [0 0 0 0 1 0] [0 0 0 0 0 1] [0 0 1 0 0 0]]
*/


alphas, betas, gammas, zx := snark.Utils.PF.R1CSToQAP(a, b, c)


ax, bx, cx, px := snark.Utils.PF.CombinePolynomials(w, alphas, betas, gammas)

hx := snark.Utils.PF.DivisorPolinomial(px, zx)

// hx==px/zx so px==hx*zx
assert.Equal(t, px, snark.Utils.PF.Mul(hx, zx))

// p(x) = a(x) * b(x) - c(x) == h(x) * z(x)
abc := snark.Utils.PF.Sub(pf.Mul(ax, bx), cx)
assert.Equal(t, abc, px)
hz := snark.Utils.PF.Mul(hx, zx)
assert.Equal(t, abc, hz)
	
div, rem := snark.Utils.PF.Div(px, zx)
assert.Equal(t, hx, div)
assert.Equal(t, rem, r1csqap.ArrayOfBigZeros(4))

// calculate trusted setup
setup, err := snark.GenerateTrustedSetup(len(w), circuit, alphas, betas, gammas, zx)
assert.Nil(t, err)
fmt.Println("t", setup.Toxic.T)

// piA = g1 * A(t), piB = g2 * B(t), piC = g1 * C(t), piH = g1 * H(t)
proof, err := snark.GenerateProofs(circuit, setup, hx, w)
assert.Nil(t, err)

assert.True(t, snark.VerifyProof(circuit, setup, proof))

Test

go test ./... -v

Thanks to @jbaylina, @bellesmarta, @adriamb for their explanations that helped to understand this a little bit. Also thanks to @vbuterin for all the published articles explaining the zkSNARKs.