mpi: Constant time mpi_inv with some conditions.

* mpi/mpi-inv.c (mpih_add_n_cond, mpih_sub_n_cond, mpih_swap_cond)
(mpih_abs_cond): New.
(mpi_invm_odd): New.
(mpi_invm_generic): Rename from _gcry_mpi_invm.
(_gcry_mpi_invm): Use mpi_invm_odd for usual odd cases.

--

GnuPG-bug-id: 4869
Signed-off-by: NIIBE Yutaka <gniibe@fsij.org>

1 files changed, 200 insertions, 7 deletions

diff --git a/mpi/mpi-inv.c b/mpi/mpi-inv.c
index ee6813b1..0114622d 100644
--- a/mpi/mpi-inv.c
+++ b/mpi/mpi-inv.c
@@ -23,13 +23,196 @@
 #include "mpi-internal.h"
 #include "g10lib.h"
 
+/*
+ *  W = U + V when OP_ENABLED=1
+ *  otherwise, W = U
+ */
+static mpi_limb_t
+mpih_add_n_cond (mpi_ptr_t wp, mpi_ptr_t up, mpi_ptr_t vp, mpi_size_t usize,
+                 unsigned long op_enable)
+{
+  mpi_size_t i;
+  mpi_limb_t cy;
+  mpi_limb_t mask = ((mpi_limb_t)0) - op_enable;
+
+  cy = 0;
+  for (i = 0; i < usize; i++)
+    {
+      mpi_limb_t x = up[i] + (vp[i] & mask);
+      mpi_limb_t cy1 = x < up[i];
+      mpi_limb_t cy2;
+
+      x = x + cy;
+      cy2 = x < cy;
+      cy = cy1 | cy2;
+      wp[i] = x;
+    }
+
+  return cy;
+}
+
+
+/*
+ *  W = U - V when OP_ENABLED=1
+ *  otherwise, W = U
+ */
+static mpi_limb_t
+mpih_sub_n_cond (mpi_ptr_t wp, mpi_ptr_t up, mpi_ptr_t vp, mpi_size_t usize,
+                 unsigned long op_enable)
+{
+  mpi_size_t i;
+  mpi_limb_t cy;
+  mpi_limb_t mask = ((mpi_limb_t)0) - op_enable;
+
+  cy = 0;
+  for (i = 0; i < usize; i++)
+    {
+      mpi_limb_t x = up[i] - (vp[i] & mask);
+      mpi_limb_t cy1 = x > up[i];
+      mpi_limb_t cy2;
+
+      cy2 = x < cy;
+      x = x - cy;
+      cy = cy1 | cy2;
+      wp[i] = x;
+    }
+
+  return cy;
+}
+
+
+/*
+ *  Swap value of U and V when OP_ENABLED=1
+ *  otherwise, no change
+ */
+static void
+mpih_swap_cond (mpi_ptr_t up, mpi_ptr_t vp, mpi_size_t usize,
+                unsigned long op_enable)
+{
+  mpi_size_t i;
+  mpi_limb_t mask = ((mpi_limb_t)0) - op_enable;
+
+  for (i = 0; i < usize; i++)
+    {
+      mpi_limb_t x = mask & (up[i] ^ vp[i]);
+
+      up[i] = up[i] ^ x;
+      vp[i] = vp[i] ^ x;
+    }
+}
+
+
+/*
+ *  W = -U when OP_ENABLED=1
+ *  otherwise, W = U
+ */
+static void
+mpih_abs_cond (mpi_limb_t *wp, const mpi_limb_t *up, mpi_size_t usize,
+               unsigned long op_enable)
+{
+  mpi_size_t i;
+  mpi_limb_t mask = ((mpi_limb_t)0) - op_enable;
+  mpi_limb_t cy = op_enable;
+
+  for (i = 0; i < usize; i++)
+    {
+      mpi_limb_t x = ~up[i] + cy;
+
+      cy = (x < ~up[i]);
+      wp[i] = up[i] ^ (mask & (x ^ up[i]));
+    }
+}
+
+
+/*
+ * This uses a modular inversion algorithm designed by Niels Möller
+ * which was implemented in Nettle.  The same algorithm was later also
+ * adapted to GMP in mpn_sec_invert.
+ *
+ * For the description of the algorithm, see Algorithm 5 in Appendix A
+ * of "Fast Software Polynomial Multiplication on ARM Processors using
+ * the NEON Engine" by Danilo Câmara, Conrado P. L. Gouvêa, Julio
+ * López, and Ricardo Dahab:
+ *   https://hal.inria.fr/hal-01506572/document
+ *
+ * Note that in the reference above, at the line 2 of Algorithm 5,
+ * initial value of V was described as V:=1 wrongly.  It must be V:=0.
+ */
+static int
+mpi_invm_odd (gcry_mpi_t x, gcry_mpi_t a_orig, gcry_mpi_t n)
+{
+  mpi_size_t nsize;
+  gcry_mpi_t a, b, n1h;
+  gcry_mpi_t u;
+  unsigned int iterations;
+  mpi_ptr_t ap, bp, n1hp;
+  mpi_ptr_t up, vp;
+  int is_gcd_one;
+
+  nsize = n->nlimbs;
+
+  a = mpi_copy (a_orig);
+  mpi_resize (a, nsize);
+  ap = a->d;
+
+  b = mpi_copy (n);
+  bp = b->d;
+
+  u = mpi_alloc_set_ui (1);
+  mpi_resize (u, nsize);
+  up = u->d;
+
+  mpi_resize (x, nsize);
+  x->nlimbs = nsize;
+  vp = x->d;
+  memset (vp, 0, nsize * BYTES_PER_MPI_LIMB);
+
+  n1h = mpi_copy (n);
+  mpi_rshift (n1h, n1h, 1);
+  mpi_add_ui (n1h, n1h, 1);
+  mpi_resize (n1h, nsize);
+
+  n1hp = n1h->d;
+
+  iterations = 2 * nsize * BITS_PER_MPI_LIMB;
+
+  while (iterations-- > 0)
+    {
+      mpi_limb_t odd_a, odd_u, underflow, borrow;
+
+      odd_a = ap[0] & 1;
+
+      underflow = mpih_sub_n_cond (ap, ap, bp, nsize, odd_a);
+      mpih_add_n_cond (bp, bp, ap, nsize, underflow);
+      mpih_abs_cond (ap, ap, nsize, underflow);
+      mpih_swap_cond (up, vp, nsize, underflow);
+
+      _gcry_mpih_rshift (ap, ap, nsize, 1);
+
+      borrow = mpih_sub_n_cond (up, up, vp, nsize, odd_a);
+      mpih_add_n_cond (up, up, n->d, nsize, borrow);
+
+      odd_u = _gcry_mpih_rshift (up, up, nsize, 1) != 0;
+      mpih_add_n_cond (up, up, n1hp, nsize, odd_u);
+    }
+
+  is_gcd_one = (mpi_cmp_ui (b, 1) == 0);
+
+  mpi_free (n1h);
+  mpi_free (u);
+  mpi_free (b);
+  mpi_free (a);
+
+  return is_gcd_one;
+}
+
 /****************
  * Calculate the multiplicative inverse X of A mod N
  * That is: Find the solution x for
  *		1 = (a*x) mod n
  */
-int
-_gcry_mpi_invm (gcry_mpi_t x, gcry_mpi_t a, gcry_mpi_t n)
+static int
+mpi_invm_generic (gcry_mpi_t x, gcry_mpi_t a, gcry_mpi_t n)
 {
 #if 0
     gcry_mpi_t u, v, u1, u2, u3, v1, v2, v3, q, t1, t2, t3;
@@ -165,11 +348,6 @@ _gcry_mpi_invm (gcry_mpi_t x, gcry_mpi_t a, gcry_mpi_t n)
     int sign;
     int odd ;
 
-    if (!mpi_cmp_ui (a, 0))
-        return 0; /* Inverse does not exists.  */
-    if (!mpi_cmp_ui (n, 1))
-        return 0; /* Inverse does not exists.  */
-
     u = mpi_copy(a);
     v = mpi_copy(n);
 
@@ -270,3 +448,18 @@ _gcry_mpi_invm (gcry_mpi_t x, gcry_mpi_t a, gcry_mpi_t n)
 #endif
     return 1;
 }
+
+
+int
+_gcry_mpi_invm (gcry_mpi_t x, gcry_mpi_t a, gcry_mpi_t n)
+{
+  if (!mpi_cmp_ui (a, 0))
+    return 0; /* Inverse does not exists.  */
+  if (!mpi_cmp_ui (n, 1))
+    return 0; /* Inverse does not exists.  */
+
+  if (mpi_test_bit (n, 0) && mpi_cmp (a, n) < 0)
+    return mpi_invm_odd (x, a, n);
+  else
+    return mpi_invm_generic (x, a, n);
+}