--- /dev/null
+/*
+ * 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 */
+};
+
+
+
+#ifdef USE_64BIT
+
+/* 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
+
+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);
+ }
+}
projects / zpackage / blobdiff