/src/relic/src/eb/relic_eb_mul.c

Source
/*
 * RELIC is an Efficient LIbrary for Cryptography
 * Copyright (c) 2009 RELIC Authors
 *
 * This file is part of RELIC. RELIC is legal property of its developers,
 * whose names are not listed here. Please refer to the COPYRIGHT file
 * for contact information.
 *
 * RELIC is free software; you can redistribute it and/or modify it under the
 * terms of the version 2.1 (or later) of the GNU Lesser General Public License
 * as published by the Free Software Foundation; or version 2.0 of the Apache
 * License as published by the Apache Software Foundation. See the LICENSE files
 * for more details.
 *
 * RELIC is distributed in the hope that it will be useful, but WITHOUT ANY
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 * A PARTICULAR PURPOSE. See the LICENSE files for more details.
 *
 * You should have received a copy of the GNU Lesser General Public or the
 * Apache License along with RELIC. If not, see <https://www.gnu.org/licenses/>
 * or <https://www.apache.org/licenses/>.
 */

/**
 * @file
 *
 * Implementation of point multiplication on binary elliptic curves.
 *
 * @ingroup eb
 */

#include "relic_core.h"
#include "relic_fb_low.h"
#include "relic_ep_mul_tmpl.h"

/*============================================================================*/
/* Private definitions                                                        */
/*============================================================================*/

/**
 * Combines the accumulators of a right-to-left w-NAF method, computing
 * R = sum_i (2i + 1) * T[i] as T[0] + 2 * sum_{i >= 1} S[i], where
 * S[i] = sum_{j >= i} T[j]. The accumulators are overwritten.
 *
 * @param[out] r        - the result.
 * @param[in,out] t       - the accumulators.
 */
static void eb_mul_odd_sum(eb_t r, eb_t *t) {
  int i;

  eb_set_infty(r);
  for (i = (1 << (RLC_WIDTH - 2)) - 1; i > 0; i--) {
    eb_add(t[i - 1], t[i - 1], t[i]);
    eb_add(r, r, t[i]);
  }
  eb_dbl(r, r);
  eb_add(r, r, t[0]);
}

/**
 * Multiplies a point on a Koblitz curve by a small element
 * alpha = b + g * tau of Z[tau], using the tau-adic NAF of alpha.
 *
 * @param[out] r        - the result.
 * @param[in] p         - the point to multiply.
 * @param[in] b         - the integer part of alpha.
 * @param[in] g         - the coefficient of tau in alpha.
 * @param[in] u         - the parameter mu of the curve.
 */
static void eb_mul_alpha(eb_t r, const eb_t p, int b, int g, int u) {
  int8_t d[2 * RLC_WIDTH + 4];
  int c, l = 0;

  /* Compute the tau-adic NAF of alpha, least significant digit first. */
  while (b != 0 || g != 0) {
    d[l] = 0;
    if (b & 1) {
      d[l] = 2 - ((b - 2 * g) & 3);
      b -= d[l];
    }
    /* Divide by tau, using tau^2 = mu * tau - 2. */
    c = b / 2;
    b = g + u * c;
    g = -c;
    l++;
  }

  /* Evaluate the expansion with Horner's rule. */
  if (d[l - 1] > 0) {
    eb_copy(r, p);
  } else {
    eb_neg(r, p);
  }
  for (l -= 2; l >= 0; l--) {
    eb_frb(r, r);
    if (d[l] > 0) {
      eb_add(r, r, p);
    }
    if (d[l] < 0) {
      eb_sub(r, r, p);
    }
  }
}

/*============================================================================*/
/* Public definitions                                                         */
/*============================================================================*/

TMPL_EC_MUL_BASIC(eb);

void eb_mul_lodah(eb_t r, const eb_t p, const bn_t k) {
  dv_t x1, z1, x2, z2, r1, r2, r3, r4, r5;
  const dig_t *b;
  bn_t t, n;
  size_t bits;
  eb_t q;

  if (bn_is_zero(k) || eb_is_infty(p)) {
    eb_set_infty(r);
    return;
  }

  eb_null(q);
  bn_null_all(n, t);
  dv_null_all(x1, z1, x2, z2, r1, r2, r3, r4, r5);

  RLC_TRY {
    /* The ladder needs the input point in affine coordinates. */
    eb_new(q);
    eb_norm(q, p);
    bn_new_all(n, t);
    dv_new_all(x1, z1, x2, z2, r1, r2, r3, r4, r5);

    fb_sqr(z2, q->x);
    fb_sqr(x2, z2);
    dv_zero(r5, 2 * RLC_FB_DIGS);

    b = eb_curve_get_b();
    eb_curve_get_ord(n);
    bits = bn_bits(n);

    bn_mod(t, k, n);
    bn_add(t, t, n);
    bn_add(n, t, n);
    dv_swap_sec(t->dp, n->dp, RLC_MAX(t->used, n->used),
      bn_get_bit(t, bits) == 0);
    t->used = RLC_SEL(t->used, n->used, bn_get_bit(t, bits) == 0);

    switch (eb_curve_opt_b()) {
      case RLC_ZERO:
        break;
      case RLC_ONE:
        fb_add_dig(x2, x2, (dig_t)1);
        break;
      case RLC_TINY:
        fb_add_dig(x2, x2, b[0]);
        break;
      default:
        fb_addn_low(x2, x2, b);
        break;
    }

    /* Blind both points indebendently. */
    fb_rand(z1);
    fb_mul(x1, z1, q->x);
    fb_rand(r1);
    fb_mul(z2, z2, r1);
    fb_mul(x2, x2, r1);

    for (int i = bits - 1; i >= 0; i--) {
      int j = bn_get_bit(t, i);
      fb_mul(r1, x1, z2);
      fb_mul(r2, x2, z1);
      fb_add(r3, r1, r2);
      fb_muln_low(r4, r1, r2);
      dv_swap_sec(x1, x2, RLC_FB_DIGS, j ^ 1);
      dv_swap_sec(z1, z2, RLC_FB_DIGS, j ^ 1);
      fb_sqr(z1, r3);
      fb_muln_low(r1, z1, q->x);
      fb_addd_low(x1, r1, r4, 2 * RLC_FB_DIGS);
      fb_rdcn_low(x1, x1);
      fb_sqr(r1, z2);
      fb_sqr(r2, x2);
      fb_mul(z2, r1, r2);
      switch (eb_curve_opt_b()) {
        case RLC_ZERO:
          fb_sqr(x2, r2);
          break;
        case RLC_ONE:
          fb_add(r1, r1, r2);
          fb_sqr(x2, r1);
          break;
        case RLC_TINY:
          fb_sqr(r1, r1);
          fb_sqrl_low(x2, r2);
          fb_mul1_low(r5, r1, b[0]);
          fb_addd_low(x2, x2, r5, RLC_FB_DIGS + 1);
          fb_rdcn_low(x2, x2);
          break;
        default:
          fb_sqr(r1, r1);
          fb_sqrl_low(x2, r2);
          fb_muln_low(r5, r1, b);
          fb_addd_low(x2, x2, r5, 2 * RLC_FB_DIGS);
          fb_rdcn_low(x2, x2);
          break;
      }
      dv_swap_sec(x1, x2, RLC_FB_DIGS, j ^ 1);
      dv_swap_sec(z1, z2, RLC_FB_DIGS, j ^ 1);
    }

    if (fb_is_zero(z1)) {
      /* The point q is at infinity. */
      eb_set_infty(r);
    } else {
      if (fb_is_zero(z2)) {
        fb_copy(r->x, q->x);
        fb_add(r->y, q->x, q->y);
        fb_set_dig(r->z, 1);
      } else {
        /* r3 = z1 * z2. */
        fb_mul(r3, z1, z2);
        /* z1 = (x1 + x * z1). */
        fb_mul(z1, z1, q->x);
        fb_add(z1, z1, x1);
        /* z2 = x * z2. */
        fb_mul(z2, z2, q->x);
        /* x1 = x1 * z2. */
        fb_mul(x1, x1, z2);
        /* z2 = (x2 + x * z2)(x1 + x * z1). */
        fb_add(z2, z2, x2);
        fb_mul(z2, z2, z1);

        /* r4 = (x^2 + y) * z1 * z2 + (x2 + x * z2)(x1 + x * z1). */
        fb_sqr(r4, q->x);
        fb_add(r4, r4, q->y);
        fb_mul(r4, r4, r3);
        fb_add(r4, r4, z2);

        /* r3 = (z1 * z2 * x)^{-1}. */
        fb_mul(r3, r3, q->x);
        fb_inv(r3, r3);
        /* r4 = (x^2 + y) * z1 * z2 + (x2 + x * z2)(x1 + x * z1) * r3. */
        fb_mul(r4, r4, r3);
        /* x2 = x1 * x * z2 * (z1 * z2 * x)^{-1} = x1/z1. */
        fb_mul(x2, x1, r3);
        /* z2 = x + x1/z1. */
        fb_add(z2, x2, q->x);

        /* z2 = z2 * r4 + y. */
        fb_mul(z2, z2, r4);
        fb_add(z2, z2, q->y);

        fb_copy(r->x, x2);
        fb_copy(r->y, z2);
        fb_set_dig(r->z, 1);
      }
    }

    r->coord = BASIC;
  }
  RLC_CATCH_ANY {
    RLC_THROW(ERR_CAUGHT);
  }
  RLC_FINALLY {
    eb_free(q);
    bn_free_all(n, t);
    dv_free_all(x1, z1, x2, z2, r1, r2, r3, r4, r5);
  }
}

void eb_mul_lwnaf(eb_t r, const eb_t p, const bn_t k) {
  int i, d, kbltz = eb_curve_is_kbltz();
  int8_t u, naf[RLC_FB_BITS + 8];
  eb_t t[1 << (RLC_WIDTH - 2)];
  bn_t n, _k;
  size_t l;

  if (bn_is_zero(k) || eb_is_infty(p)) {
    eb_set_infty(r);
    return;
  }

  bn_null_all(n, _k);
  for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
    eb_null(t[i]);
  }

  RLC_TRY {
    bn_new_all(n, _k);
    for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
      eb_new(t[i]);
    }

    /* Reduce the scalar to bound the length of its recoding. */
    eb_curve_get_ord(n);
    bn_mod(_k, k, n);

    /* Compute the precomputation table. */
    eb_tab(t, p, RLC_WIDTH);

    /* Compute the w-TNAF or w-NAF representation of k. */
    l = sizeof(naf);
    if (kbltz) {
      u = (eb_curve_opt_a() == RLC_ZERO ? -1 : 1);
      bn_rec_tnaf(naf, &l, _k, u, RLC_FB_BITS, RLC_WIDTH);
    } else {
      bn_rec_naf(naf, &l, _k, RLC_WIDTH);
    }

    eb_set_infty(r);
    for (i = l - 1; i >= 0; i--) {
      /* The Frobenius map replaces doubling on Koblitz curves. */
      if (kbltz) {
        eb_frb(r, r);
      } else {
        eb_dbl(r, r);
      }

      d = naf[i];
      if (d > 0) {
        eb_add(r, r, t[d / 2]);
      }
      if (d < 0) {
        eb_sub(r, r, t[-d / 2]);
      }
    }
    /* Convert r to affine coordinates. */
    eb_norm(r, r);
  }
  RLC_CATCH_ANY {
    RLC_THROW(ERR_CAUGHT);
  }
  RLC_FINALLY {
    bn_free_all(n, _k);
    for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
      eb_free(t[i]);
    }
  }
}

void eb_mul_rwnaf(eb_t r, const eb_t p, const bn_t k) {
  size_t i;
  int d, kbltz = eb_curve_is_kbltz();
  int8_t u, naf[RLC_FB_BITS + 8];
  int8_t beta[1 << (RLC_WIDTH - 2)], gama[1 << (RLC_WIDTH - 2)];
  uint8_t t_w;
  eb_t t[1 << (RLC_WIDTH - 2)];
  bn_t n, _k;
  size_t l;

  if (bn_is_zero(k) || eb_is_infty(p)) {
    eb_set_infty(r);
    return;
  }

  bn_null_all(n, _k);
  for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
    eb_null(t[i]);
  }

  RLC_TRY {
    bn_new_all(n, _k);
    /* Prepare the accumulator table. */
    for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
      eb_new(t[i]);
      eb_set_infty(t[i]);
    }

    /* Reduce the scalar to bound the length of its recoding. */
    eb_curve_get_ord(n);
    bn_mod(_k, k, n);

    /* Compute the w-TNAF or w-NAF representation of k. */
    u = (eb_curve_opt_a() == RLC_ZERO ? -1 : 1);
    l = sizeof(naf);
    if (kbltz) {
      bn_rec_tnaf(naf, &l, _k, u, RLC_FB_BITS, RLC_WIDTH);
    } else {
      bn_rec_naf(naf, &l, _k, RLC_WIDTH);
    }

    /* The method needs the input point in affine coordinates. */
    eb_norm(r, p);
    for (i = 0; i < l; i++) {
      d = naf[i];
      if (d > 0) {
        eb_add(t[d / 2], t[d / 2], r);
      }
      if (d < 0) {
        eb_sub(t[-d / 2], t[-d / 2], r);
      }

      /* The Frobenius map replaces doubling on Koblitz curves, and
       * only squares the coordinates of an affine point. */
      if (kbltz) {
        fb_sqr(r->x, r->x);
        fb_sqr(r->y, r->y);
      } else {
        eb_dbl(r, r);
      }
    }

    if (kbltz) {
      /* Multiply each accumulator by its representative alpha. */
      bn_rec_tnaf_get(&t_w, beta, gama, u, RLC_WIDTH);
      for (i = 1; i < (1 << (RLC_WIDTH - 2)); i++) {
        eb_copy(r, t[i]);
        eb_mul_alpha(t[i], r, beta[i], gama[i], u);
      }
      /* Add the accumulators. */
      eb_copy(r, t[0]);
      for (i = 1; i < (1 << (RLC_WIDTH - 2)); i++) {
        eb_add(r, r, t[i]);
      }
    } else {
      eb_mul_odd_sum(r, t);
    }
    /* Convert r to affine coordinates. */
    eb_norm(r, r);
  }
  RLC_CATCH_ANY {
    RLC_THROW(ERR_CAUGHT);
  }
  RLC_FINALLY {
    bn_free_all(n, _k);
    for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
      eb_free(t[i]);
    }
  }
}

void eb_mul_halve(eb_t r, const eb_t p, const bn_t k) {
  int i, j, cof;
  dig_t trc;
  int8_t naf[RLC_FB_BITS + 1], *_k;
  eb_t a, q, s, t[1 << (RLC_WIDTH - 2)];
  bn_t n, m;
  fb_t u, v, w, z;
  size_t l;

  if (bn_is_zero(k) || eb_is_infty(p)) {
    eb_set_infty(r);
    return;
  }

  bn_null_all(m, n);
  eb_null_all(a, q, s);
  for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
    eb_null(t[i]);
  }
  fb_null_all(u, v, w, z);

  RLC_TRY {
    /* The method needs the input point in affine coordinates. */
    eb_new(a);
    eb_norm(a, p);
    bn_new_all(n, m);
    eb_new_all(q, s);
    fb_new_all(u, v, w, z);

    /* Prepare the precomputation table. */
    for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
      eb_new(t[i]);
      eb_set_infty(t[i]);
    }

    /* Convert k to alternate representation k' = (2^{t-1}k mod n). */
    eb_curve_get_ord(n);
    bn_lsh(m, k, bn_bits(n) - 1);
    bn_mod(m, m, n);

    /* Compute the w-NAF representation of k'. */
    l = sizeof(naf);
    bn_rec_naf(naf, &l, m, RLC_WIDTH);

    if (naf[bn_bits(n)] == 1) {
      eb_dbl(t[0], a);
    }
    l = bn_bits(n);
    _k = naf + l - 1;

    eb_copy(q, a);
    eb_curve_get_cof(n);

    /* Test if curve has a cofactor bigger than 2. */
    if (bn_cmp_dig(n, 2) == RLC_GT) {
      cof = 1;
    } else {
      cof = 0;
    }

    trc = fb_trc(eb_curve_get_a());

    if (cof) {
      /* Curves with cofactor > 2, u = sqrt(a), v = Solve(u). */
      fb_srt(u, eb_curve_get_a());
      fb_slv(v, u);

      bn_rand(n, RLC_POS, l);

      for (i = l - 1; i >= 0; i--, _k--) {
        j = *_k;
        if (j > 0) {
          eb_norm(s, q);
          eb_add(t[j / 2], t[j / 2], s);
        }
        if (j < 0) {
          eb_norm(s, q);
          eb_sub(t[-j / 2], t[-j / 2], s);
        }

        /* T = 1/2(Q). */
        eb_hlv(s, q);

        /* If Tr(x_T) != Tr(a). */
        if (fb_trc(s->x) != 0) {
          /* z = l_t, w = sqrt(l_Q), l_T = l_T + sqrt(l_Q) + v. */
          fb_copy(z, s->y);
          fb_srt(w, q->y);
          fb_add(s->y, s->y, w);
          fb_add(s->y, s->y, v);
          /* z = (z + x_Q + v + sqrt(a)). */
          fb_add(z, z, q->x);
          fb_add(z, z, v);
          fb_add(z, z, u);
          /* w = sqrt(w + x_Q + l_Q + sqrt(a)). */
          fb_add(w, w, q->x);
          fb_add(w, w, q->y);
          fb_add(w, w, u);
          /* x_T = sqrt(w * z), . */
          fb_mul(w, w, z);
          fb_srt(s->x, w);
          fb_set_dig(s->z, 1);
          s->coord = HALVE;
        }
        eb_copy(q, s);
      }
    } else {
      for (i = l - 1; i >= 0; i--, _k--) {
        j = *_k;
        if (j > 0) {
          eb_norm(q, q);
          eb_add(t[j / 2], t[j / 2], q);
        }
        if (j < 0) {
          eb_norm(q, q);
          eb_sub(t[-j / 2], t[-j / 2], q);
        }
        eb_hlv(q, q);
      }
    }

    eb_mul_odd_sum(r, t);
    eb_norm(r, r);

    /* We may need to fix an error of a 2-torsion point if the curve has a
     * 4-cofactor. */
    if (cof) {
      eb_hlv(s, r);
      if (fb_trc(s->x) != trc) {
        fb_zero(s->x);
        fb_srt(s->y, eb_curve_get_b());
        fb_set_dig(s->z, 1);
        eb_add(r, r, s);
        eb_norm(r, r);
      }
    }
  }
  RLC_CATCH_ANY {
    RLC_THROW(ERR_CAUGHT);
  }
  RLC_FINALLY {
    /* Free the precomputation table. */
    for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
      eb_free(t[i]);
    }
    bn_free_all(n, m);
    eb_free_all(a, q, s);
    fb_free_all(u, v, w, z);
  }
}

#ifdef EB_PRECO
TMPL_EC_MUL_GEN(eb, preco);
#else
TMPL_EC_MUL_GEN(eb, plain);
#endif

TMPL_EC_MUL_DIG(eb);

Coverage Report

Created: 2026-10-17 06:43

Line	Count	Source
1		/*
2		* RELIC is an Efficient LIbrary for Cryptography
3		* Copyright (c) 2009 RELIC Authors
4		*
5		* This file is part of RELIC. RELIC is legal property of its developers,
6		* whose names are not listed here. Please refer to the COPYRIGHT file
7		* for contact information.
8		*
9		* RELIC is free software; you can redistribute it and/or modify it under the
10		* terms of the version 2.1 (or later) of the GNU Lesser General Public License
11		* as published by the Free Software Foundation; or version 2.0 of the Apache
12		* License as published by the Apache Software Foundation. See the LICENSE files
13		* for more details.
14		*
15		* RELIC is distributed in the hope that it will be useful, but WITHOUT ANY
16		* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
17		* A PARTICULAR PURPOSE. See the LICENSE files for more details.
18		*
19		* You should have received a copy of the GNU Lesser General Public or the
20		* Apache License along with RELIC. If not, see <https://www.gnu.org/licenses/>
21		* or <https://www.apache.org/licenses/>.
22		*/
23
24		/**
25		* @file
26		*
27		* Implementation of point multiplication on binary elliptic curves.
28		*
29		* @ingroup eb
30		*/
31
32		#include "relic_core.h"
33		#include "relic_fb_low.h"
34		#include "relic_ep_mul_tmpl.h"
35
36		/============================================================================/
37		/* Private definitions */
38		/============================================================================/
39
40		/**
41		* Combines the accumulators of a right-to-left w-NAF method, computing
42		* R = sum_i (2i + 1) * T[i] as T[0] + 2 * sum_{i >= 1} S[i], where
43		* S[i] = sum_{j >= i} T[j]. The accumulators are overwritten.
44		*
45		* @param[out] r - the result.
46		* @param[in,out] t - the accumulators.
47		*/
48	0	static void eb_mul_odd_sum(eb_t r, eb_t *t) {
49	0	int i;
50
51	0	eb_set_infty(r);
52	0	for (i = (1 << (RLC_WIDTH - 2)) - 1; i > 0; i--) {
53	0	eb_add(t[i - 1], t[i - 1], t[i]);
54	0	eb_add(r, r, t[i]);
55	0	}
56	0	eb_dbl(r, r);
57	0	eb_add(r, r, t[0]);
58	0	}
59
60		/**
61		* Multiplies a point on a Koblitz curve by a small element
62		* alpha = b + g * tau of Z[tau], using the tau-adic NAF of alpha.
63		*
64		* @param[out] r - the result.
65		* @param[in] p - the point to multiply.
66		* @param[in] b - the integer part of alpha.
67		* @param[in] g - the coefficient of tau in alpha.
68		* @param[in] u - the parameter mu of the curve.
69		*/
70	0	static void eb_mul_alpha(eb_t r, const eb_t p, int b, int g, int u) {
71	0	int8_t d[2 * RLC_WIDTH + 4];
72	0	int c, l = 0;
73
74		/* Compute the tau-adic NAF of alpha, least significant digit first. */
75	0	while (b != 0 \|\| g != 0) {
76	0	d[l] = 0;
77	0	if (b & 1) {
78	0	d[l] = 2 - ((b - 2 * g) & 3);
79	0	b -= d[l];
80	0	}
81		/* Divide by tau, using tau^2 = mu * tau - 2. */
82	0	c = b / 2;
83	0	b = g + u * c;
84	0	g = -c;
85	0	l++;
86	0	}
87
88		/* Evaluate the expansion with Horner's rule. */
89	0	if (d[l - 1] > 0) {
90	0	eb_copy(r, p);
91	0	} else {
92	0	eb_neg(r, p);
93	0	}
94	0	for (l -= 2; l >= 0; l--) {
95	0	eb_frb(r, r);
96	0	if (d[l] > 0) {
97	0	eb_add(r, r, p);
98	0	}
99	0	if (d[l] < 0) {
100	0	eb_sub(r, r, p);
101	0	}
102	0	}
103	0	}
104
105		/============================================================================/
106		/* Public definitions */
107		/============================================================================/
108
109	0	TMPL_EC_MUL_BASIC(eb);
110
111	0	void eb_mul_lodah(eb_t r, const eb_t p, const bn_t k) {
112	0	dv_t x1, z1, x2, z2, r1, r2, r3, r4, r5;
113	0	const dig_t *b;
114	0	bn_t t, n;
115	0	size_t bits;
116	0	eb_t q;
117
118	0	if (bn_is_zero(k) \|\| eb_is_infty(p)) {
119	0	eb_set_infty(r);
120	0	return;
121	0	}
122
123	0	eb_null(q);
124	0	bn_null_all(n, t);
125	0	dv_null_all(x1, z1, x2, z2, r1, r2, r3, r4, r5);
126
127	0	RLC_TRY {
128		/* The ladder needs the input point in affine coordinates. */
129	0	eb_new(q);
130	0	eb_norm(q, p);
131	0	bn_new_all(n, t);
132	0	dv_new_all(x1, z1, x2, z2, r1, r2, r3, r4, r5);
133
134	0	fb_sqr(z2, q->x);
135	0	fb_sqr(x2, z2);
136	0	dv_zero(r5, 2 * RLC_FB_DIGS);
137
138	0	b = eb_curve_get_b();
139	0	eb_curve_get_ord(n);
140	0	bits = bn_bits(n);
141
142	0	bn_mod(t, k, n);
143	0	bn_add(t, t, n);
144	0	bn_add(n, t, n);
145	0	dv_swap_sec(t->dp, n->dp, RLC_MAX(t->used, n->used),
146	0	bn_get_bit(t, bits) == 0);
147	0	t->used = RLC_SEL(t->used, n->used, bn_get_bit(t, bits) == 0);
148
149	0	switch (eb_curve_opt_b()) {
150	0	case RLC_ZERO:
151	0	break;
152	0	case RLC_ONE:
153	0	fb_add_dig(x2, x2, (dig_t)1);
154	0	break;
155	0	case RLC_TINY:
156	0	fb_add_dig(x2, x2, b[0]);
157	0	break;
158	0	default:
159	0	fb_addn_low(x2, x2, b);
160	0	break;
161	0	}
162
163		/* Blind both points indebendently. */
164	0	fb_rand(z1);
165	0	fb_mul(x1, z1, q->x);
166	0	fb_rand(r1);
167	0	fb_mul(z2, z2, r1);
168	0	fb_mul(x2, x2, r1);
169
170	0	for (int i = bits - 1; i >= 0; i--) {
171	0	int j = bn_get_bit(t, i);
172	0	fb_mul(r1, x1, z2);
173	0	fb_mul(r2, x2, z1);
174	0	fb_add(r3, r1, r2);
175	0	fb_muln_low(r4, r1, r2);
176	0	dv_swap_sec(x1, x2, RLC_FB_DIGS, j ^ 1);
177	0	dv_swap_sec(z1, z2, RLC_FB_DIGS, j ^ 1);
178	0	fb_sqr(z1, r3);
179	0	fb_muln_low(r1, z1, q->x);
180	0	fb_addd_low(x1, r1, r4, 2 * RLC_FB_DIGS);
181	0	fb_rdcn_low(x1, x1);
182	0	fb_sqr(r1, z2);
183	0	fb_sqr(r2, x2);
184	0	fb_mul(z2, r1, r2);
185	0	switch (eb_curve_opt_b()) {
186	0	case RLC_ZERO:
187	0	fb_sqr(x2, r2);
188	0	break;
189	0	case RLC_ONE:
190	0	fb_add(r1, r1, r2);
191	0	fb_sqr(x2, r1);
192	0	break;
193	0	case RLC_TINY:
194	0	fb_sqr(r1, r1);
195	0	fb_sqrl_low(x2, r2);
196	0	fb_mul1_low(r5, r1, b[0]);
197	0	fb_addd_low(x2, x2, r5, RLC_FB_DIGS + 1);
198	0	fb_rdcn_low(x2, x2);
199	0	break;
200	0	default:
201	0	fb_sqr(r1, r1);
202	0	fb_sqrl_low(x2, r2);
203	0	fb_muln_low(r5, r1, b);
204	0	fb_addd_low(x2, x2, r5, 2 * RLC_FB_DIGS);
205	0	fb_rdcn_low(x2, x2);
206	0	break;
207	0	}
208	0	dv_swap_sec(x1, x2, RLC_FB_DIGS, j ^ 1);
209	0	dv_swap_sec(z1, z2, RLC_FB_DIGS, j ^ 1);
210	0	}
211
212	0	if (fb_is_zero(z1)) {
213		/* The point q is at infinity. */
214	0	eb_set_infty(r);
215	0	} else {
216	0	if (fb_is_zero(z2)) {
217	0	fb_copy(r->x, q->x);
218	0	fb_add(r->y, q->x, q->y);
219	0	fb_set_dig(r->z, 1);
220	0	} else {
221		/* r3 = z1 * z2. */
222	0	fb_mul(r3, z1, z2);
223		/* z1 = (x1 + x * z1). */
224	0	fb_mul(z1, z1, q->x);
225	0	fb_add(z1, z1, x1);
226		/* z2 = x * z2. */
227	0	fb_mul(z2, z2, q->x);
228		/* x1 = x1 * z2. */
229	0	fb_mul(x1, x1, z2);
230		/* z2 = (x2 + x * z2)(x1 + x * z1). */
231	0	fb_add(z2, z2, x2);
232	0	fb_mul(z2, z2, z1);
233
234		/* r4 = (x^2 + y) * z1 * z2 + (x2 + x * z2)(x1 + x * z1). */
235	0	fb_sqr(r4, q->x);
236	0	fb_add(r4, r4, q->y);
237	0	fb_mul(r4, r4, r3);
238	0	fb_add(r4, r4, z2);
239
240		/* r3 = (z1 * z2 * x)^{-1}. */
241	0	fb_mul(r3, r3, q->x);
242	0	fb_inv(r3, r3);
243		/* r4 = (x^2 + y) * z1 * z2 + (x2 + x * z2)(x1 + x * z1) * r3. */
244	0	fb_mul(r4, r4, r3);
245		/* x2 = x1 * x * z2 * (z1 * z2 * x)^{-1} = x1/z1. */
246	0	fb_mul(x2, x1, r3);
247		/* z2 = x + x1/z1. */
248	0	fb_add(z2, x2, q->x);
249
250		/* z2 = z2 * r4 + y. */
251	0	fb_mul(z2, z2, r4);
252	0	fb_add(z2, z2, q->y);
253
254	0	fb_copy(r->x, x2);
255	0	fb_copy(r->y, z2);
256	0	fb_set_dig(r->z, 1);
257	0	}
258	0	}
259
260	0	r->coord = BASIC;
261	0	}
262	0	RLC_CATCH_ANY {
263	0	RLC_THROW(ERR_CAUGHT);
264	0	}
265	0	RLC_FINALLY {
266	0	eb_free(q);
267	0	bn_free_all(n, t);
268	0	dv_free_all(x1, z1, x2, z2, r1, r2, r3, r4, r5);
269	0	}
270	0	}
271
272	0	void eb_mul_lwnaf(eb_t r, const eb_t p, const bn_t k) {
273	0	int i, d, kbltz = eb_curve_is_kbltz();
274	0	int8_t u, naf[RLC_FB_BITS + 8];
275	0	eb_t t[1 << (RLC_WIDTH - 2)];
276	0	bn_t n, _k;
277	0	size_t l;
278
279	0	if (bn_is_zero(k) \|\| eb_is_infty(p)) {
280	0	eb_set_infty(r);
281	0	return;
282	0	}
283
284	0	bn_null_all(n, _k);
285	0	for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
286	0	eb_null(t[i]);
287	0	}
288
289	0	RLC_TRY {
290	0	bn_new_all(n, _k);
291	0	for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
292	0	eb_new(t[i]);
293	0	}
294
295		/* Reduce the scalar to bound the length of its recoding. */
296	0	eb_curve_get_ord(n);
297	0	bn_mod(_k, k, n);
298
299		/* Compute the precomputation table. */
300	0	eb_tab(t, p, RLC_WIDTH);
301
302		/* Compute the w-TNAF or w-NAF representation of k. */
303	0	l = sizeof(naf);
304	0	if (kbltz) {
305	0	u = (eb_curve_opt_a() == RLC_ZERO ? -1 : 1);
306	0	bn_rec_tnaf(naf, &l, _k, u, RLC_FB_BITS, RLC_WIDTH);
307	0	} else {
308	0	bn_rec_naf(naf, &l, _k, RLC_WIDTH);
309	0	}
310
311	0	eb_set_infty(r);
312	0	for (i = l - 1; i >= 0; i--) {
313		/* The Frobenius map replaces doubling on Koblitz curves. */
314	0	if (kbltz) {
315	0	eb_frb(r, r);
316	0	} else {
317	0	eb_dbl(r, r);
318	0	}
319
320	0	d = naf[i];
321	0	if (d > 0) {
322	0	eb_add(r, r, t[d / 2]);
323	0	}
324	0	if (d < 0) {
325	0	eb_sub(r, r, t[-d / 2]);
326	0	}
327	0	}
328		/* Convert r to affine coordinates. */
329	0	eb_norm(r, r);
330	0	}
331	0	RLC_CATCH_ANY {
332	0	RLC_THROW(ERR_CAUGHT);
333	0	}
334	0	RLC_FINALLY {
335	0	bn_free_all(n, _k);
336	0	for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
337	0	eb_free(t[i]);
338	0	}
339	0	}
340	0	}
341
342	0	void eb_mul_rwnaf(eb_t r, const eb_t p, const bn_t k) {
343	0	size_t i;
344	0	int d, kbltz = eb_curve_is_kbltz();
345	0	int8_t u, naf[RLC_FB_BITS + 8];
346	0	int8_t beta[1 << (RLC_WIDTH - 2)], gama[1 << (RLC_WIDTH - 2)];
347	0	uint8_t t_w;
348	0	eb_t t[1 << (RLC_WIDTH - 2)];
349	0	bn_t n, _k;
350	0	size_t l;
351
352	0	if (bn_is_zero(k) \|\| eb_is_infty(p)) {
353	0	eb_set_infty(r);
354	0	return;
355	0	}
356
357	0	bn_null_all(n, _k);
358	0	for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
359	0	eb_null(t[i]);
360	0	}
361
362	0	RLC_TRY {
363	0	bn_new_all(n, _k);
364		/* Prepare the accumulator table. */
365	0	for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
366	0	eb_new(t[i]);
367	0	eb_set_infty(t[i]);
368	0	}
369
370		/* Reduce the scalar to bound the length of its recoding. */
371	0	eb_curve_get_ord(n);
372	0	bn_mod(_k, k, n);
373
374		/* Compute the w-TNAF or w-NAF representation of k. */
375	0	u = (eb_curve_opt_a() == RLC_ZERO ? -1 : 1);
376	0	l = sizeof(naf);
377	0	if (kbltz) {
378	0	bn_rec_tnaf(naf, &l, _k, u, RLC_FB_BITS, RLC_WIDTH);
379	0	} else {
380	0	bn_rec_naf(naf, &l, _k, RLC_WIDTH);
381	0	}
382
383		/* The method needs the input point in affine coordinates. */
384	0	eb_norm(r, p);
385	0	for (i = 0; i < l; i++) {
386	0	d = naf[i];
387	0	if (d > 0) {
388	0	eb_add(t[d / 2], t[d / 2], r);
389	0	}
390	0	if (d < 0) {
391	0	eb_sub(t[-d / 2], t[-d / 2], r);
392	0	}
393
394		/* The Frobenius map replaces doubling on Koblitz curves, and
395		* only squares the coordinates of an affine point. */
396	0	if (kbltz) {
397	0	fb_sqr(r->x, r->x);
398	0	fb_sqr(r->y, r->y);
399	0	} else {
400	0	eb_dbl(r, r);
401	0	}
402	0	}
403
404	0	if (kbltz) {
405		/* Multiply each accumulator by its representative alpha. */
406	0	bn_rec_tnaf_get(&t_w, beta, gama, u, RLC_WIDTH);
407	0	for (i = 1; i < (1 << (RLC_WIDTH - 2)); i++) {
408	0	eb_copy(r, t[i]);
409	0	eb_mul_alpha(t[i], r, beta[i], gama[i], u);
410	0	}
411		/* Add the accumulators. */
412	0	eb_copy(r, t[0]);
413	0	for (i = 1; i < (1 << (RLC_WIDTH - 2)); i++) {
414	0	eb_add(r, r, t[i]);
415	0	}
416	0	} else {
417	0	eb_mul_odd_sum(r, t);
418	0	}
419		/* Convert r to affine coordinates. */
420	0	eb_norm(r, r);
421	0	}
422	0	RLC_CATCH_ANY {
423	0	RLC_THROW(ERR_CAUGHT);
424	0	}
425	0	RLC_FINALLY {
426	0	bn_free_all(n, _k);
427	0	for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
428	0	eb_free(t[i]);
429	0	}
430	0	}
431	0	}
432
433	0	void eb_mul_halve(eb_t r, const eb_t p, const bn_t k) {
434	0	int i, j, cof;
435	0	dig_t trc;
436	0	int8_t naf[RLC_FB_BITS + 1], *_k;
437	0	eb_t a, q, s, t[1 << (RLC_WIDTH - 2)];
438	0	bn_t n, m;
439	0	fb_t u, v, w, z;
440	0	size_t l;
441
442	0	if (bn_is_zero(k) \|\| eb_is_infty(p)) {
443	0	eb_set_infty(r);
444	0	return;
445	0	}
446
447	0	bn_null_all(m, n);
448	0	eb_null_all(a, q, s);
449	0	for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
450	0	eb_null(t[i]);
451	0	}
452	0	fb_null_all(u, v, w, z);
453
454	0	RLC_TRY {
455		/* The method needs the input point in affine coordinates. */
456	0	eb_new(a);
457	0	eb_norm(a, p);
458	0	bn_new_all(n, m);
459	0	eb_new_all(q, s);
460	0	fb_new_all(u, v, w, z);
461
462		/* Prepare the precomputation table. */
463	0	for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
464	0	eb_new(t[i]);
465	0	eb_set_infty(t[i]);
466	0	}
467
468		/* Convert k to alternate representation k' = (2^{t-1}k mod n). */
469	0	eb_curve_get_ord(n);
470	0	bn_lsh(m, k, bn_bits(n) - 1);
471	0	bn_mod(m, m, n);
472
473		/* Compute the w-NAF representation of k'. */
474	0	l = sizeof(naf);
475	0	bn_rec_naf(naf, &l, m, RLC_WIDTH);
476
477	0	if (naf[bn_bits(n)] == 1) {
478	0	eb_dbl(t[0], a);
479	0	}
480	0	l = bn_bits(n);
481	0	_k = naf + l - 1;
482
483	0	eb_copy(q, a);
484	0	eb_curve_get_cof(n);
485
486		/* Test if curve has a cofactor bigger than 2. */
487	0	if (bn_cmp_dig(n, 2) == RLC_GT) {
488	0	cof = 1;
489	0	} else {
490	0	cof = 0;
491	0	}
492
493	0	trc = fb_trc(eb_curve_get_a());
494
495	0	if (cof) {
496		/* Curves with cofactor > 2, u = sqrt(a), v = Solve(u). */
497	0	fb_srt(u, eb_curve_get_a());
498	0	fb_slv(v, u);
499
500	0	bn_rand(n, RLC_POS, l);
501
502	0	for (i = l - 1; i >= 0; i--, _k--) {
503	0	j = *_k;
504	0	if (j > 0) {
505	0	eb_norm(s, q);
506	0	eb_add(t[j / 2], t[j / 2], s);
507	0	}
508	0	if (j < 0) {
509	0	eb_norm(s, q);
510	0	eb_sub(t[-j / 2], t[-j / 2], s);
511	0	}
512
513		/* T = 1/2(Q). */
514	0	eb_hlv(s, q);
515
516		/* If Tr(x_T) != Tr(a). */
517	0	if (fb_trc(s->x) != 0) {
518		/* z = l_t, w = sqrt(l_Q), l_T = l_T + sqrt(l_Q) + v. */
519	0	fb_copy(z, s->y);
520	0	fb_srt(w, q->y);
521	0	fb_add(s->y, s->y, w);
522	0	fb_add(s->y, s->y, v);
523		/* z = (z + x_Q + v + sqrt(a)). */
524	0	fb_add(z, z, q->x);
525	0	fb_add(z, z, v);
526	0	fb_add(z, z, u);
527		/* w = sqrt(w + x_Q + l_Q + sqrt(a)). */
528	0	fb_add(w, w, q->x);
529	0	fb_add(w, w, q->y);
530	0	fb_add(w, w, u);
531		/* x_T = sqrt(w * z), . */
532	0	fb_mul(w, w, z);
533	0	fb_srt(s->x, w);
534	0	fb_set_dig(s->z, 1);
535	0	s->coord = HALVE;
536	0	}
537	0	eb_copy(q, s);
538	0	}
539	0	} else {
540	0	for (i = l - 1; i >= 0; i--, _k--) {
541	0	j = *_k;
542	0	if (j > 0) {
543	0	eb_norm(q, q);
544	0	eb_add(t[j / 2], t[j / 2], q);
545	0	}
546	0	if (j < 0) {
547	0	eb_norm(q, q);
548	0	eb_sub(t[-j / 2], t[-j / 2], q);
549	0	}
550	0	eb_hlv(q, q);
551	0	}
552	0	}
553
554	0	eb_mul_odd_sum(r, t);
555	0	eb_norm(r, r);
556
557		/* We may need to fix an error of a 2-torsion point if the curve has a
558		* 4-cofactor. */
559	0	if (cof) {
560	0	eb_hlv(s, r);
561	0	if (fb_trc(s->x) != trc) {
562	0	fb_zero(s->x);
563	0	fb_srt(s->y, eb_curve_get_b());
564	0	fb_set_dig(s->z, 1);
565	0	eb_add(r, r, s);
566	0	eb_norm(r, r);
567	0	}
568	0	}
569	0	}
570	0	RLC_CATCH_ANY {
571	0	RLC_THROW(ERR_CAUGHT);
572	0	}
573	0	RLC_FINALLY {
574		/* Free the precomputation table. */
575	0	for (i = 0; i < (1 << (RLC_WIDTH - 2)); i++) {
576	0	eb_free(t[i]);
577	0	}
578	0	bn_free_all(n, m);
579	0	eb_free_all(a, q, s);
580	0	fb_free_all(u, v, w, z);
581	0	}
582	0	}
583
584		#ifdef EB_PRECO
585		TMPL_EC_MUL_GEN(eb, preco);
586		#else
587		TMPL_EC_MUL_GEN(eb, plain);
588		#endif
589
590		TMPL_EC_MUL_DIG(eb);