1 #include <assert.h>
2
3 #include "ssh.h"
4 #include "mpint.h"
5 #include "ecc.h"
6
7 /* ----------------------------------------------------------------------
8 * Weierstrass curves.
9 */
10
11 struct WeierstrassPoint {
12 /*
13 * Internally, we represent a point using 'Jacobian coordinates',
14 * which are three values X,Y,Z whose relation to the affine
15 * coordinates x,y is that x = X/Z^2 and y = Y/Z^3.
16 *
17 * This allows us to do most of our calculations without having to
18 * take an inverse mod p: every time the obvious affine formulae
19 * would need you to divide by something, you instead multiply it
20 * into the 'denominator' coordinate Z. You only have to actually
21 * take the inverse of Z when you need to get the affine
22 * coordinates back out, which means you do it once after your
23 * entire computation instead of at every intermediate step.
24 *
25 * The point at infinity is represented by setting all three
26 * coordinates to zero.
27 *
28 * These values are also stored in the Montgomery-multiplication
29 * transformed representation.
30 */
31 mp_int *X, *Y, *Z;
32
33 WeierstrassCurve *wc;
34 };
35
36 struct WeierstrassCurve {
37 /* Prime modulus of the finite field. */
38 mp_int *p;
39
40 /* Persistent Montgomery context for doing arithmetic mod p. */
41 MontyContext *mc;
42
43 /* Modsqrt context for point decompression. NULL if this curve was
44 * constructed without providing nonsquare_mod_p. */
45 ModsqrtContext *sc;
46
47 /* Parameters of the curve, in Montgomery-multiplication
48 * transformed form. */
49 mp_int *a, *b;
50 };
51
ecc_weierstrass_curve(mp_int * p,mp_int * a,mp_int * b,mp_int * nonsquare_mod_p)52 WeierstrassCurve *ecc_weierstrass_curve(
53 mp_int *p, mp_int *a, mp_int *b, mp_int *nonsquare_mod_p)
54 {
55 WeierstrassCurve *wc = snew(WeierstrassCurve);
56 wc->p = mp_copy(p);
57 wc->mc = monty_new(p);
58 wc->a = monty_import(wc->mc, a);
59 wc->b = monty_import(wc->mc, b);
60
61 if (nonsquare_mod_p)
62 wc->sc = modsqrt_new(p, nonsquare_mod_p);
63 else
64 wc->sc = NULL;
65
66 return wc;
67 }
68
ecc_weierstrass_curve_free(WeierstrassCurve * wc)69 void ecc_weierstrass_curve_free(WeierstrassCurve *wc)
70 {
71 mp_free(wc->p);
72 mp_free(wc->a);
73 mp_free(wc->b);
74 monty_free(wc->mc);
75 if (wc->sc)
76 modsqrt_free(wc->sc);
77 sfree(wc);
78 }
79
ecc_weierstrass_point_new_empty(WeierstrassCurve * wc)80 static WeierstrassPoint *ecc_weierstrass_point_new_empty(WeierstrassCurve *wc)
81 {
82 WeierstrassPoint *wp = snew(WeierstrassPoint);
83 wp->wc = wc;
84 wp->X = wp->Y = wp->Z = NULL;
85 return wp;
86 }
87
ecc_weierstrass_point_new_imported(WeierstrassCurve * wc,mp_int * monty_x,mp_int * monty_y)88 static WeierstrassPoint *ecc_weierstrass_point_new_imported(
89 WeierstrassCurve *wc, mp_int *monty_x, mp_int *monty_y)
90 {
91 WeierstrassPoint *wp = ecc_weierstrass_point_new_empty(wc);
92 wp->X = monty_x;
93 wp->Y = monty_y;
94 wp->Z = mp_copy(monty_identity(wc->mc));
95 return wp;
96 }
97
ecc_weierstrass_point_new(WeierstrassCurve * wc,mp_int * x,mp_int * y)98 WeierstrassPoint *ecc_weierstrass_point_new(
99 WeierstrassCurve *wc, mp_int *x, mp_int *y)
100 {
101 return ecc_weierstrass_point_new_imported(
102 wc, monty_import(wc->mc, x), monty_import(wc->mc, y));
103 }
104
ecc_weierstrass_point_new_identity(WeierstrassCurve * wc)105 WeierstrassPoint *ecc_weierstrass_point_new_identity(WeierstrassCurve *wc)
106 {
107 WeierstrassPoint *wp = ecc_weierstrass_point_new_empty(wc);
108 size_t bits = mp_max_bits(wc->p);
109 wp->X = mp_new(bits);
110 wp->Y = mp_new(bits);
111 wp->Z = mp_new(bits);
112 return wp;
113 }
114
ecc_weierstrass_point_copy_into(WeierstrassPoint * dest,WeierstrassPoint * src)115 void ecc_weierstrass_point_copy_into(
116 WeierstrassPoint *dest, WeierstrassPoint *src)
117 {
118 mp_copy_into(dest->X, src->X);
119 mp_copy_into(dest->Y, src->Y);
120 mp_copy_into(dest->Z, src->Z);
121 }
122
ecc_weierstrass_point_copy(WeierstrassPoint * orig)123 WeierstrassPoint *ecc_weierstrass_point_copy(WeierstrassPoint *orig)
124 {
125 WeierstrassPoint *wp = ecc_weierstrass_point_new_empty(orig->wc);
126 wp->X = mp_copy(orig->X);
127 wp->Y = mp_copy(orig->Y);
128 wp->Z = mp_copy(orig->Z);
129 return wp;
130 }
131
ecc_weierstrass_point_free(WeierstrassPoint * wp)132 void ecc_weierstrass_point_free(WeierstrassPoint *wp)
133 {
134 mp_free(wp->X);
135 mp_free(wp->Y);
136 mp_free(wp->Z);
137 smemclr(wp, sizeof(*wp));
138 sfree(wp);
139 }
140
ecc_weierstrass_point_new_from_x(WeierstrassCurve * wc,mp_int * xorig,unsigned desired_y_parity)141 WeierstrassPoint *ecc_weierstrass_point_new_from_x(
142 WeierstrassCurve *wc, mp_int *xorig, unsigned desired_y_parity)
143 {
144 assert(wc->sc);
145
146 /*
147 * The curve equation is y^2 = x^3 + ax + b, which is already
148 * conveniently in a form where we can compute the RHS and take
149 * the square root of it to get y.
150 */
151 unsigned success;
152
153 mp_int *x = monty_import(wc->mc, xorig);
154
155 /*
156 * Compute the RHS of the curve equation. We don't need to take
157 * account of z here, because we're constructing the point from
158 * scratch. So it really is just x^3 + ax + b.
159 */
160 mp_int *x2 = monty_mul(wc->mc, x, x);
161 mp_int *x2_plus_a = monty_add(wc->mc, x2, wc->a);
162 mp_int *x3_plus_ax = monty_mul(wc->mc, x2_plus_a, x);
163 mp_int *rhs = monty_add(wc->mc, x3_plus_ax, wc->b);
164 mp_free(x2);
165 mp_free(x2_plus_a);
166 mp_free(x3_plus_ax);
167
168 mp_int *y = monty_modsqrt(wc->sc, rhs, &success);
169 mp_free(rhs);
170
171 if (!success) {
172 /* Failure! x^3+ax+b worked out to be a number that has no
173 * square root mod p. In this situation there's no point in
174 * trying to be time-constant, since the protocol sequence is
175 * going to diverge anyway when we complain to whoever gave us
176 * this bogus value. */
177 mp_free(x);
178 mp_free(y);
179 return NULL;
180 }
181
182 /*
183 * Choose whichever of y and p-y has the specified parity (of its
184 * lowest positive residue mod p).
185 */
186 mp_int *tmp = monty_export(wc->mc, y);
187 unsigned flip = (mp_get_bit(tmp, 0) ^ desired_y_parity) & 1;
188 mp_sub_into(tmp, wc->p, y);
189 mp_select_into(y, y, tmp, flip);
190 mp_free(tmp);
191
192 return ecc_weierstrass_point_new_imported(wc, x, y);
193 }
194
ecc_weierstrass_cond_overwrite(WeierstrassPoint * dest,WeierstrassPoint * src,unsigned overwrite)195 static void ecc_weierstrass_cond_overwrite(
196 WeierstrassPoint *dest, WeierstrassPoint *src, unsigned overwrite)
197 {
198 mp_select_into(dest->X, dest->X, src->X, overwrite);
199 mp_select_into(dest->Y, dest->Y, src->Y, overwrite);
200 mp_select_into(dest->Z, dest->Z, src->Z, overwrite);
201 }
202
ecc_weierstrass_cond_swap(WeierstrassPoint * P,WeierstrassPoint * Q,unsigned swap)203 static void ecc_weierstrass_cond_swap(
204 WeierstrassPoint *P, WeierstrassPoint *Q, unsigned swap)
205 {
206 mp_cond_swap(P->X, Q->X, swap);
207 mp_cond_swap(P->Y, Q->Y, swap);
208 mp_cond_swap(P->Z, Q->Z, swap);
209 }
210
211 /*
212 * Shared code between all three of the basic arithmetic functions:
213 * once we've determined the slope of the line that we're intersecting
214 * the curve with, this takes care of finding the coordinates of the
215 * third intersection point (given the two input x-coordinates and one
216 * of the y-coords) and negating it to generate the output.
217 */
ecc_weierstrass_epilogue(mp_int * Px,mp_int * Qx,mp_int * Py,mp_int * common_Z,mp_int * lambda_n,mp_int * lambda_d,WeierstrassPoint * out)218 static inline void ecc_weierstrass_epilogue(
219 mp_int *Px, mp_int *Qx, mp_int *Py, mp_int *common_Z,
220 mp_int *lambda_n, mp_int *lambda_d, WeierstrassPoint *out)
221 {
222 WeierstrassCurve *wc = out->wc;
223
224 /* Powers of the numerator and denominator of the slope lambda */
225 mp_int *lambda_n2 = monty_mul(wc->mc, lambda_n, lambda_n);
226 mp_int *lambda_d2 = monty_mul(wc->mc, lambda_d, lambda_d);
227 mp_int *lambda_d3 = monty_mul(wc->mc, lambda_d, lambda_d2);
228
229 /* Make the output x-coordinate */
230 mp_int *xsum = monty_add(wc->mc, Px, Qx);
231 mp_int *lambda_d2_xsum = monty_mul(wc->mc, lambda_d2, xsum);
232 out->X = monty_sub(wc->mc, lambda_n2, lambda_d2_xsum);
233
234 /* Make the output y-coordinate */
235 mp_int *lambda_d2_Px = monty_mul(wc->mc, lambda_d2, Px);
236 mp_int *xdiff = monty_sub(wc->mc, lambda_d2_Px, out->X);
237 mp_int *lambda_n_xdiff = monty_mul(wc->mc, lambda_n, xdiff);
238 mp_int *lambda_d3_Py = monty_mul(wc->mc, lambda_d3, Py);
239 out->Y = monty_sub(wc->mc, lambda_n_xdiff, lambda_d3_Py);
240
241 /* Make the output z-coordinate */
242 out->Z = monty_mul(wc->mc, common_Z, lambda_d);
243
244 mp_free(lambda_n2);
245 mp_free(lambda_d2);
246 mp_free(lambda_d3);
247 mp_free(xsum);
248 mp_free(xdiff);
249 mp_free(lambda_d2_xsum);
250 mp_free(lambda_n_xdiff);
251 mp_free(lambda_d2_Px);
252 mp_free(lambda_d3_Py);
253 }
254
255 /*
256 * Shared code between add and add_general: put the two input points
257 * over a common denominator, and determine the slope lambda of the
258 * line through both of them. If the points have the same
259 * x-coordinate, then the slope will be returned with a zero
260 * denominator.
261 */
ecc_weierstrass_add_prologue(WeierstrassPoint * P,WeierstrassPoint * Q,mp_int ** Px,mp_int ** Py,mp_int ** Qx,mp_int ** denom,mp_int ** lambda_n,mp_int ** lambda_d)262 static inline void ecc_weierstrass_add_prologue(
263 WeierstrassPoint *P, WeierstrassPoint *Q,
264 mp_int **Px, mp_int **Py, mp_int **Qx, mp_int **denom,
265 mp_int **lambda_n, mp_int **lambda_d)
266 {
267 WeierstrassCurve *wc = P->wc;
268
269 /* Powers of the points' denominators */
270 mp_int *Pz2 = monty_mul(wc->mc, P->Z, P->Z);
271 mp_int *Pz3 = monty_mul(wc->mc, Pz2, P->Z);
272 mp_int *Qz2 = monty_mul(wc->mc, Q->Z, Q->Z);
273 mp_int *Qz3 = monty_mul(wc->mc, Qz2, Q->Z);
274
275 /* Points' x,y coordinates scaled by the other one's denominator
276 * (raised to the appropriate power) */
277 *Px = monty_mul(wc->mc, P->X, Qz2);
278 *Py = monty_mul(wc->mc, P->Y, Qz3);
279 *Qx = monty_mul(wc->mc, Q->X, Pz2);
280 mp_int *Qy = monty_mul(wc->mc, Q->Y, Pz3);
281
282 /* Common denominator */
283 *denom = monty_mul(wc->mc, P->Z, Q->Z);
284
285 /* Slope of the line through the two points, if P != Q */
286 *lambda_n = monty_sub(wc->mc, Qy, *Py);
287 *lambda_d = monty_sub(wc->mc, *Qx, *Px);
288
289 mp_free(Pz2);
290 mp_free(Pz3);
291 mp_free(Qz2);
292 mp_free(Qz3);
293 mp_free(Qy);
294 }
295
ecc_weierstrass_add(WeierstrassPoint * P,WeierstrassPoint * Q)296 WeierstrassPoint *ecc_weierstrass_add(WeierstrassPoint *P, WeierstrassPoint *Q)
297 {
298 WeierstrassCurve *wc = P->wc;
299 assert(Q->wc == wc);
300
301 WeierstrassPoint *S = ecc_weierstrass_point_new_empty(wc);
302
303 mp_int *Px, *Py, *Qx, *denom, *lambda_n, *lambda_d;
304 ecc_weierstrass_add_prologue(
305 P, Q, &Px, &Py, &Qx, &denom, &lambda_n, &lambda_d);
306
307 /* Never expect to have received two mutually inverse inputs, or
308 * two identical ones (which would make this a doubling). In other
309 * words, the two input x-coordinates (after putting over a common
310 * denominator) should never have been equal. */
311 assert(!mp_eq_integer(lambda_n, 0));
312
313 /* Now go to the common epilogue code. */
314 ecc_weierstrass_epilogue(Px, Qx, Py, denom, lambda_n, lambda_d, S);
315
316 mp_free(Px);
317 mp_free(Py);
318 mp_free(Qx);
319 mp_free(denom);
320 mp_free(lambda_n);
321 mp_free(lambda_d);
322
323 return S;
324 }
325
326 /*
327 * Code to determine the slope of the line you need to intersect with
328 * the curve in the case where you're adding a point to itself. In
329 * this situation you can't just say "the line through both input
330 * points" because that's under-determined; instead, you have to take
331 * the _tangent_ to the curve at the given point, by differentiating
332 * the curve equation y^2=x^3+ax+b to get 2y dy/dx = 3x^2+a.
333 */
ecc_weierstrass_tangent_slope(WeierstrassPoint * P,mp_int ** lambda_n,mp_int ** lambda_d)334 static inline void ecc_weierstrass_tangent_slope(
335 WeierstrassPoint *P, mp_int **lambda_n, mp_int **lambda_d)
336 {
337 WeierstrassCurve *wc = P->wc;
338
339 mp_int *X2 = monty_mul(wc->mc, P->X, P->X);
340 mp_int *twoX2 = monty_add(wc->mc, X2, X2);
341 mp_int *threeX2 = monty_add(wc->mc, twoX2, X2);
342 mp_int *Z2 = monty_mul(wc->mc, P->Z, P->Z);
343 mp_int *Z4 = monty_mul(wc->mc, Z2, Z2);
344 mp_int *aZ4 = monty_mul(wc->mc, wc->a, Z4);
345
346 *lambda_n = monty_add(wc->mc, threeX2, aZ4);
347 *lambda_d = monty_add(wc->mc, P->Y, P->Y);
348
349 mp_free(X2);
350 mp_free(twoX2);
351 mp_free(threeX2);
352 mp_free(Z2);
353 mp_free(Z4);
354 mp_free(aZ4);
355 }
356
ecc_weierstrass_double(WeierstrassPoint * P)357 WeierstrassPoint *ecc_weierstrass_double(WeierstrassPoint *P)
358 {
359 WeierstrassCurve *wc = P->wc;
360 WeierstrassPoint *D = ecc_weierstrass_point_new_empty(wc);
361
362 mp_int *lambda_n, *lambda_d;
363 ecc_weierstrass_tangent_slope(P, &lambda_n, &lambda_d);
364 ecc_weierstrass_epilogue(P->X, P->X, P->Y, P->Z, lambda_n, lambda_d, D);
365 mp_free(lambda_n);
366 mp_free(lambda_d);
367
368 return D;
369 }
370
ecc_weierstrass_select_into(WeierstrassPoint * dest,WeierstrassPoint * P,WeierstrassPoint * Q,unsigned choose_Q)371 static inline void ecc_weierstrass_select_into(
372 WeierstrassPoint *dest, WeierstrassPoint *P, WeierstrassPoint *Q,
373 unsigned choose_Q)
374 {
375 mp_select_into(dest->X, P->X, Q->X, choose_Q);
376 mp_select_into(dest->Y, P->Y, Q->Y, choose_Q);
377 mp_select_into(dest->Z, P->Z, Q->Z, choose_Q);
378 }
379
ecc_weierstrass_add_general(WeierstrassPoint * P,WeierstrassPoint * Q)380 WeierstrassPoint *ecc_weierstrass_add_general(
381 WeierstrassPoint *P, WeierstrassPoint *Q)
382 {
383 WeierstrassCurve *wc = P->wc;
384 assert(Q->wc == wc);
385
386 WeierstrassPoint *S = ecc_weierstrass_point_new_empty(wc);
387
388 /* Parameters for the epilogue, and slope of the line if P != Q */
389 mp_int *Px, *Py, *Qx, *denom, *lambda_n, *lambda_d;
390 ecc_weierstrass_add_prologue(
391 P, Q, &Px, &Py, &Qx, &denom, &lambda_n, &lambda_d);
392
393 /* Slope if P == Q */
394 mp_int *lambda_n_tangent, *lambda_d_tangent;
395 ecc_weierstrass_tangent_slope(P, &lambda_n_tangent, &lambda_d_tangent);
396
397 /* Select between those slopes depending on whether P == Q */
398 unsigned same_x_coord = mp_eq_integer(lambda_d, 0);
399 unsigned same_y_coord = mp_eq_integer(lambda_n, 0);
400 unsigned equality = same_x_coord & same_y_coord;
401 mp_select_into(lambda_n, lambda_n, lambda_n_tangent, equality);
402 mp_select_into(lambda_d, lambda_d, lambda_d_tangent, equality);
403
404 /* Now go to the common code between addition and doubling */
405 ecc_weierstrass_epilogue(Px, Qx, Py, denom, lambda_n, lambda_d, S);
406
407 /* Check for the input identity cases, and overwrite the output if
408 * necessary. */
409 ecc_weierstrass_select_into(S, S, Q, mp_eq_integer(P->Z, 0));
410 ecc_weierstrass_select_into(S, S, P, mp_eq_integer(Q->Z, 0));
411
412 /*
413 * In the case where P == -Q and so the output is the identity,
414 * we'll have calculated lambda_d = 0 and so the output will have
415 * z==0 already. Detect that and use it to normalise the other two
416 * coordinates to zero.
417 */
418 unsigned output_id = mp_eq_integer(S->Z, 0);
419 mp_cond_clear(S->X, output_id);
420 mp_cond_clear(S->Y, output_id);
421
422 mp_free(Px);
423 mp_free(Py);
424 mp_free(Qx);
425 mp_free(denom);
426 mp_free(lambda_n);
427 mp_free(lambda_d);
428 mp_free(lambda_n_tangent);
429 mp_free(lambda_d_tangent);
430
431 return S;
432 }
433
ecc_weierstrass_multiply(WeierstrassPoint * B,mp_int * n)434 WeierstrassPoint *ecc_weierstrass_multiply(WeierstrassPoint *B, mp_int *n)
435 {
436 WeierstrassPoint *two_B = ecc_weierstrass_double(B);
437 WeierstrassPoint *k_B = ecc_weierstrass_point_copy(B);
438 WeierstrassPoint *kplus1_B = ecc_weierstrass_point_copy(two_B);
439
440 /*
441 * This multiply routine more or less follows the shape of the
442 * 'Montgomery ladder' technique that you have to use under the
443 * extra constraint on addition in Montgomery curves, because it
444 * was fresh in my mind and easier to just do it the same way. See
445 * the comment in ecc_montgomery_multiply.
446 */
447
448 unsigned not_started_yet = 1;
449 for (size_t bitindex = mp_max_bits(n); bitindex-- > 0 ;) {
450 unsigned nbit = mp_get_bit(n, bitindex);
451
452 WeierstrassPoint *sum = ecc_weierstrass_add(k_B, kplus1_B);
453 ecc_weierstrass_cond_swap(k_B, kplus1_B, nbit);
454 WeierstrassPoint *other = ecc_weierstrass_double(k_B);
455 ecc_weierstrass_point_free(k_B);
456 ecc_weierstrass_point_free(kplus1_B);
457 k_B = other;
458 kplus1_B = sum;
459 ecc_weierstrass_cond_swap(k_B, kplus1_B, nbit);
460
461 ecc_weierstrass_cond_overwrite(k_B, B, not_started_yet);
462 ecc_weierstrass_cond_overwrite(kplus1_B, two_B, not_started_yet);
463 not_started_yet &= ~nbit;
464 }
465
466 ecc_weierstrass_point_free(two_B);
467 ecc_weierstrass_point_free(kplus1_B);
468 return k_B;
469 }
470
ecc_weierstrass_is_identity(WeierstrassPoint * wp)471 unsigned ecc_weierstrass_is_identity(WeierstrassPoint *wp)
472 {
473 return mp_eq_integer(wp->Z, 0);
474 }
475
476 /*
477 * Normalise a point by scaling its Jacobian coordinates so that Z=1.
478 * This doesn't change what point is represented by the triple, but it
479 * means the affine x,y can now be easily recovered from X and Y.
480 */
ecc_weierstrass_normalise(WeierstrassPoint * wp)481 static void ecc_weierstrass_normalise(WeierstrassPoint *wp)
482 {
483 WeierstrassCurve *wc = wp->wc;
484 mp_int *zinv = monty_invert(wc->mc, wp->Z);
485 mp_int *zinv2 = monty_mul(wc->mc, zinv, zinv);
486 mp_int *zinv3 = monty_mul(wc->mc, zinv2, zinv);
487 monty_mul_into(wc->mc, wp->X, wp->X, zinv2);
488 monty_mul_into(wc->mc, wp->Y, wp->Y, zinv3);
489 monty_mul_into(wc->mc, wp->Z, wp->Z, zinv);
490 mp_free(zinv);
491 mp_free(zinv2);
492 mp_free(zinv3);
493 }
494
ecc_weierstrass_get_affine(WeierstrassPoint * wp,mp_int ** x,mp_int ** y)495 void ecc_weierstrass_get_affine(
496 WeierstrassPoint *wp, mp_int **x, mp_int **y)
497 {
498 WeierstrassCurve *wc = wp->wc;
499
500 ecc_weierstrass_normalise(wp);
501
502 if (x)
503 *x = monty_export(wc->mc, wp->X);
504 if (y)
505 *y = monty_export(wc->mc, wp->Y);
506 }
507
ecc_weierstrass_point_valid(WeierstrassPoint * P)508 unsigned ecc_weierstrass_point_valid(WeierstrassPoint *P)
509 {
510 WeierstrassCurve *wc = P->wc;
511
512 /*
513 * The projective version of the curve equation is
514 * Y^2 = X^3 + a X Z^4 + b Z^6
515 */
516 mp_int *lhs = monty_mul(P->wc->mc, P->Y, P->Y);
517 mp_int *x2 = monty_mul(wc->mc, P->X, P->X);
518 mp_int *x3 = monty_mul(wc->mc, x2, P->X);
519 mp_int *z2 = monty_mul(wc->mc, P->Z, P->Z);
520 mp_int *z4 = monty_mul(wc->mc, z2, z2);
521 mp_int *az4 = monty_mul(wc->mc, wc->a, z4);
522 mp_int *axz4 = monty_mul(wc->mc, az4, P->X);
523 mp_int *x3_plus_axz4 = monty_add(wc->mc, x3, axz4);
524 mp_int *z6 = monty_mul(wc->mc, z2, z4);
525 mp_int *bz6 = monty_mul(wc->mc, wc->b, z6);
526 mp_int *rhs = monty_add(wc->mc, x3_plus_axz4, bz6);
527
528 unsigned valid = mp_cmp_eq(lhs, rhs);
529
530 mp_free(lhs);
531 mp_free(x2);
532 mp_free(x3);
533 mp_free(z2);
534 mp_free(z4);
535 mp_free(az4);
536 mp_free(axz4);
537 mp_free(x3_plus_axz4);
538 mp_free(z6);
539 mp_free(bz6);
540 mp_free(rhs);
541
542 return valid;
543 }
544
545 /* ----------------------------------------------------------------------
546 * Montgomery curves.
547 */
548
549 struct MontgomeryPoint {
550 /* XZ coordinates. These represent the affine x coordinate by the
551 * relationship x = X/Z. */
552 mp_int *X, *Z;
553
554 MontgomeryCurve *mc;
555 };
556
557 struct MontgomeryCurve {
558 /* Prime modulus of the finite field. */
559 mp_int *p;
560
561 /* Montgomery context for arithmetic mod p. */
562 MontyContext *mc;
563
564 /* Parameters of the curve, in Montgomery-multiplication
565 * transformed form. */
566 mp_int *a, *b;
567
568 /* (a+2)/4, also in Montgomery-multiplication form. */
569 mp_int *aplus2over4;
570 };
571
ecc_montgomery_curve(mp_int * p,mp_int * a,mp_int * b)572 MontgomeryCurve *ecc_montgomery_curve(
573 mp_int *p, mp_int *a, mp_int *b)
574 {
575 MontgomeryCurve *mc = snew(MontgomeryCurve);
576 mc->p = mp_copy(p);
577 mc->mc = monty_new(p);
578 mc->a = monty_import(mc->mc, a);
579 mc->b = monty_import(mc->mc, b);
580
581 mp_int *four = mp_from_integer(4);
582 mp_int *fourinverse = mp_invert(four, mc->p);
583 mp_int *aplus2 = mp_copy(a);
584 mp_add_integer_into(aplus2, aplus2, 2);
585 mp_int *aplus2over4 = mp_modmul(aplus2, fourinverse, mc->p);
586 mc->aplus2over4 = monty_import(mc->mc, aplus2over4);
587 mp_free(four);
588 mp_free(fourinverse);
589 mp_free(aplus2);
590 mp_free(aplus2over4);
591
592 return mc;
593 }
594
ecc_montgomery_curve_free(MontgomeryCurve * mc)595 void ecc_montgomery_curve_free(MontgomeryCurve *mc)
596 {
597 mp_free(mc->p);
598 mp_free(mc->a);
599 mp_free(mc->b);
600 mp_free(mc->aplus2over4);
601 monty_free(mc->mc);
602 sfree(mc);
603 }
604
ecc_montgomery_point_new_empty(MontgomeryCurve * mc)605 static MontgomeryPoint *ecc_montgomery_point_new_empty(MontgomeryCurve *mc)
606 {
607 MontgomeryPoint *mp = snew(MontgomeryPoint);
608 mp->mc = mc;
609 mp->X = mp->Z = NULL;
610 return mp;
611 }
612
ecc_montgomery_point_new(MontgomeryCurve * mc,mp_int * x)613 MontgomeryPoint *ecc_montgomery_point_new(MontgomeryCurve *mc, mp_int *x)
614 {
615 MontgomeryPoint *mp = ecc_montgomery_point_new_empty(mc);
616 mp->X = monty_import(mc->mc, x);
617 mp->Z = mp_copy(monty_identity(mc->mc));
618 return mp;
619 }
620
ecc_montgomery_point_copy_into(MontgomeryPoint * dest,MontgomeryPoint * src)621 void ecc_montgomery_point_copy_into(
622 MontgomeryPoint *dest, MontgomeryPoint *src)
623 {
624 mp_copy_into(dest->X, src->X);
625 mp_copy_into(dest->Z, src->Z);
626 }
627
ecc_montgomery_point_copy(MontgomeryPoint * orig)628 MontgomeryPoint *ecc_montgomery_point_copy(MontgomeryPoint *orig)
629 {
630 MontgomeryPoint *mp = ecc_montgomery_point_new_empty(orig->mc);
631 mp->X = mp_copy(orig->X);
632 mp->Z = mp_copy(orig->Z);
633 return mp;
634 }
635
ecc_montgomery_point_free(MontgomeryPoint * mp)636 void ecc_montgomery_point_free(MontgomeryPoint *mp)
637 {
638 mp_free(mp->X);
639 mp_free(mp->Z);
640 smemclr(mp, sizeof(*mp));
641 sfree(mp);
642 }
643
ecc_montgomery_cond_overwrite(MontgomeryPoint * dest,MontgomeryPoint * src,unsigned overwrite)644 static void ecc_montgomery_cond_overwrite(
645 MontgomeryPoint *dest, MontgomeryPoint *src, unsigned overwrite)
646 {
647 mp_select_into(dest->X, dest->X, src->X, overwrite);
648 mp_select_into(dest->Z, dest->Z, src->Z, overwrite);
649 }
650
ecc_montgomery_cond_swap(MontgomeryPoint * P,MontgomeryPoint * Q,unsigned swap)651 static void ecc_montgomery_cond_swap(
652 MontgomeryPoint *P, MontgomeryPoint *Q, unsigned swap)
653 {
654 mp_cond_swap(P->X, Q->X, swap);
655 mp_cond_swap(P->Z, Q->Z, swap);
656 }
657
ecc_montgomery_diff_add(MontgomeryPoint * P,MontgomeryPoint * Q,MontgomeryPoint * PminusQ)658 MontgomeryPoint *ecc_montgomery_diff_add(
659 MontgomeryPoint *P, MontgomeryPoint *Q, MontgomeryPoint *PminusQ)
660 {
661 MontgomeryCurve *mc = P->mc;
662 assert(Q->mc == mc);
663 assert(PminusQ->mc == mc);
664
665 /*
666 * Differential addition is achieved using the following formula
667 * that relates the affine x-coordinates of P, Q, P+Q and P-Q:
668 *
669 * x(P+Q) x(P-Q) (x(Q)-x(P))^2 = (x(P)x(Q) - 1)^2
670 *
671 * As with the Weierstrass coordinates, the code below transforms
672 * that affine relation into a projective one to avoid having to
673 * do a division during the main arithmetic.
674 */
675
676 MontgomeryPoint *S = ecc_montgomery_point_new_empty(mc);
677
678 mp_int *Px_m_Pz = monty_sub(mc->mc, P->X, P->Z);
679 mp_int *Px_p_Pz = monty_add(mc->mc, P->X, P->Z);
680 mp_int *Qx_m_Qz = monty_sub(mc->mc, Q->X, Q->Z);
681 mp_int *Qx_p_Qz = monty_add(mc->mc, Q->X, Q->Z);
682 mp_int *PmQp = monty_mul(mc->mc, Px_m_Pz, Qx_p_Qz);
683 mp_int *PpQm = monty_mul(mc->mc, Px_p_Pz, Qx_m_Qz);
684 mp_int *Xpre = monty_add(mc->mc, PmQp, PpQm);
685 mp_int *Zpre = monty_sub(mc->mc, PmQp, PpQm);
686 mp_int *Xpre2 = monty_mul(mc->mc, Xpre, Xpre);
687 mp_int *Zpre2 = monty_mul(mc->mc, Zpre, Zpre);
688 S->X = monty_mul(mc->mc, Xpre2, PminusQ->Z);
689 S->Z = monty_mul(mc->mc, Zpre2, PminusQ->X);
690
691 mp_free(Px_m_Pz);
692 mp_free(Px_p_Pz);
693 mp_free(Qx_m_Qz);
694 mp_free(Qx_p_Qz);
695 mp_free(PmQp);
696 mp_free(PpQm);
697 mp_free(Xpre);
698 mp_free(Zpre);
699 mp_free(Xpre2);
700 mp_free(Zpre2);
701
702 return S;
703 }
704
ecc_montgomery_double(MontgomeryPoint * P)705 MontgomeryPoint *ecc_montgomery_double(MontgomeryPoint *P)
706 {
707 MontgomeryCurve *mc = P->mc;
708 MontgomeryPoint *D = ecc_montgomery_point_new_empty(mc);
709
710 /*
711 * To double a point in affine coordinates, in principle you can
712 * use the same technique as for Weierstrass: differentiate the
713 * curve equation to get the tangent line at the input point, use
714 * that to get an expression for y which you substitute back into
715 * the curve equation, and subtract the known two roots (in this
716 * case both the same) from the x^2 coefficient of the resulting
717 * cubic.
718 *
719 * In this case, we don't have an input y-coordinate, so you have
720 * to do a bit of extra transformation to find a formula that can
721 * work without it. The tangent formula is (3x^2 + 2ax + 1)/(2y),
722 * and when that appears in the final formula it will be squared -
723 * so we can substitute the y^2 in the denominator for the RHS of
724 * the curve equation. Put together, that gives
725 *
726 * x_out = (x+1)^2 (x-1)^2 / 4(x^3+ax^2+x)
727 *
728 * and, as usual, the code below transforms that into projective
729 * form to avoid the division.
730 */
731
732 mp_int *Px_m_Pz = monty_sub(mc->mc, P->X, P->Z);
733 mp_int *Px_p_Pz = monty_add(mc->mc, P->X, P->Z);
734 mp_int *Px_m_Pz_2 = monty_mul(mc->mc, Px_m_Pz, Px_m_Pz);
735 mp_int *Px_p_Pz_2 = monty_mul(mc->mc, Px_p_Pz, Px_p_Pz);
736 D->X = monty_mul(mc->mc, Px_m_Pz_2, Px_p_Pz_2);
737 mp_int *XZ = monty_mul(mc->mc, P->X, P->Z);
738 mp_int *twoXZ = monty_add(mc->mc, XZ, XZ);
739 mp_int *fourXZ = monty_add(mc->mc, twoXZ, twoXZ);
740 mp_int *fourXZ_scaled = monty_mul(mc->mc, fourXZ, mc->aplus2over4);
741 mp_int *Zpre = monty_add(mc->mc, Px_m_Pz_2, fourXZ_scaled);
742 D->Z = monty_mul(mc->mc, fourXZ, Zpre);
743
744 mp_free(Px_m_Pz);
745 mp_free(Px_p_Pz);
746 mp_free(Px_m_Pz_2);
747 mp_free(Px_p_Pz_2);
748 mp_free(XZ);
749 mp_free(twoXZ);
750 mp_free(fourXZ);
751 mp_free(fourXZ_scaled);
752 mp_free(Zpre);
753
754 return D;
755 }
756
ecc_montgomery_normalise(MontgomeryPoint * mp)757 static void ecc_montgomery_normalise(MontgomeryPoint *mp)
758 {
759 MontgomeryCurve *mc = mp->mc;
760 mp_int *zinv = monty_invert(mc->mc, mp->Z);
761 monty_mul_into(mc->mc, mp->X, mp->X, zinv);
762 monty_mul_into(mc->mc, mp->Z, mp->Z, zinv);
763 mp_free(zinv);
764 }
765
ecc_montgomery_multiply(MontgomeryPoint * B,mp_int * n)766 MontgomeryPoint *ecc_montgomery_multiply(MontgomeryPoint *B, mp_int *n)
767 {
768 /*
769 * 'Montgomery ladder' technique, to compute an arbitrary integer
770 * multiple of B under the constraint that you can only add two
771 * unequal points if you also know their difference.
772 *
773 * The setup is that you maintain two curve points one of which is
774 * always the other one plus B. Call them kB and (k+1)B, where k
775 * is some integer that evolves as we go along. We begin by
776 * doubling the input B, to initialise those points to B and 2B,
777 * so that k=1.
778 *
779 * At each stage, we add kB and (k+1)B together - which we can do
780 * under the differential-addition constraint because we know
781 * their difference is always just B - to give us (2k+1)B. Then we
782 * double one of kB or (k+1)B, and depending on which one we
783 * choose, we end up with (2k)B or (2k+2)B. Either way, that
784 * differs by B from the other value we've just computed. So in
785 * each iteration, we do one diff-add and one doubling, plus a
786 * couple of conditional swaps to choose which value we double and
787 * which way round we put the output points, and the effect is to
788 * replace k with either 2k or 2k+1, which we choose based on the
789 * appropriate bit of the desired exponent.
790 *
791 * This routine doesn't assume we know the exact location of the
792 * topmost set bit of the exponent. So to maintain constant time
793 * it does an iteration for every _potential_ bit, starting from
794 * the top downwards; after each iteration in which we haven't
795 * seen a set exponent bit yet, we just overwrite the two points
796 * with B and 2B again,
797 */
798
799 MontgomeryPoint *two_B = ecc_montgomery_double(B);
800 MontgomeryPoint *k_B = ecc_montgomery_point_copy(B);
801 MontgomeryPoint *kplus1_B = ecc_montgomery_point_copy(two_B);
802
803 unsigned not_started_yet = 1;
804 for (size_t bitindex = mp_max_bits(n); bitindex-- > 0 ;) {
805 unsigned nbit = mp_get_bit(n, bitindex);
806
807 MontgomeryPoint *sum = ecc_montgomery_diff_add(k_B, kplus1_B, B);
808 ecc_montgomery_cond_swap(k_B, kplus1_B, nbit);
809 MontgomeryPoint *other = ecc_montgomery_double(k_B);
810 ecc_montgomery_point_free(k_B);
811 ecc_montgomery_point_free(kplus1_B);
812 k_B = other;
813 kplus1_B = sum;
814 ecc_montgomery_cond_swap(k_B, kplus1_B, nbit);
815
816 ecc_montgomery_cond_overwrite(k_B, B, not_started_yet);
817 ecc_montgomery_cond_overwrite(kplus1_B, two_B, not_started_yet);
818 not_started_yet &= ~nbit;
819 }
820
821 ecc_montgomery_point_free(two_B);
822 ecc_montgomery_point_free(kplus1_B);
823 return k_B;
824 }
825
ecc_montgomery_get_affine(MontgomeryPoint * mp,mp_int ** x)826 void ecc_montgomery_get_affine(MontgomeryPoint *mp, mp_int **x)
827 {
828 MontgomeryCurve *mc = mp->mc;
829
830 ecc_montgomery_normalise(mp);
831
832 if (x)
833 *x = monty_export(mc->mc, mp->X);
834 }
835
ecc_montgomery_is_identity(MontgomeryPoint * mp)836 unsigned ecc_montgomery_is_identity(MontgomeryPoint *mp)
837 {
838 return mp_eq_integer(mp->Z, 0);
839 }
840
841 /* ----------------------------------------------------------------------
842 * Twisted Edwards curves.
843 */
844
845 struct EdwardsPoint {
846 /*
847 * We represent an Edwards curve point in 'extended coordinates'.
848 * There's more than one coordinate system going by that name,
849 * unfortunately. These ones have the semantics that X,Y,Z are
850 * ordinary projective coordinates (so x=X/Z and y=Y/Z), but also,
851 * we store the extra value T = xyZ = XY/Z.
852 */
853 mp_int *X, *Y, *Z, *T;
854
855 EdwardsCurve *ec;
856 };
857
858 struct EdwardsCurve {
859 /* Prime modulus of the finite field. */
860 mp_int *p;
861
862 /* Montgomery context for arithmetic mod p. */
863 MontyContext *mc;
864
865 /* Modsqrt context for point decompression. */
866 ModsqrtContext *sc;
867
868 /* Parameters of the curve, in Montgomery-multiplication
869 * transformed form. */
870 mp_int *d, *a;
871 };
872
ecc_edwards_curve(mp_int * p,mp_int * d,mp_int * a,mp_int * nonsquare_mod_p)873 EdwardsCurve *ecc_edwards_curve(mp_int *p, mp_int *d, mp_int *a,
874 mp_int *nonsquare_mod_p)
875 {
876 EdwardsCurve *ec = snew(EdwardsCurve);
877 ec->p = mp_copy(p);
878 ec->mc = monty_new(p);
879 ec->d = monty_import(ec->mc, d);
880 ec->a = monty_import(ec->mc, a);
881
882 if (nonsquare_mod_p)
883 ec->sc = modsqrt_new(p, nonsquare_mod_p);
884 else
885 ec->sc = NULL;
886
887 return ec;
888 }
889
ecc_edwards_curve_free(EdwardsCurve * ec)890 void ecc_edwards_curve_free(EdwardsCurve *ec)
891 {
892 mp_free(ec->p);
893 mp_free(ec->d);
894 mp_free(ec->a);
895 monty_free(ec->mc);
896 if (ec->sc)
897 modsqrt_free(ec->sc);
898 sfree(ec);
899 }
900
ecc_edwards_point_new_empty(EdwardsCurve * ec)901 static EdwardsPoint *ecc_edwards_point_new_empty(EdwardsCurve *ec)
902 {
903 EdwardsPoint *ep = snew(EdwardsPoint);
904 ep->ec = ec;
905 ep->X = ep->Y = ep->Z = ep->T = NULL;
906 return ep;
907 }
908
ecc_edwards_point_new_imported(EdwardsCurve * ec,mp_int * monty_x,mp_int * monty_y)909 static EdwardsPoint *ecc_edwards_point_new_imported(
910 EdwardsCurve *ec, mp_int *monty_x, mp_int *monty_y)
911 {
912 EdwardsPoint *ep = ecc_edwards_point_new_empty(ec);
913 ep->X = monty_x;
914 ep->Y = monty_y;
915 ep->T = monty_mul(ec->mc, ep->X, ep->Y);
916 ep->Z = mp_copy(monty_identity(ec->mc));
917 return ep;
918 }
919
ecc_edwards_point_new(EdwardsCurve * ec,mp_int * x,mp_int * y)920 EdwardsPoint *ecc_edwards_point_new(
921 EdwardsCurve *ec, mp_int *x, mp_int *y)
922 {
923 return ecc_edwards_point_new_imported(
924 ec, monty_import(ec->mc, x), monty_import(ec->mc, y));
925 }
926
ecc_edwards_point_copy_into(EdwardsPoint * dest,EdwardsPoint * src)927 void ecc_edwards_point_copy_into(EdwardsPoint *dest, EdwardsPoint *src)
928 {
929 mp_copy_into(dest->X, src->X);
930 mp_copy_into(dest->Y, src->Y);
931 mp_copy_into(dest->Z, src->Z);
932 mp_copy_into(dest->T, src->T);
933 }
934
ecc_edwards_point_copy(EdwardsPoint * orig)935 EdwardsPoint *ecc_edwards_point_copy(EdwardsPoint *orig)
936 {
937 EdwardsPoint *ep = ecc_edwards_point_new_empty(orig->ec);
938 ep->X = mp_copy(orig->X);
939 ep->Y = mp_copy(orig->Y);
940 ep->Z = mp_copy(orig->Z);
941 ep->T = mp_copy(orig->T);
942 return ep;
943 }
944
ecc_edwards_point_free(EdwardsPoint * ep)945 void ecc_edwards_point_free(EdwardsPoint *ep)
946 {
947 mp_free(ep->X);
948 mp_free(ep->Y);
949 mp_free(ep->Z);
950 mp_free(ep->T);
951 smemclr(ep, sizeof(*ep));
952 sfree(ep);
953 }
954
ecc_edwards_point_new_from_y(EdwardsCurve * ec,mp_int * yorig,unsigned desired_x_parity)955 EdwardsPoint *ecc_edwards_point_new_from_y(
956 EdwardsCurve *ec, mp_int *yorig, unsigned desired_x_parity)
957 {
958 assert(ec->sc);
959
960 /*
961 * The curve equation is ax^2 + y^2 = 1 + dx^2y^2, which
962 * rearranges to x^2(dy^2-a) = y^2-1. So we compute
963 * (y^2-1)/(dy^2-a) and take its square root.
964 */
965 unsigned success;
966
967 mp_int *y = monty_import(ec->mc, yorig);
968 mp_int *y2 = monty_mul(ec->mc, y, y);
969 mp_int *dy2 = monty_mul(ec->mc, ec->d, y2);
970 mp_int *dy2ma = monty_sub(ec->mc, dy2, ec->a);
971 mp_int *y2m1 = monty_sub(ec->mc, y2, monty_identity(ec->mc));
972 mp_int *recip_denominator = monty_invert(ec->mc, dy2ma);
973 mp_int *radicand = monty_mul(ec->mc, y2m1, recip_denominator);
974 mp_int *x = monty_modsqrt(ec->sc, radicand, &success);
975 mp_free(y2);
976 mp_free(dy2);
977 mp_free(dy2ma);
978 mp_free(y2m1);
979 mp_free(recip_denominator);
980 mp_free(radicand);
981
982 if (!success) {
983 /* Failure! x^2 worked out to be a number that has no square
984 * root mod p. In this situation there's no point in trying to
985 * be time-constant, since the protocol sequence is going to
986 * diverge anyway when we complain to whoever gave us this
987 * bogus value. */
988 mp_free(x);
989 mp_free(y);
990 return NULL;
991 }
992
993 /*
994 * Choose whichever of x and p-x has the specified parity (of its
995 * lowest positive residue mod p).
996 */
997 mp_int *tmp = monty_export(ec->mc, x);
998 unsigned flip = (mp_get_bit(tmp, 0) ^ desired_x_parity) & 1;
999 mp_sub_into(tmp, ec->p, x);
1000 mp_select_into(x, x, tmp, flip);
1001 mp_free(tmp);
1002
1003 return ecc_edwards_point_new_imported(ec, x, y);
1004 }
1005
ecc_edwards_cond_overwrite(EdwardsPoint * dest,EdwardsPoint * src,unsigned overwrite)1006 static void ecc_edwards_cond_overwrite(
1007 EdwardsPoint *dest, EdwardsPoint *src, unsigned overwrite)
1008 {
1009 mp_select_into(dest->X, dest->X, src->X, overwrite);
1010 mp_select_into(dest->Y, dest->Y, src->Y, overwrite);
1011 mp_select_into(dest->Z, dest->Z, src->Z, overwrite);
1012 mp_select_into(dest->T, dest->T, src->T, overwrite);
1013 }
1014
ecc_edwards_cond_swap(EdwardsPoint * P,EdwardsPoint * Q,unsigned swap)1015 static void ecc_edwards_cond_swap(
1016 EdwardsPoint *P, EdwardsPoint *Q, unsigned swap)
1017 {
1018 mp_cond_swap(P->X, Q->X, swap);
1019 mp_cond_swap(P->Y, Q->Y, swap);
1020 mp_cond_swap(P->Z, Q->Z, swap);
1021 mp_cond_swap(P->T, Q->T, swap);
1022 }
1023
ecc_edwards_add(EdwardsPoint * P,EdwardsPoint * Q)1024 EdwardsPoint *ecc_edwards_add(EdwardsPoint *P, EdwardsPoint *Q)
1025 {
1026 EdwardsCurve *ec = P->ec;
1027 assert(Q->ec == ec);
1028
1029 EdwardsPoint *S = ecc_edwards_point_new_empty(ec);
1030
1031 /*
1032 * The affine rule for Edwards addition of (x1,y1) and (x2,y2) is
1033 *
1034 * x_out = (x1 y2 + y1 x2) / (1 + d x1 x2 y1 y2)
1035 * y_out = (y1 y2 - a x1 x2) / (1 - d x1 x2 y1 y2)
1036 *
1037 * The formulae below are listed as 'add-2008-hwcd' in
1038 * https://hyperelliptic.org/EFD/g1p/auto-twisted-extended.html
1039 *
1040 * and if you undo the careful optimisation to find out what
1041 * they're actually computing, it comes out to
1042 *
1043 * X_out = (X1 Y2 + Y1 X2) (Z1 Z2 - d T1 T2)
1044 * Y_out = (Y1 Y2 - a X1 X2) (Z1 Z2 + d T1 T2)
1045 * Z_out = (Z1 Z2 - d T1 T2) (Z1 Z2 + d T1 T2)
1046 * T_out = (X1 Y2 + Y1 X2) (Y1 Y2 - a X1 X2)
1047 */
1048 mp_int *PxQx = monty_mul(ec->mc, P->X, Q->X);
1049 mp_int *PyQy = monty_mul(ec->mc, P->Y, Q->Y);
1050 mp_int *PtQt = monty_mul(ec->mc, P->T, Q->T);
1051 mp_int *PzQz = monty_mul(ec->mc, P->Z, Q->Z);
1052 mp_int *Psum = monty_add(ec->mc, P->X, P->Y);
1053 mp_int *Qsum = monty_add(ec->mc, Q->X, Q->Y);
1054 mp_int *aPxQx = monty_mul(ec->mc, ec->a, PxQx);
1055 mp_int *dPtQt = monty_mul(ec->mc, ec->d, PtQt);
1056 mp_int *sumprod = monty_mul(ec->mc, Psum, Qsum);
1057 mp_int *xx_p_yy = monty_add(ec->mc, PxQx, PyQy);
1058 mp_int *E = monty_sub(ec->mc, sumprod, xx_p_yy);
1059 mp_int *F = monty_sub(ec->mc, PzQz, dPtQt);
1060 mp_int *G = monty_add(ec->mc, PzQz, dPtQt);
1061 mp_int *H = monty_sub(ec->mc, PyQy, aPxQx);
1062 S->X = monty_mul(ec->mc, E, F);
1063 S->Z = monty_mul(ec->mc, F, G);
1064 S->Y = monty_mul(ec->mc, G, H);
1065 S->T = monty_mul(ec->mc, H, E);
1066
1067 mp_free(PxQx);
1068 mp_free(PyQy);
1069 mp_free(PtQt);
1070 mp_free(PzQz);
1071 mp_free(Psum);
1072 mp_free(Qsum);
1073 mp_free(aPxQx);
1074 mp_free(dPtQt);
1075 mp_free(sumprod);
1076 mp_free(xx_p_yy);
1077 mp_free(E);
1078 mp_free(F);
1079 mp_free(G);
1080 mp_free(H);
1081
1082 return S;
1083 }
1084
ecc_edwards_normalise(EdwardsPoint * ep)1085 static void ecc_edwards_normalise(EdwardsPoint *ep)
1086 {
1087 EdwardsCurve *ec = ep->ec;
1088 mp_int *zinv = monty_invert(ec->mc, ep->Z);
1089 monty_mul_into(ec->mc, ep->X, ep->X, zinv);
1090 monty_mul_into(ec->mc, ep->Y, ep->Y, zinv);
1091 monty_mul_into(ec->mc, ep->Z, ep->Z, zinv);
1092 mp_free(zinv);
1093 monty_mul_into(ec->mc, ep->T, ep->X, ep->Y);
1094 }
1095
ecc_edwards_multiply(EdwardsPoint * B,mp_int * n)1096 EdwardsPoint *ecc_edwards_multiply(EdwardsPoint *B, mp_int *n)
1097 {
1098 EdwardsPoint *two_B = ecc_edwards_add(B, B);
1099 EdwardsPoint *k_B = ecc_edwards_point_copy(B);
1100 EdwardsPoint *kplus1_B = ecc_edwards_point_copy(two_B);
1101
1102 /*
1103 * Another copy of the same exponentiation routine following the
1104 * pattern of the Montgomery ladder, because it works as well as
1105 * any other technique and this way I didn't have to debug two of
1106 * them.
1107 */
1108
1109 unsigned not_started_yet = 1;
1110 for (size_t bitindex = mp_max_bits(n); bitindex-- > 0 ;) {
1111 unsigned nbit = mp_get_bit(n, bitindex);
1112
1113 EdwardsPoint *sum = ecc_edwards_add(k_B, kplus1_B);
1114 ecc_edwards_cond_swap(k_B, kplus1_B, nbit);
1115 EdwardsPoint *other = ecc_edwards_add(k_B, k_B);
1116 ecc_edwards_point_free(k_B);
1117 ecc_edwards_point_free(kplus1_B);
1118 k_B = other;
1119 kplus1_B = sum;
1120 ecc_edwards_cond_swap(k_B, kplus1_B, nbit);
1121
1122 ecc_edwards_cond_overwrite(k_B, B, not_started_yet);
1123 ecc_edwards_cond_overwrite(kplus1_B, two_B, not_started_yet);
1124 not_started_yet &= ~nbit;
1125 }
1126
1127 ecc_edwards_point_free(two_B);
1128 ecc_edwards_point_free(kplus1_B);
1129 return k_B;
1130 }
1131
1132 /*
1133 * Helper routine to determine whether two values each given as a pair
1134 * of projective coordinates represent the same affine value.
1135 */
projective_eq(MontyContext * mc,mp_int * An,mp_int * Ad,mp_int * Bn,mp_int * Bd)1136 static inline unsigned projective_eq(
1137 MontyContext *mc, mp_int *An, mp_int *Ad,
1138 mp_int *Bn, mp_int *Bd)
1139 {
1140 mp_int *AnBd = monty_mul(mc, An, Bd);
1141 mp_int *BnAd = monty_mul(mc, Bn, Ad);
1142 unsigned toret = mp_cmp_eq(AnBd, BnAd);
1143 mp_free(AnBd);
1144 mp_free(BnAd);
1145 return toret;
1146 }
1147
ecc_edwards_eq(EdwardsPoint * P,EdwardsPoint * Q)1148 unsigned ecc_edwards_eq(EdwardsPoint *P, EdwardsPoint *Q)
1149 {
1150 EdwardsCurve *ec = P->ec;
1151 assert(Q->ec == ec);
1152
1153 return (projective_eq(ec->mc, P->X, P->Z, Q->X, Q->Z) &
1154 projective_eq(ec->mc, P->Y, P->Z, Q->Y, Q->Z));
1155 }
1156
ecc_edwards_get_affine(EdwardsPoint * ep,mp_int ** x,mp_int ** y)1157 void ecc_edwards_get_affine(EdwardsPoint *ep, mp_int **x, mp_int **y)
1158 {
1159 EdwardsCurve *ec = ep->ec;
1160
1161 ecc_edwards_normalise(ep);
1162
1163 if (x)
1164 *x = monty_export(ec->mc, ep->X);
1165 if (y)
1166 *y = monty_export(ec->mc, ep->Y);
1167 }
1168