github.com/hellobchain/newcryptosm@v0.0.0-20221019060107-edb949a317e9/sm9/rate.go (about)

     1  package sm9
     2  
     3  func lineFunctionAdd(r, p *twistPoint, q *curvePoint, r2 *gfP2) (a, b, c *gfP2, rOut *twistPoint) {
     4  	// See the mixed addition algorithm from "Faster Computation of the
     5  	// Tate Pairing", http://arxiv.org/pdf/0904.0854v3.pdf
     6  	B := (&gfP2{}).Mul(&p.x, &r.t)
     7  	D := (&gfP2{}).Add(&p.y, &r.z)
     8  	D.Square(D)
     9  	D.Sub(D, r2)
    10  	D.Sub(D, &r.t)
    11  	D.Mul(D, &r.t)
    12  
    13  	H := (&gfP2{}).Sub(B, &r.x)
    14  	I := (&gfP2{}).Square(H)
    15  
    16  	E := (&gfP2{}).Add(I, I)
    17  	E.Add(E, E)
    18  
    19  	J := (&gfP2{}).Mul(H, E)
    20  
    21  	L1 := (&gfP2{}).Sub(D, &r.y)
    22  	L1.Sub(L1, &r.y)
    23  
    24  	V := (&gfP2{}).Mul(&r.x, E)
    25  
    26  	rOut = &twistPoint{} //siomn:2R part
    27  	rOut.x.Square(L1)
    28  	rOut.x.Sub(&rOut.x, J)
    29  	rOut.x.Sub(&rOut.x, V)
    30  	rOut.x.Sub(&rOut.x, V)
    31  
    32  	rOut.z.Add(&r.z, H)
    33  	rOut.z.Square(&rOut.z)
    34  	rOut.z.Sub(&rOut.z, &r.t)
    35  	rOut.z.Sub(&rOut.z, I)
    36  
    37  	t := (&gfP2{}).Sub(V, &rOut.x)
    38  	t.Mul(t, L1)
    39  	t2 := (&gfP2{}).Mul(&r.y, J)
    40  	t2.Add(t2, t2)
    41  	rOut.y.Sub(t, t2)
    42  	rOut.t.Square(&rOut.z) //R=P+Q end;
    43  
    44  	//cf version
    45  	t.Add(&p.y, &rOut.z).Square(t).Sub(t, r2).Sub(t, &rOut.t)
    46  
    47  	t2.Mul(L1, &p.x)
    48  	t2.Add(t2, t2)
    49  	a = (&gfP2{}).Sub(t2, t)
    50  
    51  	c = (&gfP2{}).MulScalar(&rOut.z, &q.y)
    52  	c.Add(c, c)
    53  	c.MulXi(c) //simon add1030;first test success;
    54  
    55  	b = (&gfP2{}).Neg(L1)
    56  	b.MulScalar(b, &q.x).Add(b, b)
    57  
    58  	return
    59  }
    60  
    61  func lineFunctionDouble(r *twistPoint, q *curvePoint) (a, b, c *gfP2, rOut *twistPoint) {
    62  	// See the doubling algorithm for a=0 from "Faster Computation of the
    63  	// Tate Pairing", http://arxiv.org/pdf/0904.0854v3.pdf
    64  
    65  	A := (&gfP2{}).Square(&r.x)
    66  	B := (&gfP2{}).Square(&r.y)
    67  	C := (&gfP2{}).Square(B)
    68  
    69  	D := (&gfP2{}).Add(&r.x, B)
    70  	D.Square(D)
    71  	D.Sub(D, A)
    72  	D.Sub(D, C)
    73  	D.Add(D, D)
    74  
    75  	E := (&gfP2{}).Add(A, A)
    76  	E.Add(E, A) //E = 3A
    77  
    78  	G := (&gfP2{}).Square(E) //G = E^2
    79  
    80  	rOut = &twistPoint{} //X3 = G-2D
    81  	rOut.x.Sub(G, D)
    82  	rOut.x.Sub(&rOut.x, D)
    83  
    84  	rOut.z.Add(&r.y, &r.z) //Z3 = (Y1+Z1)^2-B-T1;
    85  	rOut.z.Square(&rOut.z)
    86  	rOut.z.Sub(&rOut.z, B)
    87  	rOut.z.Sub(&rOut.z, &r.t)
    88  
    89  	rOut.y.Sub(D, &rOut.x) //Y3 = E(D-X3)-8C
    90  	rOut.y.Mul(&rOut.y, E)
    91  	t := (&gfP2{}).Add(C, C)
    92  	t.Add(t, t)
    93  	t.Add(t, t)
    94  	rOut.y.Sub(&rOut.y, t) //t =8C
    95  
    96  	rOut.t.Square(&rOut.z) // T3 = Z3^2; R = 2R end
    97  
    98  	//bn256 version//1030first test success;speed up 25%
    99  	t.Mul(E, &r.t) //
   100  	t.Add(t, t)
   101  	b = &gfP2{}
   102  	b.SetZero()
   103  	b.Sub(b, t)
   104  	b.MulScalar(b, &q.x) //b = -2E*T1*q.x;
   105  
   106  	a = &gfP2{}
   107  	a.Add(&r.x, E) //a = (X1 +E)^2 -A-G-4B; // a(W^3)
   108  	a.Square(a)
   109  	a.Sub(a, A)
   110  	a.Sub(a, G)
   111  	t.Add(B, B)
   112  	t.Add(t, t)
   113  	a.Sub(a, t)
   114  
   115  	c = &gfP2{}
   116  	c.Mul(&rOut.z, &r.t)
   117  	c.Add(c, c)
   118  	c.MulScalar(c, &q.y) //c = 2*Z3*T1*q.y
   119  	c.MulXi(c)
   120  
   121  	return
   122  }
   123  
   124  //simon: 0928change this mulLine;
   125  //0929version
   126  func mulLine(ret *gfP12, a, b, c *gfP2) {
   127  	a2 := &gfP6{}
   128  	a3 := &gfP6{}
   129  	a12 := &gfP12{}
   130  
   131  	a2.z.SetZero()
   132  	a2.y.Set(a)
   133  	a2.x.Set(b)
   134  
   135  	//c.MulXi(c, pool) //(0,0,c)
   136  	a3.z.Set(c)
   137  	a3.y.SetZero()
   138  	a3.x.SetZero()
   139  
   140  	a12.x.Set(a2) //test
   141  	a12.y.Set(a3)
   142  	ret.Mul(ret, a12)
   143  
   144  }
   145  
   146  // sixuPlus2NAF is 6u+2 in non-adjacent form.//NAF is checked again in 0918;but the pair result is wrong!
   147  //var sixuPlus2NAF1 = []int8{0, -1, 0, 0, 0, 0, 1, 0, 1, 0, 0, -1, 0, -1, 0,
   148  //	0, 0, -1, 0, -1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   149  //	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}
   150  
   151  var sixuPlus2NAF = [66]int8{0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1,
   152  	1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
   153  	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   154  	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}
   155  
   156  // miller implements the Miller loop for calculating the Optimal Ate pairing.
   157  // See algorithm 1 from http://cryptojedi.org/papers/dclxvi-20100714.pdf
   158  func miller(q *twistPoint, p *curvePoint) *gfP12 {
   159  	ret := &gfP12{} //ret = f
   160  	ret.SetOne()
   161  	//fmt.Println("In miller:")
   162  
   163  	aAffine := &twistPoint{}
   164  	aAffine.Set(q)
   165  	aAffine.MakeAffine()
   166  
   167  	bAffine := &curvePoint{}
   168  	bAffine.Set(p)
   169  	bAffine.MakeAffine()
   170  
   171  	minusA := &twistPoint{}
   172  	minusA.Neg(aAffine)
   173  
   174  	r := &twistPoint{} // R<--Q
   175  	r.Set(aAffine)
   176  
   177  	r2 := (&gfP2{})
   178  	r2.Square(&aAffine.y) //r2 <-- Q.y^2
   179  	//fmt.Println("loop:")
   180  	for i := len(sixuPlus2NAF) - 2; i >= 0; i-- {
   181  		a, b, c, newR := lineFunctionDouble(r, bAffine) // lR,R(P); R <--2R
   182  
   183  		if i != len(sixuPlus2NAF)-1 {
   184  			ret.Square(ret)
   185  		}
   186  
   187  		mulLine(ret, a, b, c)
   188  
   189  		r = newR
   190  
   191  		switch sixuPlus2NAF[i] {
   192  		case 1:
   193  			a, b, c, newR = lineFunctionAdd(r, aAffine, bAffine, r2)
   194  		case -1:
   195  			a, b, c, newR = lineFunctionAdd(r, minusA, bAffine, r2)
   196  		default:
   197  			continue
   198  		}
   199  
   200  		mulLine(ret, a, b, c)
   201  		r = newR
   202  	}
   203  
   204  	// In order to calculate Q1 we have to convert q from the sextic twist
   205  	// to the full GF(p^12) group, apply the Frobenius there, and convert
   206  	// back.
   207  	//
   208  	// The twist isomorphism is (x', y') -> (xω², yω³). If we consider just
   209  	// x for a moment, then after applying the Frobenius, we have x̄ω^(2p)
   210  	// where x̄ is the conjugate of x. If we are going to apply the inverse
   211  	// isomorphism we need a value with a single coefficient of ω² so we
   212  	// rewrite this as x̄ω^(2p-2)ω². ξ⁶ = ω and, due to the construction of
   213  	// p, 2p-2 is a multiple of six. Therefore we can rewrite as
   214  	// x̄ξ^((p-1)/3)ω² and applying the inverse isomorphism eliminates the
   215  	// ω².
   216  	//
   217  	// A similar argument can be made for the y value.
   218  
   219  	q1 := &twistPoint{}
   220  	q1.x.Conjugate(&aAffine.x) //q1.x = 共轭(q.x)
   221  	q1.x.Mul(&q1.x, xiToPMinus1Over3i)
   222  	q1.y.Conjugate(&aAffine.y)
   223  	q1.y.Mul(&q1.y, xiToPMinus1Over2i)
   224  	q1.z.SetOne()
   225  	q1.t.SetOne() //simon:0926 test pass;
   226  
   227  	// For Q2 we are applying the p² Frobenius. The two conjugations cancel
   228  	// out and we are left only with the factors from the isomorphism. In
   229  	// the case of x, we end up with a pure number which is why
   230  	// xiToPSquaredMinus1Over3 is ∈ GF(p). With y we get a factor of -1. We
   231  	// ignore this to end up with -Q2.
   232  
   233  	minusQ2 := &twistPoint{}
   234  	minusQ2.x.MulScalar(&aAffine.x, xiToPSquaredMinus1Over3i)
   235  	minusQ2.y.Set(&aAffine.y)
   236  	minusQ2.z.SetOne()
   237  	minusQ2.t.SetOne()
   238  
   239  	r2.Square(&q1.y)
   240  	a, b, c, newR := lineFunctionAdd(r, q1, bAffine, r2) //step11
   241  	mulLine(ret, a, b, c)
   242  	r = newR
   243  
   244  	r2.Square(&minusQ2.y)
   245  	a, b, c, newR = lineFunctionAdd(r, minusQ2, bAffine, r2) //step12
   246  	mulLine(ret, a, b, c)
   247  
   248  	r = newR
   249  
   250  	return ret
   251  }
   252  
   253  // finalExponentiation computes the (p¹²-1)/Order-th power of an element of
   254  // GF(p¹²) to obtain an element of GT (steps 13-15 of algorithm 1 from
   255  // http://cryptojedi.org/papers/dclxvi-20100714.pdf)
   256  func finalExponentiation(in *gfP12) *gfP12 {
   257  	t1 := &gfP12{}
   258  
   259  	// This is the p^6-Frobenius
   260  	t1.x.Neg(&in.x)
   261  	t1.y.Set(&in.y)
   262  
   263  	inv := &gfP12{}
   264  	inv.Invert(in)
   265  	t1.Mul(t1, inv) //f^(p^6-1);//simon check ok;0930
   266  
   267  	t2 := (&gfP12{}).FrobeniusP2(t1)
   268  	t1.Mul(t1, t2) //simon:(f^(p^6-1))^(p^2+1) 0920
   269  
   270  	//simon1031version:
   271  	st0 := (&gfP12{}).Frobenius(t1)
   272  	sx0 := (&gfP12{}).Frobenius(st0)
   273  	sx1 := (&gfP12{}).Mul(t1, st0)
   274  
   275  	sx0 = sx0.Mul(sx0, sx1).Frobenius(sx0)
   276  	sx1 = sx1.Invert(t1) //simon check ok 1031;
   277  
   278  	sx4 := (&gfP12{})
   279  	//sx4.Exp(t1, tneg)
   280  	sx4 = sx4.Exp(t1, t)
   281  	sx4.Invert(sx4)
   282  
   283  	sx3 := (&gfP12{}).Frobenius(sx4)
   284  	//sx2 := (&gfP12{}).Exp(sx4, tneg)
   285  	sx2 := (&gfP12{}).Exp(sx4, t)
   286  	sx2.Invert(sx2)
   287  	sx5 := (&gfP12{}).Invert(sx2)
   288  
   289  	//st0.Exp(sx2, tneg)
   290  	st0.Exp(sx2, t)
   291  	st0.Invert(st0)
   292  
   293  	sx2.Frobenius(sx2)
   294  	sx2inv := (&gfP12{}).Invert(sx2)
   295  	sx4.Mul(sx4, sx2inv)
   296  	sx2.Frobenius(sx2)
   297  
   298  	t1.Frobenius(st0)
   299  	st0.Mul(st0, t1)
   300  
   301  	st0.Mul(st0, st0)
   302  	st0.Mul(st0, sx4)
   303  	st0.Mul(st0, sx5)
   304  	t1.Mul(sx3, sx5)
   305  	t1.Mul(t1, st0)
   306  
   307  	st0.Mul(st0, sx2)
   308  	t1.Mul(t1, t1)
   309  	t1.Mul(t1, st0)
   310  
   311  	t1.Mul(t1, t1)
   312  	st0.Mul(t1, sx1)
   313  	t1.Mul(t1, sx0)
   314  
   315  	st0.Mul(st0, st0)
   316  	st0.Mul(st0, t1)
   317  	//1031version end;
   318  
   319  	return st0
   320  }
   321  
   322  func optimalAte(a *twistPoint, b *curvePoint) *gfP12 {
   323  	e := miller(a, b)
   324  	//t1 := time.Now()
   325  	ret := finalExponentiation(e)
   326  	//elapsed := time.Since(t1)
   327  	//fmt.Println("final spend time: ", elapsed)
   328  
   329  	if a.IsInfinity() || b.IsInfinity() {
   330  		ret.SetOne()
   331  	}
   332  
   333  	return ret
   334  }