Reimplement zz_powm()

mpn.c

@@ -361,6 +361,49 @@ mpn_import (mp_ptr zp, mp_size_t *zsize, size_t count, int order,
     MPN_NORMALIZE (zp, *zsize);
 }
 
+extern void __gmpn_powm(mp_limb_t *rp, const mp_limb_t *bp, mp_size_t bn,
+                        const mp_limb_t *ep, mp_size_t en, const mp_limb_t *mp,
+                        mp_size_t n, mp_limb_t *tp);
+extern mp_size_t __gmpn_binvert_itch(mp_size_t n);
+extern void __gmpn_powlo(mp_limb_t *rp, const mp_limb_t *bp,
+                         const mp_limb_t *ep, mp_size_t en, mp_size_t n,
+                         mp_limb_t *tp);
+extern void __gmpn_binvert(mp_limb_t *rp, const mp_limb_t *up, mp_size_t n,
+                           mp_limb_t *scratch);
+extern void __gmpn_mullo_n(mp_limb_t *rp, const mp_limb_t *xp,
+                           const mp_limb_t *yp, mp_size_t n);
+
+void mpn_powm(mp_limb_t *rp, const mp_limb_t *bp, mp_size_t bn,
+              const mp_limb_t *ep, mp_size_t en, const mp_limb_t *mp,
+              mp_size_t n, mp_limb_t *tp)
+{
+    __gmpn_powm(rp, bp, bn, ep, en, mp, n, tp);
+}
+
+mp_size_t mpn_binvert_itch(mp_size_t n)
+{
+    return __gmpn_binvert_itch(n);
+}
+
+void mpn_powlo(mp_limb_t *rp, const mp_limb_t *bp,
+               const mp_limb_t *ep, mp_size_t en, mp_size_t n,
+               mp_limb_t *tp)
+{
+    __gmpn_powlo(rp, bp, ep, en, n, tp);
+}
+
+void mpn_binvert(mp_limb_t *rp, const mp_limb_t *up, mp_size_t n,
+                 mp_limb_t *scratch)
+{
+    __gmpn_binvert(rp, up, n, scratch);
+}
+
+void mpn_mullo_n(mp_limb_t *rp, const mp_limb_t *xp,
+                 const mp_limb_t *yp, mp_size_t n)
+{
+    __gmpn_mullo_n(rp, xp, yp, n);
+}
+
 #ifdef __GNUC__
 #  pragma GCC diagnostic pop
 #endif

mpn.h

@@ -23,4 +23,25 @@ void * mpn_export(void *data, size_t *countp, int order,
 void mpn_import(mp_ptr zp, mp_size_t *zsize, size_t count, int order,
                 size_t size, int endian, size_t nail, const void *data);
 
+/* Compute r = b^e mod m.  Requires that m is odd and e > 1.
+   Uses scratch space at tp of MAX(mpn_binvert_itch(n), 2n) limbs. */
+void mpn_powm(mp_limb_t *rp, const mp_limb_t *bp, mp_size_t bn,
+              const mp_limb_t *ep, mp_size_t en, const mp_limb_t *mp,
+              mp_size_t n, mp_limb_t *tp);
+mp_size_t mpn_binvert_itch(mp_size_t n);
+
+/* Compute r = b^e mod B^n, B is the limb base.
+   Requires normalized e.  Uses scratch space of 3n words in tp. */
+void mpn_powlo(mp_limb_t *rp, const mp_limb_t *bp,
+               const mp_limb_t *ep, mp_size_t en, mp_size_t n,
+               mp_limb_t *tp);
+
+/* Compute r = u^(-1) mod B^n, B is the limb base. */
+void mpn_binvert(mp_limb_t *rp, const mp_limb_t *up, mp_size_t n,
+                 mp_limb_t *scratch);
+
+/* Multiply two n-limb numbers and return the low n limbs of their products. */
+void mpn_mullo_n(mp_limb_t *rp, const mp_limb_t *xp,
+                 const mp_limb_t *yp, mp_size_t n);
+
 #endif /* MPN_H */

tests/t-pow.c

@@ -23,12 +23,16 @@ zz_ref_powm(const zz_t *u, const zz_t *v, const zz_t *w, zz_t *r)
         return ZZ_MEM;
     }
     mpz_init(z);
+    mpz_abs(mw, mw);
     mpz_powm(z, mu, mv, mw);
     if (zz_set_mpz_t(z, r)) {
         mpz_clear(z);
         return ZZ_MEM;
     }
     mpz_clear(z);
+    if (zz_isneg(w) && !zz_iszero(r) && !zz_isneg(r)) {
+        return zz_add(w, r, r);
+    }
     return ZZ_OK;
 }
 
@@ -62,10 +66,10 @@ check_powm_bulk(void)
         if (zz_init(&u) || zz_random(bs, true, &u)) {
             abort();
         }
-        if (zz_init(&v) || zz_random(32, true, &v)) {
+        if (zz_init(&v) || zz_random(128, true, &v)) {
             abort();
         }
-        if (zz_init(&w) || zz_random(bs, false, &w)) {
+        if (zz_init(&w) || zz_random(bs, true, &w)) {
             abort();
         }
         if (zz_init(&z)) {
@@ -82,6 +86,12 @@ check_powm_bulk(void)
             zz_clear(&r);
         }
         else {
+            if (zz_iszero(&w)) {
+                if (zz_powm(&u, &v, &w, &z) != ZZ_VAL) {
+                    abort();
+                }
+                goto clear;
+            }
             if (zz_ref_gcd(&u, &w, &z) || zz_cmp(&z, 1) == ZZ_EQ) {
                 abort();
             }
@@ -131,6 +141,7 @@ check_powm_bulk(void)
                 abort();
             }
         }
+clear:
         zz_clear(&u);
         zz_clear(&v);
         zz_clear(&w);
@@ -168,6 +179,31 @@ check_powm_examples(void)
     if (zz_set(0, &w) || zz_powm(&u, &v, &w, &w) != ZZ_VAL) {
         abort();
     }
+    if (zz_set(123, &u) || zz_set(321, &v) || zz_set(1, &w)
+        || zz_powm(&u, &v, &w, &w) || zz_cmp(&w, 0))
+    {
+        abort();
+    }
+    if (zz_set(123, &u) || zz_set(0, &v) || zz_set(321, &w)
+        || zz_powm(&u, &v, &w, &w) || zz_cmp(&w, 1))
+    {
+        abort();
+    }
+    if (zz_set(0, &u) || zz_set(321, &v) || zz_set(123, &w)
+        || zz_powm(&u, &v, &w, &w) || zz_cmp(&w, 0))
+    {
+        abort();
+    }
+    if (zz_set(1, &u) || zz_set(321, &v) || zz_set(123, &w)
+        || zz_powm(&u, &v, &w, &w) || zz_cmp(&w, 1))
+    {
+        abort();
+    }
+    if (zz_set(321, &u) || zz_set(1, &v) || zz_set(123, &w)
+        || zz_powm(&u, &v, &w, &w) || zz_cmp(&w, 75))
+    {
+        abort();
+    }
     zz_clear(&u);
     zz_clear(&v);
     zz_clear(&w);