diff --git a/sm9/bn256/bn_pair_b6.go b/sm9/bn256/bn_pair_b6.go
new file mode 100644
index 0000000..c016783
--- /dev/null
+++ b/sm9/bn256/bn_pair_b6.go
@@ -0,0 +1,66 @@
+package bn256
+
+func pairingB6(a *twistPoint, b *curvePoint) *gfP12 {
+	e := miller(a, b)
+	eb6 := (&gfP12b6{}).SetGfP12(e)
+	ret := finalExponentiationB6(eb6)
+
+	if a.IsInfinity() || b.IsInfinity() {
+		ret.SetOne()
+	}
+	return ret.ToGfP12()
+}
+
+// finalExponentiation computes the (p¹²-1)/Order-th power of an element of
+// GF(p¹²) to obtain an element of GT. https://eprint.iacr.org/2007/390.pdf
+// http://cryptojedi.org/papers/dclxvi-20100714.pdf
+func finalExponentiationB6(in *gfP12b6) *gfP12b6 {
+	t1 := &gfP12b6{}
+
+	// This is the p^6-Frobenius
+	t1.x.Neg(&in.x)
+	t1.y.Set(&in.y)
+
+	inv := &gfP12b6{}
+	inv.Invert(in)
+	t1.Mul(t1, inv)
+
+	t2 := (&gfP12b6{}).FrobeniusP2(t1)
+	t1.Mul(t1, t2)
+
+	fp := (&gfP12b6{}).Frobenius(t1)
+	fp2 := (&gfP12b6{}).FrobeniusP2(t1)
+	fp3 := (&gfP12b6{}).Frobenius(fp2)
+
+	fu := (&gfP12b6{}).Exp(t1, u)
+	fu2 := (&gfP12b6{}).Exp(fu, u)
+	fu3 := (&gfP12b6{}).Exp(fu2, u)
+
+	y3 := (&gfP12b6{}).Frobenius(fu)
+	fu2p := (&gfP12b6{}).Frobenius(fu2)
+	fu3p := (&gfP12b6{}).Frobenius(fu3)
+	y2 := (&gfP12b6{}).FrobeniusP2(fu2)
+
+	y0 := &gfP12b6{}
+	y0.Mul(fp, fp2).Mul(y0, fp3)
+
+	y1 := (&gfP12b6{}).Conjugate(t1)
+	y5 := (&gfP12b6{}).Conjugate(fu2)
+	y3.Conjugate(y3)
+	y4 := (&gfP12b6{}).Mul(fu, fu2p)
+	y4.Conjugate(y4)
+
+	y6 := (&gfP12b6{}).Mul(fu3, fu3p)
+	y6.Conjugate(y6)
+
+	t0 := (&gfP12b6{}).Square(y6)
+	t0.Mul(t0, y4).Mul(t0, y5)
+	t1.Mul(y3, y5).Mul(t1, t0)
+	t0.Mul(t0, y2)
+	t1.Square(t1).Mul(t1, t0).Square(t1)
+	t0.Mul(t1, y1)
+	t1.Mul(t1, y0)
+	t0.Square(t0).Mul(t0, t1)
+
+	return t0
+}
diff --git a/sm9/bn256/bn_pair_b6_test.go b/sm9/bn256/bn_pair_b6_test.go
new file mode 100644
index 0000000..6427114
--- /dev/null
+++ b/sm9/bn256/bn_pair_b6_test.go
@@ -0,0 +1,91 @@
+package bn256
+
+import (
+	"math/big"
+	"testing"
+)
+
+func Test_finalExponentiationB6(t *testing.T) {
+	x := &gfP12b6{
+		p6,
+		p6,
+	}
+	got := finalExponentiationB6(x)
+
+	exp := new(big.Int).Exp(p, big.NewInt(12), nil)
+	exp.Sub(exp, big.NewInt(1))
+	exp.Div(exp, Order)
+	expected := (&gfP12b6{}).Exp(x, exp)
+
+	if *got != *expected {
+		t.Errorf("got %v, expected %v\n", got, expected)
+	}
+}
+
+func Test_PairingB6_A2(t *testing.T) {
+	pk := bigFromHex("0130E78459D78545CB54C587E02CF480CE0B66340F319F348A1D5B1F2DC5F4")
+	g2 := &G2{}
+	_, err := g2.ScalarBaseMult(NormalizeScalar(pk.Bytes()))
+	if err != nil {
+		t.Fatal(err)
+	}
+	ret := pairingB6(g2.p, curveGen)
+	if *ret != *expected1 {
+		t.Errorf("not expected")
+	}
+}
+
+func Test_PairingB6_B2(t *testing.T) {
+	deB := &twistPoint{}
+	deB.x.x = *fromBigInt(bigFromHex("74CCC3AC9C383C60AF083972B96D05C75F12C8907D128A17ADAFBAB8C5A4ACF7"))
+	deB.x.y = *fromBigInt(bigFromHex("01092FF4DE89362670C21711B6DBE52DCD5F8E40C6654B3DECE573C2AB3D29B2"))
+	deB.y.x = *fromBigInt(bigFromHex("44B0294AA04290E1524FF3E3DA8CFD432BB64DE3A8040B5B88D1B5FC86A4EBC1"))
+	deB.y.y = *fromBigInt(bigFromHex("8CFC48FB4FF37F1E27727464F3C34E2153861AD08E972D1625FC1A7BD18D5539"))
+	deB.z.SetOne()
+	deB.t.SetOne()
+
+	rA := &curvePoint{}
+	rA.x = *fromBigInt(bigFromHex("7CBA5B19069EE66AA79D490413D11846B9BA76DD22567F809CF23B6D964BB265"))
+	rA.y = *fromBigInt(bigFromHex("A9760C99CB6F706343FED05637085864958D6C90902ABA7D405FBEDF7B781599"))
+	rA.z = *one
+	rA.t = *one
+
+	ret := pairingB6(deB, rA)
+	if ret.x != expected_b2.x || ret.y != expected_b2.y || ret.z != expected_b2.z {
+		t.Errorf("not expected")
+	}
+}
+
+func Test_PairingB6_B2_2(t *testing.T) {
+	pubE := &curvePoint{}
+	pubE.x = *fromBigInt(bigFromHex("9174542668E8F14AB273C0945C3690C66E5DD09678B86F734C4350567ED06283"))
+	pubE.y = *fromBigInt(bigFromHex("54E598C6BF749A3DACC9FFFEDD9DB6866C50457CFC7AA2A4AD65C3168FF74210"))
+	pubE.z = *one
+	pubE.t = *one
+
+	ret := pairingB6(twistGen, pubE)
+	ret.Exp(ret, bigFromHex("00018B98C44BEF9F8537FB7D071B2C928B3BC65BD3D69E1EEE213564905634FE"))
+	if ret.x != expected_b2_2.x || ret.y != expected_b2_2.y || ret.z != expected_b2_2.z {
+		t.Errorf("not expected")
+	}
+}
+
+func BenchmarkFinalExponentiationB6(b *testing.B) {
+	x := &gfP12b6{
+		p6,
+		p6,
+	}
+	exp := new(big.Int).Exp(p, big.NewInt(12), nil)
+	exp.Sub(exp, big.NewInt(1))
+	exp.Div(exp, Order)
+	expected := (&gfP12b6{}).Exp(x, exp)
+
+	b.ReportAllocs()
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		got := finalExponentiationB6(x)
+		if *got != *expected {
+			b.Errorf("got %v, expected %v\n", got, expected)
+		}
+	}
+}
diff --git a/sm9/bn256/constants.go b/sm9/bn256/constants.go
index 1a81df6..cef6ebb 100644
--- a/sm9/bn256/constants.go
+++ b/sm9/bn256/constants.go
@@ -85,3 +85,15 @@ var betaToNegP2Plus1Over3 = fromBigInt(bigFromHex("b640000002a3a6f0e303ab4ff2eb2
 
 // betaToNegP2Plus1Over2 = i^(-(p^2-1)/2)
 var betaToNegP2Plus1Over2 = fromBigInt(bigFromHex("b640000002a3a6f1d603ab4ff58ec74521f2934b1a7aeedbe56f9b27e351457c"))
+
+var sToPMinus1 = w2ToPMinus1
+
+var sTo2PMinus2 = w2ToP2Minus1
+
+var sToPSquaredMinus1 = w2ToP2Minus1
+
+var sTo2PSquaredMinus2 = betaToNegP2Plus1Over3
+
+var sToPMinus1Over2 = frobConstant
+
+var sToPSquaredMinus1Over2 = sToPMinus1
diff --git a/sm9/bn256/gfp12_b6.go b/sm9/bn256/gfp12_b6.go
new file mode 100644
index 0000000..762d444
--- /dev/null
+++ b/sm9/bn256/gfp12_b6.go
@@ -0,0 +1,259 @@
+package bn256
+
+import "math/big"
+
+// gfP12b6 implements the field of size p¹² as a quadratic extension of gfP6
+// where t²=s.
+type gfP12b6 struct {
+	x, y gfP6 // value is xt + y
+}
+
+func gfP12b6Decode(in *gfP12b6) *gfP12b6 {
+	out := &gfP12b6{}
+	out.x = *gfP6Decode(&in.x)
+	out.y = *gfP6Decode(&in.y)
+	return out
+}
+
+var gfP12b6Gen *gfP12b6 = &gfP12b6{
+	x: gfP6{
+		x: gfP2{
+			x: *fromBigInt(bigFromHex("5f23b1ce3ac438e768411e843fea2a9be192c39d7c3e6440eb2aeaa2823d010c")),
+			y: *fromBigInt(bigFromHex("61ebca2110d736bf0d4eba5c8c017781e2447d6d5edfdda6065c1ad6d376db4f")),
+		},
+		y: gfP2{
+			x: *fromBigInt(bigFromHex("3b0fc03d5711c93dff31b23d5bc78184c4e4ad66027f8f55c219536a54552cae")),
+			y: *fromBigInt(bigFromHex("254eb32dea84e64dfa196a2583564700074e1694c800c130290e1c8bdb9441aa")),
+		},
+		z: gfP2{
+			x: *fromBigInt(bigFromHex("a3eec3cd6a795be8671d686fd9c9271dd32d71f71d7bd3de24fb5abe38626c9c")),
+			y: *fromBigInt(bigFromHex("b101d668bfbf8ac8e546ccb8d6e1f9b89b988c0c238fb05e7b9c733c1f964b52")),
+		},
+	},
+	y: gfP6{
+		x: gfP2{
+			x: *fromBigInt(bigFromHex("2efe33f18332bb77282c24f00c10930f7e2a3e36c6c822c7487ab1a6229d91f3")),
+			y: *fromBigInt(bigFromHex("a6db7142e0ca24ae9ba7630e295a5ce7ed43ed38c0ce33e6346965f4dc5b5813")),
+		},
+		y: gfP2{
+			x: *fromBigInt(bigFromHex("9ee43c7e3740bcd8e9d6067a4cf3c571441e074b4573390cfea0bce10965b32b")),
+			y: *fromBigInt(bigFromHex("aa07010f9d42787cb0ebd9852fc780efb01ab631f2f10a180e06727b47ee6118")),
+		},
+		z: gfP2{
+			x: *fromBigInt(bigFromHex("6a5fed210720de5844e199bee3498d4d2a72158dbf514e31be7381e2bce90a00")),
+			y: *fromBigInt(bigFromHex("a0f422c35d7b6262796c802ec3f1370b9ef5d413e3176666b55d63ee8d7a8468")),
+		},
+	},
+}
+
+func (e *gfP12b6) String() string {
+	return "(" + e.x.String() + "," + e.y.String() + ")"
+}
+
+func (e *gfP12b6) ToGfP12() *gfP12 {
+	ret := &gfP12{}
+
+	ret.z.y.Set(&e.y.z)
+	ret.x.y.Set(&e.y.y)
+	ret.y.x.Set(&e.y.x)
+	ret.y.y.Set(&e.x.z)
+	ret.z.x.Set(&e.x.y)
+	ret.x.x.Set(&e.x.x)
+
+	return ret
+}
+
+func (e *gfP12b6) SetGfP12(a *gfP12) *gfP12b6 {
+	e.y.z.Set(&a.z.y)
+	e.y.y.Set(&a.x.y)
+	e.y.x.Set(&a.y.x)
+	e.x.z.Set(&a.y.y)
+	e.x.y.Set(&a.z.x)
+	e.x.x.Set(&a.x.x)
+
+	return e
+}
+
+func (e *gfP12b6) Set(a *gfP12b6) *gfP12b6 {
+	e.x.Set(&a.x)
+	e.y.Set(&a.y)
+	return e
+}
+
+func (e *gfP12b6) SetZero() *gfP12b6 {
+	e.x.SetZero()
+	e.y.SetZero()
+	return e
+}
+
+func (e *gfP12b6) SetOne() *gfP12b6 {
+	e.x.SetZero()
+	e.y.SetOne()
+	return e
+}
+
+func (e *gfP12b6) IsZero() bool {
+	return e.x.IsZero() && e.y.IsZero()
+}
+
+func (e *gfP12b6) IsOne() bool {
+	return e.x.IsZero() && e.y.IsOne()
+}
+
+func (e *gfP12b6) Neg(a *gfP12b6) *gfP12b6 {
+	e.x.Neg(&a.x)
+	e.y.Neg(&a.y)
+	return e
+}
+
+func (e *gfP12b6) Conjugate(a *gfP12b6) *gfP12b6 {
+	e.x.Neg(&a.x)
+	e.y.Set(&a.y)
+	return e
+}
+
+func (e *gfP12b6) Add(a, b *gfP12b6) *gfP12b6 {
+	e.x.Add(&a.x, &b.x)
+	e.y.Add(&a.y, &b.y)
+	return e
+}
+
+func (e *gfP12b6) Sub(a, b *gfP12b6) *gfP12b6 {
+	e.x.Sub(&a.x, &b.x)
+	e.y.Sub(&a.y, &b.y)
+	return e
+}
+
+func (e *gfP12b6) Mul(a, b *gfP12b6) *gfP12b6 {
+	// "Multiplication and Squaring on Pairing-Friendly Fields"
+	// Section 4, Karatsuba method.
+	// http://eprint.iacr.org/2006/471.pdf
+	//(a0+a1*t)(b0+b1*t)=c0+c1*t, where
+	//c0 = a0*b0 +a1*b1*s
+	//c1 = (a0 + a1)(b0 + b1) - a0*b0 - a1*b1 = a0*b1 + a1*b0
+	tx, ty, v0, v1 := &gfP6{}, &gfP6{}, &gfP6{}, &gfP6{}
+	v0.Mul(&a.y, &b.y)
+	v1.Mul(&a.x, &b.x)
+
+	tx.Add(&a.x, &a.y)
+	ty.Add(&b.x, &b.y)
+	tx.Mul(tx, ty)
+	tx.Sub(tx, v0)
+	tx.Sub(tx, v1)
+
+	ty.MulS(v1)
+	ty.Add(ty, v0)
+
+	e.x.Set(tx)
+	e.y.Set(ty)
+	return e
+}
+
+func (e *gfP12b6) MulScalar(a *gfP12b6, b *gfP6) *gfP12b6 {
+	e.x.Mul(&a.x, b)
+	e.y.Mul(&a.y, b)
+	return e
+}
+
+func (e *gfP12b6) MulGfP(a *gfP12b6, b *gfP) *gfP12b6 {
+	e.x.MulGfP(&a.x, b)
+	e.y.MulGfP(&a.y, b)
+	return e
+}
+
+func (e *gfP12b6) MulGfP2(a *gfP12b6, b *gfP2) *gfP12b6 {
+	e.x.MulScalar(&a.x, b)
+	e.y.MulScalar(&a.y, b)
+	return e
+}
+
+func (e *gfP12b6) Square(a *gfP12b6) *gfP12b6 {
+	// Complex squaring algorithm
+	// (xt+y)² = (x^2*s + y^2) + 2*x*y*t
+	tx, ty := &gfP6{}, &gfP6{}
+	tx.Square(&a.x).MulS(tx)
+	ty.Square(&a.y)
+	ty.Add(tx, ty)
+
+	tx.Mul(&a.x, &a.y)
+	tx.Add(tx, tx)
+
+	e.x.Set(tx)
+	e.y.Set(ty)
+	return e
+}
+
+func (c *gfP12b6) Exp(a *gfP12b6, power *big.Int) *gfP12b6 {
+	sum := (&gfP12b6{}).SetOne()
+	t := &gfP12b6{}
+
+	for i := power.BitLen() - 1; i >= 0; i-- {
+		t.Square(sum)
+		if power.Bit(i) != 0 {
+			sum.Mul(t, a)
+		} else {
+			sum.Set(t)
+		}
+	}
+
+	c.Set(sum)
+	return c
+}
+
+func (e *gfP12b6) Invert(a *gfP12b6) *gfP12b6 {
+	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
+	// ftp://136.206.11.249/pub/crypto/pairings.pdf
+
+	t0, t1 := &gfP6{}, &gfP6{}
+
+	t0.Mul(&a.y, &a.y)
+	t1.Mul(&a.x, &a.x).MulS(t1)
+	t0.Sub(t0, t1)
+	t0.Invert(t0)
+
+	e.x.Neg(&a.x)
+	e.y.Set(&a.y)
+	e.MulScalar(e, t0)
+
+	return e
+}
+
+// Frobenius computes (xt+y)^p
+// = x^p t^p + y^p
+// = x^p t^(p-1) t + y^p
+// = x^p s^((p-1)/2) t + y^p
+// sToPMinus1Over2
+func (e *gfP12b6) Frobenius(a *gfP12b6) *gfP12b6 {
+	e.x.Frobenius(&a.x)
+	e.y.Frobenius(&a.y)
+	e.x.MulGfP(&e.x, sToPMinus1Over2)
+	return e
+}
+
+// FrobeniusP2 computes (xt+y)^p² = x^p² t ·s^((p²-1)/2) + y^p²
+func (e *gfP12b6) FrobeniusP2(a *gfP12b6) *gfP12b6 {
+	e.x.FrobeniusP2(&a.x)
+	e.y.FrobeniusP2(&a.y)
+	e.x.MulGfP(&e.x, sToPSquaredMinus1Over2)
+	return e
+}
+
+func (e *gfP12b6) FrobeniusP4(a *gfP12b6) *gfP12b6 {
+	e.x.FrobeniusP4(&a.x)
+	e.y.FrobeniusP4(&a.y)
+	e.x.MulGfP(&e.x, sToPSquaredMinus1)
+	return e
+}
+
+func (e *gfP12b6) FrobeniusP6(a *gfP12b6) *gfP12b6 {
+	e.x.Neg(&a.x)
+	e.y.Set(&a.y)
+	return e
+}
+
+// Select sets q to p1 if cond == 1, and to p2 if cond == 0.
+func (q *gfP12b6) Select(p1, p2 *gfP12b6, cond int) *gfP12b6 {
+	q.x.Select(&p1.x, &p2.x, cond)
+	q.y.Select(&p1.y, &p2.y, cond)
+	return q
+}
diff --git a/sm9/bn256/gfp12_b6_test.go b/sm9/bn256/gfp12_b6_test.go
new file mode 100644
index 0000000..3b7a1dd
--- /dev/null
+++ b/sm9/bn256/gfp12_b6_test.go
@@ -0,0 +1,313 @@
+package bn256
+
+import (
+	"math/big"
+	"testing"
+)
+
+var p6 = gfP6{
+	gfP2{
+		*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+		*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+	},
+	gfP2{
+		*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+		*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+	},
+	gfP2{
+		*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+		*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+	},
+}
+
+func testGfP12b6Square(t *testing.T, x *gfP12b6) {
+	xmulx := &gfP12b6{}
+	xmulx.Mul(x, x)
+	xmulx = gfP12b6Decode(xmulx)
+
+	x2 := &gfP12b6{}
+	x2.Square(x)
+	x2 = gfP12b6Decode(x2)
+
+	if *xmulx != *x2 {
+		t.Errorf("xmulx=%v, x2=%v", xmulx, x2)
+	}
+}
+
+func Test_gfP12b6Square(t *testing.T) {
+	x := &gfP12b6{
+		p6,
+		p6,
+	}
+	testGfP12b6Square(t, x)
+
+	x = &gfP12b6{
+		p6,
+		*(&gfP6{}).SetOne(),
+	}
+	testGfP12b6Square(t, x)
+
+	x = &gfP12b6{
+		*(&gfP6{}).SetOne(),
+		p6,
+	}
+	testGfP12b6Square(t, x)
+
+	x = &gfP12b6{
+		*(&gfP6{}).SetZero(),
+		p6,
+	}
+	testGfP12b6Square(t, x)
+
+	x = &gfP12b6{
+		p6,
+		*(&gfP6{}).SetZero(),
+	}
+	testGfP12b6Square(t, x)
+}
+
+func testGfP12b6Invert(t *testing.T, x *gfP12b6) {
+	xInv := &gfP12b6{}
+	xInv.Invert(x)
+
+	y := &gfP12b6{}
+	y.Mul(x, xInv)
+	if !y.IsOne() {
+		t.Fail()
+	}
+}
+
+func TestToGfP12(t *testing.T) {
+	x := &gfP12b6{
+		p6,
+		p6,
+	}
+
+	p12 := x.ToGfP12()
+
+	x1 := &gfP12b6{}
+	x1.SetGfP12(p12)
+
+	if *x1 != *x {
+		t.Errorf("not same")
+	}
+
+	// after add
+	x2 := (&gfP12b6{}).Add(x, x)
+	p12_1 := (&gfP12{}).Add(p12, p12)
+	x3 := (&gfP12b6{}).SetGfP12(p12_1)
+	if *x2 != *x3 {
+		x3 = gfP12b6Decode(x3)
+		x2 = gfP12b6Decode(x2)
+		t.Errorf("not same after double, %v\n%v\n", x3, x2)
+	}
+
+	// after sub
+	x2 = (&gfP12b6{}).Sub(x, x)
+	p12_1 = (&gfP12{}).Sub(p12, p12)
+	x3 = (&gfP12b6{}).SetGfP12(p12_1)
+	if *x2 != *x3 {
+		x3 = gfP12b6Decode(x3)
+		x2 = gfP12b6Decode(x2)
+		t.Errorf("not same after double, %v\n%v\n", x3, x2)
+	}
+
+	// after neg
+	x2 = (&gfP12b6{}).Neg(x)
+	p12_1 = (&gfP12{}).Neg(p12)
+	x3 = (&gfP12b6{}).SetGfP12(p12_1)
+	if *x2 != *x3 {
+		x3 = gfP12b6Decode(x3)
+		x2 = gfP12b6Decode(x2)
+		t.Errorf("not same after double, %v\n%v\n", x3, x2)
+	}
+
+	// after mul gfp
+	x2.MulGfP(x, fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")))
+	p12_1.MulGFP(p12, fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")))
+	x3.SetGfP12(p12_1)
+	if *x2 != *x3 {
+		x3 = gfP12b6Decode(x3)
+		x2 = gfP12b6Decode(x2)
+		t.Errorf("not same after mul gfp, %v\n%v\n", x3, x2)
+	}
+
+	// after mul gfp2
+	gfp2 := &gfP2{
+		*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+		*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+	}
+
+	x2.MulGfP2(x, gfp2)
+	p12_1.MulGFP2(p12, gfp2)
+	x3.SetGfP12(p12_1)
+	if *x2 != *x3 {
+		x3 = gfP12b6Decode(x3)
+		x2 = gfP12b6Decode(x2)
+		t.Errorf("not same after mul gfp2, %v\n%v\n", x3, x2)
+	}
+
+	// after squre
+	x2.Square(x)
+	p12_1.Square(p12)
+	x3.SetGfP12(p12_1)
+	if *x2 != *x3 {
+		x3 = gfP12b6Decode(x3)
+		x2 = gfP12b6Decode(x2)
+		t.Errorf("not same after square, %v\n%v\n", x3, x2)
+	}
+
+	// after mul
+	x2.Mul(x, x)
+	p12_1.Mul(p12, p12)
+	x3.SetGfP12(p12_1)
+	if *x2 != *x3 {
+		x3 = gfP12b6Decode(x3)
+		x2 = gfP12b6Decode(x2)
+		t.Errorf("not same after mul, %v\n%v\n", x3, x2)
+	}
+}
+
+func Test_gfP12b6Invert(t *testing.T) {
+	x := &gfP12b6{
+		*(&gfP6{}).SetZero(),
+		p6,
+	}
+	testGfP12b6Invert(t, x)
+	x = &gfP12b6{
+		*(&gfP6{}).SetOne(),
+		p6,
+	}
+	testGfP12b6Invert(t, x)
+}
+
+func TestSToPMinus1Over2(t *testing.T) {
+	expected := &gfP2{}
+	expected.y.Set(fromBigInt(bigFromHex("3f23ea58e5720bdb843c6cfa9c08674947c5c86e0ddd04eda91d8354377b698b")))
+	expected.x.Set(zero)
+
+	s := &gfP6{}
+	s.SetS()
+	s.Exp(s, pMinus1Over2Big)
+	if !(s.x.IsZero() && s.y.IsZero() && s.z == *expected) {
+		s = gfP6Decode(s)
+		t.Errorf("not same as expected %v\n", s)
+	}
+}
+
+func Test_gfP12b6Frobenius(t *testing.T) {
+	x := &gfP12b6{
+		p6,
+		p6,
+	}
+	expected := &gfP12b6{}
+	expected.Exp(x, p)
+	got := &gfP12b6{}
+	got.Frobenius(x)
+	if *expected != *got {
+		t.Errorf("got %v, expected %v", got, expected)
+	}
+}
+
+func TestSToPSquaredMinus1Over2(t *testing.T) {
+	s := &gfP6{}
+	s.SetS()
+	p2 := new(big.Int).Mul(p, p)
+	p2 = new(big.Int).Sub(p2, big.NewInt(1))
+	p2.Rsh(p2, 1)
+	s.Exp(s, p2)
+
+	expected := &gfP2{}
+	expected.y.Set(fromBigInt(bigFromHex("0000000000000000f300000002a3a6f2780272354f8b78f4d5fc11967be65334")))
+	expected.x.Set(zero)
+
+	if !(s.x.IsZero() && s.y.IsZero() && s.z == *expected) {
+		s = gfP6Decode(s)
+		t.Errorf("not same as expected %v\n", s)
+	}
+}
+
+func Test_gfP12b6FrobeniusP2(t *testing.T) {
+	x := &gfP12b6{
+		p6,
+		p6,
+	}
+	expected := &gfP12b6{}
+	p2 := new(big.Int).Mul(p, p)
+	expected.Exp(x, p2)
+	got := &gfP12b6{}
+	got.FrobeniusP2(x)
+	if *expected != *got {
+		t.Errorf("got %v, expected %v", got, expected)
+	}
+}
+
+func TestSToP4Minus1Over2(t *testing.T) {
+	s := &gfP6{}
+	s.SetS()
+	p4 := new(big.Int).Mul(p, p)
+	p4.Mul(p4, p4)
+	p4 = new(big.Int).Sub(p4, big.NewInt(1))
+	p4.Rsh(p4, 1)
+	s.Exp(s, p4)
+
+	expected := &gfP2{}
+	expected.y.Set(fromBigInt(bigFromHex("0000000000000000f300000002a3a6f2780272354f8b78f4d5fc11967be65333")))
+	expected.x.Set(zero)
+
+	if !(s.x.IsZero() && s.y.IsZero() && s.z == *expected) {
+		s = gfP6Decode(s)
+		t.Errorf("not same as expected %v\n", s)
+	}
+}
+
+func Test_gfP12b6FrobeniusP4(t *testing.T) {
+	x := &gfP12b6{
+		p6,
+		p6,
+	}
+	expected := &gfP12b6{}
+	p4 := new(big.Int).Mul(p, p)
+	p4.Mul(p4, p4)
+	expected.Exp(x, p4)
+	got := &gfP12b6{}
+	got.FrobeniusP4(x)
+	if *expected != *got {
+		t.Errorf("got %v, expected %v", got, expected)
+	}
+}
+
+func Test_gfP12b6FrobeniusP6(t *testing.T) {
+	x := &gfP12b6{
+		p6,
+		p6,
+	}
+	expected := &gfP12b6{}
+	p6 := new(big.Int).Mul(p, p)
+	p6.Mul(p6, p)
+	p6.Mul(p6, p6)
+	expected.Exp(x, p6)
+	got := &gfP12b6{}
+	got.FrobeniusP6(x)
+	if *expected != *got {
+		t.Errorf("got %v, expected %v", got, expected)
+	}
+}
+
+func BenchmarkGfP12b6Frobenius(b *testing.B) {
+	x := &gfP12b6{
+		p6,
+		p6,
+	}
+	expected := &gfP12b6{}
+	expected.Exp(x, p)
+	got := &gfP12b6{}
+	b.ReportAllocs()
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		got.Frobenius(x)
+		if *expected != *got {
+			b.Errorf("got %v, expected %v", got, expected)
+		}
+	}
+}
diff --git a/sm9/bn256/gfp6.go b/sm9/bn256/gfp6.go
new file mode 100644
index 0000000..9ff3cdf
--- /dev/null
+++ b/sm9/bn256/gfp6.go
@@ -0,0 +1,285 @@
+package bn256
+
+import "math/big"
+
+// For details of the algorithms used, see "Multiplication and Squaring on
+// Pairing-Friendly Fields, Devegili et al.
+// http://eprint.iacr.org/2006/471.pdf.
+//
+
+// gfP6 implements the field of size p^6 as a cubic extension of gfP2
+// where s³=ξ and ξ=u
+type gfP6 struct {
+	x, y, z gfP2 // value is xs² + ys + z
+}
+
+func gfP6Decode(in *gfP6) *gfP6 {
+	out := &gfP6{}
+	out.x = *gfP2Decode(&in.x)
+	out.y = *gfP2Decode(&in.y)
+	out.z = *gfP2Decode(&in.z)
+	return out
+}
+
+func (e *gfP6) String() string {
+	return "(" + e.x.String() + ", " + e.y.String() + ", " + e.z.String() + ")"
+}
+
+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) SetS() *gfP6 {
+	e.x.SetZero()
+	e.y.SetOne()
+	e.z.SetZero()
+	return e
+}
+
+func (e *gfP6) SetS2() *gfP6 {
+	e.x.SetOne()
+	e.y.SetZero()
+	e.z.SetZero()
+	return e
+}
+
+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) Neg(a *gfP6) *gfP6 {
+	e.x.Neg(&a.x)
+	e.y.Neg(&a.y)
+	e.z.Neg(&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) MulScalar(a *gfP6, b *gfP2) *gfP6 {
+	e.x.Mul(&a.x, b)
+	e.y.Mul(&a.y, b)
+	e.z.Mul(&a.z, b)
+	return e
+}
+
+func (e *gfP6) MulGfP(a *gfP6, b *gfP) *gfP6 {
+	e.x.MulScalar(&a.x, b)
+	e.y.MulScalar(&a.y, b)
+	e.z.MulScalar(&a.z, b)
+	return e
+}
+
+func (e *gfP6) Mul(a, b *gfP6) *gfP6 {
+	// (z0 + y0*s + x0*s²)* (z1 + y1*s + x1*s²)
+	//  z0*z1 + z0*y1*s + z0*x1*s²
+	// +y0*z1*s + y0*y1*s² + y0*x1*u
+	// +x0*z1*s² + x0*y1*u + x0*x1*s*u
+	//=(z0*z1+y0*x1*u+x0*y1*u) + (z0*y1+y0*z1+x0*x1*u)s + (z0*x1 + y0*y1 + x0*z1)*s²
+	tx, ty, tz, t, v0, v1, v2 := &gfP2{}, &gfP2{}, &gfP2{}, &gfP2{}, &gfP2{}, &gfP2{}, &gfP2{}
+	v0.Mul(&a.z, &b.z)
+	v1.Mul(&a.y, &b.y)
+	v2.Mul(&a.x, &b.x)
+
+	t.Add(&a.y, &a.x)
+	tz.Add(&b.y, &b.x)
+	t.Mul(t, tz)
+	t.Sub(t, v1)
+	t.Sub(t, v2)
+	t.MulU1(t)
+	tz.Add(t, v0)
+
+	t.Add(&a.z, &a.y)
+	ty.Add(&b.z, &b.y)
+	ty.Mul(t, ty)
+	ty.Sub(ty, v0)
+	ty.Sub(ty, v1)
+	t.MulU1(v2)
+	ty.Add(ty, t)
+
+	t.Add(&a.z, &a.x)
+	tx.Add(&b.z, &b.x)
+	tx.Mul(tx, t)
+	tx.Sub(tx, v0)
+	tx.Add(tx, v1)
+	tx.Sub(tx, v2)
+
+	e.x.Set(tx)
+	e.y.Set(ty)
+	e.z.Set(tz)
+	return e
+}
+
+// MulS returns (z + y*s + x*s²)*s
+// = ys² + zs + xu
+func (e *gfP6) MulS(a *gfP6) *gfP6 {
+	ty := (&gfP2{}).Set(&a.y)
+	tz := &gfP2{}
+
+	tz.x.Set(&a.x.y)
+	gfpAdd(&tz.y, &a.x.x, &a.x.x)
+	gfpNeg(&tz.y, &tz.y)
+
+	e.y.Set(&a.z)
+	e.x.Set(ty)
+	e.z.Set(tz)
+
+	return e
+}
+
+func (e *gfP6) Square(a *gfP6) *gfP6 {
+	// (z + y*s + x*s²)* (z + y*s + x*s²)
+	// z^2 + z*y*s + z*x*s² + y*z*s + y^2*s² + y*x*u + x*z*s² + x*y*u + x^2 *u *s
+	// (z^2 + y*x*s + x*y*u) + (z*y + y*z + u * x^2)s + (z*x + y^2 + x*z)*s²
+	// (z^2 + 2*x*y*u) + (u*x^2 + 2*y*z) * s + (y^2 + 2*x*z) * s²
+	tx, ty, tz, t := &gfP2{}, &gfP2{}, &gfP2{}, &gfP2{}
+
+	tz.Square(&a.z)
+	t.MulU(&a.x, &a.y)
+	t.Add(t, t)
+	tz.Add(tz, t)
+
+	ty.SquareU(&a.x)
+	t.Mul(&a.y, &a.z)
+	t.Add(t, t)
+	ty.Add(ty, t)
+
+	tx.Square(&a.y)
+	t.Mul(&a.x, &a.z)
+	t.Add(t, t)
+	tx.Add(tx, t)
+
+	e.x.Set(tx)
+	e.y.Set(ty)
+	e.z.Set(tz)
+	return e
+}
+
+func (e *gfP6) Exp(f *gfP6, power *big.Int) *gfP6 {
+	sum := (&gfP6{}).SetOne()
+	t := &gfP6{}
+
+	for i := power.BitLen() - 1; i >= 0; i-- {
+		t.Square(sum)
+		if power.Bit(i) != 0 {
+			sum.Mul(t, f)
+		} else {
+			sum.Set(t)
+		}
+	}
+
+	e.Set(sum)
+	return e
+}
+
+func (e *gfP6) Invert(a *gfP6) *gfP6 {
+	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
+	// ftp://136.206.11.249/pub/crypto/pairings.pdf
+
+	t1 := (&gfP2{}).MulU(&a.x, &a.y)
+	A := (&gfP2{}).Square(&a.z)
+	A.Sub(A, t1)
+
+	B := (&gfP2{}).SquareU(&a.x)
+	t1.Mul(&a.y, &a.z)
+	B.Sub(B, t1)
+
+	C := (&gfP2{}).Square(&a.y)
+	t1.Mul(&a.x, &a.z)
+	C.Sub(C, t1)
+
+	F := (&gfP2{}).MulU(C, &a.y)
+	t1.Mul(A, &a.z)
+	F.Add(F, t1)
+	t1.MulU(B, &a.x)
+	F.Add(F, t1)
+
+	F.Invert(F)
+
+	e.x.Mul(C, F)
+	e.y.Mul(B, F)
+	e.z.Mul(A, F)
+	return e
+}
+
+// (z + y*s + x*s²)^p
+//= z^p + y^p*s*s^(p-1)+x^p*s²*(s²)^(p-1)
+//= z^p + (s^(p-1)*y^p)*s+((s²)^(p-1)*x^p)*s²
+//= f(z) + (s^(p-1)*f(y))*s+((s²)^(p-1)*f(x))*s²
+// sToPMinus1^3 = p - 1
+// sTo2PMinus2 = sToPMinus1 ^ 2
+func (e *gfP6) Frobenius(a *gfP6) *gfP6 {
+	e.z.Conjugate(&a.z)
+	e.y.Conjugate(&a.y)
+	e.x.Conjugate(&a.x)
+	e.y.MulScalar(&e.y, sToPMinus1)
+	e.x.MulScalar(&e.x, sTo2PMinus2)
+
+	return e
+}
+
+// FrobeniusP2 computes (xs²+ys+z)^(p²)
+// = z^p² + y^p²*s^p² + x^p²*s²^p²
+// = z + y*s^p² + x * s^(2p²)
+// = z + y*s*s^(p²-1) + x * s² * s^(2p² - 2)
+// = z + y*s*u^((p²-1)/3) + x * s² * u^((2p² - 2)/3)
+// uToPSquaredMinus1Over3 = sToPSquaredMinus1
+// s^(p²-1) = s^((p-1)(p+1)) = (s^(p-1))^(p+1) = sTo2PMinus2
+//
+// uToPSquaredMinus1Over3 = sTo2PSquaredMinus2
+// = sTo2PMinus2^2
+func (e *gfP6) FrobeniusP2(a *gfP6) *gfP6 {
+	e.x.MulScalar(&a.x, sTo2PSquaredMinus2)
+	e.y.MulScalar(&a.y, sToPSquaredMinus1)
+	e.z.Set(&a.z)
+	return e
+}
+
+// FrobeniusP4 computes (xs²+ys+z)^(p^4)
+func (e *gfP6) FrobeniusP4(a *gfP6) *gfP6 {
+	e.x.MulScalar(&a.x, sToPSquaredMinus1)
+	e.y.MulScalar(&a.y, sTo2PSquaredMinus2)
+	e.z.Set(&a.z)
+	return e
+}
+
+// Select sets q to p1 if cond == 1, and to p2 if cond == 0.
+func (q *gfP6) Select(p1, p2 *gfP6, cond int) *gfP6 {
+	q.x.Select(&p1.x, &p2.x, cond)
+	q.y.Select(&p1.y, &p2.y, cond)
+	q.z.Select(&p1.z, &p2.z, cond)
+	return q
+}
diff --git a/sm9/bn256/gfp6_test.go b/sm9/bn256/gfp6_test.go
new file mode 100644
index 0000000..d7a9199
--- /dev/null
+++ b/sm9/bn256/gfp6_test.go
@@ -0,0 +1,299 @@
+package bn256
+
+import (
+	"math/big"
+	"testing"
+)
+
+func TestMulS(t *testing.T) {
+	x := &gfP6{
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+	}
+	s := &gfP6{}
+	s.SetS()
+
+	xmuls1, xmuls2 := &gfP6{}, &gfP6{}
+	xmuls1.MulS(x)
+	xmuls1 = gfP6Decode(xmuls1)
+	xmuls2.Mul(x, s)
+	xmuls2 = gfP6Decode(xmuls2)
+
+	if *xmuls1 != *xmuls2 {
+		t.Errorf("xmulx=%v, x2=%v", xmuls1, xmuls2)
+	}
+}
+
+func testGfP6Square(t *testing.T, x *gfP6) {
+	xmulx := &gfP6{}
+	xmulx.Mul(x, x)
+	xmulx = gfP6Decode(xmulx)
+
+	x2 := &gfP6{}
+	x2.Square(x)
+	x2 = gfP6Decode(x2)
+
+	if xmulx.x != x2.x || xmulx.y != x2.y {
+		t.Errorf("xmulx=%v, x2=%v", xmulx, x2)
+	}
+}
+
+func Test_gfP6Square(t *testing.T) {
+	gfp2Zero := (&gfP2{}).SetZero()
+	x := &gfP6{
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+	}
+	testGfP6Square(t, x)
+	x = &gfP6{
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+		*gfp2Zero,
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+	}
+	testGfP6Square(t, x)
+}
+
+func testGfP6Invert(t *testing.T, x *gfP6) {
+	xInv := &gfP6{}
+	xInv.Invert(x)
+
+	y := &gfP6{}
+	y.Mul(x, xInv)
+	if !y.IsOne() {
+		t.Fail()
+	}
+}
+
+func Test_gfP6Invert(t *testing.T) {
+	gfp2Zero := (&gfP2{}).SetZero()
+	x := &gfP6{
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+	}
+	testGfP6Invert(t, x)
+
+	x = &gfP6{
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+		*gfp2Zero,
+	}
+	testGfP6Invert(t, x)
+
+	x = &gfP6{
+		*gfp2Zero,
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+	}
+	testGfP6Invert(t, x)
+
+	
+	x = &gfP6{
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+		*gfp2Zero,
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+	}
+	testGfP6Invert(t, x)
+}
+
+// sToPMinus1 = s^(p-1) = u ^ ((p-1) / 3)
+// sToPMinus1 ^ 3 = -1
+// sToPMinus1 = 0000000000000000f300000002a3a6f2780272354f8b78f4d5fc11967be65334
+func TestSToPMinus1(t *testing.T) {
+	expected := &gfP2{}
+	expected.y.Set(fromBigInt(bigFromHex("0000000000000000f300000002a3a6f2780272354f8b78f4d5fc11967be65334")))
+	expected.x.Set(zero)
+
+	s := &gfP6{}
+	s.SetS()
+	s.Exp(s, bigFromHex("b640000002a3a6f1d603ab4ff58ec74521f2934b1a7aeedbe56f9b27e351457c"))
+	if !(s.x.IsZero() && s.y.IsZero() && s.z == *expected) {
+		t.Error("not same as expected")
+	}
+}
+
+// s2ToPMinus1 = (s^2)^(p-1) = sToPMinus1 ^ 2
+// s2ToPMinus1 = sToPMinus1^2
+// s2ToPMinus1 = 0000000000000000f300000002a3a6f2780272354f8b78f4d5fc11967be65333
+func TestS2ToPMinus1(t *testing.T) {
+	expected := &gfP2{}
+	expected.y.Set(fromBigInt(bigFromHex("0000000000000000f300000002a3a6f2780272354f8b78f4d5fc11967be65333")))
+	expected.x.Set(zero)
+
+	s := &gfP6{}
+	s.SetS()
+	s.Square(s)
+	s.Exp(s, bigFromHex("b640000002a3a6f1d603ab4ff58ec74521f2934b1a7aeedbe56f9b27e351457c"))
+	if !(s.x.IsZero() && s.y.IsZero() && s.z == *expected) {
+		t.Error("not same as expected")
+	}
+
+	s2 := &gfP2{}
+	s2.y.Set(fromBigInt(bigFromHex("0000000000000000f300000002a3a6f2780272354f8b78f4d5fc11967be65334")))
+	s2.x.Set(zero)
+	s2.Square(s2)
+
+	if *s2 != *expected {
+		t.Errorf("not same as expected")
+	}
+}
+
+func Test_gfP6Frobenius(t *testing.T) {
+	x := &gfP6{
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+	}
+	expected := &gfP6{}
+	expected.Exp(x, p)
+	got := &gfP6{}
+	got.Frobenius(x)
+	if expected.x != got.x || expected.y != got.y || expected.z != got.z {
+		t.Errorf("got %v, expected %v", got, expected)
+	}
+}
+
+func TestSToPSquaredMinus1(t *testing.T) {
+	s := &gfP6{}
+	s.SetS()
+	p2 := new(big.Int).Mul(p, p)
+	p2 = new(big.Int).Sub(p2, big.NewInt(1))
+	s.Exp(s, p2)
+
+	expected := &gfP2{}
+	expected.y.Set(fromBigInt(bigFromHex("0000000000000000f300000002a3a6f2780272354f8b78f4d5fc11967be65333")))
+	expected.x.Set(zero)
+
+	if !(s.x.IsZero() && s.y.IsZero() && s.z == *expected) {
+		t.Error("not same as expected")
+	}
+}
+
+func TestSTo2PSquaredMinus2(t *testing.T) {
+	expected := &gfP2{}
+	expected.y.Set(fromBigInt(bigFromHex("b640000002a3a6f0e303ab4ff2eb2052a9f02115caef75e70f738991676af249")))
+	expected.x.Set(zero)
+
+	s2 := &gfP2{}
+	s2.y.Set(fromBigInt(bigFromHex("0000000000000000f300000002a3a6f2780272354f8b78f4d5fc11967be65333")))
+	s2.x.Set(zero)
+	s2.Square(s2)
+
+	if *s2 != *expected {
+		s2 = gfP2Decode(s2)
+		t.Errorf("not same as expected: %v", s2)
+	}
+}
+
+func Test_gfP6FrobeniusP2(t *testing.T) {
+	x := &gfP6{
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+	}
+	expected := &gfP6{}
+	p2 := new(big.Int).Mul(p, p)
+	expected.Exp(x, p2)
+	got := &gfP6{}
+	got.FrobeniusP2(x)
+	if expected.x != got.x || expected.y != got.y || expected.z != got.z {
+		t.Errorf("got %v, expected %v", got, expected)
+	}
+}
+
+func Test_gfP6FrobeniusP4(t *testing.T) {
+	x := &gfP6{
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+	}
+	expected := &gfP6{}
+	p4 := new(big.Int).Mul(p, p)
+	p4.Mul(p4, p4)
+	expected.Exp(x, p4)
+	got := &gfP6{}
+	got.FrobeniusP4(x)
+	if expected.x != got.x || expected.y != got.y || expected.z != got.z {
+		t.Errorf("got %v, expected %v", got, expected)
+	}
+}