Refactor bigint functions and add remaining helpers

A second variant of the mpz-based iterator is added, and both functions
are moved into a new 'mpz' module, removing 'iterations' and the common
pszip.h header.

These iterators have been updated to have all (de)allocation performed
externally as to reduce overhead / boost performance.  Related helper
functions have also been added to this module.

Signed-off-by: Malfurious <m@lfurio.us>

4 files changed, 118 insertions, 45 deletions

diff --git a/iterations.c b/iterations.c
deleted file mode 100644
index 517bcaa..0000000
--- a/iterations.c
+++ /dev/null
@@ -1,40 +0,0 @@
-#include "pszip.h"
-
-/*
- * Bit hack source:
- * https://graphics.stanford.edu/~seander/bithacks.html#NextBitPermutation
- */
-
-void pszip_iterate_bithack01(mpz_t w, const mpz_t v)
-{
-    mpz_t vmo, t, tpo, ct, nct, act, actmo;
-    mp_limb_t *actmo_l;
-    size_t actmo_s;
-
-    mpz_inits(vmo, t, tpo, ct, nct, act, actmo, 0);
-
-    /* Gets v's least significant 0 bits set to 1.
-     * t = v | (v - 1); */
-    mpz_sub_ui(vmo, v, 1);
-    mpz_ior(t, v, vmo);
-
-    /* Next set to 1 the most significant bit to change,
-     * set to 0 the least significant ones, and add the necessary 1 bits.
-     * w = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(v) + 1));
-     * __builtin_ctz(v) -> mpz_scan1(v, 0) */
-    mpz_add_ui(tpo, t, 1);
-    mpz_com(ct, t);
-    mpz_neg(nct, ct);
-    mpz_and(act, ct, nct);
-    mpz_sub_ui(actmo, act, 1);
-
-    /* the bitshift requires obtaining direct memory access */
-    actmo_s = mpz_size(actmo);
-    actmo_l = mpz_limbs_modify(actmo, actmo_s);
-    mpn_rshift(actmo_l, actmo_l, actmo_s, (mpz_scan1(v, 0) + 1));
-    mpz_limbs_finish(actmo, actmo_s);
-
-    mpz_ior(w, tpo, actmo);
-
-    mpz_clears(vmo, t, tpo, ct, nct, act, actmo, 0);
-}
diff --git a/mpz.c b/mpz.c
new file mode 100644
index 0000000..c264a44
--- /dev/null
+++ b/mpz.c
@@ -0,0 +1,96 @@
+#include "mpz.h"
+#include <string.h>
+
+static size_t __ceil(double f)
+{
+    /* broken for negative values, which are outside of our need here */
+    return (size_t)f + (f != (size_t)f);
+}
+
+static size_t sizeinitmpz(mpz_t n, size_t size)
+{
+    size_t limbs = __ceil((double)size / sizeof(mp_limb_t));
+    size_t bytes = limbs * sizeof(mp_limb_t);
+
+    mpz_init(n);
+
+    mp_limb_t *ptr = mpz_limbs_write(n, limbs);
+    memset(ptr, 0, bytes);
+    mpz_limbs_finish(n, 0);
+
+    return limbs;
+}
+
+void pszip_mpz_state_init(struct pszip_mpz_state *mpzs, size_t popcnt, size_t size)
+{
+    sizeinitmpz(mpzs->current, size);
+    mpz_inits(mpzs->a, mpzs->b, mpzs->c, mpzs->d, 0);
+
+    mpz_ui_pow_ui(mpzs->current, 2, popcnt);
+    mpz_sub_ui(mpzs->current, mpzs->current, 1);
+}
+
+void pszip_mpz_state_clear(struct pszip_mpz_state *mpzs)
+{
+    mpz_clears(mpzs->current, mpzs->a, mpzs->b, mpzs->c, mpzs->d, 0);
+}
+
+void pszip_mpz_init_data(mpz_t n, const void *data, size_t size)
+{
+    size_t limbs = sizeinitmpz(n, size);
+
+    mp_limb_t *ptr = mpz_limbs_modify(n, limbs);
+    memcpy(ptr, data, size);
+    mpz_limbs_finish(n, limbs);
+}
+
+/*
+ * Bit hack source:
+ * https://graphics.stanford.edu/~seander/bithacks.html#NextBitPermutation
+ */
+
+void pszip_mpz_iterate_bithack01(struct pszip_mpz_state *mpzs)
+{
+    /* Gets current's least significant 0 bits set to 1
+     * a = cur | (cur - 1); */
+    mpz_sub_ui(mpzs->b, mpzs->current, 1);
+    mpz_ior(mpzs->a, mpzs->current, mpzs->b);
+
+    /* Next set to 1 the most significant bit to change,
+     * set to 0 the least significant ones, and add the necessary 1 bits.
+     * new = (a + 1) | (((~a & -~a) - 1) >> (__builtin_ctz(cur) + 1));
+     * __builtin_ctz(cur) -> mpz_scan1(cur, 0) */
+    mpz_add_ui(mpzs->b, mpzs->a, 1);    // (a + 1)
+    mpz_com(mpzs->c, mpzs->a);          // (~a)
+    mpz_neg(mpzs->d, mpzs->c);          // (-~a)
+    mpz_and(mpzs->c, mpzs->c, mpzs->d); // (~a & -~a)
+    mpz_sub_ui(mpzs->c, mpzs->c, 1);    // ((~a & -~a) - 1)
+
+    /* bitshift requires direct memory access */
+    size_t limbs = mpz_size(mpzs->c);
+    mp_limb_t *ptr = mpz_limbs_modify(mpzs->c, limbs);
+    mpn_rshift(ptr, ptr, limbs, (mpz_scan1(mpzs->current, 0) + 1));
+    mpz_limbs_finish(mpzs->c, limbs);
+
+    mpz_ior(mpzs->current, mpzs->b, mpzs->c);
+}
+
+void pszip_mpz_iterate_bithack02(struct pszip_mpz_state *mpzs)
+{
+    /* See comments from bithack01.
+     * a = (cur | (cur - 1)) + 1; */
+    mpz_sub_ui(mpzs->a, mpzs->current, 1);
+    mpz_ior(mpzs->a, mpzs->current, mpzs->a);
+    mpz_add_ui(mpzs->a, mpzs->a, 1);
+
+    /* new = a | ((((a & -a) / (cur & -cur)) >> 1) - 1);
+     * The right-shift by 1 is implemented as a divide by 2. */
+    mpz_neg(mpzs->b, mpzs->a);                // (-a)
+    mpz_and(mpzs->b, mpzs->a, mpzs->b);       // (a & -a)
+    mpz_neg(mpzs->c, mpzs->current);          // (-cur)
+    mpz_and(mpzs->c, mpzs->current, mpzs->c); // (cur & -cur)
+    mpz_mul_ui(mpzs->c, mpzs->c, 2);          // (... / (cur & -cur)) >> 1
+    mpz_tdiv_q(mpzs->b, mpzs->b, mpzs->c);    // (((a & -a) / (cur & -cur)) >> 1)
+    mpz_sub_ui(mpzs->b, mpzs->b, 1);          // (ior) right hand side
+    mpz_ior(mpzs->current, mpzs->a, mpzs->b);
+}
diff --git a/mpz.h b/mpz.h
new file mode 100644
index 0000000..ef41ec6
--- /dev/null
+++ b/mpz.h
@@ -0,0 +1,22 @@
+#pragma once
+
+#include <gmp.h>
+#include <stddef.h>
+
+/* This struct contains the instance data used by the iterators.  Allocation is
+ * performed at a single point to reduce run-time overhead.
+ *
+ * 'current' contains the iterator's current value.  The rest of the members are
+ * temporary variables.  See the appropriate memory management functions for
+ * this type. */
+struct pszip_mpz_state
+{
+    mpz_t current, a, b, c, d;
+};
+
+void pszip_mpz_state_init(struct pszip_mpz_state *mpzs, size_t popcnt, size_t size);
+void pszip_mpz_state_clear(struct pszip_mpz_state *mpzs);
+void pszip_mpz_init_data(mpz_t n, const void *data, size_t size);
+
+void pszip_mpz_iterate_bithack01(struct pszip_mpz_state *mpzs);
+void pszip_mpz_iterate_bithack02(struct pszip_mpz_state *mpzs);
diff --git a/pszip.h b/pszip.h
deleted file mode 100644
index 9d07e49..0000000
--- a/pszip.h
+++ /dev/null
@@ -1,5 +0,0 @@
-#pragma once
-
-#include <gmp.h>
-
-void pszip_iterate_bithack01(mpz_t, const mpz_t);