ecc: Fix ec_mulm_25519.
[libgcrypt.git] / mpi / ec.c
index 9e007cd..b0eed97 100644 (file)
--- a/mpi/ec.c
+++ b/mpi/ec.c
@@ -35,7 +35,7 @@
 #define point_free(a)  _gcry_mpi_point_free_parts ((a))
 
 
-/* Print a point using the log fucntions.  If CTX is not NULL affine
+/* Print a point using the log functions.  If CTX is not NULL affine
    coordinates will be printed.  */
 void
 _gcry_mpi_point_log (const char *name, mpi_point_t point, mpi_ec_t ctx)
@@ -139,6 +139,60 @@ point_set (mpi_point_t d, mpi_point_t s)
 }
 
 
+/* Return a copy of POINT. */
+gcry_mpi_point_t
+_gcry_mpi_point_copy (gcry_mpi_point_t point)
+{
+  mpi_point_t newpoint;
+
+  newpoint = _gcry_mpi_point_new (0);
+  if (point)
+    point_set (newpoint, point);
+
+  return newpoint;
+}
+
+
+static void
+point_resize (mpi_point_t p, mpi_ec_t ctx)
+{
+  size_t nlimbs;
+
+  if (ctx->model == MPI_EC_MONTGOMERY)
+    {
+      nlimbs = ctx->p->nlimbs;
+
+      mpi_resize (p->x, nlimbs);
+      mpi_resize (p->z, nlimbs);
+      p->x->nlimbs = nlimbs;
+      p->z->nlimbs = nlimbs;
+    }
+  else
+    {
+      /*
+       * For now, we allocate enough limbs for our EC computation of ec_*.
+       * Once we will improve ec_* to be constant size (and constant
+       * time), NLIMBS can be ctx->p->nlimbs.
+       */
+      nlimbs = 2*ctx->p->nlimbs+1;
+      mpi_resize (p->x, nlimbs);
+      mpi_resize (p->y, nlimbs);
+      mpi_resize (p->z, nlimbs);
+    }
+}
+
+
+static void
+point_swap_cond (mpi_point_t d, mpi_point_t s, unsigned long swap,
+                 mpi_ec_t ctx)
+{
+  mpi_swap_cond (d->x, s->x, swap);
+  if (ctx->model != MPI_EC_MONTGOMERY)
+    mpi_swap_cond (d->y, s->y, swap);
+  mpi_swap_cond (d->z, s->z, swap);
+}
+
+
 /* Set the projective coordinates from POINT into X, Y, and Z.  If a
    coordinate is not required, X, Y, or Z may be passed as NULL.  */
 void
@@ -247,8 +301,9 @@ ec_addm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
 static void
 ec_subm (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ec)
 {
-  (void)ec;
   mpi_sub (w, u, v);
+  while (w->sign)
+    mpi_add (w, w, ec->p);
   /*ec_mod (w, ec);*/
 }
 
@@ -308,8 +363,166 @@ ec_invm (gcry_mpi_t x, gcry_mpi_t a, mpi_ec_t ctx)
       log_mpidump ("  p", ctx->p);
     }
 }
+\f
+static void
+mpih_set_cond (mpi_ptr_t wp, mpi_ptr_t up, mpi_size_t usize, unsigned long set)
+{
+  mpi_size_t i;
+  mpi_limb_t mask = ((mpi_limb_t)0) - set;
+  mpi_limb_t x;
+
+  for (i = 0; i < usize; i++)
+    {
+      x = mask & (wp[i] ^ up[i]);
+      wp[i] = wp[i] ^ x;
+    }
+}
+
+/* Routines for 2^255 - 19.  */
+
+static void
+ec_mod_25519 (gcry_mpi_t w, mpi_ec_t ec)
+{
+  _gcry_mpi_mod (w, w, ec->p);
+}
+
+#define LIMB_SIZE_25519 ((256+BITS_PER_MPI_LIMB-1)/BITS_PER_MPI_LIMB)
 
+static void
+ec_addm_25519 (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
+{
+  mpi_ptr_t wp, up, vp;
+  mpi_size_t wsize = LIMB_SIZE_25519;
+  mpi_limb_t n[LIMB_SIZE_25519];
+  mpi_limb_t borrow;
+
+  if (w->alloced != wsize || u->alloced != wsize || v->alloced != wsize)
+    log_bug ("addm_25519: different sizes\n");
+
+  memset (n, 0, sizeof n);
+  up = u->d;
+  vp = v->d;
+  wp = w->d;
+
+  _gcry_mpih_add_n (wp, up, vp, wsize);
+  borrow = _gcry_mpih_sub_n (wp, wp, ctx->p->d, wsize);
+  mpih_set_cond (n, ctx->p->d, wsize, (borrow != 0UL));
+  _gcry_mpih_add_n (wp, wp, n, wsize);
+  wp[LIMB_SIZE_25519-1] &= ~(1UL << (255 % BITS_PER_MPI_LIMB));
+}
+
+static void
+ec_subm_25519 (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
+{
+  mpi_ptr_t wp, up, vp;
+  mpi_size_t wsize = LIMB_SIZE_25519;
+  mpi_limb_t n[LIMB_SIZE_25519];
+  mpi_limb_t borrow;
+
+  if (w->alloced != wsize || u->alloced != wsize || v->alloced != wsize)
+    log_bug ("subm_25519: different sizes\n");
+
+  memset (n, 0, sizeof n);
+  up = u->d;
+  vp = v->d;
+  wp = w->d;
+
+  borrow = _gcry_mpih_sub_n (wp, up, vp, wsize);
+  mpih_set_cond (n, ctx->p->d, wsize, (borrow != 0UL));
+  _gcry_mpih_add_n (wp, wp, n, wsize);
+  wp[LIMB_SIZE_25519-1] &= ~(1UL << (255 % BITS_PER_MPI_LIMB));
+}
+
+static void
+ec_mulm_25519 (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx)
+{
+  mpi_ptr_t wp, up, vp;
+  mpi_size_t wsize = LIMB_SIZE_25519;
+  mpi_limb_t n[LIMB_SIZE_25519*2];
+  mpi_limb_t m[LIMB_SIZE_25519+1];
+  mpi_limb_t cy;
+  int msb;
 
+  (void)ctx;
+  if (w->alloced != wsize || u->alloced != wsize || v->alloced != wsize)
+    log_bug ("mulm_25519: different sizes\n");
+
+  up = u->d;
+  vp = v->d;
+  wp = w->d;
+
+  _gcry_mpih_mul_n (n, up, vp, wsize);
+  memcpy (wp, n, wsize * BYTES_PER_MPI_LIMB);
+  wp[LIMB_SIZE_25519-1] &= ~(1UL << (255 % BITS_PER_MPI_LIMB));
+
+  memcpy (m, n+LIMB_SIZE_25519-1, (wsize+1) * BYTES_PER_MPI_LIMB);
+  _gcry_mpih_rshift (m, m, LIMB_SIZE_25519+1, (255 % BITS_PER_MPI_LIMB));
+
+  memcpy (n, m, wsize * BYTES_PER_MPI_LIMB);
+  cy = _gcry_mpih_lshift (m, m, LIMB_SIZE_25519, 4);
+  m[LIMB_SIZE_25519] = cy;
+  cy = _gcry_mpih_add_n (m, m, n, wsize);
+  m[LIMB_SIZE_25519] += cy;
+  cy = _gcry_mpih_add_n (m, m, n, wsize);
+  m[LIMB_SIZE_25519] += cy;
+  cy = _gcry_mpih_add_n (m, m, n, wsize);
+  m[LIMB_SIZE_25519] += cy;
+
+  cy = _gcry_mpih_add_n (wp, wp, m, wsize);
+  m[LIMB_SIZE_25519] += cy;
+
+  memset (m, 0, wsize * BYTES_PER_MPI_LIMB);
+  m[0] = m[LIMB_SIZE_25519] * 2 * 19;
+  cy = _gcry_mpih_add_n (wp, wp, m, wsize);
+
+  msb = (wp[LIMB_SIZE_25519-1] >> (255 % BITS_PER_MPI_LIMB));
+  m[0] = (cy * 2 + msb) * 19;
+  _gcry_mpih_add_n (wp, wp, m, wsize);
+  wp[LIMB_SIZE_25519-1] &= ~(1UL << (255 % BITS_PER_MPI_LIMB));
+
+  m[0] = 0;
+  cy = _gcry_mpih_sub_n (wp, wp, ctx->p->d, wsize);
+  mpih_set_cond (m, ctx->p->d, wsize, (cy != 0UL));
+  _gcry_mpih_add_n (wp, wp, m, wsize);
+}
+
+static void
+ec_mul2_25519 (gcry_mpi_t w, gcry_mpi_t u, mpi_ec_t ctx)
+{
+  ec_addm_25519 (w, u, u, ctx);
+}
+
+static void
+ec_pow2_25519 (gcry_mpi_t w, const gcry_mpi_t b, mpi_ec_t ctx)
+{
+  ec_mulm_25519 (w, b, b, ctx);
+}
+
+struct field_table {
+  const char *p;
+
+  /* computation routines for the field.  */
+  void (* mod) (gcry_mpi_t w, mpi_ec_t ctx);
+  void (* addm) (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx);
+  void (* subm) (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx);
+  void (* mulm) (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v, mpi_ec_t ctx);
+  void (* mul2) (gcry_mpi_t w, gcry_mpi_t u, mpi_ec_t ctx);
+  void (* pow2) (gcry_mpi_t w, const gcry_mpi_t b, mpi_ec_t ctx);
+};
+
+static const struct field_table field_table[] = {
+  {
+    "0x7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFED",
+    ec_mod_25519,
+    ec_addm_25519,
+    ec_subm_25519,
+    ec_mulm_25519,
+    ec_mul2_25519,
+    ec_pow2_25519
+  },
+  { NULL, NULL, NULL, NULL, NULL, NULL, NULL },
+};
+\f
 /* Force recomputation of all helper variables.  */
 void
 _gcry_mpi_ec_get_reset (mpi_ec_t ec)
@@ -353,6 +566,29 @@ ec_get_two_inv_p (mpi_ec_t ec)
 }
 
 
+static const char *curve25519_bad_points[] = {
+  "0x0000000000000000000000000000000000000000000000000000000000000000",
+  "0x0000000000000000000000000000000000000000000000000000000000000001",
+  "0x00b8495f16056286fdb1329ceb8d09da6ac49ff1fae35616aeb8413b7c7aebe0",
+  "0x57119fd0dd4e22d8868e1c58c45c44045bef839c55b1d0b1248c50a3bc959c5f",
+  "0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffec",
+  "0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed",
+  "0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffee",
+  NULL
+};
+
+static gcry_mpi_t
+scanval (const char *string)
+{
+  gpg_err_code_t rc;
+  gcry_mpi_t val;
+
+  rc = _gcry_mpi_scan (&val, GCRYMPI_FMT_HEX, string, 0, NULL);
+  if (rc)
+    log_fatal ("scanning ECC parameter failed: %s\n", gpg_strerror (rc));
+  return val;
+}
+
 
 /* This function initialized a context for elliptic curve based on the
    field GF(p).  P is the prime specifying this field, A is the first
@@ -391,9 +627,51 @@ ec_p_init (mpi_ec_t ctx, enum gcry_mpi_ec_models model,
 
   _gcry_mpi_ec_get_reset (ctx);
 
-  /* Allocate scratch variables.  */
-  for (i=0; i< DIM(ctx->t.scratch); i++)
-    ctx->t.scratch[i] = mpi_alloc_like (ctx->p);
+  if (model == MPI_EC_MONTGOMERY)
+    {
+      for (i=0; i< DIM(ctx->t.scratch) && curve25519_bad_points[i]; i++)
+        ctx->t.scratch[i] = scanval (curve25519_bad_points[i]);
+    }
+  else
+    {
+      /* Allocate scratch variables.  */
+      for (i=0; i< DIM(ctx->t.scratch); i++)
+        ctx->t.scratch[i] = mpi_alloc_like (ctx->p);
+    }
+
+  ctx->mod = ec_mod;
+  ctx->addm = ec_addm;
+  ctx->subm = ec_subm;
+  ctx->mulm = ec_mulm;
+  ctx->mul2 = ec_mul2;
+  ctx->pow2 = ec_pow2;
+
+  for (i=0; field_table[i].p; i++)
+    {
+      gcry_mpi_t f_p;
+      gpg_err_code_t rc;
+
+      rc = _gcry_mpi_scan (&f_p, GCRYMPI_FMT_HEX, field_table[i].p, 0, NULL);
+      if (rc)
+        log_fatal ("scanning ECC parameter failed: %s\n", gpg_strerror (rc));
+
+      if (!mpi_cmp (p, f_p))
+        {
+          ctx->mod = field_table[i].mod;
+          ctx->addm = field_table[i].addm;
+          ctx->subm = field_table[i].subm;
+          ctx->mulm = field_table[i].mulm;
+          ctx->mul2 = field_table[i].mul2;
+          ctx->pow2 = field_table[i].pow2;
+          _gcry_mpi_release (f_p);
+
+          mpi_resize (ctx->a, ctx->p->nlimbs);
+          ctx->a->nlimbs = ctx->p->nlimbs;
+          break;
+        }
+
+      _gcry_mpi_release (f_p);
+    }
 
   /* Prepare for fast reduction.  */
   /* FIXME: need a test for NIST values.  However it does not gain us
@@ -429,6 +707,7 @@ ec_deinit (void *opaque)
   mpi_free (ctx->b);
   _gcry_mpi_point_release (ctx->G);
   mpi_free (ctx->n);
+  mpi_free (ctx->h);
 
   /* The key.  */
   _gcry_mpi_point_release (ctx->Q);
@@ -495,7 +774,7 @@ _gcry_mpi_ec_p_new (gcry_ctx_t *r_ctx,
   mpi_ec_t ec;
 
   *r_ctx = NULL;
-  if (!p || !a || !mpi_cmp_ui (a, 0))
+  if (!p || !a)
     return GPG_ERR_EINVAL;
 
   ctx = _gcry_ctx_alloc (CONTEXT_TYPE_EC, sizeof *ec, ec_deinit);
@@ -560,6 +839,27 @@ _gcry_mpi_ec_set_point (const char *name, gcry_mpi_point_t newvalue,
 }
 
 
+/* Given an encoded point in the MPI VALUE and a context EC, decode
+ * the point according to the context and store it in RESULT.  On
+ * error an error code is return but RESULT might have been changed.
+ * If no context is given the function tries to decode VALUE by
+ * assuming a 0x04 prefixed uncompressed encoding.  */
+gpg_err_code_t
+_gcry_mpi_ec_decode_point (mpi_point_t result, gcry_mpi_t value, mpi_ec_t ec)
+{
+  gcry_err_code_t rc;
+
+  if (ec && ec->dialect == ECC_DIALECT_ED25519)
+    rc = _gcry_ecc_eddsa_decodepoint (value, ec, result, NULL, NULL);
+  else if (ec && ec->model == MPI_EC_MONTGOMERY)
+    rc = _gcry_ecc_mont_decodepoint (value, ec, result);
+  else
+    rc = _gcry_ecc_os2ec (result, value);
+
+  return rc;
+}
+
+
 /* Compute the affine coordinates from the projective coordinates in
    POINT.  Set them into X and Y.  If one coordinate is not required,
    X or Y may be passed as NULL.  CTX is the usual context. Returns: 0
@@ -600,10 +900,17 @@ _gcry_mpi_ec_get_affine (gcry_mpi_t x, gcry_mpi_t y, mpi_point_t point,
 
     case MPI_EC_MONTGOMERY:
       {
-        log_fatal ("%s: %s not yet supported\n",
-                   "_gcry_mpi_ec_get_affine", "Montgomery");
+        if (x)
+          mpi_set (x, point->x);
+
+        if (y)
+          {
+            log_fatal ("%s: Getting Y-coordinate on %s is not supported\n",
+                       "_gcry_mpi_ec_get_affine", "Montgomery");
+            return -1;
+          }
       }
-      return -1;
+      return 0;
 
     case MPI_EC_EDWARDS:
       {
@@ -754,10 +1061,7 @@ dup_point_edwards (mpi_point_t result, mpi_point_t point, mpi_ec_t ctx)
 
   /* E = aC */
   if (ctx->dialect == ECC_DIALECT_ED25519)
-    {
-      mpi_set (E, C);
-      _gcry_mpi_neg (E, E);
-    }
+    mpi_sub (E, ctx->p, C);
   else
     ec_mulm (E, ctx->a, C, ctx);
 
@@ -1035,11 +1339,7 @@ add_points_edwards (mpi_point_t result,
   /* Y_3 = A · G · (D - aC) */
   if (ctx->dialect == ECC_DIALECT_ED25519)
     {
-      /* Using ec_addm (Y3, D, C, ctx) is possible but a litte bit
-         slower because a subm does currently skip the mod step.  */
-      mpi_set (Y3, C);
-      _gcry_mpi_neg (Y3, Y3);
-      ec_subm (Y3, D, Y3, ctx);
+      ec_addm (Y3, D, C, ctx);
     }
   else
     {
@@ -1073,6 +1373,35 @@ add_points_edwards (mpi_point_t result,
 }
 
 
+/* Compute a step of Montgomery Ladder (only use X and Z in the point).
+   Inputs:  P1, P2, and x-coordinate of DIF = P1 - P1.
+   Outputs: PRD = 2 * P1 and  SUM = P1 + P2. */
+static void
+montgomery_ladder (mpi_point_t prd, mpi_point_t sum,
+                   mpi_point_t p1, mpi_point_t p2, gcry_mpi_t dif_x,
+                   mpi_ec_t ctx)
+{
+  ctx->addm (sum->x, p2->x, p2->z, ctx);
+  ctx->subm (p2->z, p2->x, p2->z, ctx);
+  ctx->addm (prd->x, p1->x, p1->z, ctx);
+  ctx->subm (p1->z, p1->x, p1->z, ctx);
+  ctx->mulm (p2->x, p1->z, sum->x, ctx);
+  ctx->mulm (p2->z, prd->x, p2->z, ctx);
+  ctx->pow2 (p1->x, prd->x, ctx);
+  ctx->pow2 (p1->z, p1->z, ctx);
+  ctx->addm (sum->x, p2->x, p2->z, ctx);
+  ctx->subm (p2->z, p2->x, p2->z, ctx);
+  ctx->mulm (prd->x, p1->x, p1->z, ctx);
+  ctx->subm (p1->z, p1->x, p1->z, ctx);
+  ctx->pow2 (sum->x, sum->x, ctx);
+  ctx->pow2 (sum->z, p2->z, ctx);
+  ctx->mulm (prd->z, p1->z, ctx->a, ctx); /* CTX->A: (a-2)/4 */
+  ctx->mulm (sum->z, sum->z, dif_x, ctx);
+  ctx->addm (prd->z, p1->x, prd->z, ctx);
+  ctx->mulm (prd->z, prd->z, p1->z, ctx);
+}
+
+
 /* RESULT = P1 + P2 */
 void
 _gcry_mpi_ec_add_points (mpi_point_t result,
@@ -1094,6 +1423,71 @@ _gcry_mpi_ec_add_points (mpi_point_t result,
 }
 
 
+/* RESULT = P1 - P2  (Weierstrass version).*/
+static void
+sub_points_weierstrass (mpi_point_t result,
+                        mpi_point_t p1, mpi_point_t p2,
+                        mpi_ec_t ctx)
+{
+  (void)result;
+  (void)p1;
+  (void)p2;
+  (void)ctx;
+  log_fatal ("%s: %s not yet supported\n",
+             "_gcry_mpi_ec_sub_points", "Weierstrass");
+}
+
+
+/* RESULT = P1 - P2  (Montgomery version).*/
+static void
+sub_points_montgomery (mpi_point_t result,
+                       mpi_point_t p1, mpi_point_t p2,
+                       mpi_ec_t ctx)
+{
+  (void)result;
+  (void)p1;
+  (void)p2;
+  (void)ctx;
+  log_fatal ("%s: %s not yet supported\n",
+             "_gcry_mpi_ec_sub_points", "Montgomery");
+}
+
+
+/* RESULT = P1 - P2  (Twisted Edwards version).*/
+static void
+sub_points_edwards (mpi_point_t result,
+                    mpi_point_t p1, mpi_point_t p2,
+                    mpi_ec_t ctx)
+{
+  mpi_point_t p2i = _gcry_mpi_point_new (0);
+  point_set (p2i, p2);
+  mpi_sub (p2i->x, ctx->p, p2i->x);
+  add_points_edwards (result, p1, p2i, ctx);
+  _gcry_mpi_point_release (p2i);
+}
+
+
+/* RESULT = P1 - P2 */
+void
+_gcry_mpi_ec_sub_points (mpi_point_t result,
+                         mpi_point_t p1, mpi_point_t p2,
+                         mpi_ec_t ctx)
+{
+  switch (ctx->model)
+    {
+    case MPI_EC_WEIERSTRASS:
+      sub_points_weierstrass (result, p1, p2, ctx);
+      break;
+    case MPI_EC_MONTGOMERY:
+      sub_points_montgomery (result, p1, p2, ctx);
+      break;
+    case MPI_EC_EDWARDS:
+      sub_points_edwards (result, p1, p2, ctx);
+      break;
+    }
+}
+
+
 /* Scalar point multiplication - the main function for ECC.  If takes
    an integer SCALAR and a POINT as well as the usual context CTX.
    RESULT will be set to the resulting point. */
@@ -1106,16 +1500,32 @@ _gcry_mpi_ec_mul_point (mpi_point_t result,
   unsigned int i, loops;
   mpi_point_struct p1, p2, p1inv;
 
-  if (ctx->model == MPI_EC_EDWARDS)
+  if (ctx->model == MPI_EC_EDWARDS
+      || (ctx->model == MPI_EC_WEIERSTRASS
+          && mpi_is_secure (scalar)))
     {
-      /* Simple left to right binary method.  GECC Algorithm 3.27 */
+      /* Simple left to right binary method.  Algorithm 3.27 from
+       * {author={Hankerson, Darrel and Menezes, Alfred J. and Vanstone, Scott},
+       *  title = {Guide to Elliptic Curve Cryptography},
+       *  year = {2003}, isbn = {038795273X},
+       *  url = {http://www.cacr.math.uwaterloo.ca/ecc/},
+       *  publisher = {Springer-Verlag New York, Inc.}} */
       unsigned int nbits;
       int j;
 
       nbits = mpi_get_nbits (scalar);
-      mpi_set_ui (result->x, 0);
-      mpi_set_ui (result->y, 1);
-      mpi_set_ui (result->z, 1);
+      if (ctx->model == MPI_EC_WEIERSTRASS)
+        {
+          mpi_set_ui (result->x, 1);
+          mpi_set_ui (result->y, 1);
+          mpi_set_ui (result->z, 0);
+        }
+      else
+        {
+          mpi_set_ui (result->x, 0);
+          mpi_set_ui (result->y, 1);
+          mpi_set_ui (result->z, 1);
+        }
 
       if (mpi_is_secure (scalar))
         {
@@ -1124,12 +1534,13 @@ _gcry_mpi_ec_mul_point (mpi_point_t result,
           mpi_point_struct tmppnt;
 
           point_init (&tmppnt);
+          point_resize (result, ctx);
+          point_resize (&tmppnt, ctx);
           for (j=nbits-1; j >= 0; j--)
             {
               _gcry_mpi_ec_dup_point (result, result, ctx);
               _gcry_mpi_ec_add_points (&tmppnt, result, point, ctx);
-              if (mpi_test_bit (scalar, j))
-                point_set (result, &tmppnt);
+              point_swap_cond (result, &tmppnt, mpi_test_bit (scalar, j), ctx);
             }
           point_free (&tmppnt);
         }
@@ -1144,6 +1555,80 @@ _gcry_mpi_ec_mul_point (mpi_point_t result,
         }
       return;
     }
+  else if (ctx->model == MPI_EC_MONTGOMERY)
+    {
+      unsigned int nbits;
+      int j;
+      mpi_point_struct p1_, p2_;
+      mpi_point_t q1, q2, prd, sum;
+      unsigned long sw;
+      mpi_size_t rsize;
+
+      /* Compute scalar point multiplication with Montgomery Ladder.
+         Note that we don't use Y-coordinate in the points at all.
+         RESULT->Y will be filled by zero.  */
+
+      nbits = mpi_get_nbits (scalar);
+      point_init (&p1);
+      point_init (&p2);
+      point_init (&p1_);
+      point_init (&p2_);
+      mpi_set_ui (p1.x, 1);
+      mpi_free (p2.x);
+      p2.x  = mpi_copy (point->x);
+      mpi_set_ui (p2.z, 1);
+
+      point_resize (&p1, ctx);
+      point_resize (&p2, ctx);
+      point_resize (&p1_, ctx);
+      point_resize (&p2_, ctx);
+
+      mpi_resize (point->x, ctx->p->nlimbs);
+      point->x->nlimbs = ctx->p->nlimbs;
+
+      q1 = &p1;
+      q2 = &p2;
+      prd = &p1_;
+      sum = &p2_;
+
+      for (j=nbits-1; j >= 0; j--)
+        {
+          mpi_point_t t;
+
+          sw = mpi_test_bit (scalar, j);
+          point_swap_cond (q1, q2, sw, ctx);
+          montgomery_ladder (prd, sum, q1, q2, point->x, ctx);
+          point_swap_cond (prd, sum, sw, ctx);
+          t = q1;  q1 = prd;  prd = t;
+          t = q2;  q2 = sum;  sum = t;
+        }
+
+      mpi_clear (result->y);
+      sw = (nbits & 1);
+      point_swap_cond (&p1, &p1_, sw, ctx);
+
+      rsize = p1.z->nlimbs;
+      MPN_NORMALIZE (p1.z->d, rsize);
+      if (rsize == 0)
+        {
+          mpi_set_ui (result->x, 1);
+          mpi_set_ui (result->z, 0);
+        }
+      else
+        {
+          z1 = mpi_new (0);
+          ec_invm (z1, p1.z, ctx);
+          ec_mulm (result->x, p1.x, z1, ctx);
+          mpi_set_ui (result->z, 1);
+          mpi_free (z1);
+        }
+
+      point_free (&p1);
+      point_free (&p2);
+      point_free (&p1_);
+      point_free (&p2_);
+      return;
+    }
 
   x1 = mpi_alloc_like (ctx->p);
   y1 = mpi_alloc_like (ctx->p);
@@ -1205,6 +1690,10 @@ _gcry_mpi_ec_mul_point (mpi_point_t result,
   point_init (&p2);
   point_init (&p1inv);
 
+  /* Invert point: y = p - y mod p  */
+  point_set (&p1inv, &p1);
+  ec_subm (p1inv.y, ctx->p, p1inv.y, ctx);
+
   for (i=loops-2; i > 0; i--)
     {
       _gcry_mpi_ec_dup_point (result, result, ctx);
@@ -1216,9 +1705,6 @@ _gcry_mpi_ec_mul_point (mpi_point_t result,
       if (mpi_test_bit (h, i) == 0 && mpi_test_bit (k, i) == 1)
         {
           point_set (&p2, result);
-          /* Invert point: y = p - y mod p  */
-          point_set (&p1inv, &p1);
-          ec_subm (p1inv.y, ctx->p, p1inv.y, ctx);
           _gcry_mpi_ec_add_points (result, &p2, &p1inv, ctx);
         }
     }
@@ -1242,14 +1728,16 @@ _gcry_mpi_ec_curve_point (gcry_mpi_point_t point, mpi_ec_t ctx)
   y = mpi_new (0);
   w = mpi_new (0);
 
-  if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
-    return 0;
-
   switch (ctx->model)
     {
     case MPI_EC_WEIERSTRASS:
       {
-        gcry_mpi_t xxx = mpi_new (0);
+        gcry_mpi_t xxx;
+
+        if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
+          goto leave;
+
+        xxx = mpi_new (0);
 
         /* y^2 == x^3 + a·x + b */
         ec_pow2 (y, y, ctx);
@@ -1266,19 +1754,45 @@ _gcry_mpi_ec_curve_point (gcry_mpi_point_t point, mpi_ec_t ctx)
       }
       break;
     case MPI_EC_MONTGOMERY:
-      log_fatal ("%s: %s not yet supported\n",
-                 "_gcry_mpi_ec_curve_point", "Montgomery");
+      {
+#define xx y
+        /* With Montgomery curve, only X-coordinate is valid.  */
+        if (_gcry_mpi_ec_get_affine (x, NULL, point, ctx))
+          goto leave;
+
+        /* The equation is: b * y^2 == x^3 + a · x^2 + x */
+        /* We check if right hand is quadratic residue or not by
+           Euler's criterion.  */
+        /* CTX->A has (a-2)/4 and CTX->B has b^-1 */
+        ec_mulm (w, ctx->a, mpi_const (MPI_C_FOUR), ctx);
+        ec_addm (w, w, mpi_const (MPI_C_TWO), ctx);
+        ec_mulm (w, w, x, ctx);
+        ec_pow2 (xx, x, ctx);
+        ec_addm (w, w, xx, ctx);
+        ec_addm (w, w, mpi_const (MPI_C_ONE), ctx);
+        ec_mulm (w, w, x, ctx);
+        ec_mulm (w, w, ctx->b, ctx);
+#undef xx
+        /* Compute Euler's criterion: w^(p-1)/2 */
+#define p_minus1 y
+        ec_subm (p_minus1, ctx->p, mpi_const (MPI_C_ONE), ctx);
+        mpi_rshift (p_minus1, p_minus1, 1);
+        ec_powm (w, w, p_minus1, ctx);
+
+        res = !mpi_cmp_ui (w, 1);
+#undef p_minus1
+      }
       break;
     case MPI_EC_EDWARDS:
       {
+        if (_gcry_mpi_ec_get_affine (x, y, point, ctx))
+          goto leave;
+
         /* a · x^2 + y^2 - 1 - b · x^2 · y^2 == 0 */
         ec_pow2 (x, x, ctx);
         ec_pow2 (y, y, ctx);
         if (ctx->dialect == ECC_DIALECT_ED25519)
-          {
-            mpi_set (w, x);
-            _gcry_mpi_neg (w, w);
-          }
+          mpi_sub (w, ctx->p, x);
         else
           ec_mulm (w, ctx->a, x, ctx);
         ec_addm (w, w, y, ctx);
@@ -1292,9 +1806,24 @@ _gcry_mpi_ec_curve_point (gcry_mpi_point_t point, mpi_ec_t ctx)
       break;
     }
 
+ leave:
   _gcry_mpi_release (w);
   _gcry_mpi_release (x);
   _gcry_mpi_release (y);
 
   return res;
 }
+
+
+int
+_gcry_mpi_ec_bad_point (gcry_mpi_point_t point, mpi_ec_t ctx)
+{
+  int i;
+  gcry_mpi_t x_bad;
+
+  for (i = 0; (x_bad = ctx->t.scratch[i]); i++)
+    if (!mpi_cmp (point->x, x_bad))
+      return 1;
+
+  return 0;
+}