1 #include <stdarg.h>
2 #include <stdio.h>
3 #include <stdint.h> // for intptr_t
4 #include <stdlib.h>
5 #include <gmp.h>
6 #include "pbc_utils.h"
7 #include "pbc_field.h"
8 #include "pbc_curve.h"
9 #include "pbc_param.h"
10 #include "pbc_pairing.h"
11 #include "pbc_fp.h"
12 #include "pbc_memory.h"
13
14 //TODO: Store as integer mod ring instead and convert at last minute?
15 struct point_s {
16 int inf_flag;
17 element_t x;
18 element_t y;
19 };
20 typedef struct point_s *point_ptr;
21 typedef struct point_s point_t[1];
22
sn_init(element_ptr e)23 static void sn_init(element_ptr e) {
24 field_ptr f = e->field->data;
25 e->data = pbc_malloc(sizeof(point_t));
26 point_ptr p = e->data;
27 element_init(p->x, f);
28 element_init(p->y, f);
29 p->inf_flag = 1;
30 }
31
sn_clear(element_ptr e)32 static void sn_clear(element_ptr e) {
33 point_ptr p = e->data;
34 element_clear(p->x);
35 element_clear(p->y);
36 pbc_free(e->data);
37 }
38
sn_set0(element_ptr x)39 static void sn_set0(element_ptr x) {
40 point_ptr p = x->data;
41 p->inf_flag = 1;
42 }
43
sn_is0(element_ptr x)44 static int sn_is0(element_ptr x) {
45 point_ptr p = x->data;
46 return p->inf_flag;
47 }
48
49 //singular with node: y^2 = x^3 + x^2
sn_random(element_t a)50 static void sn_random(element_t a) {
51 point_ptr p = a->data;
52 element_t t;
53
54 element_init(t, p->x->field);
55 p->inf_flag = 0;
56 do {
57 element_random(p->x);
58 if (element_is0(p->x)) continue;
59 element_square(t, p->x);
60 element_add(t, t, p->x);
61 element_mul(t, t, p->x);
62 } while (!element_is_sqr(t));
63 element_sqrt(p->y, t);
64
65 element_clear(t);
66 }
67
sn_double_no_check(point_ptr r,point_ptr p)68 static inline void sn_double_no_check(point_ptr r, point_ptr p) {
69 element_t lambda, e0, e1;
70
71 element_init(lambda, p->x->field);
72 element_init(e0, p->x->field);
73 element_init(e1, p->x->field);
74 //same point: double them
75
76 //lambda = (3x^2 + 2x) / 2y
77 element_mul_si(lambda, p->x, 3);
78 element_set_si(e0, 2);
79 element_add(lambda, lambda, e0);
80 element_mul(lambda, lambda, p->x);
81 element_add(e0, p->y, p->y);
82 element_invert(e0, e0);
83 element_mul(lambda, lambda, e0);
84 //x1 = lambda^2 - 2x - 1
85 element_add(e1, p->x, p->x);
86 element_square(e0, lambda);
87 element_sub(e0, e0, e1);
88 element_set_si(e1, 1);
89 element_sub(e0, e0, e1);
90 //y1 = (x - x1)lambda - y
91 element_sub(e1, p->x, e0);
92 element_mul(e1, e1, lambda);
93 element_sub(e1, e1, p->y);
94
95 element_set(r->x, e0);
96 element_set(r->y, e1);
97 r->inf_flag = 0;
98
99 element_clear(lambda);
100 element_clear(e0);
101 element_clear(e1);
102 return;
103 }
104
sn_double(element_t c,element_t a)105 static void sn_double(element_t c, element_t a) {
106 point_ptr r = c->data;
107 point_ptr p = a->data;
108 if (p->inf_flag) {
109 r->inf_flag = 1;
110 return;
111 }
112 if (element_is0(p->y)) {
113 r->inf_flag = 1;
114 return;
115 }
116 sn_double_no_check(r, p);
117 }
118
sn_set(element_ptr c,element_ptr a)119 static void sn_set(element_ptr c, element_ptr a) {
120 point_ptr r = c->data, p = a->data;
121 if (p->inf_flag) {
122 r->inf_flag = 1;
123 return;
124 }
125 r->inf_flag = 0;
126 element_set(r->x, p->x);
127 element_set(r->y, p->y);
128 }
129
sn_add(element_t c,element_t a,element_t b)130 static void sn_add(element_t c, element_t a, element_t b) {
131 point_ptr r = c->data;
132 point_ptr p = a->data;
133 point_ptr q = b->data;
134 if (p->inf_flag) {
135 sn_set(c, b);
136 return;
137 }
138 if (q->inf_flag) {
139 sn_set(c, a);
140 return;
141 }
142 if (!element_cmp(p->x, q->x)) {
143 if (!element_cmp(p->y, q->y)) {
144 if (element_is0(p->y)) {
145 r->inf_flag = 1;
146 return;
147 } else {
148 sn_double_no_check(r, p);
149 return;
150 }
151 }
152 //points are inverses of each other
153 r->inf_flag = 1;
154 return;
155 } else {
156 element_t lambda, e0, e1;
157
158 element_init(lambda, p->x->field);
159 element_init(e0, p->x->field);
160 element_init(e1, p->x->field);
161
162 //lambda = (y2-y1)/(x2-x1)
163 element_sub(e0, q->x, p->x);
164 element_invert(e0, e0);
165 element_sub(lambda, q->y, p->y);
166 element_mul(lambda, lambda, e0);
167 //x3 = lambda^2 - x1 - x2 - 1
168 element_square(e0, lambda);
169 element_sub(e0, e0, p->x);
170 element_sub(e0, e0, q->x);
171 element_set1(e1);
172 element_sub(e0, e0, e1);
173 //y3 = (x1-x3)lambda - y1
174 element_sub(e1, p->x, e0);
175 element_mul(e1, e1, lambda);
176 element_sub(e1, e1, p->y);
177
178 element_set(r->x, e0);
179 element_set(r->y, e1);
180 r->inf_flag = 0;
181
182 element_clear(lambda);
183 element_clear(e0);
184 element_clear(e1);
185 }
186 }
187
sn_invert(element_ptr c,element_ptr a)188 static void sn_invert(element_ptr c, element_ptr a) {
189 point_ptr r = c->data, p = a->data;
190
191 if (p->inf_flag) {
192 r->inf_flag = 1;
193 return;
194 }
195 r->inf_flag = 0;
196 element_set(r->x, p->x);
197 element_neg(r->y, p->y);
198 }
199
sn_field_clear(field_ptr c)200 static void sn_field_clear(field_ptr c) {
201 UNUSED_VAR(c);
202 }
203
204 /* TODO: Write a test program that uses these functions.
205
206 // Nonsingular points on sn curves map to finite field elements via
207 // (x, y) --> (y + x)/(y - x)
208 // The reverse map is
209 // a --> (4a/(a-1)^2, 4a(a+1)/(a-1)^3)
210
211 void sn_point_to_field(element_t out, point_ptr P) {
212 element_t e0, e1;
213 if (P->inf_flag) {
214 element_set1(out);
215 return;
216 }
217 element_init(e0, out->field);
218 element_init(e1, out->field);
219 element_add(e0, P->y, P->x);
220 element_sub(e1, P->y, P->x);
221 element_invert(e1, e1);
222 element_mul(out, e0, e1);
223 element_clear(e0);
224 element_clear(e1);
225 }
226
227 static void sn_field_to_point(point_ptr P, element_t in) {
228 element_t e0, e1, e2;
229
230 if (element_is1(in)) {
231 P->inf_flag = 1;
232 return;
233 }
234 element_init(e0, in->field);
235 element_init(e1, in->field);
236 element_init(e2, in->field);
237
238 element_set1(e1);
239 element_sub(e0, in, e1);
240 element_invert(e0, e0);
241
242 element_mul_si(e2, in, 4);
243
244 element_add(P->y, in, e1);
245
246 element_mul(e1, e0, e0);
247 element_mul(P->x, e1, e2);
248 element_mul(P->y, P->y, e2);
249 element_mul(P->y, P->y, e0);
250 element_mul(P->y, P->y, e1);
251 P->inf_flag = 0;
252
253 element_clear(e0);
254 element_clear(e1);
255 element_clear(e2);
256 }
257 */
258
sn_out_str(FILE * stream,int base,element_ptr a)259 static size_t sn_out_str(FILE *stream, int base, element_ptr a) {
260 point_ptr p = a->data;
261 size_t result, status;
262 if (p->inf_flag) {
263 if (EOF == fputc('O', stream)) return 0;
264 return 1;
265 }
266 result = element_out_str(stream, base, p->x);
267 if (!result) return 0;
268 if (EOF == fputc(' ', stream)) return 0;
269 status = element_out_str(stream, base, p->y);
270 if (!status) return 0;
271 return result + status + 1;
272 }
273
274 void naive_generic_pow_mpz(element_ptr x, element_ptr a, mpz_ptr n);
field_init_curve_singular_with_node(field_t c,field_t field)275 void field_init_curve_singular_with_node(field_t c, field_t field) {
276 mpz_set(c->order, field->order);
277 c->data = (void *) field;
278 c->init = sn_init;
279 c->clear = sn_clear;
280 c->random = sn_random;
281 //c->from_x = cc_from_x;
282 //c->from_hash = cc_from_hash;
283 c->set = sn_set;
284 c->invert = c->neg = sn_invert;
285 c->square = c->doub = sn_double;
286 c->mul = c->add = sn_add;
287 c->set1 = c->set0 = sn_set0;
288 c->is1 = c->is0 = sn_is0;
289 c->mul_mpz = element_pow_mpz;
290 c->out_str = sn_out_str;
291 c->field_clear = sn_field_clear;
292 }
293
294 //TODO: the following code is useless as the Tate pairing is degenerate on singular curves
sn_miller(element_t res,mpz_t q,element_t P,element_ptr Qx,element_ptr Qy)295 static void sn_miller(element_t res, mpz_t q, element_t P,
296 element_ptr Qx, element_ptr Qy) {
297 //collate divisions
298 int m;
299 element_t v, vd;
300 element_t Z;
301 element_t a, b, c;
302 element_t e0, e1;
303 element_ptr Zx;
304 element_ptr Zy;
305 const element_ptr Px = curve_x_coord(P);
306 const element_ptr Py = curve_y_coord(P);
307
308 #define do_vertical(e) \
309 element_sub(e0, Qx, Zx); \
310 element_mul(e, e, e0);
311
312 //a = -slope_tangent(Z.x, Z.y);
313 //b = 1;
314 //c = -(Z.y + a * Z.x);
315 //but we multiply by 2*Z.y to avoid division
316 //a = -Zx * (Zx + Zx + Zx + 2)
317 //b = 2 * Zy
318 //c = -(2 Zy^2 + a Zx);
319 #define do_tangent(e) \
320 element_double(e0, Zx); \
321 element_add(a, Zx, e0); \
322 element_set_si(e0, 2); \
323 element_add(a, a, e0); \
324 element_mul(a, a, Zx); \
325 element_neg(a, a); \
326 element_add(b, Zy, Zy); \
327 element_mul(e0, b, Zy); \
328 element_mul(c, a, Zx); \
329 element_add(c, c, e0); \
330 element_neg(c, c); \
331 element_mul(e0, a, Qx); \
332 element_mul(e1, b, Qy); \
333 element_add(e0, e0, e1); \
334 element_add(e0, e0, c); \
335 element_mul(e, e, e0);
336
337 //a = -(B.y - A.y) / (B.x - A.x);
338 //b = 1;
339 //c = -(A.y + a * A.x);
340 //but we'll multiply by B.x - A.x to avoid division
341 #define do_line(e) \
342 element_sub(b, Px, Zx); \
343 element_sub(a, Zy, Py); \
344 element_mul(e0, b, Zy); \
345 element_mul(c, a, Zx); \
346 element_add(c, c, e0); \
347 element_neg(c, c); \
348 element_mul(e0, a, Qx); \
349 element_mul(e1, b, Qy); \
350 element_add(e0, e0, e1); \
351 element_add(e0, e0, c); \
352 element_mul(e, e, e0);
353
354 element_init(a, Px->field);
355 element_init(b, Px->field);
356 element_init(c, Px->field);
357 element_init(e0, res->field);
358 element_init(e1, res->field);
359
360 element_init(v, res->field);
361 element_init(vd, res->field);
362 element_init(Z, P->field);
363
364 element_set(Z, P);
365 Zx = curve_x_coord(Z);
366 Zy = curve_y_coord(Z);
367
368 element_set1(v);
369 element_set1(vd);
370 m = mpz_sizeinbase(q, 2) - 2;
371
372 while(m >= 0) {
373 element_mul(v, v, v);
374 element_mul(vd, vd, vd);
375 do_tangent(v);
376 element_double(Z, Z);
377 do_vertical(vd);
378 if (mpz_tstbit(q, m)) {
379 do_line(v);
380 element_add(Z, Z, P);
381 do_vertical(vd);
382 }
383 m--;
384 }
385 #undef do_tangent
386 #undef do_vertical
387 #undef do_line
388
389 element_invert(vd, vd);
390 element_mul(res, v, vd);
391
392 element_clear(v);
393 element_clear(vd);
394 element_clear(Z);
395 element_clear(a);
396 element_clear(b);
397 element_clear(c);
398 element_clear(e0);
399 element_clear(e1);
400 }
401
402 struct sn_pairing_data_s {
403 field_t Fq, Eq;
404 };
405 typedef struct sn_pairing_data_s sn_pairing_data_t[1];
406 typedef struct sn_pairing_data_s *sn_pairing_data_ptr;
407
sn_pairing(element_ptr out,element_ptr in1,element_ptr in2,pairing_t pairing)408 static void sn_pairing(element_ptr out, element_ptr in1, element_ptr in2,
409 pairing_t pairing) {
410 sn_pairing_data_ptr p = pairing->data;
411 element_ptr Q = in2;
412 element_t e0;
413 element_t R, QR;
414 element_init(R, p->Eq);
415 element_init(QR, p->Eq);
416 element_random(R);
417 element_init(e0, out->field);
418 element_add(QR, Q, R);
419 sn_miller(out, pairing->r, in1, curve_x_coord(QR), curve_y_coord(QR));
420 sn_miller(e0, pairing->r, in1, curve_x_coord(R), curve_y_coord(R));
421 element_invert(e0, e0);
422 element_mul(out, out, e0);
423 //element_pow_mpz(out, out, p->tateexp);
424 element_clear(R);
425 element_clear(QR);
426 }
427
pairing_init_singular_with_node(pairing_t pairing,mpz_t q)428 void pairing_init_singular_with_node(pairing_t pairing, mpz_t q) {
429 sn_pairing_data_ptr p;
430
431 mpz_init(pairing->r);
432 mpz_sub_ui(pairing->r, q, 1);
433 field_init_fp(pairing->Zr, pairing->r);
434 pairing->map = sn_pairing;
435
436 p = pairing->data = pbc_malloc(sizeof(sn_pairing_data_t));
437 field_init_fp(p->Fq, q);
438 field_init_curve_singular_with_node(p->Eq, p->Fq);
439
440 //mpz_init(p->tateexp);
441 //mpz_sub_ui(p->tateexp, p->Fq->order, 1);
442 //mpz_divexact(p->tateexp, p->tateexp, pairing->r);
443
444 pairing->G2 = pairing->G1 = p->Eq;
445
446 pairing_GT_init(pairing, p->Fq);
447 }
448