Merge pull request #14336 from obscuren/metropolis-preparation

consensus, core/*, params: metropolis preparation refactor

11 files changed, 2538 insertions, 0 deletions

diff --git a/crypto/bn256/bn256.go b/crypto/bn256/bn256.go
new file mode 100644
index 000000000..92418369b
--- /dev/null
+++ b/crypto/bn256/bn256.go
@@ -0,0 +1,428 @@
+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package bn256 implements a particular bilinear group at the 128-bit security level.
+//
+// Bilinear groups are the basis of many of the new cryptographic protocols
+// that have been proposed over the past decade. They consist of a triplet of
+// groups (G₁, G₂ and GT) such that there exists a function e(g₁ˣ,g₂ʸ)=gTˣʸ
+// (where gₓ is a generator of the respective group). That function is called
+// a pairing function.
+//
+// This package specifically implements the Optimal Ate pairing over a 256-bit
+// Barreto-Naehrig curve as described in
+// http://cryptojedi.org/papers/dclxvi-20100714.pdf. Its output is compatible
+// with the implementation described in that paper.
+package bn256
+
+import (
+	"crypto/rand"
+	"io"
+	"math/big"
+)
+
+// BUG(agl): this implementation is not constant time.
+// TODO(agl): keep GF(p²) elements in Mongomery form.
+
+// G1 is an abstract cyclic group. The zero value is suitable for use as the
+// output of an operation, but cannot be used as an input.
+type G1 struct {
+	p *curvePoint
+}
+
+// RandomG1 returns x and g₁ˣ where x is a random, non-zero number read from r.
+func RandomG1(r io.Reader) (*big.Int, *G1, error) {
+	var k *big.Int
+	var err error
+
+	for {
+		k, err = rand.Int(r, Order)
+		if err != nil {
+			return nil, nil, err
+		}
+		if k.Sign() > 0 {
+			break
+		}
+	}
+
+	return k, new(G1).ScalarBaseMult(k), nil
+}
+
+func (g *G1) String() string {
+	return "bn256.G1" + g.p.String()
+}
+
+// CurvePoints returns p's curve points in big integer
+func (e *G1) CurvePoints() (*big.Int, *big.Int, *big.Int, *big.Int) {
+	return e.p.x, e.p.y, e.p.z, e.p.t
+}
+
+// ScalarBaseMult sets e to g*k where g is the generator of the group and
+// then returns e.
+func (e *G1) ScalarBaseMult(k *big.Int) *G1 {
+	if e.p == nil {
+		e.p = newCurvePoint(nil)
+	}
+	e.p.Mul(curveGen, k, new(bnPool))
+	return e
+}
+
+// ScalarMult sets e to a*k and then returns e.
+func (e *G1) ScalarMult(a *G1, k *big.Int) *G1 {
+	if e.p == nil {
+		e.p = newCurvePoint(nil)
+	}
+	e.p.Mul(a.p, k, new(bnPool))
+	return e
+}
+
+// Add sets e to a+b and then returns e.
+// BUG(agl): this function is not complete: a==b fails.
+func (e *G1) Add(a, b *G1) *G1 {
+	if e.p == nil {
+		e.p = newCurvePoint(nil)
+	}
+	e.p.Add(a.p, b.p, new(bnPool))
+	return e
+}
+
+// Neg sets e to -a and then returns e.
+func (e *G1) Neg(a *G1) *G1 {
+	if e.p == nil {
+		e.p = newCurvePoint(nil)
+	}
+	e.p.Negative(a.p)
+	return e
+}
+
+// Marshal converts n to a byte slice.
+func (n *G1) Marshal() []byte {
+	n.p.MakeAffine(nil)
+
+	xBytes := new(big.Int).Mod(n.p.x, P).Bytes()
+	yBytes := new(big.Int).Mod(n.p.y, P).Bytes()
+
+	// Each value is a 256-bit number.
+	const numBytes = 256 / 8
+
+	ret := make([]byte, numBytes*2)
+	copy(ret[1*numBytes-len(xBytes):], xBytes)
+	copy(ret[2*numBytes-len(yBytes):], yBytes)
+
+	return ret
+}
+
+// Unmarshal sets e to the result of converting the output of Marshal back into
+// a group element and then returns e.
+func (e *G1) Unmarshal(m []byte) (*G1, bool) {
+	// Each value is a 256-bit number.
+	const numBytes = 256 / 8
+
+	if len(m) != 2*numBytes {
+		return nil, false
+	}
+
+	if e.p == nil {
+		e.p = newCurvePoint(nil)
+	}
+
+	e.p.x.SetBytes(m[0*numBytes : 1*numBytes])
+	e.p.y.SetBytes(m[1*numBytes : 2*numBytes])
+
+	if e.p.x.Sign() == 0 && e.p.y.Sign() == 0 {
+		// This is the point at infinity.
+		e.p.y.SetInt64(1)
+		e.p.z.SetInt64(0)
+		e.p.t.SetInt64(0)
+	} else {
+		e.p.z.SetInt64(1)
+		e.p.t.SetInt64(1)
+
+		if !e.p.IsOnCurve() {
+			return nil, false
+		}
+	}
+
+	return e, true
+}
+
+// G2 is an abstract cyclic group. The zero value is suitable for use as the
+// output of an operation, but cannot be used as an input.
+type G2 struct {
+	p *twistPoint
+}
+
+// RandomG1 returns x and g₂ˣ where x is a random, non-zero number read from r.
+func RandomG2(r io.Reader) (*big.Int, *G2, error) {
+	var k *big.Int
+	var err error
+
+	for {
+		k, err = rand.Int(r, Order)
+		if err != nil {
+			return nil, nil, err
+		}
+		if k.Sign() > 0 {
+			break
+		}
+	}
+
+	return k, new(G2).ScalarBaseMult(k), nil
+}
+
+func (g *G2) String() string {
+	return "bn256.G2" + g.p.String()
+}
+
+// CurvePoints returns the curve points of p which includes the real
+// and imaginary parts of the curve point.
+func (e *G2) CurvePoints() (*gfP2, *gfP2, *gfP2, *gfP2) {
+	return e.p.x, e.p.y, e.p.z, e.p.t
+}
+
+// ScalarBaseMult sets e to g*k where g is the generator of the group and
+// then returns out.
+func (e *G2) ScalarBaseMult(k *big.Int) *G2 {
+	if e.p == nil {
+		e.p = newTwistPoint(nil)
+	}
+	e.p.Mul(twistGen, k, new(bnPool))
+	return e
+}
+
+// ScalarMult sets e to a*k and then returns e.
+func (e *G2) ScalarMult(a *G2, k *big.Int) *G2 {
+	if e.p == nil {
+		e.p = newTwistPoint(nil)
+	}
+	e.p.Mul(a.p, k, new(bnPool))
+	return e
+}
+
+// Add sets e to a+b and then returns e.
+// BUG(agl): this function is not complete: a==b fails.
+func (e *G2) Add(a, b *G2) *G2 {
+	if e.p == nil {
+		e.p = newTwistPoint(nil)
+	}
+	e.p.Add(a.p, b.p, new(bnPool))
+	return e
+}
+
+// Marshal converts n into a byte slice.
+func (n *G2) Marshal() []byte {
+	n.p.MakeAffine(nil)
+
+	xxBytes := new(big.Int).Mod(n.p.x.x, P).Bytes()
+	xyBytes := new(big.Int).Mod(n.p.x.y, P).Bytes()
+	yxBytes := new(big.Int).Mod(n.p.y.x, P).Bytes()
+	yyBytes := new(big.Int).Mod(n.p.y.y, P).Bytes()
+
+	// Each value is a 256-bit number.
+	const numBytes = 256 / 8
+
+	ret := make([]byte, numBytes*4)
+	copy(ret[1*numBytes-len(xxBytes):], xxBytes)
+	copy(ret[2*numBytes-len(xyBytes):], xyBytes)
+	copy(ret[3*numBytes-len(yxBytes):], yxBytes)
+	copy(ret[4*numBytes-len(yyBytes):], yyBytes)
+
+	return ret
+}
+
+// Unmarshal sets e to the result of converting the output of Marshal back into
+// a group element and then returns e.
+func (e *G2) Unmarshal(m []byte) (*G2, bool) {
+	// Each value is a 256-bit number.
+	const numBytes = 256 / 8
+
+	if len(m) != 4*numBytes {
+		return nil, false
+	}
+
+	if e.p == nil {
+		e.p = newTwistPoint(nil)
+	}
+
+	e.p.x.x.SetBytes(m[0*numBytes : 1*numBytes])
+	e.p.x.y.SetBytes(m[1*numBytes : 2*numBytes])
+	e.p.y.x.SetBytes(m[2*numBytes : 3*numBytes])
+	e.p.y.y.SetBytes(m[3*numBytes : 4*numBytes])
+
+	if e.p.x.x.Sign() == 0 &&
+		e.p.x.y.Sign() == 0 &&
+		e.p.y.x.Sign() == 0 &&
+		e.p.y.y.Sign() == 0 {
+		// This is the point at infinity.
+		e.p.y.SetOne()
+		e.p.z.SetZero()
+		e.p.t.SetZero()
+	} else {
+		e.p.z.SetOne()
+		e.p.t.SetOne()
+
+		if !e.p.IsOnCurve() {
+			return nil, false
+		}
+	}
+
+	return e, true
+}
+
+// GT is an abstract cyclic group. The zero value is suitable for use as the
+// output of an operation, but cannot be used as an input.
+type GT struct {
+	p *gfP12
+}
+
+func (g *GT) String() string {
+	return "bn256.GT" + g.p.String()
+}
+
+// ScalarMult sets e to a*k and then returns e.
+func (e *GT) ScalarMult(a *GT, k *big.Int) *GT {
+	if e.p == nil {
+		e.p = newGFp12(nil)
+	}
+	e.p.Exp(a.p, k, new(bnPool))
+	return e
+}
+
+// Add sets e to a+b and then returns e.
+func (e *GT) Add(a, b *GT) *GT {
+	if e.p == nil {
+		e.p = newGFp12(nil)
+	}
+	e.p.Mul(a.p, b.p, new(bnPool))
+	return e
+}
+
+// Neg sets e to -a and then returns e.
+func (e *GT) Neg(a *GT) *GT {
+	if e.p == nil {
+		e.p = newGFp12(nil)
+	}
+	e.p.Invert(a.p, new(bnPool))
+	return e
+}
+
+// Marshal converts n into a byte slice.
+func (n *GT) Marshal() []byte {
+	n.p.Minimal()
+
+	xxxBytes := n.p.x.x.x.Bytes()
+	xxyBytes := n.p.x.x.y.Bytes()
+	xyxBytes := n.p.x.y.x.Bytes()
+	xyyBytes := n.p.x.y.y.Bytes()
+	xzxBytes := n.p.x.z.x.Bytes()
+	xzyBytes := n.p.x.z.y.Bytes()
+	yxxBytes := n.p.y.x.x.Bytes()
+	yxyBytes := n.p.y.x.y.Bytes()
+	yyxBytes := n.p.y.y.x.Bytes()
+	yyyBytes := n.p.y.y.y.Bytes()
+	yzxBytes := n.p.y.z.x.Bytes()
+	yzyBytes := n.p.y.z.y.Bytes()
+
+	// Each value is a 256-bit number.
+	const numBytes = 256 / 8
+
+	ret := make([]byte, numBytes*12)
+	copy(ret[1*numBytes-len(xxxBytes):], xxxBytes)
+	copy(ret[2*numBytes-len(xxyBytes):], xxyBytes)
+	copy(ret[3*numBytes-len(xyxBytes):], xyxBytes)
+	copy(ret[4*numBytes-len(xyyBytes):], xyyBytes)
+	copy(ret[5*numBytes-len(xzxBytes):], xzxBytes)
+	copy(ret[6*numBytes-len(xzyBytes):], xzyBytes)
+	copy(ret[7*numBytes-len(yxxBytes):], yxxBytes)
+	copy(ret[8*numBytes-len(yxyBytes):], yxyBytes)
+	copy(ret[9*numBytes-len(yyxBytes):], yyxBytes)
+	copy(ret[10*numBytes-len(yyyBytes):], yyyBytes)
+	copy(ret[11*numBytes-len(yzxBytes):], yzxBytes)
+	copy(ret[12*numBytes-len(yzyBytes):], yzyBytes)
+
+	return ret
+}
+
+// Unmarshal sets e to the result of converting the output of Marshal back into
+// a group element and then returns e.
+func (e *GT) Unmarshal(m []byte) (*GT, bool) {
+	// Each value is a 256-bit number.
+	const numBytes = 256 / 8
+
+	if len(m) != 12*numBytes {
+		return nil, false
+	}
+
+	if e.p == nil {
+		e.p = newGFp12(nil)
+	}
+
+	e.p.x.x.x.SetBytes(m[0*numBytes : 1*numBytes])
+	e.p.x.x.y.SetBytes(m[1*numBytes : 2*numBytes])
+	e.p.x.y.x.SetBytes(m[2*numBytes : 3*numBytes])
+	e.p.x.y.y.SetBytes(m[3*numBytes : 4*numBytes])
+	e.p.x.z.x.SetBytes(m[4*numBytes : 5*numBytes])
+	e.p.x.z.y.SetBytes(m[5*numBytes : 6*numBytes])
+	e.p.y.x.x.SetBytes(m[6*numBytes : 7*numBytes])
+	e.p.y.x.y.SetBytes(m[7*numBytes : 8*numBytes])
+	e.p.y.y.x.SetBytes(m[8*numBytes : 9*numBytes])
+	e.p.y.y.y.SetBytes(m[9*numBytes : 10*numBytes])
+	e.p.y.z.x.SetBytes(m[10*numBytes : 11*numBytes])
+	e.p.y.z.y.SetBytes(m[11*numBytes : 12*numBytes])
+
+	return e, true
+}
+
+// Pair calculates an Optimal Ate pairing.
+func Pair(g1 *G1, g2 *G2) *GT {
+	return &GT{optimalAte(g2.p, g1.p, new(bnPool))}
+}
+
+func PairingCheck(a []*G1, b []*G2) bool {
+	pool := new(bnPool)
+	e := newGFp12(pool)
+	e.SetOne()
+	for i := 0; i < len(a); i++ {
+		new_e := miller(b[i].p, a[i].p, pool)
+		e.Mul(e, new_e, pool)
+	}
+	ret := finalExponentiation(e, pool)
+	e.Put(pool)
+	return ret.IsOne()
+}
+
+// bnPool implements a tiny cache of *big.Int objects that's used to reduce the
+// number of allocations made during processing.
+type bnPool struct {
+	bns   []*big.Int
+	count int
+}
+
+func (pool *bnPool) Get() *big.Int {
+	if pool == nil {
+		return new(big.Int)
+	}
+
+	pool.count++
+	l := len(pool.bns)
+	if l == 0 {
+		return new(big.Int)
+	}
+
+	bn := pool.bns[l-1]
+	pool.bns = pool.bns[:l-1]
+	return bn
+}
+
+func (pool *bnPool) Put(bn *big.Int) {
+	if pool == nil {
+		return
+	}
+	pool.bns = append(pool.bns, bn)
+	pool.count--
+}
+
+func (pool *bnPool) Count() int {
+	return pool.count
+}
diff --git a/crypto/bn256/bn256_test.go b/crypto/bn256/bn256_test.go
new file mode 100644
index 000000000..866065d0c
--- /dev/null
+++ b/crypto/bn256/bn256_test.go
@@ -0,0 +1,304 @@
+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package bn256
+
+import (
+	"bytes"
+	"crypto/rand"
+	"math/big"
+	"testing"
+)
+
+func TestGFp2Invert(t *testing.T) {
+	pool := new(bnPool)
+
+	a := newGFp2(pool)
+	a.x.SetString("23423492374", 10)
+	a.y.SetString("12934872398472394827398470", 10)
+
+	inv := newGFp2(pool)
+	inv.Invert(a, pool)
+
+	b := newGFp2(pool).Mul(inv, a, pool)
+	if b.x.Int64() != 0 || b.y.Int64() != 1 {
+		t.Fatalf("bad result for a^-1*a: %s %s", b.x, b.y)
+	}
+
+	a.Put(pool)
+	b.Put(pool)
+	inv.Put(pool)
+
+	if c := pool.Count(); c > 0 {
+		t.Errorf("Pool count non-zero: %d\n", c)
+	}
+}
+
+func isZero(n *big.Int) bool {
+	return new(big.Int).Mod(n, P).Int64() == 0
+}
+
+func isOne(n *big.Int) bool {
+	return new(big.Int).Mod(n, P).Int64() == 1
+}
+
+func TestGFp6Invert(t *testing.T) {
+	pool := new(bnPool)
+
+	a := newGFp6(pool)
+	a.x.x.SetString("239487238491", 10)
+	a.x.y.SetString("2356249827341", 10)
+	a.y.x.SetString("082659782", 10)
+	a.y.y.SetString("182703523765", 10)
+	a.z.x.SetString("978236549263", 10)
+	a.z.y.SetString("64893242", 10)
+
+	inv := newGFp6(pool)
+	inv.Invert(a, pool)
+
+	b := newGFp6(pool).Mul(inv, a, pool)
+	if !isZero(b.x.x) ||
+		!isZero(b.x.y) ||
+		!isZero(b.y.x) ||
+		!isZero(b.y.y) ||
+		!isZero(b.z.x) ||
+		!isOne(b.z.y) {
+		t.Fatalf("bad result for a^-1*a: %s", b)
+	}
+
+	a.Put(pool)
+	b.Put(pool)
+	inv.Put(pool)
+
+	if c := pool.Count(); c > 0 {
+		t.Errorf("Pool count non-zero: %d\n", c)
+	}
+}
+
+func TestGFp12Invert(t *testing.T) {
+	pool := new(bnPool)
+
+	a := newGFp12(pool)
+	a.x.x.x.SetString("239846234862342323958623", 10)
+	a.x.x.y.SetString("2359862352529835623", 10)
+	a.x.y.x.SetString("928836523", 10)
+	a.x.y.y.SetString("9856234", 10)
+	a.x.z.x.SetString("235635286", 10)
+	a.x.z.y.SetString("5628392833", 10)
+	a.y.x.x.SetString("252936598265329856238956532167968", 10)
+	a.y.x.y.SetString("23596239865236954178968", 10)
+	a.y.y.x.SetString("95421692834", 10)
+	a.y.y.y.SetString("236548", 10)
+	a.y.z.x.SetString("924523", 10)
+	a.y.z.y.SetString("12954623", 10)
+
+	inv := newGFp12(pool)
+	inv.Invert(a, pool)
+
+	b := newGFp12(pool).Mul(inv, a, pool)
+	if !isZero(b.x.x.x) ||
+		!isZero(b.x.x.y) ||
+		!isZero(b.x.y.x) ||
+		!isZero(b.x.y.y) ||
+		!isZero(b.x.z.x) ||
+		!isZero(b.x.z.y) ||
+		!isZero(b.y.x.x) ||
+		!isZero(b.y.x.y) ||
+		!isZero(b.y.y.x) ||
+		!isZero(b.y.y.y) ||
+		!isZero(b.y.z.x) ||
+		!isOne(b.y.z.y) {
+		t.Fatalf("bad result for a^-1*a: %s", b)
+	}
+
+	a.Put(pool)
+	b.Put(pool)
+	inv.Put(pool)
+
+	if c := pool.Count(); c > 0 {
+		t.Errorf("Pool count non-zero: %d\n", c)
+	}
+}
+
+func TestCurveImpl(t *testing.T) {
+	pool := new(bnPool)
+
+	g := &curvePoint{
+		pool.Get().SetInt64(1),
+		pool.Get().SetInt64(-2),
+		pool.Get().SetInt64(1),
+		pool.Get().SetInt64(0),
+	}
+
+	x := pool.Get().SetInt64(32498273234)
+	X := newCurvePoint(pool).Mul(g, x, pool)
+
+	y := pool.Get().SetInt64(98732423523)
+	Y := newCurvePoint(pool).Mul(g, y, pool)
+
+	s1 := newCurvePoint(pool).Mul(X, y, pool).MakeAffine(pool)
+	s2 := newCurvePoint(pool).Mul(Y, x, pool).MakeAffine(pool)
+
+	if s1.x.Cmp(s2.x) != 0 ||
+		s2.x.Cmp(s1.x) != 0 {
+		t.Errorf("DH points don't match: (%s, %s) (%s, %s)", s1.x, s1.y, s2.x, s2.y)
+	}
+
+	pool.Put(x)
+	X.Put(pool)
+	pool.Put(y)
+	Y.Put(pool)
+	s1.Put(pool)
+	s2.Put(pool)
+	g.Put(pool)
+
+	if c := pool.Count(); c > 0 {
+		t.Errorf("Pool count non-zero: %d\n", c)
+	}
+}
+
+func TestOrderG1(t *testing.T) {
+	g := new(G1).ScalarBaseMult(Order)
+	if !g.p.IsInfinity() {
+		t.Error("G1 has incorrect order")
+	}
+
+	one := new(G1).ScalarBaseMult(new(big.Int).SetInt64(1))
+	g.Add(g, one)
+	g.p.MakeAffine(nil)
+	if g.p.x.Cmp(one.p.x) != 0 || g.p.y.Cmp(one.p.y) != 0 {
+		t.Errorf("1+0 != 1 in G1")
+	}
+}
+
+func TestOrderG2(t *testing.T) {
+	g := new(G2).ScalarBaseMult(Order)
+	if !g.p.IsInfinity() {
+		t.Error("G2 has incorrect order")
+	}
+
+	one := new(G2).ScalarBaseMult(new(big.Int).SetInt64(1))
+	g.Add(g, one)
+	g.p.MakeAffine(nil)
+	if g.p.x.x.Cmp(one.p.x.x) != 0 ||
+		g.p.x.y.Cmp(one.p.x.y) != 0 ||
+		g.p.y.x.Cmp(one.p.y.x) != 0 ||
+		g.p.y.y.Cmp(one.p.y.y) != 0 {
+		t.Errorf("1+0 != 1 in G2")
+	}
+}
+
+func TestOrderGT(t *testing.T) {
+	gt := Pair(&G1{curveGen}, &G2{twistGen})
+	g := new(GT).ScalarMult(gt, Order)
+	if !g.p.IsOne() {
+		t.Error("GT has incorrect order")
+	}
+}
+
+func TestBilinearity(t *testing.T) {
+	for i := 0; i < 2; i++ {
+		a, p1, _ := RandomG1(rand.Reader)
+		b, p2, _ := RandomG2(rand.Reader)
+		e1 := Pair(p1, p2)
+
+		e2 := Pair(&G1{curveGen}, &G2{twistGen})
+		e2.ScalarMult(e2, a)
+		e2.ScalarMult(e2, b)
+
+		minusE2 := new(GT).Neg(e2)
+		e1.Add(e1, minusE2)
+
+		if !e1.p.IsOne() {
+			t.Fatalf("bad pairing result: %s", e1)
+		}
+	}
+}
+
+func TestG1Marshal(t *testing.T) {
+	g := new(G1).ScalarBaseMult(new(big.Int).SetInt64(1))
+	form := g.Marshal()
+	_, ok := new(G1).Unmarshal(form)
+	if !ok {
+		t.Fatalf("failed to unmarshal")
+	}
+
+	g.ScalarBaseMult(Order)
+	form = g.Marshal()
+	g2, ok := new(G1).Unmarshal(form)
+	if !ok {
+		t.Fatalf("failed to unmarshal ∞")
+	}
+	if !g2.p.IsInfinity() {
+		t.Fatalf("∞ unmarshaled incorrectly")
+	}
+}
+
+func TestG2Marshal(t *testing.T) {
+	g := new(G2).ScalarBaseMult(new(big.Int).SetInt64(1))
+	form := g.Marshal()
+	_, ok := new(G2).Unmarshal(form)
+	if !ok {
+		t.Fatalf("failed to unmarshal")
+	}
+
+	g.ScalarBaseMult(Order)
+	form = g.Marshal()
+	g2, ok := new(G2).Unmarshal(form)
+	if !ok {
+		t.Fatalf("failed to unmarshal ∞")
+	}
+	if !g2.p.IsInfinity() {
+		t.Fatalf("∞ unmarshaled incorrectly")
+	}
+}
+
+func TestG1Identity(t *testing.T) {
+	g := new(G1).ScalarBaseMult(new(big.Int).SetInt64(0))
+	if !g.p.IsInfinity() {
+		t.Error("failure")
+	}
+}
+
+func TestG2Identity(t *testing.T) {
+	g := new(G2).ScalarBaseMult(new(big.Int).SetInt64(0))
+	if !g.p.IsInfinity() {
+		t.Error("failure")
+	}
+}
+
+func TestTripartiteDiffieHellman(t *testing.T) {
+	a, _ := rand.Int(rand.Reader, Order)
+	b, _ := rand.Int(rand.Reader, Order)
+	c, _ := rand.Int(rand.Reader, Order)
+
+	pa, _ := new(G1).Unmarshal(new(G1).ScalarBaseMult(a).Marshal())
+	qa, _ := new(G2).Unmarshal(new(G2).ScalarBaseMult(a).Marshal())
+	pb, _ := new(G1).Unmarshal(new(G1).ScalarBaseMult(b).Marshal())
+	qb, _ := new(G2).Unmarshal(new(G2).ScalarBaseMult(b).Marshal())
+	pc, _ := new(G1).Unmarshal(new(G1).ScalarBaseMult(c).Marshal())
+	qc, _ := new(G2).Unmarshal(new(G2).ScalarBaseMult(c).Marshal())
+
+	k1 := Pair(pb, qc)
+	k1.ScalarMult(k1, a)
+	k1Bytes := k1.Marshal()
+
+	k2 := Pair(pc, qa)
+	k2.ScalarMult(k2, b)
+	k2Bytes := k2.Marshal()
+
+	k3 := Pair(pa, qb)
+	k3.ScalarMult(k3, c)
+	k3Bytes := k3.Marshal()
+
+	if !bytes.Equal(k1Bytes, k2Bytes) || !bytes.Equal(k2Bytes, k3Bytes) {
+		t.Errorf("keys didn't agree")
+	}
+}
+
+func BenchmarkPairing(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		Pair(&G1{curveGen}, &G2{twistGen})
+	}
+}
diff --git a/crypto/bn256/constants.go b/crypto/bn256/constants.go
new file mode 100644
index 000000000..ab649d7f3
--- /dev/null
+++ b/crypto/bn256/constants.go
@@ -0,0 +1,44 @@
+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package bn256
+
+import (
+	"math/big"
+)
+
+func bigFromBase10(s string) *big.Int {
+	n, _ := new(big.Int).SetString(s, 10)
+	return n
+}
+
+// u is the BN parameter that determines the prime: 1868033³.
+var u = bigFromBase10("4965661367192848881")
+
+// p is a prime over which we form a basic field: 36u⁴+36u³+24u²+6u+1.
+var P = bigFromBase10("21888242871839275222246405745257275088696311157297823662689037894645226208583")
+
+// Order is the number of elements in both G₁ and G₂: 36u⁴+36u³+18u²+6u+1.
+var Order = bigFromBase10("21888242871839275222246405745257275088548364400416034343698204186575808495617")
+
+// xiToPMinus1Over6 is ξ^((p-1)/6) where ξ = i+9.
+var xiToPMinus1Over6 = &gfP2{bigFromBase10("16469823323077808223889137241176536799009286646108169935659301613961712198316"), bigFromBase10("8376118865763821496583973867626364092589906065868298776909617916018768340080")}
+
+// xiToPMinus1Over3 is ξ^((p-1)/3) where ξ = i+9.
+var xiToPMinus1Over3 = &gfP2{bigFromBase10("10307601595873709700152284273816112264069230130616436755625194854815875713954"), bigFromBase10("21575463638280843010398324269430826099269044274347216827212613867836435027261")}
+
+// xiToPMinus1Over2 is ξ^((p-1)/2) where ξ = i+9.
+var xiToPMinus1Over2 = &gfP2{bigFromBase10("3505843767911556378687030309984248845540243509899259641013678093033130930403"), bigFromBase10("2821565182194536844548159561693502659359617185244120367078079554186484126554")}
+
+// xiToPSquaredMinus1Over3 is ξ^((p²-1)/3) where ξ = i+9.
+var xiToPSquaredMinus1Over3 = bigFromBase10("21888242871839275220042445260109153167277707414472061641714758635765020556616")
+
+// xiTo2PSquaredMinus2Over3 is ξ^((2p²-2)/3) where ξ = i+9 (a cubic root of unity, mod p).
+var xiTo2PSquaredMinus2Over3 = bigFromBase10("2203960485148121921418603742825762020974279258880205651966")
+
+// xiToPSquaredMinus1Over6 is ξ^((1p²-1)/6) where ξ = i+9 (a cubic root of -1, mod p).
+var xiToPSquaredMinus1Over6 = bigFromBase10("21888242871839275220042445260109153167277707414472061641714758635765020556617")
+
+// xiTo2PMinus2Over3 is ξ^((2p-2)/3) where ξ = i+9.
+var xiTo2PMinus2Over3 = &gfP2{bigFromBase10("19937756971775647987995932169929341994314640652964949448313374472400716661030"), bigFromBase10("2581911344467009335267311115468803099551665605076196740867805258568234346338")}
diff --git a/crypto/bn256/curve.go b/crypto/bn256/curve.go
new file mode 100644
index 000000000..233b1f252
--- /dev/null
+++ b/crypto/bn256/curve.go
@@ -0,0 +1,278 @@
+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package bn256
+
+import (
+	"math/big"
+)
+
+// curvePoint implements the elliptic curve y²=x³+3. Points are kept in
+// Jacobian form and t=z² when valid. G₁ is the set of points of this curve on
+// GF(p).
+type curvePoint struct {
+	x, y, z, t *big.Int
+}
+
+var curveB = new(big.Int).SetInt64(3)
+
+// curveGen is the generator of G₁.
+var curveGen = &curvePoint{
+	new(big.Int).SetInt64(1),
+	new(big.Int).SetInt64(-2),
+	new(big.Int).SetInt64(1),
+	new(big.Int).SetInt64(1),
+}
+
+func newCurvePoint(pool *bnPool) *curvePoint {
+	return &curvePoint{
+		pool.Get(),
+		pool.Get(),
+		pool.Get(),
+		pool.Get(),
+	}
+}
+
+func (c *curvePoint) String() string {
+	c.MakeAffine(new(bnPool))
+	return "(" + c.x.String() + ", " + c.y.String() + ")"
+}
+
+func (c *curvePoint) Put(pool *bnPool) {
+	pool.Put(c.x)
+	pool.Put(c.y)
+	pool.Put(c.z)
+	pool.Put(c.t)
+}
+
+func (c *curvePoint) Set(a *curvePoint) {
+	c.x.Set(a.x)
+	c.y.Set(a.y)
+	c.z.Set(a.z)
+	c.t.Set(a.t)
+}
+
+// IsOnCurve returns true iff c is on the curve where c must be in affine form.
+func (c *curvePoint) IsOnCurve() bool {
+	yy := new(big.Int).Mul(c.y, c.y)
+	xxx := new(big.Int).Mul(c.x, c.x)
+	xxx.Mul(xxx, c.x)
+	yy.Sub(yy, xxx)
+	yy.Sub(yy, curveB)
+	if yy.Sign() < 0 || yy.Cmp(P) >= 0 {
+		yy.Mod(yy, P)
+	}
+	return yy.Sign() == 0
+}
+
+func (c *curvePoint) SetInfinity() {
+	c.z.SetInt64(0)
+}
+
+func (c *curvePoint) IsInfinity() bool {
+	return c.z.Sign() == 0
+}
+
+func (c *curvePoint) Add(a, b *curvePoint, pool *bnPool) {
+	if a.IsInfinity() {
+		c.Set(b)
+		return
+	}
+	if b.IsInfinity() {
+		c.Set(a)
+		return
+	}
+
+	// See http://hyperelliptic.org/EFD/g1p/auto-code/shortw/jacobian-0/addition/add-2007-bl.op3
+
+	// Normalize the points by replacing a = [x1:y1:z1] and b = [x2:y2:z2]
+	// by [u1:s1:z1·z2] and [u2:s2:z1·z2]
+	// where u1 = x1·z2², s1 = y1·z2³ and u1 = x2·z1², s2 = y2·z1³
+	z1z1 := pool.Get().Mul(a.z, a.z)
+	z1z1.Mod(z1z1, P)
+	z2z2 := pool.Get().Mul(b.z, b.z)
+	z2z2.Mod(z2z2, P)
+	u1 := pool.Get().Mul(a.x, z2z2)
+	u1.Mod(u1, P)
+	u2 := pool.Get().Mul(b.x, z1z1)
+	u2.Mod(u2, P)
+
+	t := pool.Get().Mul(b.z, z2z2)
+	t.Mod(t, P)
+	s1 := pool.Get().Mul(a.y, t)
+	s1.Mod(s1, P)
+
+	t.Mul(a.z, z1z1)
+	t.Mod(t, P)
+	s2 := pool.Get().Mul(b.y, t)
+	s2.Mod(s2, P)
+
+	// Compute x = (2h)²(s²-u1-u2)
+	// where s = (s2-s1)/(u2-u1) is the slope of the line through
+	// (u1,s1) and (u2,s2). The extra factor 2h = 2(u2-u1) comes from the value of z below.
+	// This is also:
+	// 4(s2-s1)² - 4h²(u1+u2) = 4(s2-s1)² - 4h³ - 4h²(2u1)
+	//                        = r² - j - 2v
+	// with the notations below.
+	h := pool.Get().Sub(u2, u1)
+	xEqual := h.Sign() == 0
+
+	t.Add(h, h)
+	// i = 4h²
+	i := pool.Get().Mul(t, t)
+	i.Mod(i, P)
+	// j = 4h³
+	j := pool.Get().Mul(h, i)
+	j.Mod(j, P)
+
+	t.Sub(s2, s1)
+	yEqual := t.Sign() == 0
+	if xEqual && yEqual {
+		c.Double(a, pool)
+		return
+	}
+	r := pool.Get().Add(t, t)
+
+	v := pool.Get().Mul(u1, i)
+	v.Mod(v, P)
+
+	// t4 = 4(s2-s1)²
+	t4 := pool.Get().Mul(r, r)
+	t4.Mod(t4, P)
+	t.Add(v, v)
+	t6 := pool.Get().Sub(t4, j)
+	c.x.Sub(t6, t)
+
+	// Set y = -(2h)³(s1 + s*(x/4h²-u1))
+	// This is also
+	// y = - 2·s1·j - (s2-s1)(2x - 2i·u1) = r(v-x) - 2·s1·j
+	t.Sub(v, c.x) // t7
+	t4.Mul(s1, j) // t8
+	t4.Mod(t4, P)
+	t6.Add(t4, t4) // t9
+	t4.Mul(r, t)   // t10
+	t4.Mod(t4, P)
+	c.y.Sub(t4, t6)
+
+	// Set z = 2(u2-u1)·z1·z2 = 2h·z1·z2
+	t.Add(a.z, b.z) // t11
+	t4.Mul(t, t)    // t12
+	t4.Mod(t4, P)
+	t.Sub(t4, z1z1) // t13
+	t4.Sub(t, z2z2) // t14
+	c.z.Mul(t4, h)
+	c.z.Mod(c.z, P)
+
+	pool.Put(z1z1)
+	pool.Put(z2z2)
+	pool.Put(u1)
+	pool.Put(u2)
+	pool.Put(t)
+	pool.Put(s1)
+	pool.Put(s2)
+	pool.Put(h)
+	pool.Put(i)
+	pool.Put(j)
+	pool.Put(r)
+	pool.Put(v)
+	pool.Put(t4)
+	pool.Put(t6)
+}
+
+func (c *curvePoint) Double(a *curvePoint, pool *bnPool) {
+	// See http://hyperelliptic.org/EFD/g1p/auto-code/shortw/jacobian-0/doubling/dbl-2009-l.op3
+	A := pool.Get().Mul(a.x, a.x)
+	A.Mod(A, P)
+	B := pool.Get().Mul(a.y, a.y)
+	B.Mod(B, P)
+	C_ := pool.Get().Mul(B, B)
+	C_.Mod(C_, P)
+
+	t := pool.Get().Add(a.x, B)
+	t2 := pool.Get().Mul(t, t)
+	t2.Mod(t2, P)
+	t.Sub(t2, A)
+	t2.Sub(t, C_)
+	d := pool.Get().Add(t2, t2)
+	t.Add(A, A)
+	e := pool.Get().Add(t, A)
+	f := pool.Get().Mul(e, e)
+	f.Mod(f, P)
+
+	t.Add(d, d)
+	c.x.Sub(f, t)
+
+	t.Add(C_, C_)
+	t2.Add(t, t)
+	t.Add(t2, t2)
+	c.y.Sub(d, c.x)
+	t2.Mul(e, c.y)
+	t2.Mod(t2, P)
+	c.y.Sub(t2, t)
+
+	t.Mul(a.y, a.z)
+	t.Mod(t, P)
+	c.z.Add(t, t)
+
+	pool.Put(A)
+	pool.Put(B)
+	pool.Put(C_)
+	pool.Put(t)
+	pool.Put(t2)
+	pool.Put(d)
+	pool.Put(e)
+	pool.Put(f)
+}
+
+func (c *curvePoint) Mul(a *curvePoint, scalar *big.Int, pool *bnPool) *curvePoint {
+	sum := newCurvePoint(pool)
+	sum.SetInfinity()
+	t := newCurvePoint(pool)
+
+	for i := scalar.BitLen(); i >= 0; i-- {
+		t.Double(sum, pool)
+		if scalar.Bit(i) != 0 {
+			sum.Add(t, a, pool)
+		} else {
+			sum.Set(t)
+		}
+	}
+
+	c.Set(sum)
+	sum.Put(pool)
+	t.Put(pool)
+	return c
+}
+
+func (c *curvePoint) MakeAffine(pool *bnPool) *curvePoint {
+	if words := c.z.Bits(); len(words) == 1 && words[0] == 1 {
+		return c
+	}
+
+	zInv := pool.Get().ModInverse(c.z, P)
+	t := pool.Get().Mul(c.y, zInv)
+	t.Mod(t, P)
+	zInv2 := pool.Get().Mul(zInv, zInv)
+	zInv2.Mod(zInv2, P)
+	c.y.Mul(t, zInv2)
+	c.y.Mod(c.y, P)
+	t.Mul(c.x, zInv2)
+	t.Mod(t, P)
+	c.x.Set(t)
+	c.z.SetInt64(1)
+	c.t.SetInt64(1)
+
+	pool.Put(zInv)
+	pool.Put(t)
+	pool.Put(zInv2)
+
+	return c
+}
+
+func (c *curvePoint) Negative(a *curvePoint) {
+	c.x.Set(a.x)
+	c.y.Neg(a.y)
+	c.z.Set(a.z)
+	c.t.SetInt64(0)
+}
diff --git a/crypto/bn256/example_test.go b/crypto/bn256/example_test.go
new file mode 100644
index 000000000..b2d19807a
--- /dev/null
+++ b/crypto/bn256/example_test.go
@@ -0,0 +1,43 @@
+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package bn256
+
+import (
+	"crypto/rand"
+)
+
+func ExamplePair() {
+	// This implements the tripartite Diffie-Hellman algorithm from "A One
+	// Round Protocol for Tripartite Diffie-Hellman", A. Joux.
+	// http://www.springerlink.com/content/cddc57yyva0hburb/fulltext.pdf
+
+	// Each of three parties, a, b and c, generate a private value.
+	a, _ := rand.Int(rand.Reader, Order)
+	b, _ := rand.Int(rand.Reader, Order)
+	c, _ := rand.Int(rand.Reader, Order)
+
+	// Then each party calculates g₁ and g₂ times their private value.
+	pa := new(G1).ScalarBaseMult(a)
+	qa := new(G2).ScalarBaseMult(a)
+
+	pb := new(G1).ScalarBaseMult(b)
+	qb := new(G2).ScalarBaseMult(b)
+
+	pc := new(G1).ScalarBaseMult(c)
+	qc := new(G2).ScalarBaseMult(c)
+
+	// Now each party exchanges its public values with the other two and
+	// all parties can calculate the shared key.
+	k1 := Pair(pb, qc)
+	k1.ScalarMult(k1, a)
+
+	k2 := Pair(pc, qa)
+	k2.ScalarMult(k2, b)
+
+	k3 := Pair(pa, qb)
+	k3.ScalarMult(k3, c)
+
+	// k1, k2 and k3 will all be equal.
+}
diff --git a/crypto/bn256/gfp12.go b/crypto/bn256/gfp12.go
new file mode 100644
index 000000000..f084eddf2
--- /dev/null
+++ b/crypto/bn256/gfp12.go
@@ -0,0 +1,200 @@
+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package bn256
+
+// For details of the algorithms used, see "Multiplication and Squaring on
+// Pairing-Friendly Fields, Devegili et al.
+// http://eprint.iacr.org/2006/471.pdf.
+
+import (
+	"math/big"
+)
+
+// gfP12 implements the field of size p¹² as a quadratic extension of gfP6
+// where ω²=τ.
+type gfP12 struct {
+	x, y *gfP6 // value is xω + y
+}
+
+func newGFp12(pool *bnPool) *gfP12 {
+	return &gfP12{newGFp6(pool), newGFp6(pool)}
+}
+
+func (e *gfP12) String() string {
+	return "(" + e.x.String() + "," + e.y.String() + ")"
+}
+
+func (e *gfP12) Put(pool *bnPool) {
+	e.x.Put(pool)
+	e.y.Put(pool)
+}
+
+func (e *gfP12) Set(a *gfP12) *gfP12 {
+	e.x.Set(a.x)
+	e.y.Set(a.y)
+	return e
+}
+
+func (e *gfP12) SetZero() *gfP12 {
+	e.x.SetZero()
+	e.y.SetZero()
+	return e
+}
+
+func (e *gfP12) SetOne() *gfP12 {
+	e.x.SetZero()
+	e.y.SetOne()
+	return e
+}
+
+func (e *gfP12) Minimal() {
+	e.x.Minimal()
+	e.y.Minimal()
+}
+
+func (e *gfP12) IsZero() bool {
+	e.Minimal()
+	return e.x.IsZero() && e.y.IsZero()
+}
+
+func (e *gfP12) IsOne() bool {
+	e.Minimal()
+	return e.x.IsZero() && e.y.IsOne()
+}
+
+func (e *gfP12) Conjugate(a *gfP12) *gfP12 {
+	e.x.Negative(a.x)
+	e.y.Set(a.y)
+	return a
+}
+
+func (e *gfP12) Negative(a *gfP12) *gfP12 {
+	e.x.Negative(a.x)
+	e.y.Negative(a.y)
+	return e
+}
+
+// Frobenius computes (xω+y)^p = x^p ω·ξ^((p-1)/6) + y^p
+func (e *gfP12) Frobenius(a *gfP12, pool *bnPool) *gfP12 {
+	e.x.Frobenius(a.x, pool)
+	e.y.Frobenius(a.y, pool)
+	e.x.MulScalar(e.x, xiToPMinus1Over6, pool)
+	return e
+}
+
+// FrobeniusP2 computes (xω+y)^p² = x^p² ω·ξ^((p²-1)/6) + y^p²
+func (e *gfP12) FrobeniusP2(a *gfP12, pool *bnPool) *gfP12 {
+	e.x.FrobeniusP2(a.x)
+	e.x.MulGFP(e.x, xiToPSquaredMinus1Over6)
+	e.y.FrobeniusP2(a.y)
+	return e
+}
+
+func (e *gfP12) Add(a, b *gfP12) *gfP12 {
+	e.x.Add(a.x, b.x)
+	e.y.Add(a.y, b.y)
+	return e
+}
+
+func (e *gfP12) Sub(a, b *gfP12) *gfP12 {
+	e.x.Sub(a.x, b.x)
+	e.y.Sub(a.y, b.y)
+	return e
+}
+
+func (e *gfP12) Mul(a, b *gfP12, pool *bnPool) *gfP12 {
+	tx := newGFp6(pool)
+	tx.Mul(a.x, b.y, pool)
+	t := newGFp6(pool)
+	t.Mul(b.x, a.y, pool)
+	tx.Add(tx, t)
+
+	ty := newGFp6(pool)
+	ty.Mul(a.y, b.y, pool)
+	t.Mul(a.x, b.x, pool)
+	t.MulTau(t, pool)
+	e.y.Add(ty, t)
+	e.x.Set(tx)
+
+	tx.Put(pool)
+	ty.Put(pool)
+	t.Put(pool)
+	return e
+}
+
+func (e *gfP12) MulScalar(a *gfP12, b *gfP6, pool *bnPool) *gfP12 {
+	e.x.Mul(e.x, b, pool)
+	e.y.Mul(e.y, b, pool)
+	return e
+}
+
+func (c *gfP12) Exp(a *gfP12, power *big.Int, pool *bnPool) *gfP12 {
+	sum := newGFp12(pool)
+	sum.SetOne()
+	t := newGFp12(pool)
+
+	for i := power.BitLen() - 1; i >= 0; i-- {
+		t.Square(sum, pool)
+		if power.Bit(i) != 0 {
+			sum.Mul(t, a, pool)
+		} else {
+			sum.Set(t)
+		}
+	}
+
+	c.Set(sum)
+
+	sum.Put(pool)
+	t.Put(pool)
+
+	return c
+}
+
+func (e *gfP12) Square(a *gfP12, pool *bnPool) *gfP12 {
+	// Complex squaring algorithm
+	v0 := newGFp6(pool)
+	v0.Mul(a.x, a.y, pool)
+
+	t := newGFp6(pool)
+	t.MulTau(a.x, pool)
+	t.Add(a.y, t)
+	ty := newGFp6(pool)
+	ty.Add(a.x, a.y)
+	ty.Mul(ty, t, pool)
+	ty.Sub(ty, v0)
+	t.MulTau(v0, pool)
+	ty.Sub(ty, t)
+
+	e.y.Set(ty)
+	e.x.Double(v0)
+
+	v0.Put(pool)
+	t.Put(pool)
+	ty.Put(pool)
+
+	return e
+}
+
+func (e *gfP12) Invert(a *gfP12, pool *bnPool) *gfP12 {
+	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
+	// ftp://136.206.11.249/pub/crypto/pairings.pdf
+	t1 := newGFp6(pool)
+	t2 := newGFp6(pool)
+
+	t1.Square(a.x, pool)
+	t2.Square(a.y, pool)
+	t1.MulTau(t1, pool)
+	t1.Sub(t2, t1)
+	t2.Invert(t1, pool)
+
+	e.x.Negative(a.x)
+	e.y.Set(a.y)
+	e.MulScalar(e, t2, pool)
+
+	t1.Put(pool)
+	t2.Put(pool)
+
+	return e
+}
diff --git a/crypto/bn256/gfp2.go b/crypto/bn256/gfp2.go
new file mode 100644
index 000000000..3981f6cb4
--- /dev/null
+++ b/crypto/bn256/gfp2.go
@@ -0,0 +1,227 @@
+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package bn256
+
+// For details of the algorithms used, see "Multiplication and Squaring on
+// Pairing-Friendly Fields, Devegili et al.
+// http://eprint.iacr.org/2006/471.pdf.
+
+import (
+	"math/big"
+)
+
+// gfP2 implements a field of size p² as a quadratic extension of the base
+// field where i²=-1.
+type gfP2 struct {
+	x, y *big.Int // value is xi+y.
+}
+
+func newGFp2(pool *bnPool) *gfP2 {
+	return &gfP2{pool.Get(), pool.Get()}
+}
+
+func (e *gfP2) String() string {
+	x := new(big.Int).Mod(e.x, P)
+	y := new(big.Int).Mod(e.y, P)
+	return "(" + x.String() + "," + y.String() + ")"
+}
+
+func (e *gfP2) Put(pool *bnPool) {
+	pool.Put(e.x)
+	pool.Put(e.y)
+}
+
+func (e *gfP2) Set(a *gfP2) *gfP2 {
+	e.x.Set(a.x)
+	e.y.Set(a.y)
+	return e
+}
+
+func (e *gfP2) SetZero() *gfP2 {
+	e.x.SetInt64(0)
+	e.y.SetInt64(0)
+	return e
+}
+
+func (e *gfP2) SetOne() *gfP2 {
+	e.x.SetInt64(0)
+	e.y.SetInt64(1)
+	return e
+}
+
+func (e *gfP2) Minimal() {
+	if e.x.Sign() < 0 || e.x.Cmp(P) >= 0 {
+		e.x.Mod(e.x, P)
+	}
+	if e.y.Sign() < 0 || e.y.Cmp(P) >= 0 {
+		e.y.Mod(e.y, P)
+	}
+}
+
+func (e *gfP2) IsZero() bool {
+	return e.x.Sign() == 0 && e.y.Sign() == 0
+}
+
+func (e *gfP2) IsOne() bool {
+	if e.x.Sign() != 0 {
+		return false
+	}
+	words := e.y.Bits()
+	return len(words) == 1 && words[0] == 1
+}
+
+func (e *gfP2) Conjugate(a *gfP2) *gfP2 {
+	e.y.Set(a.y)
+	e.x.Neg(a.x)
+	return e
+}
+
+func (e *gfP2) Negative(a *gfP2) *gfP2 {
+	e.x.Neg(a.x)
+	e.y.Neg(a.y)
+	return e
+}
+
+func (e *gfP2) Add(a, b *gfP2) *gfP2 {
+	e.x.Add(a.x, b.x)
+	e.y.Add(a.y, b.y)
+	return e
+}
+
+func (e *gfP2) Sub(a, b *gfP2) *gfP2 {
+	e.x.Sub(a.x, b.x)
+	e.y.Sub(a.y, b.y)
+	return e
+}
+
+func (e *gfP2) Double(a *gfP2) *gfP2 {
+	e.x.Lsh(a.x, 1)
+	e.y.Lsh(a.y, 1)
+	return e
+}
+
+func (c *gfP2) Exp(a *gfP2, power *big.Int, pool *bnPool) *gfP2 {
+	sum := newGFp2(pool)
+	sum.SetOne()
+	t := newGFp2(pool)
+
+	for i := power.BitLen() - 1; i >= 0; i-- {
+		t.Square(sum, pool)
+		if power.Bit(i) != 0 {
+			sum.Mul(t, a, pool)
+		} else {
+			sum.Set(t)
+		}
+	}
+
+	c.Set(sum)
+
+	sum.Put(pool)
+	t.Put(pool)
+
+	return c
+}
+
+// See "Multiplication and Squaring in Pairing-Friendly Fields",
+// http://eprint.iacr.org/2006/471.pdf
+func (e *gfP2) Mul(a, b *gfP2, pool *bnPool) *gfP2 {
+	tx := pool.Get().Mul(a.x, b.y)
+	t := pool.Get().Mul(b.x, a.y)
+	tx.Add(tx, t)
+	tx.Mod(tx, P)
+
+	ty := pool.Get().Mul(a.y, b.y)
+	t.Mul(a.x, b.x)
+	ty.Sub(ty, t)
+	e.y.Mod(ty, P)
+	e.x.Set(tx)
+
+	pool.Put(tx)
+	pool.Put(ty)
+	pool.Put(t)
+
+	return e
+}
+
+func (e *gfP2) MulScalar(a *gfP2, b *big.Int) *gfP2 {
+	e.x.Mul(a.x, b)
+	e.y.Mul(a.y, b)
+	return e
+}
+
+// MulXi sets e=ξa where ξ=i+9 and then returns e.
+func (e *gfP2) MulXi(a *gfP2, pool *bnPool) *gfP2 {
+	// (xi+y)(i+3) = (9x+y)i+(9y-x)
+	tx := pool.Get().Lsh(a.x, 3)
+	tx.Add(tx, a.x)
+	tx.Add(tx, a.y)
+
+	ty := pool.Get().Lsh(a.y, 3)
+	ty.Add(ty, a.y)
+	ty.Sub(ty, a.x)
+
+	e.x.Set(tx)
+	e.y.Set(ty)
+
+	pool.Put(tx)
+	pool.Put(ty)
+
+	return e
+}
+
+func (e *gfP2) Square(a *gfP2, pool *bnPool) *gfP2 {
+	// Complex squaring algorithm:
+	// (xi+b)² = (x+y)(y-x) + 2*i*x*y
+	t1 := pool.Get().Sub(a.y, a.x)
+	t2 := pool.Get().Add(a.x, a.y)
+	ty := pool.Get().Mul(t1, t2)
+	ty.Mod(ty, P)
+
+	t1.Mul(a.x, a.y)
+	t1.Lsh(t1, 1)
+
+	e.x.Mod(t1, P)
+	e.y.Set(ty)
+
+	pool.Put(t1)
+	pool.Put(t2)
+	pool.Put(ty)
+
+	return e
+}
+
+func (e *gfP2) Invert(a *gfP2, pool *bnPool) *gfP2 {
+	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
+	// ftp://136.206.11.249/pub/crypto/pairings.pdf
+	t := pool.Get()
+	t.Mul(a.y, a.y)
+	t2 := pool.Get()
+	t2.Mul(a.x, a.x)
+	t.Add(t, t2)
+
+	inv := pool.Get()
+	inv.ModInverse(t, P)
+
+	e.x.Neg(a.x)
+	e.x.Mul(e.x, inv)
+	e.x.Mod(e.x, P)
+
+	e.y.Mul(a.y, inv)
+	e.y.Mod(e.y, P)
+
+	pool.Put(t)
+	pool.Put(t2)
+	pool.Put(inv)
+
+	return e
+}
+
+func (e *gfP2) Real() *big.Int {
+	return e.x
+}
+
+func (e *gfP2) Imag() *big.Int {
+	return e.y
+}
diff --git a/crypto/bn256/gfp6.go b/crypto/bn256/gfp6.go
new file mode 100644
index 000000000..218856617
--- /dev/null
+++ b/crypto/bn256/gfp6.go
@@ -0,0 +1,296 @@
+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package bn256
+
+// For details of the algorithms used, see "Multiplication and Squaring on
+// Pairing-Friendly Fields, Devegili et al.
+// http://eprint.iacr.org/2006/471.pdf.
+
+import (
+	"math/big"
+)
+
+// gfP6 implements the field of size p⁶ as a cubic extension of gfP2 where τ³=ξ
+// and ξ=i+9.
+type gfP6 struct {
+	x, y, z *gfP2 // value is xτ² + yτ + z
+}
+
+func newGFp6(pool *bnPool) *gfP6 {
+	return &gfP6{newGFp2(pool), newGFp2(pool), newGFp2(pool)}
+}
+
+func (e *gfP6) String() string {
+	return "(" + e.x.String() + "," + e.y.String() + "," + e.z.String() + ")"
+}
+
+func (e *gfP6) Put(pool *bnPool) {
+	e.x.Put(pool)
+	e.y.Put(pool)
+	e.z.Put(pool)
+}
+
+func (e *gfP6) Set(a *gfP6) *gfP6 {
+	e.x.Set(a.x)
+	e.y.Set(a.y)
+	e.z.Set(a.z)
+	return e
+}
+
+func (e *gfP6) SetZero() *gfP6 {
+	e.x.SetZero()
+	e.y.SetZero()
+	e.z.SetZero()
+	return e
+}
+
+func (e *gfP6) SetOne() *gfP6 {
+	e.x.SetZero()
+	e.y.SetZero()
+	e.z.SetOne()
+	return e
+}
+
+func (e *gfP6) Minimal() {
+	e.x.Minimal()
+	e.y.Minimal()
+	e.z.Minimal()
+}
+
+func (e *gfP6) IsZero() bool {
+	return e.x.IsZero() && e.y.IsZero() && e.z.IsZero()
+}
+
+func (e *gfP6) IsOne() bool {
+	return e.x.IsZero() && e.y.IsZero() && e.z.IsOne()
+}
+
+func (e *gfP6) Negative(a *gfP6) *gfP6 {
+	e.x.Negative(a.x)
+	e.y.Negative(a.y)
+	e.z.Negative(a.z)
+	return e
+}
+
+func (e *gfP6) Frobenius(a *gfP6, pool *bnPool) *gfP6 {
+	e.x.Conjugate(a.x)
+	e.y.Conjugate(a.y)
+	e.z.Conjugate(a.z)
+
+	e.x.Mul(e.x, xiTo2PMinus2Over3, pool)
+	e.y.Mul(e.y, xiToPMinus1Over3, pool)
+	return e
+}
+
+// FrobeniusP2 computes (xτ²+yτ+z)^(p²) = xτ^(2p²) + yτ^(p²) + z
+func (e *gfP6) FrobeniusP2(a *gfP6) *gfP6 {
+	// τ^(2p²) = τ²τ^(2p²-2) = τ²ξ^((2p²-2)/3)
+	e.x.MulScalar(a.x, xiTo2PSquaredMinus2Over3)
+	// τ^(p²) = ττ^(p²-1) = τξ^((p²-1)/3)
+	e.y.MulScalar(a.y, xiToPSquaredMinus1Over3)
+	e.z.Set(a.z)
+	return e
+}
+
+func (e *gfP6) Add(a, b *gfP6) *gfP6 {
+	e.x.Add(a.x, b.x)
+	e.y.Add(a.y, b.y)
+	e.z.Add(a.z, b.z)
+	return e
+}
+
+func (e *gfP6) Sub(a, b *gfP6) *gfP6 {
+	e.x.Sub(a.x, b.x)
+	e.y.Sub(a.y, b.y)
+	e.z.Sub(a.z, b.z)
+	return e
+}
+
+func (e *gfP6) Double(a *gfP6) *gfP6 {
+	e.x.Double(a.x)
+	e.y.Double(a.y)
+	e.z.Double(a.z)
+	return e
+}
+
+func (e *gfP6) Mul(a, b *gfP6, pool *bnPool) *gfP6 {
+	// "Multiplication and Squaring on Pairing-Friendly Fields"
+	// Section 4, Karatsuba method.
+	// http://eprint.iacr.org/2006/471.pdf
+
+	v0 := newGFp2(pool)
+	v0.Mul(a.z, b.z, pool)
+	v1 := newGFp2(pool)
+	v1.Mul(a.y, b.y, pool)
+	v2 := newGFp2(pool)
+	v2.Mul(a.x, b.x, pool)
+
+	t0 := newGFp2(pool)
+	t0.Add(a.x, a.y)
+	t1 := newGFp2(pool)
+	t1.Add(b.x, b.y)
+	tz := newGFp2(pool)
+	tz.Mul(t0, t1, pool)
+
+	tz.Sub(tz, v1)
+	tz.Sub(tz, v2)
+	tz.MulXi(tz, pool)
+	tz.Add(tz, v0)
+
+	t0.Add(a.y, a.z)
+	t1.Add(b.y, b.z)
+	ty := newGFp2(pool)
+	ty.Mul(t0, t1, pool)
+	ty.Sub(ty, v0)
+	ty.Sub(ty, v1)
+	t0.MulXi(v2, pool)
+	ty.Add(ty, t0)
+
+	t0.Add(a.x, a.z)
+	t1.Add(b.x, b.z)
+	tx := newGFp2(pool)
+	tx.Mul(t0, t1, pool)
+	tx.Sub(tx, v0)
+	tx.Add(tx, v1)
+	tx.Sub(tx, v2)
+
+	e.x.Set(tx)
+	e.y.Set(ty)
+	e.z.Set(tz)
+
+	t0.Put(pool)
+	t1.Put(pool)
+	tx.Put(pool)
+	ty.Put(pool)
+	tz.Put(pool)
+	v0.Put(pool)
+	v1.Put(pool)
+	v2.Put(pool)
+	return e
+}
+
+func (e *gfP6) MulScalar(a *gfP6, b *gfP2, pool *bnPool) *gfP6 {
+	e.x.Mul(a.x, b, pool)
+	e.y.Mul(a.y, b, pool)
+	e.z.Mul(a.z, b, pool)
+	return e
+}
+
+func (e *gfP6) MulGFP(a *gfP6, b *big.Int) *gfP6 {
+	e.x.MulScalar(a.x, b)
+	e.y.MulScalar(a.y, b)
+	e.z.MulScalar(a.z, b)
+	return e
+}
+
+// MulTau computes τ·(aτ²+bτ+c) = bτ²+cτ+aξ
+func (e *gfP6) MulTau(a *gfP6, pool *bnPool) {
+	tz := newGFp2(pool)
+	tz.MulXi(a.x, pool)
+	ty := newGFp2(pool)
+	ty.Set(a.y)
+	e.y.Set(a.z)
+	e.x.Set(ty)
+	e.z.Set(tz)
+	tz.Put(pool)
+	ty.Put(pool)
+}
+
+func (e *gfP6) Square(a *gfP6, pool *bnPool) *gfP6 {
+	v0 := newGFp2(pool).Square(a.z, pool)
+	v1 := newGFp2(pool).Square(a.y, pool)
+	v2 := newGFp2(pool).Square(a.x, pool)
+
+	c0 := newGFp2(pool).Add(a.x, a.y)
+	c0.Square(c0, pool)
+	c0.Sub(c0, v1)
+	c0.Sub(c0, v2)
+	c0.MulXi(c0, pool)
+	c0.Add(c0, v0)
+
+	c1 := newGFp2(pool).Add(a.y, a.z)
+	c1.Square(c1, pool)
+	c1.Sub(c1, v0)
+	c1.Sub(c1, v1)
+	xiV2 := newGFp2(pool).MulXi(v2, pool)
+	c1.Add(c1, xiV2)
+
+	c2 := newGFp2(pool).Add(a.x, a.z)
+	c2.Square(c2, pool)
+	c2.Sub(c2, v0)
+	c2.Add(c2, v1)
+	c2.Sub(c2, v2)
+
+	e.x.Set(c2)
+	e.y.Set(c1)
+	e.z.Set(c0)
+
+	v0.Put(pool)
+	v1.Put(pool)
+	v2.Put(pool)
+	c0.Put(pool)
+	c1.Put(pool)
+	c2.Put(pool)
+	xiV2.Put(pool)
+
+	return e
+}
+
+func (e *gfP6) Invert(a *gfP6, pool *bnPool) *gfP6 {
+	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
+	// ftp://136.206.11.249/pub/crypto/pairings.pdf
+
+	// Here we can give a short explanation of how it works: let j be a cubic root of
+	// unity in GF(p²) so that 1+j+j²=0.
+	// Then (xτ² + yτ + z)(xj²τ² + yjτ + z)(xjτ² + yj²τ + z)
+	// = (xτ² + yτ + z)(Cτ²+Bτ+A)
+	// = (x³ξ²+y³ξ+z³-3ξxyz) = F is an element of the base field (the norm).
+	//
+	// On the other hand (xj²τ² + yjτ + z)(xjτ² + yj²τ + z)
+	// = τ²(y²-ξxz) + τ(ξx²-yz) + (z²-ξxy)
+	//
+	// So that's why A = (z²-ξxy), B = (ξx²-yz), C = (y²-ξxz)
+	t1 := newGFp2(pool)
+
+	A := newGFp2(pool)
+	A.Square(a.z, pool)
+	t1.Mul(a.x, a.y, pool)
+	t1.MulXi(t1, pool)
+	A.Sub(A, t1)
+
+	B := newGFp2(pool)
+	B.Square(a.x, pool)
+	B.MulXi(B, pool)
+	t1.Mul(a.y, a.z, pool)
+	B.Sub(B, t1)
+
+	C_ := newGFp2(pool)
+	C_.Square(a.y, pool)
+	t1.Mul(a.x, a.z, pool)
+	C_.Sub(C_, t1)
+
+	F := newGFp2(pool)
+	F.Mul(C_, a.y, pool)
+	F.MulXi(F, pool)
+	t1.Mul(A, a.z, pool)
+	F.Add(F, t1)
+	t1.Mul(B, a.x, pool)
+	t1.MulXi(t1, pool)
+	F.Add(F, t1)
+
+	F.Invert(F, pool)
+
+	e.x.Mul(C_, F, pool)
+	e.y.Mul(B, F, pool)
+	e.z.Mul(A, F, pool)
+
+	t1.Put(pool)
+	A.Put(pool)
+	B.Put(pool)
+	C_.Put(pool)
+	F.Put(pool)
+
+	return e
+}
diff --git a/crypto/bn256/main_test.go b/crypto/bn256/main_test.go
new file mode 100644
index 000000000..0230f1b19
--- /dev/null
+++ b/crypto/bn256/main_test.go
@@ -0,0 +1,71 @@
+package bn256
+
+import (
+	"testing"
+
+	"crypto/rand"
+)
+
+func TestRandomG2Marshal(t *testing.T) {
+	for i := 0; i < 10; i++ {
+		n, g2, err := RandomG2(rand.Reader)
+		if err != nil {
+			t.Error(err)
+			continue
+		}
+		t.Logf("%d: %x\n", n, g2.Marshal())
+	}
+}
+
+func TestPairings(t *testing.T) {
+	a1 := new(G1).ScalarBaseMult(bigFromBase10("1"))
+	a2 := new(G1).ScalarBaseMult(bigFromBase10("2"))
+	a37 := new(G1).ScalarBaseMult(bigFromBase10("37"))
+	an1 := new(G1).ScalarBaseMult(bigFromBase10("21888242871839275222246405745257275088548364400416034343698204186575808495616"))
+
+	b0 := new(G2).ScalarBaseMult(bigFromBase10("0"))
+	b1 := new(G2).ScalarBaseMult(bigFromBase10("1"))
+	b2 := new(G2).ScalarBaseMult(bigFromBase10("2"))
+	b27 := new(G2).ScalarBaseMult(bigFromBase10("27"))
+	b999 := new(G2).ScalarBaseMult(bigFromBase10("999"))
+	bn1 := new(G2).ScalarBaseMult(bigFromBase10("21888242871839275222246405745257275088548364400416034343698204186575808495616"))
+
+	p1 := Pair(a1, b1)
+	pn1 := Pair(a1, bn1)
+	np1 := Pair(an1, b1)
+	if pn1.String() != np1.String() {
+		t.Error("Pairing mismatch: e(a, -b) != e(-a, b)")
+	}
+	if !PairingCheck([]*G1{a1, an1}, []*G2{b1, b1}) {
+		t.Error("MultiAte check gave false negative!")
+	}
+	p0 := new(GT).Add(p1, pn1)
+	p0_2 := Pair(a1, b0)
+	if p0.String() != p0_2.String() {
+		t.Error("Pairing mismatch: e(a, b) * e(a, -b) != 1")
+	}
+	p0_3 := new(GT).ScalarMult(p1, bigFromBase10("21888242871839275222246405745257275088548364400416034343698204186575808495617"))
+	if p0.String() != p0_3.String() {
+		t.Error("Pairing mismatch: e(a, b) has wrong order")
+	}
+	p2 := Pair(a2, b1)
+	p2_2 := Pair(a1, b2)
+	p2_3 := new(GT).ScalarMult(p1, bigFromBase10("2"))
+	if p2.String() != p2_2.String() {
+		t.Error("Pairing mismatch: e(a, b * 2) != e(a * 2, b)")
+	}
+	if p2.String() != p2_3.String() {
+		t.Error("Pairing mismatch: e(a, b * 2) != e(a, b) ** 2")
+	}
+	if p2.String() == p1.String() {
+		t.Error("Pairing is degenerate!")
+	}
+	if PairingCheck([]*G1{a1, a1}, []*G2{b1, b1}) {
+		t.Error("MultiAte check gave false positive!")
+	}
+	p999 := Pair(a37, b27)
+	p999_2 := Pair(a1, b999)
+	if p999.String() != p999_2.String() {
+		t.Error("Pairing mismatch: e(a * 37, b * 27) != e(a, b * 999)")
+	}
+}
diff --git a/crypto/bn256/optate.go b/crypto/bn256/optate.go
new file mode 100644
index 000000000..68716b62b
--- /dev/null
+++ b/crypto/bn256/optate.go
@@ -0,0 +1,398 @@
+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package bn256
+
+func lineFunctionAdd(r, p *twistPoint, q *curvePoint, r2 *gfP2, pool *bnPool) (a, b, c *gfP2, rOut *twistPoint) {
+	// See the mixed addition algorithm from "Faster Computation of the
+	// Tate Pairing", http://arxiv.org/pdf/0904.0854v3.pdf
+
+	B := newGFp2(pool).Mul(p.x, r.t, pool)
+
+	D := newGFp2(pool).Add(p.y, r.z)
+	D.Square(D, pool)
+	D.Sub(D, r2)
+	D.Sub(D, r.t)
+	D.Mul(D, r.t, pool)
+
+	H := newGFp2(pool).Sub(B, r.x)
+	I := newGFp2(pool).Square(H, pool)
+
+	E := newGFp2(pool).Add(I, I)
+	E.Add(E, E)
+
+	J := newGFp2(pool).Mul(H, E, pool)
+
+	L1 := newGFp2(pool).Sub(D, r.y)
+	L1.Sub(L1, r.y)
+
+	V := newGFp2(pool).Mul(r.x, E, pool)
+
+	rOut = newTwistPoint(pool)
+	rOut.x.Square(L1, pool)
+	rOut.x.Sub(rOut.x, J)
+	rOut.x.Sub(rOut.x, V)
+	rOut.x.Sub(rOut.x, V)
+
+	rOut.z.Add(r.z, H)
+	rOut.z.Square(rOut.z, pool)
+	rOut.z.Sub(rOut.z, r.t)
+	rOut.z.Sub(rOut.z, I)
+
+	t := newGFp2(pool).Sub(V, rOut.x)
+	t.Mul(t, L1, pool)
+	t2 := newGFp2(pool).Mul(r.y, J, pool)
+	t2.Add(t2, t2)
+	rOut.y.Sub(t, t2)
+
+	rOut.t.Square(rOut.z, pool)
+
+	t.Add(p.y, rOut.z)
+	t.Square(t, pool)
+	t.Sub(t, r2)
+	t.Sub(t, rOut.t)
+
+	t2.Mul(L1, p.x, pool)
+	t2.Add(t2, t2)
+	a = newGFp2(pool)
+	a.Sub(t2, t)
+
+	c = newGFp2(pool)
+	c.MulScalar(rOut.z, q.y)
+	c.Add(c, c)
+
+	b = newGFp2(pool)
+	b.SetZero()
+	b.Sub(b, L1)
+	b.MulScalar(b, q.x)
+	b.Add(b, b)
+
+	B.Put(pool)
+	D.Put(pool)
+	H.Put(pool)
+	I.Put(pool)
+	E.Put(pool)
+	J.Put(pool)
+	L1.Put(pool)
+	V.Put(pool)
+	t.Put(pool)
+	t2.Put(pool)
+
+	return
+}
+
+func lineFunctionDouble(r *twistPoint, q *curvePoint, pool *bnPool) (a, b, c *gfP2, rOut *twistPoint) {
+	// See the doubling algorithm for a=0 from "Faster Computation of the
+	// Tate Pairing", http://arxiv.org/pdf/0904.0854v3.pdf
+
+	A := newGFp2(pool).Square(r.x, pool)
+	B := newGFp2(pool).Square(r.y, pool)
+	C_ := newGFp2(pool).Square(B, pool)
+
+	D := newGFp2(pool).Add(r.x, B)
+	D.Square(D, pool)
+	D.Sub(D, A)
+	D.Sub(D, C_)
+	D.Add(D, D)
+
+	E := newGFp2(pool).Add(A, A)
+	E.Add(E, A)
+
+	G := newGFp2(pool).Square(E, pool)
+
+	rOut = newTwistPoint(pool)
+	rOut.x.Sub(G, D)
+	rOut.x.Sub(rOut.x, D)
+
+	rOut.z.Add(r.y, r.z)
+	rOut.z.Square(rOut.z, pool)
+	rOut.z.Sub(rOut.z, B)
+	rOut.z.Sub(rOut.z, r.t)
+
+	rOut.y.Sub(D, rOut.x)
+	rOut.y.Mul(rOut.y, E, pool)
+	t := newGFp2(pool).Add(C_, C_)
+	t.Add(t, t)
+	t.Add(t, t)
+	rOut.y.Sub(rOut.y, t)
+
+	rOut.t.Square(rOut.z, pool)
+
+	t.Mul(E, r.t, pool)
+	t.Add(t, t)
+	b = newGFp2(pool)
+	b.SetZero()
+	b.Sub(b, t)
+	b.MulScalar(b, q.x)
+
+	a = newGFp2(pool)
+	a.Add(r.x, E)
+	a.Square(a, pool)
+	a.Sub(a, A)
+	a.Sub(a, G)
+	t.Add(B, B)
+	t.Add(t, t)
+	a.Sub(a, t)
+
+	c = newGFp2(pool)
+	c.Mul(rOut.z, r.t, pool)
+	c.Add(c, c)
+	c.MulScalar(c, q.y)
+
+	A.Put(pool)
+	B.Put(pool)
+	C_.Put(pool)
+	D.Put(pool)
+	E.Put(pool)
+	G.Put(pool)
+	t.Put(pool)
+
+	return
+}
+
+func mulLine(ret *gfP12, a, b, c *gfP2, pool *bnPool) {
+	a2 := newGFp6(pool)
+	a2.x.SetZero()
+	a2.y.Set(a)
+	a2.z.Set(b)
+	a2.Mul(a2, ret.x, pool)
+	t3 := newGFp6(pool).MulScalar(ret.y, c, pool)
+
+	t := newGFp2(pool)
+	t.Add(b, c)
+	t2 := newGFp6(pool)
+	t2.x.SetZero()
+	t2.y.Set(a)
+	t2.z.Set(t)
+	ret.x.Add(ret.x, ret.y)
+
+	ret.y.Set(t3)
+
+	ret.x.Mul(ret.x, t2, pool)
+	ret.x.Sub(ret.x, a2)
+	ret.x.Sub(ret.x, ret.y)
+	a2.MulTau(a2, pool)
+	ret.y.Add(ret.y, a2)
+
+	a2.Put(pool)
+	t3.Put(pool)
+	t2.Put(pool)
+	t.Put(pool)
+}
+
+// sixuPlus2NAF is 6u+2 in non-adjacent form.
+var sixuPlus2NAF = []int8{0, 0, 0, 1, 0, 1, 0, -1, 0, 0, 1, -1, 0, 0, 1, 0,
+	0, 1, 1, 0, -1, 0, 0, 1, 0, -1, 0, 0, 0, 0, 1, 1,
+	1, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 1,
+	1, 0, 0, -1, 0, 0, 0, 1, 1, 0, -1, 0, 0, 1, 0, 1, 1}
+
+// miller implements the Miller loop for calculating the Optimal Ate pairing.
+// See algorithm 1 from http://cryptojedi.org/papers/dclxvi-20100714.pdf
+func miller(q *twistPoint, p *curvePoint, pool *bnPool) *gfP12 {
+	ret := newGFp12(pool)
+	ret.SetOne()
+
+	aAffine := newTwistPoint(pool)
+	aAffine.Set(q)
+	aAffine.MakeAffine(pool)
+
+	bAffine := newCurvePoint(pool)
+	bAffine.Set(p)
+	bAffine.MakeAffine(pool)
+
+	minusA := newTwistPoint(pool)
+	minusA.Negative(aAffine, pool)
+
+	r := newTwistPoint(pool)
+	r.Set(aAffine)
+
+	r2 := newGFp2(pool)
+	r2.Square(aAffine.y, pool)
+
+	for i := len(sixuPlus2NAF) - 1; i > 0; i-- {
+		a, b, c, newR := lineFunctionDouble(r, bAffine, pool)
+		if i != len(sixuPlus2NAF)-1 {
+			ret.Square(ret, pool)
+		}
+
+		mulLine(ret, a, b, c, pool)
+		a.Put(pool)
+		b.Put(pool)
+		c.Put(pool)
+		r.Put(pool)
+		r = newR
+
+		switch sixuPlus2NAF[i-1] {
+		case 1:
+			a, b, c, newR = lineFunctionAdd(r, aAffine, bAffine, r2, pool)
+		case -1:
+			a, b, c, newR = lineFunctionAdd(r, minusA, bAffine, r2, pool)
+		default:
+			continue
+		}
+
+		mulLine(ret, a, b, c, pool)
+		a.Put(pool)
+		b.Put(pool)
+		c.Put(pool)
+		r.Put(pool)
+		r = newR
+	}
+
+	// In order to calculate Q1 we have to convert q from the sextic twist
+	// to the full GF(p^12) group, apply the Frobenius there, and convert
+	// back.
+	//
+	// The twist isomorphism is (x', y') -> (xω², yω³). If we consider just
+	// x for a moment, then after applying the Frobenius, we have x̄ω^(2p)
+	// where x̄ is the conjugate of x. If we are going to apply the inverse
+	// isomorphism we need a value with a single coefficient of ω² so we
+	// rewrite this as x̄ω^(2p-2)ω². ξ⁶ = ω and, due to the construction of
+	// p, 2p-2 is a multiple of six. Therefore we can rewrite as
+	// x̄ξ^((p-1)/3)ω² and applying the inverse isomorphism eliminates the
+	// ω².
+	//
+	// A similar argument can be made for the y value.
+
+	q1 := newTwistPoint(pool)
+	q1.x.Conjugate(aAffine.x)
+	q1.x.Mul(q1.x, xiToPMinus1Over3, pool)
+	q1.y.Conjugate(aAffine.y)
+	q1.y.Mul(q1.y, xiToPMinus1Over2, pool)
+	q1.z.SetOne()
+	q1.t.SetOne()
+
+	// For Q2 we are applying the p² Frobenius. The two conjugations cancel
+	// out and we are left only with the factors from the isomorphism. In
+	// the case of x, we end up with a pure number which is why
+	// xiToPSquaredMinus1Over3 is ∈ GF(p). With y we get a factor of -1. We
+	// ignore this to end up with -Q2.
+
+	minusQ2 := newTwistPoint(pool)
+	minusQ2.x.MulScalar(aAffine.x, xiToPSquaredMinus1Over3)
+	minusQ2.y.Set(aAffine.y)
+	minusQ2.z.SetOne()
+	minusQ2.t.SetOne()
+
+	r2.Square(q1.y, pool)
+	a, b, c, newR := lineFunctionAdd(r, q1, bAffine, r2, pool)
+	mulLine(ret, a, b, c, pool)
+	a.Put(pool)
+	b.Put(pool)
+	c.Put(pool)
+	r.Put(pool)
+	r = newR
+
+	r2.Square(minusQ2.y, pool)
+	a, b, c, newR = lineFunctionAdd(r, minusQ2, bAffine, r2, pool)
+	mulLine(ret, a, b, c, pool)
+	a.Put(pool)
+	b.Put(pool)
+	c.Put(pool)
+	r.Put(pool)
+	r = newR
+
+	aAffine.Put(pool)
+	bAffine.Put(pool)
+	minusA.Put(pool)
+	r.Put(pool)
+	r2.Put(pool)
+
+	return ret
+}
+
+// finalExponentiation computes the (p¹²-1)/Order-th power of an element of
+// GF(p¹²) to obtain an element of GT (steps 13-15 of algorithm 1 from
+// http://cryptojedi.org/papers/dclxvi-20100714.pdf)
+func finalExponentiation(in *gfP12, pool *bnPool) *gfP12 {
+	t1 := newGFp12(pool)
+
+	// This is the p^6-Frobenius
+	t1.x.Negative(in.x)
+	t1.y.Set(in.y)
+
+	inv := newGFp12(pool)
+	inv.Invert(in, pool)
+	t1.Mul(t1, inv, pool)
+
+	t2 := newGFp12(pool).FrobeniusP2(t1, pool)
+	t1.Mul(t1, t2, pool)
+
+	fp := newGFp12(pool).Frobenius(t1, pool)
+	fp2 := newGFp12(pool).FrobeniusP2(t1, pool)
+	fp3 := newGFp12(pool).Frobenius(fp2, pool)
+
+	fu, fu2, fu3 := newGFp12(pool), newGFp12(pool), newGFp12(pool)
+	fu.Exp(t1, u, pool)
+	fu2.Exp(fu, u, pool)
+	fu3.Exp(fu2, u, pool)
+
+	y3 := newGFp12(pool).Frobenius(fu, pool)
+	fu2p := newGFp12(pool).Frobenius(fu2, pool)
+	fu3p := newGFp12(pool).Frobenius(fu3, pool)
+	y2 := newGFp12(pool).FrobeniusP2(fu2, pool)
+
+	y0 := newGFp12(pool)
+	y0.Mul(fp, fp2, pool)
+	y0.Mul(y0, fp3, pool)
+
+	y1, y4, y5 := newGFp12(pool), newGFp12(pool), newGFp12(pool)
+	y1.Conjugate(t1)
+	y5.Conjugate(fu2)
+	y3.Conjugate(y3)
+	y4.Mul(fu, fu2p, pool)
+	y4.Conjugate(y4)
+
+	y6 := newGFp12(pool)
+	y6.Mul(fu3, fu3p, pool)
+	y6.Conjugate(y6)
+
+	t0 := newGFp12(pool)
+	t0.Square(y6, pool)
+	t0.Mul(t0, y4, pool)
+	t0.Mul(t0, y5, pool)
+	t1.Mul(y3, y5, pool)
+	t1.Mul(t1, t0, pool)
+	t0.Mul(t0, y2, pool)
+	t1.Square(t1, pool)
+	t1.Mul(t1, t0, pool)
+	t1.Square(t1, pool)
+	t0.Mul(t1, y1, pool)
+	t1.Mul(t1, y0, pool)
+	t0.Square(t0, pool)
+	t0.Mul(t0, t1, pool)
+
+	inv.Put(pool)
+	t1.Put(pool)
+	t2.Put(pool)
+	fp.Put(pool)
+	fp2.Put(pool)
+	fp3.Put(pool)
+	fu.Put(pool)
+	fu2.Put(pool)
+	fu3.Put(pool)
+	fu2p.Put(pool)
+	fu3p.Put(pool)
+	y0.Put(pool)
+	y1.Put(pool)
+	y2.Put(pool)
+	y3.Put(pool)
+	y4.Put(pool)
+	y5.Put(pool)
+	y6.Put(pool)
+
+	return t0
+}
+
+func optimalAte(a *twistPoint, b *curvePoint, pool *bnPool) *gfP12 {
+	e := miller(a, b, pool)
+	ret := finalExponentiation(e, pool)
+	e.Put(pool)
+
+	if a.IsInfinity() || b.IsInfinity() {
+		ret.SetOne()
+	}
+
+	return ret
+}
diff --git a/crypto/bn256/twist.go b/crypto/bn256/twist.go
new file mode 100644
index 000000000..95b966e2e
--- /dev/null
+++ b/crypto/bn256/twist.go
@@ -0,0 +1,249 @@
+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package bn256
+
+import (
+	"math/big"
+)
+
+// twistPoint implements the elliptic curve y²=x³+3/ξ over GF(p²). Points are
+// kept in Jacobian form and t=z² when valid. The group G₂ is the set of
+// n-torsion points of this curve over GF(p²) (where n = Order)
+type twistPoint struct {
+	x, y, z, t *gfP2
+}
+
+var twistB = &gfP2{
+	bigFromBase10("266929791119991161246907387137283842545076965332900288569378510910307636690"),
+	bigFromBase10("19485874751759354771024239261021720505790618469301721065564631296452457478373"),
+}
+
+// twistGen is the generator of group G₂.
+var twistGen = &twistPoint{
+	&gfP2{
+		bigFromBase10("11559732032986387107991004021392285783925812861821192530917403151452391805634"),
+		bigFromBase10("10857046999023057135944570762232829481370756359578518086990519993285655852781"),
+	},
+	&gfP2{
+		bigFromBase10("4082367875863433681332203403145435568316851327593401208105741076214120093531"),
+		bigFromBase10("8495653923123431417604973247489272438418190587263600148770280649306958101930"),
+	},
+	&gfP2{
+		bigFromBase10("0"),
+		bigFromBase10("1"),
+	},
+	&gfP2{
+		bigFromBase10("0"),
+		bigFromBase10("1"),
+	},
+}
+
+func newTwistPoint(pool *bnPool) *twistPoint {
+	return &twistPoint{
+		newGFp2(pool),
+		newGFp2(pool),
+		newGFp2(pool),
+		newGFp2(pool),
+	}
+}
+
+func (c *twistPoint) String() string {
+	return "(" + c.x.String() + ", " + c.y.String() + ", " + c.z.String() + ")"
+}
+
+func (c *twistPoint) Put(pool *bnPool) {
+	c.x.Put(pool)
+	c.y.Put(pool)
+	c.z.Put(pool)
+	c.t.Put(pool)
+}
+
+func (c *twistPoint) Set(a *twistPoint) {
+	c.x.Set(a.x)
+	c.y.Set(a.y)
+	c.z.Set(a.z)
+	c.t.Set(a.t)
+}
+
+// IsOnCurve returns true iff c is on the curve where c must be in affine form.
+func (c *twistPoint) IsOnCurve() bool {
+	pool := new(bnPool)
+	yy := newGFp2(pool).Square(c.y, pool)
+	xxx := newGFp2(pool).Square(c.x, pool)
+	xxx.Mul(xxx, c.x, pool)
+	yy.Sub(yy, xxx)
+	yy.Sub(yy, twistB)
+	yy.Minimal()
+	return yy.x.Sign() == 0 && yy.y.Sign() == 0
+}
+
+func (c *twistPoint) SetInfinity() {
+	c.z.SetZero()
+}
+
+func (c *twistPoint) IsInfinity() bool {
+	return c.z.IsZero()
+}
+
+func (c *twistPoint) Add(a, b *twistPoint, pool *bnPool) {
+	// For additional comments, see the same function in curve.go.
+
+	if a.IsInfinity() {
+		c.Set(b)
+		return
+	}
+	if b.IsInfinity() {
+		c.Set(a)
+		return
+	}
+
+	// See http://hyperelliptic.org/EFD/g1p/auto-code/shortw/jacobian-0/addition/add-2007-bl.op3
+	z1z1 := newGFp2(pool).Square(a.z, pool)
+	z2z2 := newGFp2(pool).Square(b.z, pool)
+	u1 := newGFp2(pool).Mul(a.x, z2z2, pool)
+	u2 := newGFp2(pool).Mul(b.x, z1z1, pool)
+
+	t := newGFp2(pool).Mul(b.z, z2z2, pool)
+	s1 := newGFp2(pool).Mul(a.y, t, pool)
+
+	t.Mul(a.z, z1z1, pool)
+	s2 := newGFp2(pool).Mul(b.y, t, pool)
+
+	h := newGFp2(pool).Sub(u2, u1)
+	xEqual := h.IsZero()
+
+	t.Add(h, h)
+	i := newGFp2(pool).Square(t, pool)
+	j := newGFp2(pool).Mul(h, i, pool)
+
+	t.Sub(s2, s1)
+	yEqual := t.IsZero()
+	if xEqual && yEqual {
+		c.Double(a, pool)
+		return
+	}
+	r := newGFp2(pool).Add(t, t)
+
+	v := newGFp2(pool).Mul(u1, i, pool)
+
+	t4 := newGFp2(pool).Square(r, pool)
+	t.Add(v, v)
+	t6 := newGFp2(pool).Sub(t4, j)
+	c.x.Sub(t6, t)
+
+	t.Sub(v, c.x)       // t7
+	t4.Mul(s1, j, pool) // t8
+	t6.Add(t4, t4)      // t9
+	t4.Mul(r, t, pool)  // t10
+	c.y.Sub(t4, t6)
+
+	t.Add(a.z, b.z)    // t11
+	t4.Square(t, pool) // t12
+	t.Sub(t4, z1z1)    // t13
+	t4.Sub(t, z2z2)    // t14
+	c.z.Mul(t4, h, pool)
+
+	z1z1.Put(pool)
+	z2z2.Put(pool)
+	u1.Put(pool)
+	u2.Put(pool)
+	t.Put(pool)
+	s1.Put(pool)
+	s2.Put(pool)
+	h.Put(pool)
+	i.Put(pool)
+	j.Put(pool)
+	r.Put(pool)
+	v.Put(pool)
+	t4.Put(pool)
+	t6.Put(pool)
+}
+
+func (c *twistPoint) Double(a *twistPoint, pool *bnPool) {
+	// See http://hyperelliptic.org/EFD/g1p/auto-code/shortw/jacobian-0/doubling/dbl-2009-l.op3
+	A := newGFp2(pool).Square(a.x, pool)
+	B := newGFp2(pool).Square(a.y, pool)
+	C_ := newGFp2(pool).Square(B, pool)
+
+	t := newGFp2(pool).Add(a.x, B)
+	t2 := newGFp2(pool).Square(t, pool)
+	t.Sub(t2, A)
+	t2.Sub(t, C_)
+	d := newGFp2(pool).Add(t2, t2)
+	t.Add(A, A)
+	e := newGFp2(pool).Add(t, A)
+	f := newGFp2(pool).Square(e, pool)
+
+	t.Add(d, d)
+	c.x.Sub(f, t)
+
+	t.Add(C_, C_)
+	t2.Add(t, t)
+	t.Add(t2, t2)
+	c.y.Sub(d, c.x)
+	t2.Mul(e, c.y, pool)
+	c.y.Sub(t2, t)
+
+	t.Mul(a.y, a.z, pool)
+	c.z.Add(t, t)
+
+	A.Put(pool)
+	B.Put(pool)
+	C_.Put(pool)
+	t.Put(pool)
+	t2.Put(pool)
+	d.Put(pool)
+	e.Put(pool)
+	f.Put(pool)
+}
+
+func (c *twistPoint) Mul(a *twistPoint, scalar *big.Int, pool *bnPool) *twistPoint {
+	sum := newTwistPoint(pool)
+	sum.SetInfinity()
+	t := newTwistPoint(pool)
+
+	for i := scalar.BitLen(); i >= 0; i-- {
+		t.Double(sum, pool)
+		if scalar.Bit(i) != 0 {
+			sum.Add(t, a, pool)
+		} else {
+			sum.Set(t)
+		}
+	}
+
+	c.Set(sum)
+	sum.Put(pool)
+	t.Put(pool)
+	return c
+}
+
+func (c *twistPoint) MakeAffine(pool *bnPool) *twistPoint {
+	if c.z.IsOne() {
+		return c
+	}
+
+	zInv := newGFp2(pool).Invert(c.z, pool)
+	t := newGFp2(pool).Mul(c.y, zInv, pool)
+	zInv2 := newGFp2(pool).Square(zInv, pool)
+	c.y.Mul(t, zInv2, pool)
+	t.Mul(c.x, zInv2, pool)
+	c.x.Set(t)
+	c.z.SetOne()
+	c.t.SetOne()
+
+	zInv.Put(pool)
+	t.Put(pool)
+	zInv2.Put(pool)
+
+	return c
+}
+
+func (c *twistPoint) Negative(a *twistPoint, pool *bnPool) {
+	c.x.Set(a.x)
+	c.y.SetZero()
+	c.y.Sub(c.y, a.y)
+	c.z.Set(a.z)
+	c.t.SetZero()
+}