diff --git a/sm9/bn256/bn_pair.go b/sm9/bn256/bn_pair.go
index de939db..3971cf6 100644
--- a/sm9/bn256/bn_pair.go
+++ b/sm9/bn256/bn_pair.go
@@ -98,19 +98,30 @@ func lineFunctionDouble(r *twistPoint, q *curvePoint) (a, b, c, d *gfP2, rOut *t
 }
 
 func mulLine(ret *gfP12, retDen *gfP4, a, b, c, d *gfP2) {
-	l := &gfP12{}
-	l.y.SetZero()
-	l.x.x.SetZero()
-	l.x.y.Set(b)
-	l.z.x.Set(c)
-	l.z.y.Set(a)
+	tx, ty, tz, t, bx, bz := &gfP4{}, &gfP4{}, &gfP4{}, &gfP4{}, &gfP4{}, &gfP4{}
+	bx.x.SetZero()
+	bx.y.Set(b)
+	bz.x.Set(c)
+	bz.y.Set(a)
 
-	ret.Mul(ret, l)
+	tz.Mul(&ret.z, bz)
+	t.MulV(&ret.y, bx)
+	tz.Add(tz, t)
 
-	lDen := &gfP4{}
-	lDen.x.Set(d)
-	lDen.y.SetZero()
-	retDen.Mul(retDen, lDen)
+	ty.Mul(&ret.y, bz)
+	t.MulV(&ret.x, bx)
+	ret.y.Add(ty, t)
+
+	tx.Mul(&ret.z, bx)
+	t.Mul(&ret.x, bz)
+	ret.x.Add(tx, t)
+
+	ret.z.Set(tz)
+
+	txD := &gfP2{}
+	txD.Mul(&retDen.y, d)
+	retDen.y.MulU(&retDen.x, d)
+	retDen.x.Set(txD)
 }
 
 //
@@ -160,6 +171,18 @@ func miller(q *twistPoint, p *curvePoint) *gfP12 {
 		mulLine(ret, retDen, a, b, c, d)
 		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*β^(-1/3), y*β^(-1/2)). If we consider just
+	// x for a moment, then after applying the Frobenius, we have x̄*β^(-p/3)
+	// 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 β^(-1/3) so we
+	// rewrite this as x̄*β^((-p+1)/3)*β^(-1/3).
+	//
+	// A similar argument can be made for the y value.
 	q1 := &twistPoint{}
 	q1.x.Conjugate(&aAffine.x)
 	q1.x.MulScalar(&q1.x, betaToNegPPlus1Over3)
@@ -185,7 +208,13 @@ func miller(q *twistPoint, p *curvePoint) *gfP12 {
 	a, b, c, d, _ = lineFunctionAdd(r, minusQ2, bAffine, r2)
 	mulLine(ret, retDen, a, b, c, d)
 
-	retDen.Invert(retDen)
+	//retDen.Invert(retDen)
+	t2, t3 := &gfP2{}, &gfP2{}
+	t3.SquareU(&retDen.x)
+	t3.Invert(t3)
+	t2.Mul(&retDen.x, t3)
+	retDen.x.Set(t2)
+
 	ret.MulScalar(ret, retDen)
 
 	return ret
diff --git a/sm9/bn256/gfp12_test.go b/sm9/bn256/gfp12_test.go
index df653a1..aa7b1f6 100644
--- a/sm9/bn256/gfp12_test.go
+++ b/sm9/bn256/gfp12_test.go
@@ -105,6 +105,7 @@ func Test_gfP12Invert(t *testing.T) {
 	}
 }
 
+// Generate wToPMinus1
 func Test_gfP12Frobenius_Case1(t *testing.T) {
 	expected := &gfP12{}
 	i := &gfP12{}
@@ -123,6 +124,7 @@ func Test_gfP12Frobenius_Case1(t *testing.T) {
 	}
 }
 
+// Generate w2ToPMinus1
 func Test_gfP12Frobenius_Case2(t *testing.T) {
 	expected := &gfP12{}
 	i := &gfP12{}
@@ -141,6 +143,7 @@ func Test_gfP12Frobenius_Case2(t *testing.T) {
 	}
 }
 
+// Generate wToP2Minus1
 func Test_gfP12FrobeniusP2_Case1(t *testing.T) {
 	expected := &gfP12{}
 	i := &gfP12{}
@@ -160,6 +163,7 @@ func Test_gfP12FrobeniusP2_Case1(t *testing.T) {
 	}
 }
 
+// Generate w2ToP2Minus1
 func Test_gfP12FrobeniusP2_Case2(t *testing.T) {
 	expected := &gfP12{}
 	i := &gfP12{}
@@ -179,6 +183,7 @@ func Test_gfP12FrobeniusP2_Case2(t *testing.T) {
 	}
 }
 
+// Generate wToP3Minus1
 func Test_gfP12FrobeniusP3_Case1(t *testing.T) {
 	expected := &gfP12{}
 	i := &gfP12{}
@@ -199,6 +204,8 @@ func Test_gfP12FrobeniusP3_Case1(t *testing.T) {
 	}
 }
 
+
+// Generate w2ToP3minus1
 func Test_gfP12FrobeniusP3_Case2(t *testing.T) {
 	expected := &gfP12{}
 	i := &gfP12{}
diff --git a/sm9/bn256/twist.go b/sm9/bn256/twist.go
index 01544fb..85046b7 100644
--- a/sm9/bn256/twist.go
+++ b/sm9/bn256/twist.go
@@ -218,7 +218,7 @@ func (c *twistPoint) Neg(a *twistPoint) {
 }
 
 // code logic is form https://github.com/guanzhi/GmSSL/blob/develop/src/sm9_alg.c
-// the value is not same as p*a
+// the value is not same as [p]a
 func (c *twistPoint) Frobenius(a *twistPoint) {
 	c.x.Conjugate(&a.x)
 	c.y.Conjugate(&a.y)