Implemented ecc_a_to_eh and ecc_mul_a_eh, for curve25519.

7 files changed, 481 insertions, 7 deletions

diff --git a/ChangeLog b/ChangeLog
index 62100ae2..6fcc42a8 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,23 @@
+2014-07-29  Niels Möller  <nisse@lysator.liu.se>
+
+	* ecc-internal.h (ECC_MUL_A_EH_WBITS): New constant.
+	(ECC_A_TO_EH_ITCH, ECC_MUL_A_EH_ITCH): New macros.
+	* ecc-a-to-eh.c (ecc_a_to_eh, ecc_a_to_eh_itch): New file, new
+	functions.
+	* ecc-mul-a-eh.c: New file.
+	(ecc_mul_a_eh): New function. The case [ECC_MUL_A_EH_WBITS > 0]
+	not yet working).
+	(ecc_mul_a_eh_itch): New function.
+	* ecc.h: Declare new functions.
+	* Makefile.in (hogweed_SOURCES): Added ecc-a-to-eh.c and
+	ecc-mul-a-eh.c.
+
+	* testsuite/curve25519-dh-test.c (curve25519_sqrt): New function.
+	(curve_25519): Use ecc_mul_a_eh.
+	(test_a): New function.
+	(test_main): Test construction of shared secret, using scalar
+	multiplication with points other than the fix generator.
+
 2014-07-26  Niels Möller  <nisse@lysator.liu.se>
 
 	* ecc-add-ehh.c (ecc_add_ehh): Reduce scratch need.
diff --git a/Makefile.in b/Makefile.in
index bcabbf44..1e6cdd8f 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -167,8 +167,9 @@ hogweed_SOURCES = sexp.c sexp-format.c \
 		  ecc-25519.c \
 		  ecc-size.c ecc-j-to-a.c ecc-a-to-j.c \
 		  ecc-dup-jj.c ecc-add-jja.c ecc-add-jjj.c \
-		  ecc-dup-eh.c ecc-add-eh.c ecc-add-ehh.c ecc-eh-to-a.c \
-		  ecc-mul-g-eh.c \
+		  ecc-a-to-eh.c ecc-eh-to-a.c \
+		  ecc-dup-eh.c ecc-add-eh.c ecc-add-ehh.c \
+		  ecc-mul-g-eh.c ecc-mul-a-eh.c \
 		  ecc-mul-g.c ecc-mul-a.c ecc-hash.c ecc-random.c \
 		  ecc-point.c ecc-scalar.c ecc-point-mul.c ecc-point-mul-g.c \
 		  ecc-ecdsa-sign.c ecdsa-sign.c \
diff --git a/ecc-a-to-eh.c b/ecc-a-to-eh.c
new file mode 100644
index 00000000..7f77394e
--- /dev/null
+++ b/ecc-a-to-eh.c
@@ -0,0 +1,77 @@
+/* ecc-a-to-eh.c
+
+   Copyright (C) 2014 Niels Möller
+
+   This file is part of GNU Nettle.
+
+   GNU Nettle is free software: you can redistribute it and/or
+   modify it under the terms of either:
+
+     * the GNU Lesser General Public License as published by the Free
+       Software Foundation; either version 3 of the License, or (at your
+       option) any later version.
+
+   or
+
+     * the GNU General Public License as published by the Free
+       Software Foundation; either version 2 of the License, or (at your
+       option) any later version.
+
+   or both in parallel, as here.
+
+   GNU Nettle 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 GNU
+   General Public License for more details.
+
+   You should have received copies of the GNU General Public License and
+   the GNU Lesser General Public License along with this program.  If
+   not, see http://www.gnu.org/licenses/.
+*/
+
+#if HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include "ecc.h"
+#include "ecc-internal.h"
+
+mp_size_t
+ecc_a_to_eh_itch (const struct ecc_curve *ecc)
+{
+  return ECC_A_TO_EH_ITCH (ecc->size);
+}
+
+/* Convert from affine coordinates to homogeneous coordinates on the
+   corresponding Edwards curve. */
+void
+ecc_a_to_eh (const struct ecc_curve *ecc,
+	     mp_limb_t *r, const mp_limb_t *p,
+	     mp_limb_t *scratch)
+{
+#define xp p
+#define yp (p + ecc->size)
+
+#define up r
+#define vp (r + ecc->size)
+#define wp (r + 2*ecc->size)
+
+  /* u = t x / y
+     v = (x-1) / (x+1)
+
+     or in homogeneous coordinates
+
+     U = t x (x+1)
+     V = (x-1) y
+     W = (x+1) y
+  */
+
+  ecc_modp_mul (ecc, scratch, xp, yp);
+  ecc_modp_add (ecc, wp, scratch, yp);
+  ecc_modp_sub (ecc, vp, scratch, yp);
+
+  ecc_modp_sqr (ecc, scratch, xp);
+  ecc_modp_add (ecc, up, scratch, xp);
+  ecc_modp_mul (ecc, scratch, up, ecc->edwards_root);
+  mpn_copyi (up, scratch, ecc->size);
+}
diff --git a/ecc-internal.h b/ecc-internal.h
index 99f74169..e233b64f 100644
--- a/ecc-internal.h
+++ b/ecc-internal.h
@@ -72,6 +72,8 @@
    up to 5 bits, but I don't think that's worth doubling the
    storage. */
 #define ECC_MUL_A_WBITS 4
+#define ECC_MUL_A_EH_WBITS 0
+
 
 /* Reduces from 2*ecc->size to ecc->size. */
 /* Required to return a result < 2q. This property is inherited by
@@ -238,6 +240,7 @@ sec_modinv (mp_limb_t *vp, mp_limb_t *ap, mp_size_t n,
 #define ECC_MODINV_ITCH(size) (3*(size))
 #define ECC_J_TO_A_ITCH(size) (5*(size))
 #define ECC_EH_TO_A_ITCH(size) (5*(size))
+#define ECC_A_TO_EH_ITCH(size) (2*(size))
 #define ECC_DUP_JJ_ITCH(size) (5*(size))
 #define ECC_DUP_EH_ITCH(size) (5*(size))
 #define ECC_ADD_JJA_ITCH(size) (6*(size))
@@ -252,6 +255,12 @@ sec_modinv (mp_limb_t *vp, mp_limb_t *ap, mp_size_t n,
 #define ECC_MUL_A_ITCH(size) \
   (((3 << ECC_MUL_A_WBITS) + 11) * (size))
 #endif
+#if ECC_MUL_A_EH_WBITS == 0
+#define ECC_MUL_A_EH_ITCH(size) (13*(size))
+#else
+#define ECC_MUL_A_EH_ITCH(size) \
+  (((3 << ECC_MUL_A_EH_WBITS) + 10) * (size))
+#endif
 #define ECC_ECDSA_SIGN_ITCH(size) (12*(size))
 #define ECC_ECDSA_VERIFY_ITCH(size) \
   (6*(size) + ECC_MUL_A_ITCH ((size)))
diff --git a/ecc-mul-a-eh.c b/ecc-mul-a-eh.c
new file mode 100644
index 00000000..ad017565
--- /dev/null
+++ b/ecc-mul-a-eh.c
@@ -0,0 +1,184 @@
+/* ecc-mul-a-eh.c
+
+   Copyright (C) 2013, 2014 Niels Möller
+
+   This file is part of GNU Nettle.
+
+   GNU Nettle is free software: you can redistribute it and/or
+   modify it under the terms of either:
+
+     * the GNU Lesser General Public License as published by the Free
+       Software Foundation; either version 3 of the License, or (at your
+       option) any later version.
+
+   or
+
+     * the GNU General Public License as published by the Free
+       Software Foundation; either version 2 of the License, or (at your
+       option) any later version.
+
+   or both in parallel, as here.
+
+   GNU Nettle 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 GNU
+   General Public License for more details.
+
+   You should have received copies of the GNU General Public License and
+   the GNU Lesser General Public License along with this program.  If
+   not, see http://www.gnu.org/licenses/.
+*/
+
+#if HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <assert.h>
+
+#include "ecc.h"
+#include "ecc-internal.h"
+
+mp_size_t
+ecc_mul_a_eh_itch (const struct ecc_curve *ecc)
+{
+  /* Binary algorithm needs 6*ecc->size + scratch for ecc_add_ehh,
+     total 13 ecc->size
+
+     Window algorithm needs (3<<w) * ecc->size for the table,
+     3*ecc->size for a temporary point, and scratch for
+     ecc_add_ehh. */
+  return ECC_MUL_A_EH_ITCH (ecc->size);
+}
+
+#if ECC_MUL_A_EH_WBITS == 0
+void
+ecc_mul_a_eh (const struct ecc_curve *ecc,
+	      mp_limb_t *r,
+	      const mp_limb_t *np, const mp_limb_t *p,
+	      mp_limb_t *scratch)
+{
+#define pe scratch
+#define tp (scratch + 3*ecc->size)
+#define scratch_out (scratch + 6*ecc->size)
+
+  unsigned i;
+
+  ecc_a_to_eh (ecc, pe, p, pe + 3*ecc->size);
+
+  /* x = 0, y = 1, z = 1 */
+  mpn_zero (r, 3*ecc->size);
+  r[ecc->size] = r[2*ecc->size] = 1;
+  
+  for (i = ecc->size; i-- > 0; )
+    {
+      mp_limb_t w = np[i];
+      mp_limb_t bit;
+
+      for (bit = (mp_limb_t) 1 << (GMP_NUMB_BITS - 1);
+	   bit > 0;
+	   bit >>= 1)
+	{
+	  int digit;
+
+	  ecc_dup_eh (ecc, r, r, scratch_out);
+	  ecc_add_ehh (ecc, tp, r, pe, scratch_out);
+
+	  digit = (w & bit) > 0;
+	  /* If we had a one-bit, use the sum. */
+	  cnd_copy (digit, r, tp, 3*ecc->size);
+	}
+    }
+}
+#else /* ECC_MUL_A_EH_WBITS > 1 */
+
+#error Not yet working
+
+#define TABLE_SIZE (1U << ECC_MUL_A_EH_WBITS)
+#define TABLE_MASK (TABLE_SIZE - 1)
+
+#define TABLE(j) (table + (j) * 3*ecc->size)
+
+static void
+table_init (const struct ecc_curve *ecc,
+	    mp_limb_t *table, unsigned bits,
+	    const mp_limb_t *p,
+	    mp_limb_t *scratch)
+{
+  unsigned size = 1 << bits;
+  unsigned j;
+
+  mpn_zero (TABLE(0), 3*ecc->size);
+  TABLE(0)[ecc->size] = TABLE(0)[2*ecc->size] = 1;
+
+  ecc_a_to_eh (ecc, TABLE(1), p, scratch);
+  mpn_copyi (TABLE(1), p, 3*ecc->size);
+
+  for (j = 2; j < size; j += 2)
+    {
+      ecc_dup_eh (ecc, TABLE(j), TABLE(j/2), scratch);
+      ecc_add_ehh (ecc, TABLE(j+1), TABLE(j), TABLE(1), scratch);
+    }
+}
+
+void
+ecc_mul_a_eh (const struct ecc_curve *ecc,
+	      mp_limb_t *r,
+	      const mp_limb_t *np, const mp_limb_t *p,
+	      mp_limb_t *scratch)
+{
+#define tp scratch
+#define table (scratch + 3*ecc->size)
+  mp_limb_t *scratch_out = table + (3*ecc->size << ECC_MUL_A_EH_WBITS);
+
+  /* Avoid the mp_bitcnt_t type for compatibility with older GMP
+     versions. */
+  unsigned blocks = (ecc->bit_size + ECC_MUL_A_EH_WBITS - 1) / ECC_MUL_A_EH_WBITS;
+  unsigned bit_index = (blocks-1) * ECC_MUL_A_EH_WBITS;
+
+  mp_size_t limb_index = bit_index / GMP_NUMB_BITS;
+  unsigned shift = bit_index % GMP_NUMB_BITS;
+  mp_limb_t w, bits;
+
+  table_init (ecc, table, ECC_MUL_A_EH_WBITS, p, scratch_out);
+
+  w = np[limb_index];
+  bits = w >> shift;
+  if (limb_index < ecc->size - 1)
+    bits |= np[limb_index + 1] << (GMP_NUMB_BITS - shift);
+
+  assert (bits < TABLE_SIZE);
+
+  sec_tabselect (r, 3*ecc->size, table, TABLE_SIZE, bits);
+
+  for (;;)
+    {
+      unsigned j;
+      if (shift >= ECC_MUL_A_EH_WBITS)
+	{
+	  shift -= ECC_MUL_A_EH_WBITS;
+	  bits = w >> shift;
+	}
+      else
+	{
+	  if (limb_index == 0)
+	    {
+	      assert (shift == 0);
+	      break;
+	    }
+	  bits = w << (ECC_MUL_A_EH_WBITS - shift);
+	  w = np[--limb_index];
+	  shift = shift + GMP_NUMB_BITS - ECC_MUL_A_EH_WBITS;
+	  bits |= w >> shift;
+	}
+      for (j = 0; j < ECC_MUL_A_EH_WBITS; j++)
+	ecc_dup_eh (ecc, r, r, scratch_out);
+
+      bits &= TABLE_MASK;
+      sec_tabselect (tp, 3*ecc->size, table, TABLE_SIZE, bits);
+      ecc_add_ehh (ecc, r, tp, r, scratch_out);
+    }
+#undef table
+#undef tp
+}
+
+#endif /* ECC_MUL_A_EH_WBITS > 1 */
diff --git a/ecc.h b/ecc.h
index e786da53..f37b8f34 100644
--- a/ecc.h
+++ b/ecc.h
@@ -63,6 +63,8 @@ extern "C" {
 #define ecc_j_to_a nettle_ecc_j_to_a
 #define ecc_eh_to_a_itch nettle_ecc_eh_to_a_itch
 #define ecc_eh_to_a nettle_ecc_eh_to_a
+#define ecc_a_to_eh_itch nettle_ecc_a_to_eh_itch
+#define ecc_a_to_eh nettle_ecc_a_to_eh
 #define ecc_dup_jj_itch nettle_ecc_dup_jj_itch
 #define ecc_dup_jj nettle_ecc_dup_jj
 #define ecc_add_jja_itch nettle_ecc_add_jja_itch
@@ -81,6 +83,8 @@ extern "C" {
 #define ecc_mul_a nettle_ecc_mul_a
 #define ecc_mul_g_eh_itch nettle_ecc_mul_g_eh_itch
 #define ecc_mul_g_eh nettle_ecc_mul_g_eh
+#define ecc_mul_eh_itch nettle_ecc_mul_eh_itch
+#define ecc_mul_eh nettle_ecc_mul_eh
 
 struct ecc_curve;
 
@@ -205,6 +209,13 @@ ecc_eh_to_a (const struct ecc_curve *ecc,
 	     mp_limb_t *r, const mp_limb_t *p,
 	     mp_limb_t *scratch);
 
+mp_size_t
+ecc_a_to_eh_itch (const struct ecc_curve *ecc);
+void
+ecc_a_to_eh (const struct ecc_curve *ecc,
+	     mp_limb_t *r, const mp_limb_t *p,
+	     mp_limb_t *scratch);
+
 /* Group operations */
 
 /* Point doubling, with jacobian input and output. Corner cases:
@@ -292,6 +303,15 @@ void
 ecc_mul_g_eh (const struct ecc_curve *ecc, mp_limb_t *r,
 	      const mp_limb_t *np, mp_limb_t *scratch);
 
+mp_size_t
+ecc_mul_eh_itch (const struct ecc_curve *ecc);
+void
+ecc_mul_eh (const struct ecc_curve *ecc,
+	    mp_limb_t *r,
+	    const mp_limb_t *np, const mp_limb_t *p,
+	    mp_limb_t *scratch);
+
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/testsuite/curve25519-dh-test.c b/testsuite/curve25519-dh-test.c
index 623be19c..d2467548 100644
--- a/testsuite/curve25519-dh-test.c
+++ b/testsuite/curve25519-dh-test.c
@@ -31,6 +31,89 @@
 
 #include "testutils.h"
 
+static
+int curve25519_sqrt (const struct ecc_curve *ecc,
+		     mp_limb_t *rp, const mp_limb_t *ap)
+{
+  /* p-1 = 2^{255} - 20 = 4 (2^{253] - 5), s = 2^{253} - 5, e = 2 */
+
+  mpz_t g;
+  mpz_t sm1h; /* (s-1)/2 */
+  mpz_t x;
+  mpz_t a;
+  mpz_t p;
+  mpz_t b;
+  mpz_t t;
+  int success;
+
+  mpz_init_set_str (g,
+		    "2b8324804fc1df0b2b4d00993dfbd7a7"
+		    "2f431806ad2fe478c4ee1b274a0ea0b0", 16);
+  mpz_init_set_str (sm1h,
+		    "fffffffffffffffffffffffffffffff"
+		    "fffffffffffffffffffffffffffffffd", 16);
+
+  mpz_init (x);
+  mpz_init (b);
+  mpz_init (t);
+  mpz_roinit_n (p, ecc->p, ecc->size);
+  mpz_roinit_n (a, ap, ecc->size);
+
+  mpz_powm (x, a, sm1h, p);
+  mpz_mul (b, x, x); /* s-1 */
+  mpz_mul (b, b, a); /* s */
+  mpz_mod (b, b, p);
+  mpz_mul (x, x, a); /* (s+1)/2 */
+  mpz_mod (x, x, p);
+
+  if (mpz_cmp_ui (b, 1) != 0)
+    {
+      mpz_t t;
+      unsigned m, e;
+      mpz_init (t);
+      e = 2;
+      do
+	{
+	  mpz_set (t, b);
+	  m = 0;
+	  do
+	    {	  
+	      m++;
+	      if (m == e)
+		{
+		  mpz_clear (t);
+		  success = 0;
+		  goto done;
+		}
+	      mpz_mul (t, t, t);
+	      mpz_mod (t, t, p);
+	    }
+	  while (mpz_cmp_ui (t, 1) != 0);
+	  ASSERT (m < e);
+	  mpz_set_ui (t, 1UL << (e - m - 1));
+	  mpz_powm (t, g, t, p);
+	  mpz_mul (x, x, t);
+	  mpz_mod (x, x, p);
+	  mpz_mul (g, t, t);
+	  mpz_mod (g, g, p);
+	  mpz_mul (b, b, g);
+	  mpz_mod (b, b, p);
+	  e = m-1;
+	}
+      while (mpz_cmp_ui (b, 1) != 0);
+
+      mpz_clear (t);
+    }
+  mpz_limbs_copy (rp, x, ecc->size);
+  success = 1;
+ done:
+  mpz_clear (g);
+  mpz_clear (sm1h);
+  mpz_clear (x);
+  mpz_clear (b);
+  return success;
+}
+
 /* Computes the x coordinate of s G, where g is a scalar, and G is the
    base point on the curve. If x is non-NULL, it gives the X
    coordinate of the point G, otherwise, G is the specified
@@ -43,13 +126,31 @@ curve_25519 (const struct ecc_curve *ecc,
   mp_limb_t *scratch;
   mp_size_t itch;
 
+  p = gmp_alloc_limbs (3*ecc->size);
+  
   if (x)
-    die ("Not yet implemented.\n");
+    {
+      itch = ECC_MUL_A_EH_ITCH (ecc->size);
+      scratch = gmp_alloc_limbs (itch);
+      mpn_copyi (p, x, ecc->size);
+      /* y^2 = x^3 + b x^2 + x = (x^2 + bx + 1) x = ((x+b)x + 1) x */
+      ecc_modp_sqr (ecc, scratch, x);
+      ecc_modp_addmul_1 (ecc, scratch, x, 0x76d06ULL);
+      ecc_modp_add (ecc, scratch, scratch, ecc->unit);
+      ecc_modp_mul (ecc, scratch + ecc->size, scratch, x);
+
+      if (!curve25519_sqrt (ecc, p + ecc->size, scratch + ecc->size))
+	die ("Point not on curve.\n");
+      mpn_copyi (p, x, ecc->size);
+      ecc_mul_a_eh (ecc, p, s, p, scratch);
+    }
+  else
+    {
+      itch = ECC_MUL_G_EH_ITCH (ecc->size);
+      scratch = gmp_alloc_limbs (itch);
+      ecc_mul_g_eh (ecc, p, s, scratch);
+    }
 
-  itch = ECC_MUL_G_EH_ITCH (ecc->size);
-  p = gmp_alloc_limbs (3*ecc->size);
-  scratch = gmp_alloc_limbs (itch);
-  ecc_mul_g_eh (ecc, p, s, scratch);
   ecc_eh_to_a (ecc, 2, r,  p, scratch);
 
   /* FIXME: Convert to little-endian here? */
@@ -93,6 +194,48 @@ test_g (const char *sz, const char *pz)
   mpz_clear (X);
 }
 
+static void
+test_a (const char *bz, const char *sz, const char *pz)
+{
+  mpz_t B, S, R, X;
+  const struct ecc_curve *ecc = &nettle_curve25519;
+
+  mpz_init (B);
+  mpz_init (S);
+  mpz_init (R);
+  mpz_init (X);
+
+  mpz_set_str (B, bz, 16);
+  mpz_set_str (S, sz, 16);
+  mpz_set_str (R, pz, 16);
+
+  ASSERT (mpz_size (S) == ecc->size);
+  ASSERT (mpz_size (B) == ecc->size);
+  
+  curve_25519 (ecc, mpz_limbs_write (X, ecc->size),
+	       mpz_limbs_read (S), mpz_limbs_read (B));
+
+  mpz_limbs_finish (X, ecc->size);
+  if (mpz_cmp (X, R) != 0)
+    {
+      fprintf (stderr, "curve25519 failure:\nB = ");
+      mpz_out_str (stderr, 16, B);
+      fprintf (stderr, "\nS = ");
+      mpz_out_str (stderr, 16, S);
+      fprintf (stderr, "\nX = ");
+      mpz_out_str (stderr, 16, X);
+      fprintf (stderr, " (bad)\nR = ");
+      mpz_out_str (stderr, 16, R);
+      fprintf (stderr, " (expected)\n");
+      abort ();
+    }
+
+  mpz_clear (B);
+  mpz_clear (S);
+  mpz_clear (R);
+  mpz_clear (X);
+}
+
 void
 test_main (void)
 {
@@ -100,11 +243,31 @@ test_main (void)
      the P values, though. */
   test_g ("6A2CB91DA5FB77B12A99C0EB872F4CDF"
 	  "4566B25172C1163C7DA518730A6D0770",
+
 	  "6A4E9BAA8EA9A4EBF41A38260D3ABF0D"
 	  "5AF73EB4DC7D8B7454A7308909F02085");
 
   test_g ("6BE088FF278B2F1CFDB6182629B13B6F"
 	  "E60E80838B7FE1794B8A4A627E08AB58",
+
 	  "4F2B886F147EFCAD4D67785BC843833F"
 	  "3735E4ECC2615BD3B4C17D7B7DDB9EDE");
+
+  test_a ("4F2B886F147EFCAD4D67785BC843833F"
+	  "3735E4ECC2615BD3B4C17D7B7DDB9EDE",
+
+	  "6A2CB91DA5FB77B12A99C0EB872F4CDF"
+	  "4566B25172C1163C7DA518730A6D0770",
+	  
+	  "4217161E3C9BF076339ED147C9217EE0"
+	  "250F3580F43B8E72E12DCEA45B9D5D4A");
+
+  test_a ("6A4E9BAA8EA9A4EBF41A38260D3ABF0D"
+	  "5AF73EB4DC7D8B7454A7308909F02085",
+	  
+	  "6BE088FF278B2F1CFDB6182629B13B6F"
+	  "E60E80838B7FE1794B8A4A627E08AB58",
+
+	  "4217161E3C9BF076339ED147C9217EE0"
+	  "250F3580F43B8E72E12DCEA45B9D5D4A");
 }