1// Copyright 2012 The Go Authors. All rights reserved.
2// Use of this source code is governed by a BSD-style
3// license that can be found in the LICENSE file.
4
5package bn256
6
7func lineFunctionAdd(r, p *twistPoint, q *curvePoint, r2 *gfP2, pool *bnPool) (a, b, c *gfP2, rOut *twistPoint) {
8	// See the mixed addition algorithm from "Faster Computation of the
9	// Tate Pairing", http://arxiv.org/pdf/0904.0854v3.pdf
10
11	B := newGFp2(pool).Mul(p.x, r.t, pool)
12
13	D := newGFp2(pool).Add(p.y, r.z)
14	D.Square(D, pool)
15	D.Sub(D, r2)
16	D.Sub(D, r.t)
17	D.Mul(D, r.t, pool)
18
19	H := newGFp2(pool).Sub(B, r.x)
20	I := newGFp2(pool).Square(H, pool)
21
22	E := newGFp2(pool).Add(I, I)
23	E.Add(E, E)
24
25	J := newGFp2(pool).Mul(H, E, pool)
26
27	L1 := newGFp2(pool).Sub(D, r.y)
28	L1.Sub(L1, r.y)
29
30	V := newGFp2(pool).Mul(r.x, E, pool)
31
32	rOut = newTwistPoint(pool)
33	rOut.x.Square(L1, pool)
34	rOut.x.Sub(rOut.x, J)
35	rOut.x.Sub(rOut.x, V)
36	rOut.x.Sub(rOut.x, V)
37
38	rOut.z.Add(r.z, H)
39	rOut.z.Square(rOut.z, pool)
40	rOut.z.Sub(rOut.z, r.t)
41	rOut.z.Sub(rOut.z, I)
42
43	t := newGFp2(pool).Sub(V, rOut.x)
44	t.Mul(t, L1, pool)
45	t2 := newGFp2(pool).Mul(r.y, J, pool)
46	t2.Add(t2, t2)
47	rOut.y.Sub(t, t2)
48
49	rOut.t.Square(rOut.z, pool)
50
51	t.Add(p.y, rOut.z)
52	t.Square(t, pool)
53	t.Sub(t, r2)
54	t.Sub(t, rOut.t)
55
56	t2.Mul(L1, p.x, pool)
57	t2.Add(t2, t2)
58	a = newGFp2(pool)
59	a.Sub(t2, t)
60
61	c = newGFp2(pool)
62	c.MulScalar(rOut.z, q.y)
63	c.Add(c, c)
64
65	b = newGFp2(pool)
66	b.SetZero()
67	b.Sub(b, L1)
68	b.MulScalar(b, q.x)
69	b.Add(b, b)
70
71	B.Put(pool)
72	D.Put(pool)
73	H.Put(pool)
74	I.Put(pool)
75	E.Put(pool)
76	J.Put(pool)
77	L1.Put(pool)
78	V.Put(pool)
79	t.Put(pool)
80	t2.Put(pool)
81
82	return
83}
84
85func lineFunctionDouble(r *twistPoint, q *curvePoint, pool *bnPool) (a, b, c *gfP2, rOut *twistPoint) {
86	// See the doubling algorithm for a=0 from "Faster Computation of the
87	// Tate Pairing", http://arxiv.org/pdf/0904.0854v3.pdf
88
89	A := newGFp2(pool).Square(r.x, pool)
90	B := newGFp2(pool).Square(r.y, pool)
91	C := newGFp2(pool).Square(B, pool)
92
93	D := newGFp2(pool).Add(r.x, B)
94	D.Square(D, pool)
95	D.Sub(D, A)
96	D.Sub(D, C)
97	D.Add(D, D)
98
99	E := newGFp2(pool).Add(A, A)
100	E.Add(E, A)
101
102	G := newGFp2(pool).Square(E, pool)
103
104	rOut = newTwistPoint(pool)
105	rOut.x.Sub(G, D)
106	rOut.x.Sub(rOut.x, D)
107
108	rOut.z.Add(r.y, r.z)
109	rOut.z.Square(rOut.z, pool)
110	rOut.z.Sub(rOut.z, B)
111	rOut.z.Sub(rOut.z, r.t)
112
113	rOut.y.Sub(D, rOut.x)
114	rOut.y.Mul(rOut.y, E, pool)
115	t := newGFp2(pool).Add(C, C)
116	t.Add(t, t)
117	t.Add(t, t)
118	rOut.y.Sub(rOut.y, t)
119
120	rOut.t.Square(rOut.z, pool)
121
122	t.Mul(E, r.t, pool)
123	t.Add(t, t)
124	b = newGFp2(pool)
125	b.SetZero()
126	b.Sub(b, t)
127	b.MulScalar(b, q.x)
128
129	a = newGFp2(pool)
130	a.Add(r.x, E)
131	a.Square(a, pool)
132	a.Sub(a, A)
133	a.Sub(a, G)
134	t.Add(B, B)
135	t.Add(t, t)
136	a.Sub(a, t)
137
138	c = newGFp2(pool)
139	c.Mul(rOut.z, r.t, pool)
140	c.Add(c, c)
141	c.MulScalar(c, q.y)
142
143	A.Put(pool)
144	B.Put(pool)
145	C.Put(pool)
146	D.Put(pool)
147	E.Put(pool)
148	G.Put(pool)
149	t.Put(pool)
150
151	return
152}
153
154func mulLine(ret *gfP12, a, b, c *gfP2, pool *bnPool) {
155	a2 := newGFp6(pool)
156	a2.x.SetZero()
157	a2.y.Set(a)
158	a2.z.Set(b)
159	a2.Mul(a2, ret.x, pool)
160	t3 := newGFp6(pool).MulScalar(ret.y, c, pool)
161
162	t := newGFp2(pool)
163	t.Add(b, c)
164	t2 := newGFp6(pool)
165	t2.x.SetZero()
166	t2.y.Set(a)
167	t2.z.Set(t)
168	ret.x.Add(ret.x, ret.y)
169
170	ret.y.Set(t3)
171
172	ret.x.Mul(ret.x, t2, pool)
173	ret.x.Sub(ret.x, a2)
174	ret.x.Sub(ret.x, ret.y)
175	a2.MulTau(a2, pool)
176	ret.y.Add(ret.y, a2)
177
178	a2.Put(pool)
179	t3.Put(pool)
180	t2.Put(pool)
181	t.Put(pool)
182}
183
184// sixuPlus2NAF is 6u+2 in non-adjacent form.
185var sixuPlus2NAF = []int8{0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, -1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, -1, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 1}
186
187// miller implements the Miller loop for calculating the Optimal Ate pairing.
188// See algorithm 1 from http://cryptojedi.org/papers/dclxvi-20100714.pdf
189func miller(q *twistPoint, p *curvePoint, pool *bnPool) *gfP12 {
190	ret := newGFp12(pool)
191	ret.SetOne()
192
193	aAffine := newTwistPoint(pool)
194	aAffine.Set(q)
195	aAffine.MakeAffine(pool)
196
197	bAffine := newCurvePoint(pool)
198	bAffine.Set(p)
199	bAffine.MakeAffine(pool)
200
201	minusA := newTwistPoint(pool)
202	minusA.Negative(aAffine, pool)
203
204	r := newTwistPoint(pool)
205	r.Set(aAffine)
206
207	r2 := newGFp2(pool)
208	r2.Square(aAffine.y, pool)
209
210	for i := len(sixuPlus2NAF) - 1; i > 0; i-- {
211		a, b, c, newR := lineFunctionDouble(r, bAffine, pool)
212		if i != len(sixuPlus2NAF)-1 {
213			ret.Square(ret, pool)
214		}
215
216		mulLine(ret, a, b, c, pool)
217		a.Put(pool)
218		b.Put(pool)
219		c.Put(pool)
220		r.Put(pool)
221		r = newR
222
223		switch sixuPlus2NAF[i-1] {
224		case 1:
225			a, b, c, newR = lineFunctionAdd(r, aAffine, bAffine, r2, pool)
226		case -1:
227			a, b, c, newR = lineFunctionAdd(r, minusA, bAffine, r2, pool)
228		default:
229			continue
230		}
231
232		mulLine(ret, a, b, c, pool)
233		a.Put(pool)
234		b.Put(pool)
235		c.Put(pool)
236		r.Put(pool)
237		r = newR
238	}
239
240	// In order to calculate Q1 we have to convert q from the sextic twist
241	// to the full GF(p^12) group, apply the Frobenius there, and convert
242	// back.
243	//
244	// The twist isomorphism is (x', y') -> (xω², yω³). If we consider just
245	// x for a moment, then after applying the Frobenius, we have x̄ω^(2p)
246	// where x̄ is the conjugate of x. If we are going to apply the inverse
247	// isomorphism we need a value with a single coefficient of ω² so we
248	// rewrite this as x̄ω^(2p-2)ω². ξ⁶ = ω and, due to the construction of
249	// p, 2p-2 is a multiple of six. Therefore we can rewrite as
250	// x̄ξ^((p-1)/3)ω² and applying the inverse isomorphism eliminates the
251	// ω².
252	//
253	// A similar argument can be made for the y value.
254
255	q1 := newTwistPoint(pool)
256	q1.x.Conjugate(aAffine.x)
257	q1.x.Mul(q1.x, xiToPMinus1Over3, pool)
258	q1.y.Conjugate(aAffine.y)
259	q1.y.Mul(q1.y, xiToPMinus1Over2, pool)
260	q1.z.SetOne()
261	q1.t.SetOne()
262
263	// For Q2 we are applying the p² Frobenius. The two conjugations cancel
264	// out and we are left only with the factors from the isomorphism. In
265	// the case of x, we end up with a pure number which is why
266	// xiToPSquaredMinus1Over3 is ∈ GF(p). With y we get a factor of -1. We
267	// ignore this to end up with -Q2.
268
269	minusQ2 := newTwistPoint(pool)
270	minusQ2.x.MulScalar(aAffine.x, xiToPSquaredMinus1Over3)
271	minusQ2.y.Set(aAffine.y)
272	minusQ2.z.SetOne()
273	minusQ2.t.SetOne()
274
275	r2.Square(q1.y, pool)
276	a, b, c, newR := lineFunctionAdd(r, q1, bAffine, r2, pool)
277	mulLine(ret, a, b, c, pool)
278	a.Put(pool)
279	b.Put(pool)
280	c.Put(pool)
281	r.Put(pool)
282	r = newR
283
284	r2.Square(minusQ2.y, pool)
285	a, b, c, newR = lineFunctionAdd(r, minusQ2, bAffine, r2, pool)
286	mulLine(ret, a, b, c, pool)
287	a.Put(pool)
288	b.Put(pool)
289	c.Put(pool)
290	r.Put(pool)
291	r = newR
292
293	aAffine.Put(pool)
294	bAffine.Put(pool)
295	minusA.Put(pool)
296	r.Put(pool)
297	r2.Put(pool)
298
299	return ret
300}
301
302// finalExponentiation computes the (p¹²-1)/Order-th power of an element of
303// GF(p¹²) to obtain an element of GT (steps 13-15 of algorithm 1 from
304// http://cryptojedi.org/papers/dclxvi-20100714.pdf)
305func finalExponentiation(in *gfP12, pool *bnPool) *gfP12 {
306	t1 := newGFp12(pool)
307
308	// This is the p^6-Frobenius
309	t1.x.Negative(in.x)
310	t1.y.Set(in.y)
311
312	inv := newGFp12(pool)
313	inv.Invert(in, pool)
314	t1.Mul(t1, inv, pool)
315
316	t2 := newGFp12(pool).FrobeniusP2(t1, pool)
317	t1.Mul(t1, t2, pool)
318
319	fp := newGFp12(pool).Frobenius(t1, pool)
320	fp2 := newGFp12(pool).FrobeniusP2(t1, pool)
321	fp3 := newGFp12(pool).Frobenius(fp2, pool)
322
323	fu, fu2, fu3 := newGFp12(pool), newGFp12(pool), newGFp12(pool)
324	fu.Exp(t1, u, pool)
325	fu2.Exp(fu, u, pool)
326	fu3.Exp(fu2, u, pool)
327
328	y3 := newGFp12(pool).Frobenius(fu, pool)
329	fu2p := newGFp12(pool).Frobenius(fu2, pool)
330	fu3p := newGFp12(pool).Frobenius(fu3, pool)
331	y2 := newGFp12(pool).FrobeniusP2(fu2, pool)
332
333	y0 := newGFp12(pool)
334	y0.Mul(fp, fp2, pool)
335	y0.Mul(y0, fp3, pool)
336
337	y1, y4, y5 := newGFp12(pool), newGFp12(pool), newGFp12(pool)
338	y1.Conjugate(t1)
339	y5.Conjugate(fu2)
340	y3.Conjugate(y3)
341	y4.Mul(fu, fu2p, pool)
342	y4.Conjugate(y4)
343
344	y6 := newGFp12(pool)
345	y6.Mul(fu3, fu3p, pool)
346	y6.Conjugate(y6)
347
348	t0 := newGFp12(pool)
349	t0.Square(y6, pool)
350	t0.Mul(t0, y4, pool)
351	t0.Mul(t0, y5, pool)
352	t1.Mul(y3, y5, pool)
353	t1.Mul(t1, t0, pool)
354	t0.Mul(t0, y2, pool)
355	t1.Square(t1, pool)
356	t1.Mul(t1, t0, pool)
357	t1.Square(t1, pool)
358	t0.Mul(t1, y1, pool)
359	t1.Mul(t1, y0, pool)
360	t0.Square(t0, pool)
361	t0.Mul(t0, t1, pool)
362
363	inv.Put(pool)
364	t1.Put(pool)
365	t2.Put(pool)
366	fp.Put(pool)
367	fp2.Put(pool)
368	fp3.Put(pool)
369	fu.Put(pool)
370	fu2.Put(pool)
371	fu3.Put(pool)
372	fu2p.Put(pool)
373	fu3p.Put(pool)
374	y0.Put(pool)
375	y1.Put(pool)
376	y2.Put(pool)
377	y3.Put(pool)
378	y4.Put(pool)
379	y5.Put(pool)
380	y6.Put(pool)
381
382	return t0
383}
384
385func optimalAte(a *twistPoint, b *curvePoint, pool *bnPool) *gfP12 {
386	e := miller(a, b, pool)
387	ret := finalExponentiation(e, pool)
388	e.Put(pool)
389
390	if a.IsInfinity() || b.IsInfinity() {
391		ret.SetOne()
392	}
393
394	return ret
395}
396