text: /*
   * Implementing twisted edwards curve
   * 	E : -x^2 + y^2 = 1 - (121665/121666) x^2 y^2
   * over GF(2^255-19) according to [1].
   *
   * This code is public domain.
   *
   * Philipp Lay <philipp.lay@illunis.net>.
   *
   * References:
   * [1] High-speed high-security signatures, 2011/09/26,
   *     Bernstein, Duif, Lange, Schwabe, Yang
   * [2] Twisted Edwards Curves Revisited, 2008,
   *     Hisil, Wong, Carter, Dawson
   */
  
  #include <stdint.h>
  #include <string.h>
  
  #include "bitness.h"
  #include "fld.h"
  #include "sc.h"
  #include "ed.h"
  
  
  /*
   * special pre-computed form of a point on the curve, used
   * for the lookup table and some optimizations.
   */
  struct pced {
  	fld_t		diff;		/* y - x */
  	fld_t		sum;		/* y + x */
  	fld_t		prod;		/* 2*d*t */
  };
  
  
  
  #if (USE_64BIT+0)
  
  /* lookup-table for ed_scale_base - 64bit version */
  static const struct pced ed_lookup[][8] = {
    #include "ed_lookup64.h"
  };
  
  /* base point B of our group in pre-computed form */
  const struct pced pced_B = {
  	{62697248952638, 204681361388450, 631292143396476,
  	 338455783676468, 1213667448819585},
  	{1288382639258501, 245678601348599, 269427782077623,
  	 1462984067271730, 137412439391563},
  	{301289933810280, 1259582250014073, 1422107436869536,
  	 796239922652654, 1953934009299142} };
  
  #else
  
  /* lookup-table for ed_scale_base - 32bit version */
  static const struct pced ed_lookup[][8] = {
    #include "ed_lookup32.h"
  };
  
  /* base point B of our group in pre-computed form */
  const struct pced pced_B = {
  	{54563134, 934261, 64385954, 3049989, 66381436,
  	 9406985, 12720692, 5043384, 19500929, 18085054},
  	{25967493, 19198397, 29566455, 3660896, 54414519,
  	 4014786, 27544626, 21800161, 61029707, 2047604},
  	{58370664, 4489569, 9688441, 18769238, 10184608,
  	 21191052, 29287918, 11864899, 42594502, 29115885} };
  
  #endif /* (USE_64BIT+0)*/
  
  const struct ed ed_zero = { { 0 }, { 1 }, { 0 }, { 1 } };
  const struct pced pced_zero = { { 1 }, { 1 }, { 0 } };
  
  
  /*
   * memselect - helper function to select one of two buffers in a
   * time-constant way.
   */
  static void
  memselect(void *out, const void *pa, const void *pb, size_t max, int flag)
  {
  	uint8_t *a = (uint8_t*)pa;
  	uint8_t *b = (uint8_t*)pb;
  	uint8_t *p = (uint8_t*)out;
  	uint8_t *endp = p+max;
  
  	uint8_t mB = (flag & 1) - 1;
  	uint8_t mA = ~mB;
  
  	while (p < endp)
  		*p++ = (mA & *a++) ^ (mB & *b++);
  }
  
  
  /*
   * ed_import - import a point P on the curve from packet 256bit encoding.
   *
   */
  void
  ed_import(struct ed *P, const uint8_t in[32])
  {
  	fld_t U, V, A, B;
  	uint8_t tmp[32];
  	int flag;
  
  	/* import y */
  	memcpy(tmp, in, 32);
  	tmp[31] &= 0x7f;
  
  	fld_import(P->y, tmp);
  
  
  	/* U <- y^2 - 1,  V <- d*y^2 + 1 */
  	fld_sq(U, P->y);
  	fld_mul(V, con_d, U);
  	U[0]--;
  	V[0]++;
  
  	/* A <- v^2 */
  	fld_sq(A, V);
  	/* B <- v^4 */
  	fld_sq(B, A);
  	/* A <- uv^3 */
  	fld_mul(A, A, U);
  	fld_mul(A, A, V);
  	/* B <- (uv^7)^((q-5)/8) */
  	fld_mul(B, B, A);
  	fld_pow2523(B, B);
  	/* B <- uv^3 * (uv^7)^((q-5)/8) */
  	fld_mul(B, B, A);
  
  	/* A <- v * B^2 */
  	fld_sq(A, B);
  	fld_mul(A, A, V);
  	/* flag <- v*B^2 == u */
  	flag = fld_eq(A, U);
  
  	/* A <- j * B */
  	fld_mul(A, con_j, B);
  
  	memselect(P->x, B, A, sizeof(fld_t), flag);
  	fld_reduce(P->x, P->x);
  	fld_tinyscale(P->x, P->x, 1 - 2*((in[31] >> 7) ^ (P->x[0] & 1)));
  
  	/* compute t and z */
  	fld_mul(P->t, P->x, P->y);
  	fld_set0(P->z, 1);
  }
  
  
  /*
   * ed_export - export point P to packed 256bit format.
   */
  void
  ed_export(uint8_t out[32], const struct ed *P)
  {
  	fld_t x, y, zinv;
  
  	/* divide x and y by z to get affine coordinates */
  	fld_inv(zinv, P->z);
  	fld_mul(x, P->x, zinv);
  	fld_mul(y, P->y, zinv);
  	fld_export(out, y);
  
  	/* lsb decides the sign of x */
  	fld_reduce(x, x);
  	out[31] |= (x[0] & 1) << 7;
  }
  
  
  /*
   * ed_add - add points P and Q
   */
  static void
  ed_add(struct ed *out, const struct ed *P, const struct ed *Q)
  {
  	fld_t a, b, c, d, e, f, g, h, t;
  
  	fld_sub(a, P->y, P->x);
  	fld_sub(t, Q->y, Q->x);
  	fld_mul(a, a, t);
  
  	fld_add(b, P->y, P->x);
  	fld_add(t, Q->y, Q->x);
  	fld_mul(b, b, t);
  
  	fld_mul(c, P->t, Q->t);
  	fld_mul(c, c, con_2d);
  
  	fld_mul(d, P->z, Q->z);
  	fld_scale2(d, d);
  
  	fld_sub(e, b, a);
  	fld_sub(f, d, c);
  	fld_add(g, d, c);
  	fld_add(h, b, a);
  
  	fld_mul(out->x, e, f);
  	fld_mul(out->y, g, h);
  	fld_mul(out->t, e, h);
  	fld_mul(out->z, f, g);
  }
  
  
  /*
   * ed_double - special case of ed_add for P=Q.
   *
   * little bit faster, since 4 multiplications are turned into squares.
   */
  static void
  ed_double(struct ed *out, const struct ed *P)
  {
  	fld_t a, b, c, d, e, f, g, h;
  
  	fld_sub(a, P->y, P->x);
  	fld_sq(a, a);
  
  	fld_add(b, P->y, P->x);
  	fld_sq(b, b);
  
  	fld_sq(c, P->t);
  	fld_mul(c, c, con_2d);
  
  	fld_sq(d, P->z);
  	fld_scale2(d, d);
  
  	fld_sub(e, b, a);
  	fld_sub(f, d, c);
  	fld_add(g, d, c);
  	fld_add(h, b, a);
  
  	fld_mul(out->x, e, f);
  	fld_mul(out->y, g, h);
  	fld_mul(out->t, e, h);
  	fld_mul(out->z, f, g);
  }
  
  /*
   * ed_sub - subtract two points P and Q.
   *
   * alternatively we could negate a point (which is cheap) and use
   * ed_add, but this is a little bit faster.
   */
  static void
  ed_sub(struct ed *out, const struct ed *P, const struct ed *Q)
  {
  	fld_t a, b, c, d, e, f, g, h, t;
  
  	fld_sub(a, P->y, P->x);
  	fld_add(t, Q->y, Q->x);
  	fld_mul(a, a, t);
  
  	fld_add(b, P->y, P->x);
  	fld_sub(t, Q->y, Q->x);
  	fld_mul(b, b, t);
  
  	fld_mul(c, P->t, Q->t);
  	fld_mul(c, c, con_m2d);
  
  	fld_mul(d, P->z, Q->z);
  	fld_scale2(d, d);
  
  	fld_sub(e, b, a);
  	fld_sub(f, d, c);
  	fld_add(g, d, c);
  	fld_add(h, b, a);
  
  	fld_mul(out->x, e, f);
  	fld_mul(out->y, g, h);
  	fld_mul(out->t, e, h);
  	fld_mul(out->z, f, g);
  }
  
  
  /*
   * ed_add_pc - add points P and Q where Q is in precomputed form
   *
   * the precomputed form is used for the lookup table and as an
   * optimization in ed_double_scale.
   */
  static void
  ed_add_pc(struct ed *out, const struct ed *P, const struct pced *Q)
  {
  	fld_t a, b, c, d, e, f, g, h;
  
  	fld_sub(a, P->y, P->x);
  	fld_mul(a, a, Q->diff);
  
  	fld_add(b, P->y, P->x);
  	fld_mul(b, b, Q->sum);
  
  	fld_mul(c, P->t, Q->prod);
  	fld_scale2(d, P->z);
  
  	fld_sub(e, b, a);
  	fld_sub(f, d, c);
  	fld_add(g, d, c);
  	fld_add(h, b, a);
  
  	fld_mul(out->x, e, f);
  	fld_mul(out->y, g, h);
  	fld_mul(out->t, e, h);
  	fld_mul(out->z, f, g);
  }
  
  /*
   * ed_sub_pc - subtract P and Q where Q is in precomputed form.
   */
  static void
  ed_sub_pc(struct ed *out, const struct ed *P, const struct pced *Q)
  {
  	fld_t a, b, c, d, e, f, g, h;
  
  	fld_sub(a, P->y, P->x);
  	fld_mul(a, a, Q->sum);
  
  	fld_add(b, P->y, P->x);
  	fld_mul(b, b, Q->diff);
  
  	fld_mul(c, P->t, Q->prod);
  	fld_neg(c, c);
  
  	fld_scale2(d, P->z);
  
  	fld_sub(e, b, a);
  	fld_sub(f, d, c);
  	fld_add(g, d, c);
  	fld_add(h, b, a);
  
  	fld_mul(out->x, e, f);
  	fld_mul(out->y, g, h);
  	fld_mul(out->t, e, h);
  	fld_mul(out->z, f, g);
  }
  
  
  /*
   * scale16 - helper function for ed_scale_base, returns x * 16^pow * base
   *
   * assumes:
   *  -8 <= x <= 7
   *   0 <= pow <= 62
   *   2 | pow
   */
  static void
  scale16(struct pced *out, int pow, int x)
  {
  	struct pced R = { { 0 }, { 0 }, { 0 } };
  	limb_t mA, mB, mask;
  	int neg, sgnx, absx;
  	int i, k;
  
  	neg = (x >> 3) & 1;
  	sgnx = 1 - 2*neg;
  	absx = sgnx * x;
  	pow >>= 1;
  
  	/* handle abs(x) == 0 */
  	mask = absx | (absx >> 2);
  	mask |= mask >> 1;
  	mask = (mask & 1) - 1;
  	for (i = 0; i < FLD_LIMB_NUM; i++) {
  		R.diff[i] ^= pced_zero.diff[i] & mask;
  		R.sum[i] ^= pced_zero.sum[i] & mask;
  		R.prod[i] ^= pced_zero.prod[i] & mask;
  	}
  
  
  	/* go through our table and look for abs(x) */
  	for (k = 0; k < 8; k++) {
  		absx--;
  		mask = absx | (absx >> 2);
  		mask |= mask >> 1;
  		mask = (mask & 1) - 1;
  		for (i = 0; i < FLD_LIMB_NUM; i++) {
  			R.diff[i] ^= ed_lookup[pow][k].diff[i] & mask;
  			R.sum[i] ^= ed_lookup[pow][k].sum[i] & mask;
  			R.prod[i] ^= ed_lookup[pow][k].prod[i] & mask;
  		}
  	}
  
  	/* conditionally negate R and write to out */
  	mA = neg-1;
  	mB = ~mA;
  	for (i = 0; i < FLD_LIMB_NUM; i++) {
  		out->diff[i] = (mA & R.diff[i]) ^ (mB & R.sum[i]);
  		out->sum[i]  = (mB & R.diff[i]) ^ (mA & R.sum[i]);
  		out->prod[i] = sgnx * R.prod[i];
  	}
  }
  
  
  /*
   * ed_scale_base - calculates x * base
   */
  void
  ed_scale_base(struct ed *out, const sc_t x)
  {
  	struct ed R0, R1;
  	struct pced P;
  	sc_t tmp;
  	uint8_t pack[32];
  	int r[64];
  	int i;
  
  	/* s <- x + 8 * (16^64 - 1) / 15 */
  	sc_add(tmp, x, con_off);
  	sc_export(pack, tmp);
  	for (i = 0; i < 32; i++) {
  		r[2*i] = (pack[i] & 0x0f) - 8;
  		r[2*i+1] = (pack[i] >> 4) - 8;
  	}
  
  	/*
  	 * R0 <- r0*B + r2*16^2*B + ... + r62*16^62*B
  	 * R1 <- r1*B + r3*16^2*B + ... + r63*16^62*B
  	 */
  	memcpy(&R0, &ed_zero, sizeof(struct ed));
  	memcpy(&R1, &ed_zero, sizeof(struct ed));
  	for (i = 0; i < 63; i += 2) {
  		scale16(&P, i, r[i]);
  		ed_add_pc(&R0, &R0, &P);
  
  		scale16(&P, i, r[i+1]);
  		ed_add_pc(&R1, &R1, &P);
  	}
  
  	/* R1 <- 16 * R1 */
  	for (i = 0; i < 4; i++)
  		ed_add(&R1, &R1, &R1);
  
  	/* out <- R0 + R1 */
  	ed_add(out, &R0, &R1);
  }
  
  
  /*
   * helper function to speed up ed_double_scale
   */
  static void
  ed_precompute(struct pced *R, const struct ed *P)
  {
  	fld_sub(R->diff, P->y, P->x);
  	fld_add(R->sum, P->y, P->x);
  	fld_mul(R->prod, P->t, con_2d);
  }
  
  
  /*
   * ed_dual_scale - calculates R = x*base + y*Q.  (vartime)
   *
   * Note: This algorithms does NOT run in constant time! Please use this
   * only for public information like in ed25519_verify().
   *
   * assumes:
   *   Q is affine, ie has z = 1
   *   x and y must be reduced
   */
  void
  ed_dual_scale(struct ed *R,
  	      const sc_t x,
  	      const sc_t y, const struct ed *Q)
  {
  	struct ed QpB, QmB;
  	struct pced pcQ;
  
  	int ux[SC_BITS+1], uy[SC_BITS+1];
  	int n, i;
  
  	memcpy(R, &ed_zero, sizeof(struct ed));
  
  	/* calculate joint sparse form of x and y */
  	n = sc_jsf(ux, uy, x, y);
  	if (n == -1)
  		return;
  
  	/* precompute Q, Q+B and Q-B */
  	ed_add_pc(&QpB, Q, &pced_B);
  	ed_sub_pc(&QmB, Q, &pced_B);
  	ed_precompute(&pcQ, Q);
  
  	/* now we calculate R = x * base_point + y * Q using fast shamir method */
  	for (i = n; ; i--) {
  		if (ux[i] == 1) {
  			if (uy[i] == 1)
  				ed_add(R, R, &QpB);
  			else if (uy[i] == -1)
  				ed_sub(R, R, &QmB);
  			else
  				ed_add_pc(R, R, &pced_B);
  
  		} else if (ux[i] == -1) {
  			if (uy[i] == 1)
  				ed_add(R, R, &QmB);
  			else if (uy[i] == -1)
  				ed_sub(R, R, &QpB);
  			else
  				ed_sub_pc(R, R, &pced_B);
  
  		} else {
  			if (uy[i] == 1)
  				ed_add_pc(R, R, &pcQ);
  			else if (uy[i] == -1)
  				ed_sub_pc(R, R, &pcQ);
  		}
  
  		if (i == 0) break;
  
  		ed_double(R, R);
  	}
  }
  