Make uint32 lookup table.

src/Makefile.am

@@ -1,4 +1,4 @@
-bin_PROGRAMS = rms ddlog_bench
+bin_PROGRAMS = rms ddlog_bench exp_bench
 #check_programs = test_ssl1
 
 DDLOG = ddlog.c ddlog.h
@@ -6,8 +6,10 @@ ELGAMAL = elgamal.c elgamal.h
 ENTROPY = entropy.c entropy.h
 GROUP = group.c group.h
 HSS = hss.c hss.h
+TIMEIT = timeit.h
 #TESTS = $(check_programs)
 
 #test_ssl1_SOURCES = test_ssl1.c
-rms_SOURCES = rms.c $(DDLOG) $(ELGAMAL) $(ENTROPY) $(GROUP) $(HSS)
-ddlog_bench_SOURCES = ddlog_bench.c $(DDLOG) $(ENTROPY) $(GROUP)
+rms_SOURCES = rms.c $(DDLOG) $(ELGAMAL) $(ENTROPY) $(GROUP) $(HSS) $(TIMEIT)
+ddlog_bench_SOURCES = ddlog_bench.c $(DDLOG) $(ENTROPY) $(GROUP) $(TIMEIT)
+exp_bench_SOURCES = exp_bench.c $(ENTROPY) $(GROUP) $(TIMEIT)

src/ddlog.c

@@ -10,8 +10,8 @@
 
 typedef __uint128_t uint128_t;
 
-uint8_t lookup[256];
-uint8_t offset[256];
+uint32_t lookup[256];
+uint32_t offset[256];
 
 static inline
 void add_1(uint64_t *b, const size_t start, const size_t len, uint64_t a)
@@ -45,7 +45,7 @@ uint32_t __attribute__((optimize("unroll-loops"))) convert(uint64_t * nn)
   for (uint32_t w2 = halfstrip_size; w2 < 64-halfstrip_size; w2 += halfstrip_size) {
     if (!(x & (topmask >> w2))) {
       const size_t previous = (x >> (64 - halfstrip_size - w2 + halfstrip_size)) & bottommask;
-      const uint8_t next =    (x >> (64 - halfstrip_size - w2 - halfstrip_size)) & bottommask;
+      const uint32_t next =    (x >> (64 - halfstrip_size - w2 - halfstrip_size)) & bottommask;
       if (next <= lookup[previous]) return w2 - offset[previous];
     }
   }
@@ -56,20 +56,20 @@ uint32_t __attribute__((optimize("unroll-loops"))) convert(uint64_t * nn)
 
     if (!(x & bottommask)) {
       const size_t previous = (x >> halfstrip_size) & bottommask;
-      const uint8_t next = y >> (64 - halfstrip_size);
+      const uint32_t next = y >> (64 - halfstrip_size);
       if (next <= lookup[previous]) return steps - halfstrip_size - offset[previous];
     }
 
     if (!(y & topmask)) {
       const size_t previous = x & bottommask;
-      const uint8_t next = (y >> (64 - 2*halfstrip_size)) & bottommask;
+      const uint32_t next = (y >> (64 - 2*halfstrip_size)) & bottommask;
       if (next <= lookup[previous]) return steps - offset[previous];
     }
 
     for (uint32_t w2 = halfstrip_size; w2 < 64-halfstrip_size; w2 += halfstrip_size) {
       if (!(y & (topmask >> w2))) {
         const size_t previous = (y >> (64 - halfstrip_size - w2 + halfstrip_size)) & bottommask;
-        const uint8_t next =    (y >> (64 - halfstrip_size - w2 - halfstrip_size)) & bottommask;
+        const uint32_t next =    (y >> (64 - halfstrip_size - w2 - halfstrip_size)) & bottommask;
         if (next <= lookup[previous]) return steps + w2 - offset[previous];
       }
     }

src/ddlog.h

@@ -6,8 +6,8 @@
 #define strip_size 16
 #define halfstrip_size ((strip_size)/2)
 
-extern uint8_t lookup[256];
-extern uint8_t offset[256];
+extern uint32_t lookup[256];
+extern uint32_t offset[256];
 
 
 uint32_t convert(uint64_t *nn);

src/group.c

@@ -26,3 +26,19 @@ void group_clear()
 {
   mpz_clears(p, q, NULL);
 }
+
+void powmp_ui(mpz_t x, const mpz_t base, uint64_t exp)
+{
+  mpz_t y;
+  mpz_init_set_ui(y, 1);
+  mpz_set(x, base);
+
+  while (exp > 1) {
+    if (exp &  1) mpz_mul_modp(y, x, y);
+    mpz_mul_modp(x, x, x);
+    exp >>= 1;
+  }
+
+  mpz_mul_modp(x, x, y);
+  mpz_clear(y);
+}

src/group.h

@@ -13,3 +13,54 @@ extern const uint64_t gg;
 
 void group_init();
 void group_clear();
+
+
+
+/* some gmp internal funcitons to speed up modulus… */
+
+#define SIZ(x) ((x)->_mp_size)
+#define PTR(x) ((x)->_mp_d)
+#define MPN_NORMALIZE(DST, NLIMBS)				\
+  do {									\
+    while (1)								\
+      {									\
+	if ((DST)[(NLIMBS) - 1] != 0)					\
+	  break;							\
+	(NLIMBS)--;							\
+      }									\
+  } while (0)
+
+
+
+static inline
+void remp(mpz_t rop)
+{
+  int32_t limbs = SIZ(rop) - 24;
+
+  while (limbs > 0) {
+    /* note: this declarations MUST happen after checking
+     * for positivity of limbs. */
+    uint64_t a[limbs+1];
+    /* copy the most significant part of rop into a,
+     * then set it to zero */
+    for (int i = 24; i < SIZ(rop); i++) {
+      a[i-24] = PTR(rop)[i];
+      PTR(rop)[i] = 0;
+    }
+    a[limbs] = 0;
+
+    mpn_addmul_1(PTR(rop), a, limbs+1, gg);
+    MPN_NORMALIZE(PTR(rop), SIZ(rop));
+    limbs = SIZ(rop) - 24;
+  }
+
+  if (mpn_cmp(PTR(rop), PTR(p), SIZ(rop)) >= 0) {
+    mpn_sub(PTR(rop), PTR(rop), 24, PTR(p), 24);
+    MPN_NORMALIZE(PTR(rop), SIZ(rop));
+  }
+
+}
+
+void powmp_ui(mpz_t rop, const mpz_t base, uint64_t exp);
+
+#define mpz_mul_modp(rop, op1, op2) mpz_mul(rop, op1, op2); remp(rop);

src/rms.c

@@ -15,32 +15,6 @@
 #include "rms.h"
 #include "timeit.h"
 
-static inline
-void remp(mpz_t rop)
-{
-  int limbs = rop->_mp_size - 24;
-
-  if (limbs < 0) return;
-  else if (limbs == 0 && mpz_cmp(rop, p) < 0) return;
-  else if (limbs == 0 && mpz_cmp(rop, p) >= 0) {
-    mpz_sub(rop, rop, p);
-    return;
-  }
-
-  /* XXX. this can be replaced to mpn_addmul1 */
-  mp_limb_t a[limbs + 1];
-  mpn_mul(a, rop->_mp_d + 24, limbs, &gg, 1);
-  for (int i = 24; i < rop->_mp_size; i++) rop->_mp_d[i] = 0;
-  mpn_add(rop->_mp_d, rop->_mp_d, 24, a, limbs+1);
-
-  /* XXX. these two lines are useless and have the sole purpose to restore
-   * the mpz_t to a consistent type. */
-  mpz_add_ui(rop, rop, 1);
-  mpz_sub_ui(rop, rop, 1);
-
-  remp(rop);
-}
-
 static inline
 void fbpowm(mpz_t rop, const mpz_t T[4][256], const uint32_t exp)
 {
@@ -56,29 +30,28 @@ void fbpowm(mpz_t rop, const mpz_t T[4][256], const uint32_t exp)
 
 
 static inline
-uint32_t mul_single(mpz_t op1,
-                      mpz_t op2,
-                      const elgamal_cipher_t c,
-                      const uint32_t x,
-                      const mpz_t cx)
+uint32_t mul_single(const elgamal_cipher_t c,
+                    const uint32_t x,
+                    const mpz_t cx)
 {
+  mpz_t op1, op2;
+  mpz_inits(op1, op2, NULL);
   //mpz_powm(op1, c1, cx, p);
   /* first block */
-  mpz_powm_ui(op1, c->c1, cx->_mp_d[0], p);
+  powmp_ui(op1, c->c1, cx->_mp_d[0]);
   /* second block */
-  mpz_powm_ui(op2, c->c1e64, cx->_mp_d[1], p);
-  mpz_mul(op1, op2, op1);
-  remp(op1);
+  powmp_ui(op2, c->c1e64, cx->_mp_d[1]);
+  mpz_mul_modp(op1, op2, op1);
   /* third block */
-  mpz_powm_ui(op2, c->c1e128, cx->_mp_d[2], p);
-  mpz_mul(op1, op2, op1);
-  remp(op1);
+  powmp_ui(op2, c->c1e128, cx->_mp_d[2]);
+  mpz_mul_modp(op1, op2, op1);
 
-  mpz_powm_ui(op2, c->c2, x, p);
-  mpz_mul(op2, op2, op1);
-  remp(op2);
+  powmp_ui(op2, c->c2, x);
+  mpz_mul_modp(op2, op2, op1);
 
-  return convert(op2->_mp_d);
+  const uint32_t converted = convert(op2->_mp_d);
+  mpz_clears(op1, op2, NULL);
+  return converted;
 }
 
 void hss_mul(ssl2_t rop, const ssl1_t sl1, const ssl2_t sl2)
@@ -87,16 +60,14 @@ void hss_mul(ssl2_t rop, const ssl1_t sl1, const ssl2_t sl2)
   mpz_t op1, op2;
   mpz_inits(op1, op2, NULL);
 
-  rop->x = mul_single(op1, op2, sl1->w, sl2->x, sl2->cx);
+  rop->x = mul_single(sl1->w, sl2->x, sl2->cx);
 
   mpz_set_ui(rop->cx, 0);
   for (size_t t = 0; t < 160; t++) {
     mpz_mul_2exp(rop->cx, rop->cx, 1);
-    converted = mul_single(op1, op2, sl1->cw[t], sl2->x, sl2->cx);
+    converted = mul_single(sl1->cw[t], sl2->x, sl2->cx);
     mpz_add_ui(rop->cx, rop->cx, converted);
   }
-
-  mpz_clears(op1, op2, NULL);
 }
 
 int main()