1 #include <stdarg.h>
2 #include <stdio.h>
3 #include <stdint.h>
4 #include <gmp.h>
5 #include "pbc_utils.h"
6 #include "pbc_field.h"
7 #include "pbc_fp.h"
8 #include "pbc_memory.h"
9 #include "pbc_param.h"
10 #include "pbc_pairing.h"
11 #include "pbc_ternary_extension_field.h"
12 #include "param.h"
13 
14 typedef struct { /* private data of $GF(3^m)$ */
15     unsigned int len; /* the number of native machine integers required to represent one GF(3^m) element */
16     int m; /* the irreducible polynomial is $x^m + x^t + 2$ */
17     int t; /* the irreducible polynomial is $x^m + x^t + 2$ */
18     element_ptr p; /* $p$ is the irreducible polynomial. */
19     mpz_t n; /* group order of $G_1$, $G_2$, $G_T$ */
20     mpz_t n2; /* order(elliptic curve points) / order(G_1) */
21 } params;
22 
23 struct pairing_data {
24     field_t gf3m, gf32m, gf36m;
25     mpz_t n2; // cofactor
26 };
27 typedef struct pairing_data *pairing_data_ptr;
28 
29 #define PARAM(e) ((params *)e->field->data)
30 #define ITEM(e,x,y) (element_item(element_item(e,x),y))
31 #define print(e) {printf(#e": "); element_out_str(stdout, 10, e); printf("\n");}
32 
33 struct point_s { // points on the elliptic curve $y^2=x^3-x+1$
34     int isinf;
35     element_t x, y;
36 };
37 typedef struct point_s *point_ptr;
38 typedef struct point_s point_t[1];
39 
40 #define FIELD(e) ((field_ptr) e->field)
41 #define BASE(e) ((field_ptr) FIELD(e)->data)
42 #define DATA(e) ((point_ptr) e->data)
43 
point_set(element_t e,element_t a)44 static void point_set(element_t e, element_t a) {
45     point_ptr r = DATA(e), p = DATA(a);
46     r->isinf = p->isinf;
47     if (!p->isinf) {
48         element_set(r->x, p->x);
49         element_set(r->y, p->y);
50     }
51 }
52 
point_init(element_t e)53 static void point_init(element_t e) {
54     field_ptr f = BASE(e);
55     e->data = pbc_malloc(sizeof(struct point_s));
56     point_ptr p = DATA(e);
57     element_init(p->x, f);
58     element_init(p->y, f);
59     p->isinf = 1;
60 }
61 
point_clear(element_t e)62 static void point_clear(element_t e) {
63     point_ptr p = DATA(e);
64     element_clear(p->x);
65     element_clear(p->y);
66     pbc_free(p);
67 }
68 
69 /* return 1 if $a!=b$, 0 otherwise. */
point_cmp(element_t a,element_t b)70 static int point_cmp(element_t a, element_t b) {
71     point_ptr pa = DATA(a), pb = DATA(b);
72     if (pa->isinf == pb->isinf) {
73         if (pa->isinf)
74             return 0;
75         else
76             return element_cmp(pa->x, pb->x) || element_cmp(pa->y, pb->y);
77     } else
78         return 1;
79 }
80 
point_set0(element_ptr e)81 static void point_set0(element_ptr e) {
82     DATA(e)->isinf = 1;
83 }
84 
point_is0(element_ptr e)85 static int point_is0(element_ptr e) {
86     return DATA(e)->isinf;
87 }
88 
point_random(element_t a)89 static void point_random(element_t a) {
90     point_ptr p = DATA(a);
91     element_ptr x = p->x, y = p->y;
92     field_ptr f = x->field;
93     p->isinf = 0;
94     element_t t, t2, e1;
95     element_init(t, f);
96     element_init(e1, f);
97     element_set1(e1);
98     element_init(t2, f);
99     do {
100         element_random(x);
101         if (element_is0(x))
102             continue;
103         element_cubic(t, x); // t == x^3
104         element_sub(t, t, x); // t == x^3 - x
105         element_add(t, t, e1); // t == x^3 - x + 1
106         element_sqrt(y, t);  // y == sqrt(x^3 - x + 1)
107         element_mul(t2, y, y); // t2 == x^3 - x + 1
108     } while (element_cmp(t2, t)); // t2 != t
109 
110     // make sure order of $a$ is order of $G_1$
111     pairing_ptr pairing = FIELD(a)->pairing;
112     pairing_data_ptr dp = pairing->data;
113     element_pow_mpz(a, a, dp->n2);
114 
115     element_clear(t);
116     element_clear(t2);
117     element_clear(e1);
118 }
119 
point_add(element_t c,element_t a,element_t b)120 static void point_add(element_t c, element_t a, element_t b) {
121     point_ptr p1 = DATA(a), p2 = DATA(b), p3 = DATA(c);
122     int inf1 = p1->isinf, inf2 = p2->isinf;
123     element_ptr x1 = p1->x, y1 = p1->y, x2 = p2->x, y2 = p2->y;
124     field_ptr f = FIELD(x1);
125     if (inf1) {
126         point_set(c, b);
127         return;
128     }
129     if (inf2) {
130         point_set(c, a);
131         return;
132     }
133     element_t v0, v1, v2, v3, v4, ny2;
134     element_init(v0, f);
135     element_init(v1, f);
136     element_init(v2, f);
137     element_init(v3, f);
138     element_init(v4, f);
139     element_init(ny2, f);
140     if (!element_cmp(x1, x2)) { // x1 == x2
141         element_neg(ny2, y2); // ny2 == -y2
142         if (!element_cmp(y1, ny2)) {
143             p3->isinf = 1;
144             goto end;
145         }
146         if (!element_cmp(y1, y2)) { // y1 == y2
147             element_invert(v0, y1); // v0 == y1^{-1}
148             element_mul(v1, v0, v0); // v1 == [y1^{-1}]^2
149             element_add(p3->x, v1, x1); // v1 == [y1^{-1}]^2 + x1
150             element_cubic(v2, v0); // v2 == [y1^{-1}]^3
151             element_add(v2, v2, y1); // v2 == [y1^{-1}]^3 + y1
152             element_neg(p3->y, v2); // p3 == -([y1^{-1}]^3 + y1)
153             p3->isinf = 0;
154             goto end;
155         }
156     }
157     // $P1 \ne \pm P2$
158     element_sub(v0, x2, x1); // v0 == x2-x1
159     element_invert(v1, v0); // v1 == (x2-x1)^{-1}
160     element_sub(v0, y2, y1); // v0 == y2-y1
161     element_mul(v2, v0, v1); // v2 == (y2-y1)/(x2-x1)
162     element_mul(v3, v2, v2); // v3 == [(y2-y1)/(x2-x1)]^2
163     element_cubic(v4, v2); // v4 == [(y2-y1)/(x2-x1)]^3
164     element_add(v0, x1, x2); // v0 == x1+x2
165     element_sub(v3, v3, v0); // v3 == [(y2-y1)/(x2-x1)]^2 - (x1+x2)
166     element_add(v0, y1, y2); // v0 == y1+y2
167     element_sub(v4, v0, v4); // v4 == (y1+y2) - [(y2-y1)/(x2-x1)]^3
168     p3->isinf = 0;
169     element_set(p3->x, v3);
170     element_set(p3->y, v4);
171     end: element_clear(v0);
172     element_clear(v1);
173     element_clear(v2);
174     element_clear(v3);
175     element_clear(v4);
176     element_clear(ny2);
177 }
178 
point_invert(element_ptr e,element_ptr a)179 static void point_invert(element_ptr e, element_ptr a) {
180     point_ptr r = DATA(e), p = DATA(a);
181     r->isinf = p->isinf;
182     if (!p->isinf) {
183         element_set(r->x, p->x);
184         element_neg(r->y, p->y);
185     }
186 }
187 
point_out_str(FILE * stream,int base,element_ptr a)188 static size_t point_out_str(FILE *stream, int base, element_ptr a) {
189     point_ptr p = DATA(a);
190     size_t size = 0;
191     if (p->isinf)
192         return fprintf(stream, "O");
193     else {
194         size += element_out_str(stream, base, p->x);
195         size += element_out_str(stream, base, p->y);
196         return size;
197     }
198 }
199 
point_field_clear(field_ptr f)200 static void point_field_clear(field_ptr f) {
201     UNUSED_VAR(f);
202 }
203 
field_init_eta_T_3(field_t f,field_t base)204 void field_init_eta_T_3(field_t f, field_t base) {
205     field_init(f);
206     f->data = (void *) base;
207     f->init = point_init;
208     f->clear = point_clear;
209     f->random = point_random;
210     f->set = point_set;
211     f->cmp = point_cmp;
212     f->invert = f->neg = point_invert;
213     f->mul = f->add = point_add;
214     f->set1 = f->set0 = point_set0;
215     f->is1 = f->is0 = point_is0;
216     f->mul_mpz = f->pow_mpz;
217     f->out_str = point_out_str;
218     f->field_clear = point_field_clear;
219     f->name = "eta_T_3 point group";
220 }
221 
222 /* computing of $(-t^2 +u*s -t*p -p^2)^3$
223  * The algorithm is by J.Beuchat et.al, in the paper of "Algorithms and Arithmetic Operators for Computing
224  * the $eta_T$ Pairing in Characteristic Three", algorithm 4 in the appendix */
algorithm4a(element_t S,element_t t,element_t u)225 static void algorithm4a(element_t S, element_t t, element_t u) {
226     field_ptr f = FIELD(t);
227     element_t e1, c0, c1, m0, v0, v2;
228     element_init(e1, f);
229     element_init(c0, f);
230     element_init(c1, f);
231     element_init(m0, f);
232     element_init(v0, f);
233     element_init(v2, f);
234     element_set1(e1);
235     element_cubic(c0, t); // c0 == t^3
236     element_cubic(c1, u);
237     element_neg(c1, c1); // c1 == -u^3
238     element_mul(m0, c0, c0); // m0 == c0^2
239     element_neg(v0, m0); // v0 == -c0^2
240     element_sub(v0, v0, c0); // v0 == -c0^2 -c0
241     element_sub(v0, v0, e1); // v0 == -c0^2 -c0 -1
242     element_set1(v2);
243     element_sub(v2, v2, c0); // v2 == 1 -c0
244     // v1 == c1
245     // S == [[v0, v1], [v2, f3m.zero()], [f3m.two(), f3m.zero()]]
246     element_set(ITEM(S,0,0), v0);
247     element_set(ITEM(S,0,1), c1);
248     element_set(ITEM(S,1,0), v2);
249     element_set0(ITEM(S,1,1));
250     element_neg(ITEM(S,2,0), e1);
251     element_set0(ITEM(S,2,1));
252     element_clear(e1);
253     element_clear(c0);
254     element_clear(c1);
255     element_clear(m0);
256     element_clear(v0);
257     element_clear(v2);
258 }
259 
algorithm5(element_t c,element_ptr xp,element_ptr yp,element_ptr xq,element_ptr yq)260 static void algorithm5(element_t c, element_ptr xp, element_ptr yp,
261         element_ptr xq, element_ptr yq) {
262     params *p = PARAM(xp);
263     unsigned int re = p->m % 12;
264     field_ptr f = FIELD(xp) /*GF(3^m)*/, f6 = FIELD(c) /*GF(3^{6*m})*/;
265     element_t e1, xpp, ypp, xqq, yqq, t, nt, nt2, v1, v2, a1, a2, R, u, nu, S, S2;
266     element_init(e1, f);
267     element_init(xpp, f);
268     element_init(ypp, f);
269     element_init(xqq, f);
270     element_init(yqq, f);
271     element_init(t, f);
272     element_init(nt, f);
273     element_init(nt2, f);
274     element_init(v1, f);
275     element_init(v2, f);
276     element_init(a1, f6);
277     element_init(a2, f6);
278     element_init(R, f6);
279     element_init(u, f);
280     element_init(nu, f);
281     element_init(S, f6);
282     element_init(S2, f6);
283     element_set1(e1);
284     element_set(xpp, xp);
285     xp = xpp; // clone
286     element_add(xp, xp, e1); // xp == xp + b
287     element_set(ypp, yp);
288     yp = ypp; // clone
289     if (re == 1 || re == 11)
290         element_neg(yp, yp); // yp == -\mu*b*yp, \mu == 1 when re==1, or 11
291     element_set(xqq, xq);
292     xq = xqq; // clone
293     element_cubic(xq, xq); // xq == xq^3
294     element_set(yqq, yq);
295     yq = yqq; // clone
296     element_cubic(yq, yq); // yq == yq^3
297     element_add(t, xp, xq); // t == xp+xq
298     element_neg(nt, t); // nt == -t
299     element_mul(nt2, t, nt); // nt2 == -t^2
300     element_mul(v2, yp, yq); // v2 == yp*yq
301     element_mul(v1, yp, t); // v1 == yp*t
302     if (re == 7 || re == 11) { // \lambda == 1
303         element_t nyp, nyq;
304         element_init(nyp, f);
305         element_init(nyq, f);
306         element_neg(nyp, yp); // nyp == -yp
307         element_neg(nyq, yq); // nyq == -yq
308         element_set(ITEM(a1,0,0), v1);
309         element_set(ITEM(a1,0,1), nyq);
310         element_set(ITEM(a1,1,0), nyp);
311         element_clear(nyp);
312         element_clear(nyq);
313     } else { // \lambda == -1
314         element_neg(v1, v1); // v1 == -yp*t
315         element_set(ITEM(a1,0,0), v1);
316         element_set(ITEM(a1,0,1), yq);
317         element_set(ITEM(a1,1,0), yp);
318     }
319     // a2 == -t^2 +yp*yq*s -t*p -p^2
320     element_set(ITEM(a2,0,0), nt2);
321     element_set(ITEM(a2,0,1), v2);
322     element_set(ITEM(a2,1,0), nt);
323     element_neg(ITEM(a2,2,0), e1);
324     element_mul(R, a1, a2);
325     int i;
326     for (i = 0; i < (p->m - 1) / 4; i++) {
327         element_cubic(R, R);
328         element_cubic(R, R); // R <= R^9
329         element_cubic(xq, xq);
330         element_cubic(xq, xq);
331         element_sub(xq, xq, e1); // xq <= xq^9-b
332         element_cubic(yq, yq);
333         element_cubic(yq, yq); // yq <= yq^9
334         element_add(t, xp, xq); // t == xp+xq
335         element_mul(u, yp, yq); // u == yp*yq
336         element_neg(nu, u); // nu == -yp*yq
337         algorithm4a(S, t, nu); // S == (-t^2 -u*s -t*p -p^2)^3
338         element_cubic(xq, xq);
339         element_cubic(xq, xq);
340         element_sub(xq, xq, e1); // xq <= xq^9-b
341         element_cubic(yq, yq);
342         element_cubic(yq, yq); // yq <= yq^9
343         element_add(t, xp, xq); // t == xp+xq
344         element_mul(u, yp, yq); // u == yp*yq
345         element_neg(nt, t); // nt == -t
346         element_mul(nt2, t, nt); // nt2 == -t^2
347         // S2 = [[nt2, u], [nt, f3m.zero()], [f3m.two(), f3m.zero()]]
348         // S2 == -t^2 +u*s -t*p -p^2
349         element_set(ITEM(S2,0,0), nt2);
350         element_set(ITEM(S2,0,1), u);
351         element_set(ITEM(S2,1,0), nt);
352         element_set0(ITEM(S2,1,1));
353         element_neg(ITEM(S2,2,0), e1);
354         element_set0(ITEM(S2,2,1));
355         element_mul(S, S, S2);
356         element_mul(R, R, S);
357     }
358     element_set(c, R);
359     element_clear(e1);
360     element_clear(xpp);
361     element_clear(ypp);
362     element_clear(xqq);
363     element_clear(yqq);
364     element_clear(t);
365     element_clear(nt);
366     element_clear(nt2);
367     element_clear(v1);
368     element_clear(v2);
369     element_clear(a1);
370     element_clear(a2);
371     element_clear(R);
372     element_clear(u);
373     element_clear(nu);
374     element_clear(S);
375     element_clear(S2);
376 }
377 
378 /* this is the algorithm 4 in the paper of J.Beuchat et.al, "Algorithms and Arithmetic Operators for Computing
379  * the $eta_T$ Pairing in Characteristic Three" */
algorithm4(element_t c,element_ptr xp,element_ptr yp,element_ptr xq,element_ptr yq)380 static void algorithm4(element_t c, element_ptr xp, element_ptr yp,
381         element_ptr xq, element_ptr yq) {
382     params *p = PARAM(xp);
383     unsigned int re = p->m % 12;
384     field_ptr f = FIELD(xp) /*GF(3^m)*/, f6 = FIELD(c) /*GF(3^{6*m})*/;
385     element_t e1, xpp, ypp, xqq, yqq, t, nt, nt2, v1, v2, a1, a2, R, u, S;
386     element_init(e1, f);
387     element_init(xpp, f);
388     element_init(ypp, f);
389     element_init(xqq, f);
390     element_init(yqq, f);
391     element_init(t, f);
392     element_init(nt, f);
393     element_init(nt2, f);
394     element_init(v1, f);
395     element_init(v2, f);
396     element_init(a1, f6);
397     element_init(a2, f6);
398     element_init(R, f6);
399     element_init(u, f);
400     element_init(S, f6);
401     element_set1(e1);
402     element_set(xpp, xp);
403     xp = xpp; // clone
404     element_add(xp, xp, e1); // xp == xp + b
405     element_set(ypp, yp);
406     yp = ypp; // clone
407     if (re == 1 || re == 11)
408         element_neg(yp, yp); // yp == -\mu*b*yp, \mu == 1 when re==1, or 11
409     element_set(xqq, xq);
410     xq = xqq; // clone
411     element_cubic(xq, xq); // xq == xq^3
412     element_set(yqq, yq);
413     yq = yqq; // clone
414     element_cubic(yq, yq); // yq == yq^3
415     element_add(t, xp, xq); // t == xp+xq
416     element_neg(nt, t); // nt == -t
417     element_mul(nt2, t, nt); // nt2 == -t^2
418     element_mul(v2, yp, yq); // v2 == yp*yq
419     element_mul(v1, yp, t); // v1 == yp*t
420     if (re == 7 || re == 11) { // \lambda == 1
421         element_t nyp, nyq;
422         element_init(nyp, f);
423         element_init(nyq, f);
424         element_neg(nyp, yp); // nyp == -yp
425         element_neg(nyq, yq); // nyq == -yq
426         element_set(ITEM(a1,0,0), v1);
427         element_set(ITEM(a1,0,1), nyq);
428         element_set(ITEM(a1,1,0), nyp);
429         element_clear(nyp);
430         element_clear(nyq);
431     } else { // \lambda == -1
432         element_neg(v1, v1); // v1 == -yp*t
433         element_set(ITEM(a1,0,0), v1);
434         element_set(ITEM(a1,0,1), yq);
435         element_set(ITEM(a1,1,0), yp);
436     }
437     // a2 == -t^2 +yp*yq*s -t*p -p^2
438     element_set(ITEM(a2,0,0), nt2);
439     element_set(ITEM(a2,0,1), v2);
440     element_set(ITEM(a2,1,0), nt);
441     element_neg(ITEM(a2,2,0), e1);
442     element_mul(R, a1, a2);
443     int i;
444     for (i = 0; i < (p->m - 1) / 2; i++) {
445         element_cubic(R, R);
446         element_cubic(xq, xq);
447         element_cubic(xq, xq);
448         element_sub(xq, xq, e1); // xq <= xq^9-b
449         element_cubic(yq, yq);
450         element_cubic(yq, yq);
451         element_neg(yq, yq); // yq <= -yq^9
452         element_add(t, xp, xq); // t == xp+xq
453         element_neg(nt, t); // nt == -t
454         element_mul(nt2, t, nt); // nt2 == -t^2
455         element_mul(u, yp, yq); // u == yp*yq
456         element_set0(S);
457         element_set(ITEM(S,0,0), nt2);
458         element_set(ITEM(S,0,1), u);
459         element_set(ITEM(S,1,0), nt);
460         element_neg(ITEM(S,2,0), e1);
461         element_mul(R, R, S);
462     }
463     element_set(c, R);
464     element_clear(e1);
465     element_clear(xpp);
466     element_clear(ypp);
467     element_clear(xqq);
468     element_clear(yqq);
469     element_clear(t);
470     element_clear(nt);
471     element_clear(nt2);
472     element_clear(v1);
473     element_clear(v2);
474     element_clear(a1);
475     element_clear(a2);
476     element_clear(R);
477     element_clear(u);
478     element_clear(S);
479 }
480 
481 /* computation of $c <- U ^ {3^{3m} - 1}$
482  * This is the algorithm 6 in the paper above. */
algorithm6(element_t c,element_t u)483 static void algorithm6(element_t c, element_t u) {
484     element_ptr u0 = ITEM(u,0,0), u1 = ITEM(u,0,1), u2 = ITEM(u,1,0), u3 =
485             ITEM(u,1,1), u4 = ITEM(u,2,0), u5 = ITEM(u,2,1);
486     field_ptr f = FIELD(u0); /*GF(3^m)*/
487     field_t f3; /*GF(3^{3*m})*/
488     field_init_gf33m(f3, f);
489     element_t v0, v1, m0, m1, m2, a0, a1, i;
490     element_init(v0, f3);
491     element_init(v1, f3);
492     element_init(m0, f3);
493     element_init(m1, f3);
494     element_init(m2, f3);
495     element_init(a0, f3);
496     element_init(a1, f3);
497     element_init(i, f3);
498     element_set(element_item(v0, 0), u0);
499     element_set(element_item(v0, 1), u2);
500     element_set(element_item(v0, 2), u4);
501     element_set(element_item(v1, 0), u1);
502     element_set(element_item(v1, 1), u3);
503     element_set(element_item(v1, 2), u5);
504     element_mul(m0, v0, v0);
505     element_mul(m1, v1, v1);
506     element_mul(m2, v0, v1);
507     element_sub(a0, m0, m1);
508     element_add(a1, m0, m1);
509     element_invert(i, a1);
510     element_mul(v0, a0, i);
511     element_mul(v1, m2, i);
512     element_set(ITEM(c,0,0), element_item(v0, 0));
513     element_set(ITEM(c,1,0), element_item(v0, 1));
514     element_set(ITEM(c,2,0), element_item(v0, 2));
515     element_set(ITEM(c,0,1), element_item(v1, 0));
516     element_set(ITEM(c,1,1), element_item(v1, 1));
517     element_set(ITEM(c,2,1), element_item(v1, 2));
518     element_clear(v0);
519     element_clear(v1);
520     element_clear(m0);
521     element_clear(m1);
522     element_clear(m2);
523     element_clear(a0);
524     element_clear(a1);
525     element_clear(i);
526     field_clear(f3);
527 }
528 
529 /* computation of $c <- U ^ {3^m+1}$, $U \in T_2(F_{3^3M})$
530  * This is the algorithm 7 in the paper above. */
algorithm7(element_t c,element_t u)531 static void algorithm7(element_t c, element_t u) {
532     element_ptr u0 = ITEM(u,0,0), u1 = ITEM(u,0,1), u2 = ITEM(u,1,0), u3 =
533             ITEM(u,1,1), u4 = ITEM(u,2,0), u5 = ITEM(u,2,1);
534     field_ptr f = FIELD(u0); /*GF(3^m)*/
535     params *p = PARAM(u0);
536     element_t a0, a1, a2, a3, a4, a5, a6, m0, m1, m2, m3, m4, m5, m6, m7, m8,
537             v0, v1, v2, v3, v4, v5, e1;
538     element_init(a0, f);
539     element_init(a1, f);
540     element_init(a2, f);
541     element_init(a3, f);
542     element_init(a4, f);
543     element_init(a5, f);
544     element_init(a6, f);
545     element_init(m0, f);
546     element_init(m1, f);
547     element_init(m2, f);
548     element_init(m3, f);
549     element_init(m4, f);
550     element_init(m5, f);
551     element_init(m6, f);
552     element_init(m7, f);
553     element_init(m8, f);
554     element_init(v0, f);
555     element_init(v1, f);
556     element_init(v2, f);
557     element_init(v3, f);
558     element_init(v4, f);
559     element_init(v5, f);
560     element_init(e1, f);
561     element_set1(e1);
562     element_add(a0, u0, u1);
563     element_add(a1, u2, u3);
564     element_sub(a2, u4, u5);
565     element_mul(m0, u0, u4);
566     element_mul(m1, u1, u5);
567     element_mul(m2, u2, u4);
568     element_mul(m3, u3, u5);
569     element_mul(m4, a0, a2);
570     element_mul(m5, u1, u2);
571     element_mul(m6, u0, u3);
572     element_mul(m7, a0, a1);
573     element_mul(m8, a1, a2);
574     element_add(a3, m5, m6);
575     element_sub(a3, a3, m7);
576     element_neg(a4, m2);
577     element_sub(a4, a4, m3);
578     element_sub(a5, m3, m2);
579     element_sub(a6, m1, m0);
580     element_add(a6, a6, m4);
581     if (p->m % 6 == 1) {
582         element_add(v0, m0, m1);
583         element_add(v0, v0, a4);
584         element_add(v0, e1, v0);
585         element_sub(v1, m5, m6);
586         element_add(v1, v1, a6);
587         element_sub(v2, a4, a3);
588         element_add(v3, m8, a5);
589         element_sub(v3, v3, a6);
590         element_add(v4, a3, a4);
591         element_neg(v4, v4);
592         element_add(v5, m8, a5);
593     } else { // p->m % 6 == 5
594         element_add(v0, m0, m1);
595         element_sub(v0, v0, a4);
596         element_add(v0, e1, v0);
597         element_sub(v1, m6, m5);
598         element_add(v1, v1, a6);
599         element_set(v2, a3);
600         element_add(v3, m8, a5);
601         element_add(v3, v3, a6);
602         element_add(v4, a3, a4);
603         element_neg(v4, v4);
604         element_add(v5, m8, a5);
605         element_neg(v5, v5);
606     }
607     element_set(ITEM(c,0,0), v0);
608     element_set(ITEM(c,0,1), v1);
609     element_set(ITEM(c,1,0), v2);
610     element_set(ITEM(c,1,1), v3);
611     element_set(ITEM(c,2,0), v4);
612     element_set(ITEM(c,2,1), v5);
613     element_clear(a0);
614     element_clear(a1);
615     element_clear(a2);
616     element_clear(a3);
617     element_clear(a4);
618     element_clear(a5);
619     element_clear(a6);
620     element_clear(m0);
621     element_clear(m1);
622     element_clear(m2);
623     element_clear(m3);
624     element_clear(m4);
625     element_clear(m5);
626     element_clear(m6);
627     element_clear(m7);
628     element_clear(m8);
629     element_clear(v0);
630     element_clear(v1);
631     element_clear(v2);
632     element_clear(v3);
633     element_clear(v4);
634     element_clear(v5);
635     element_clear(e1);
636 }
637 
638 /* computing $c <- U^M, M=(3^{3m}-1)*(3^m+1)*(3^m+1-\mu*b*3^{(m+1)//2})$
639  * This is the algorithm 8 in the paper above. */
algorithm8(element_t c,element_t u)640 static void algorithm8(element_t c, element_t u) {
641     field_ptr f6 = FIELD(u), f = FIELD(ITEM(u,0,0));
642     params *p = (params *) f->data;
643     element_t v, w;
644     element_init(v, f6);
645     element_init(w, f6);
646     algorithm6(v, u);
647     algorithm7(v, v);
648     element_set(w, v);
649     int i;
650     for (i = 0; i < (p->m + 1) / 2; i++)
651         element_cubic(w, w);
652     algorithm7(v, v);
653     if (p->m % 12 == 1 || p->m % 12 == 11) { // w <= w^{-\mu*b}
654         element_ptr e;
655         e = ITEM(w,0,1);
656         element_neg(e, e);
657         e = ITEM(w,1,1);
658         element_neg(e, e);
659         e = ITEM(w,2,1);
660         element_neg(e, e);
661     }
662     element_mul(c, v, w);
663     element_clear(v);
664     element_clear(w);
665 }
666 
667 /* computing the Eta_T bilinear pairing $c <- Eta_T pairing(P,R)$ */
eta_T_pairing(element_ptr c,element_ptr P,element_ptr R,struct pairing_s * p)668 static void eta_T_pairing(element_ptr c, element_ptr P, element_ptr R, struct pairing_s *p) {
669     UNUSED_VAR(p);
670     if (DATA(P)->isinf || DATA(R)->isinf)
671         element_set1(c);
672     else {
673         element_ptr x1 = DATA(P)->x, y1 = DATA(P)->y, x2 = DATA(R)->x, y2 =
674                 DATA(R)->y;
675         if((PARAM(x1)->m - 1) / 2 % 2 == 0)
676             algorithm5(c, x1, y1, x2, y2);
677         else
678             algorithm4(c, x1, y1, x2, y2);
679         algorithm8(c, c);
680     }
681 }
682 
eta_T_3_clear(params * p)683 static void eta_T_3_clear(params *p) {
684     mpz_clear(p->n);
685     mpz_clear(p->n2);
686     pbc_free(p);
687 }
688 
GT_random(element_ptr e)689 static void GT_random(element_ptr e) {
690     element_t a, b;
691     element_init(a, e->field->pairing->G1);
692     element_init(b, e->field->pairing->G1);
693     element_random(a);
694     element_random(b);
695     element_pairing(e, a, b);
696     element_clear(a);
697     element_clear(b);
698 }
699 
eta_T_3_pairing_clear(pairing_t pairing)700 static void eta_T_3_pairing_clear(pairing_t pairing) {
701     mpz_clear(pairing->r);
702     field_clear(pairing->Zr);
703     field_clear(pairing->GT);
704     field_clear(pairing->G1);
705     pbc_free(pairing->G1);
706     pairing_data_ptr dp = pairing->data;
707     field_clear(dp->gf3m);
708     field_clear(dp->gf32m);
709     field_clear(dp->gf36m);
710     mpz_clear(dp->n2);
711     pbc_free(dp);
712 }
713 
eta_T_3_init_pairing(pairing_t pairing,params * p)714 static void eta_T_3_init_pairing(pairing_t pairing, params *p) {
715     mpz_init(pairing->r);
716     mpz_set(pairing->r, p->n);
717     field_init_fp(pairing->Zr, pairing->r);
718 
719     pairing_data_ptr dp = pbc_malloc(sizeof(*dp));
720     mpz_init(dp->n2);
721     mpz_set(dp->n2, p->n2);
722     field_init_gf3m(dp->gf3m, p->m, p->t);
723     field_init_gf32m(dp->gf32m, dp->gf3m);
724     field_init_gf33m(dp->gf36m, dp->gf32m);
725     pairing_GT_init(pairing, dp->gf36m);
726     pairing->GT->name = "eta_T_3 group of roots of 1";
727     pairing->GT->random = GT_random;
728     pairing->G2 = pairing->G1 = pbc_malloc(sizeof(field_t));
729     field_init_eta_T_3(pairing->G1, dp->gf3m);
730     pairing->G1->pairing = pairing;
731     mpz_set(pairing->G1->order, p->n);
732     mpz_set(pairing->GT->order, p->n);
733     pairing->map = eta_T_pairing;
734     pairing->data = dp;
735     pairing->clear_func = eta_T_3_pairing_clear;
736 }
737 
eta_T_3_out_str(FILE * stream,params * p)738 static void eta_T_3_out_str(FILE *stream, params *p) {
739     param_out_type(stream, "i");
740     param_out_int(stream, "m", p->m);
741     param_out_int(stream, "t", p->t);
742     param_out_mpz(stream, "n", p->n);
743     param_out_mpz(stream, "n2", p->n2);
744 }
745 
param_init(pbc_param_ptr p)746 static void param_init(pbc_param_ptr p) {
747     static pbc_param_interface_t interface = {{
748       (void (*)(void *))eta_T_3_clear,
749       (void (*)(pairing_t, void *))eta_T_3_init_pairing,
750       (void (*)(FILE *, void *))eta_T_3_out_str,
751     }};
752     p->api = interface;
753     params *param = p->data = pbc_malloc(sizeof(*param));
754     mpz_init(param->n);
755     mpz_init(param->n2);
756 }
757 
pbc_param_init_i(pbc_param_ptr p,struct symtab_s * tab)758 int pbc_param_init_i(pbc_param_ptr p, struct symtab_s *tab) {
759     param_init(p);
760     params *param = p->data;
761     int err = 0;
762     err += lookup_int(&param->m, tab, "m");
763     err += lookup_int(&param->t, tab, "t");
764     err += lookup_mpz(param->n, tab, "n");
765     err += lookup_mpz(param->n2, tab, "n2");
766     return err;
767 }
768 
pbc_param_init_i_gen(pbc_param_ptr par,int group_size)769 void pbc_param_init_i_gen(pbc_param_ptr par, int group_size) {
770     param_init(par);
771     params *p = par->data;
772     if (group_size <= 150) {
773         p->m = 97;
774         p->t = 12;
775         mpz_set_str(p->n, "2726865189058261010774960798134976187171462721", 10);
776         mpz_set_str(p->n2, "7", 10);
777     } else if (group_size <= 206) {
778         p->m = 199;
779         p->t = 164;
780         mpz_set_str(p->n, "167725321489096000055336949742738378351010268990525380470313869", 10);
781         mpz_set_str(p->n2, "527874953560391326545598291952743", 10);
782     } else if (group_size <= 259) {
783         p->m = 235;
784         p->t = 26;
785         mpz_set_str(p->n, "1124316700897695330265827797088699345032488681307846555184025129863722718180241", 10);
786         mpz_set_str(p->n2, "11819693021332914275777073321995059", 10);
787     } else if (group_size <= 316) {
788         p->m = 385;
789         p->t = 22;
790         mpz_set_str(p->n, "140884762419712839999909157778648717913595360839856026704744558309545986970238264714753014287541", 10);
791         mpz_set_str(p->n2, "34899486997246711147841377458771182755186809219564106252058066150110543296498189654810187", 10);
792     } else if (group_size <= 376) {
793         p->m = 337;
794         p->t = 30;
795         mpz_set_str(p->n, "250796519030408069744426774377542635685621984993105288007781750196791322190409525696108840742205849171229571431053", 10);
796         mpz_set_str(p->n2, "245777055088325363697128811262733732423405120899", 10);
797     } else if (group_size <= 430) {
798         p->m = 373;
799         p->t = 198;
800         mpz_set_str(p->n, "2840685307599487500956683789051368080919805957805957356540760731597378326586402072132959867084691357708217739285576524329854284197", 10);
801         mpz_set_str(p->n2, "3256903458766749542151641063558247849550904613763", 10);
802     } else if (group_size <= 484) {
803         p->m = 395;
804         p->t = 338;
805         mpz_set_str(p->n, "80172097064154181257340545445945701478615643539554910656655431171167598268341527430200810544156625333601812351266052856520678455274751591367269291", 10);
806         mpz_set_str(p->n2, "3621365590261279902324876775553649595261567", 10);
807     } else if (group_size <= 552) {
808         p->m = 433;
809         p->t = 120;
810         mpz_set_str(p->n, "15699907553631673835088720676147779193076555382157913339177784853763686462870506492752576492212322736133645158157557950634628006965882177348385366381692092784577773463", 10);
811         mpz_set_str(p->n2, "24980791723059119877470531054938874784049", 10);
812     } else if (group_size <= 644) {
813         p->m = 467;
814         p->t = 48;
815         mpz_set_str(p->n, "108220469499363631995525712756135494735252733492048868417164002000654321383482753640072319529019505742300964525569770933946381504691909098938045089999753901375631613294579329433690943459352138231", 10);
816         mpz_set_str(p->n2, "60438898450096967424971813347", 10);
817     } else if (group_size <= 696) {
818         p->m = 503;
819         p->t = 104;
820         mpz_set_str(p->n, "545523657676112447260904563578912738373307867219686215849632469801471112426878939776725222290437653718473962733760874627315930933126581248465899651120481066111839081575164964589811985885719017214938514563804313", 10);
821         mpz_set_str(p->n2, "1799606423432800810122901025413", 10);
822     } else if (group_size <= 803) {
823         p->m = 509;
824         p->t = 358;
825         mpz_set_str(p->n, "102239946202586852409809887418093021457150612495255706614733003327526279081563687830782748305746187060264985869283524441819589592750998086186315250781067131293823177124077445718802216415539934838376431091001197641295264650596195201747790167311", 10);
826         mpz_set_str(p->n2, "7", 10);
827     } else if (group_size <= 892) {
828         p->m = 617;
829         p->t = 88;
830         mpz_set_str(p->n, "57591959284219511220590893724691916802833742568034971006633345422620650391172287893878655658086794200963521584019889327992536532560877385225451713282279597074750857647455565899702728629166541223955196002755787520206774906606158388947359746178875040401304783332742806641", 10);
831         mpz_set_str(p->n2, "42019638181715250622338241", 10);
832     } else
833         pbc_die("unsupported group size");
834 }
835 
836