mpi: use integer math for radix length.

* mpi/logtab.h: Regenerated.
(s_logv_2): Now table of scaled integers.
(MP_LOG_SCALE): New constant.

* mpi/make-logtab.txr (scale, type): New variables.
Generate integer table with log2 values scaled by the scale
factor, rounded up.

* mpi/config.h (MP_LOGTAB): Removed. We always use the table.

* mpi/mpi.c: Unconditionally include logtab.h.
(LOG_V_2): Macro removed.
(s_mp_outlen): Rewritten using scaled integer math.
Overflow is avoided by splitting the input into a part that is
an exact multiple of the scale factor, and a remaining
part. Only the remaining part need be multiplied by a value
from the table before dividing by the scale factor.

4 files changed, 42 insertions, 59 deletions

diff --git a/mpi/logtab.h b/mpi/logtab.h
index a4b2d04f..c6f468bf 100644
--- a/mpi/logtab.h
+++ b/mpi/logtab.h
@@ -1,35 +1,29 @@
 /*
- * A table of the logs of 2 for various bases (the 0 and 1 entries of
- * this table are meaningless and should not be referenced).
+ * A table of the logs of 2 for various bases scaled by a factor
+ * and converted to integer.
  *
- * This table is used to compute output lengths for the mp_toradix()
- * function.  Since a number n in radix r takes up about log_r(n)
- * digits, we estimate the output size by taking the least integer
- * greater than log_r(n), where:
- *
- * log_r(n) = log_2(n) * log_r(2)
- *
- * This table, therefore, is a table of log_r(2) for 2 <= r <= 36,
- * which are the output bases supported.
+ * This table is used to compute output lengths for the mp_toradix
+ * function.
  */
 
-const double s_logv_2[] = {
-  0.000000000, 0.000000000, 1.000000000, 0.630929754,
-  0.500000000, 0.430676558, 0.386852807, 0.356207187,
-  0.333333333, 0.315464877, 0.301029996, 0.289064826,
-  0.278942946, 0.270238154, 0.262649535, 0.255958025,
-  0.250000000, 0.244650542, 0.239812467, 0.235408913,
-  0.231378213, 0.227670249, 0.224243824, 0.221064729,
-  0.218104292, 0.215338279, 0.212746054, 0.210309918,
-  0.208014598, 0.205846832, 0.203795047, 0.201849087,
-  0.200000000, 0.198239863, 0.196561632, 0.194959022,
-  0.193426404, 0.191958720, 0.190551412, 0.189200360,
-  0.187901825, 0.186652411, 0.185449023, 0.184288833,
-  0.183169251, 0.182087900, 0.181042597, 0.180031327,
-  0.179052232, 0.178103594, 0.177183820, 0.176291434,
-  0.175425064, 0.174583430, 0.173765343, 0.172969690,
-  0.172195434, 0.171441601, 0.170707280, 0.169991616,
-  0.169293808, 0.168613099, 0.167948779, 0.167300179,
-  0.166666667
-};
+#define MP_LOG_SCALE 16384
 
+const unsigned int s_logv_2[] = {
+  0, 0, 16384, 10338,
+  8192, 7057, 6339, 5837,
+  5462, 5169, 4933, 4737,
+  4571, 4428, 4304, 4194,
+  4096, 4009, 3930, 3857,
+  3791, 3731, 3675, 3622,
+  3574, 3529, 3486, 3446,
+  3409, 3373, 3339, 3308,
+  3277, 3248, 3221, 3195,
+  3170, 3146, 3122, 3100,
+  3079, 3059, 3039, 3020,
+  3002, 2984, 2967, 2950,
+  2934, 2919, 2903, 2889,
+  2875, 2861, 2847, 2834,
+  2822, 2809, 2797, 2786,
+  2774, 2763, 2752, 2742,
+  2731
+};
diff --git a/mpi/make-logtab.txr b/mpi/make-logtab.txr
index 5566dd66..90e36f38 100644
--- a/mpi/make-logtab.txr
+++ b/mpi/make-logtab.txr
@@ -1,22 +1,21 @@
-@(bind logs @(tuples 4 ^(0.0 0.0 ,*(take 63 (mapcar* [chain log2 /] (range 2))))))
-@(do (set *pprint-flo-format* "~0,9f"))
+@(bind scale 16384)
+@(bind type "unsigned int")
+@(bind logs @(tuples 4 ^(0.0 0.0
+                         ,*(take 63
+                                 (mapcar* [chain log2 / (op * scale) ceil toint]
+                                          (range 2))))))
 @(output)
 /*
- * A table of the logs of 2 for various bases (the 0 and 1 entries of
- * this table are meaningless and should not be referenced).
+ * A table of the logs of 2 for various bases scaled by a factor
+ * and converted to integer.
  *
- * This table is used to compute output lengths for the mp_toradix()
- * function.  Since a number n in radix r takes up about log_r(n)
- * digits, we estimate the output size by taking the least integer
- * greater than log_r(n), where:
- *
- * log_r(n) = log_2(n) * log_r(2)
- *
- * This table, therefore, is a table of log_r(2) for 2 <= r <= 36,
- * which are the output bases supported.
+ * This table is used to compute output lengths for the mp_toradix
+ * function.
  */
 
-const double s_logv_2[] = {
+#define MP_LOG_SCALE @scale
+
+const @type s_logv_2[] = {
 @  (repeat)
   @(rep)@logs, @(last)@logs,@(end)
 @  (last)
diff --git a/mpi/mpi-config.h b/mpi/mpi-config.h
index ed7a87b6..801eaae7 100644
--- a/mpi/mpi-config.h
+++ b/mpi/mpi-config.h
@@ -12,10 +12,6 @@
 #define MP_NUMTH 1 /* include number theoretic functions? */
 #endif
 
-#ifndef MP_LOGTAB
-#define MP_LOGTAB 1 /* use table of logs instead of log()? */
-#endif
-
 #ifndef MP_MEMSET
 #define MP_MEMSET 1 /* use memset() to zero buffers? */
 #endif
diff --git a/mpi/mpi.c b/mpi/mpi.c
index 579a43a6..9ece8047 100644
--- a/mpi/mpi.c
+++ b/mpi/mpi.c
@@ -43,17 +43,7 @@ extern mem_t *chk_calloc(size_t n, size_t size);
 #define DIAG(T,V)
 #endif
 
-/* If MP_LOGTAB is not defined, use the math library to compute the
- * logarithms on the fly.  Otherwise, use the static table below.
- * Pick which works best for your system.
- */
-#if MP_LOGTAB
 #include "logtab.h"
-#define LOG_V_2(R) s_logv_2[(R)]
-#else
-#include <math.h>
-#define LOG_V_2(R) (log(2.0)/log(R))
-#endif
 
 /* Default precision for newly created mp_int's */
 static mp_size s_mp_defprec = MP_DEFPREC;
@@ -4147,5 +4137,9 @@ char s_mp_todigit(int val, int r, int low)
  */
 size_t s_mp_outlen(mp_size bits, int r)
 {
-  return convert(size_t, convert(double, bits) * LOG_V_2(r) + 0.5);
+  mp_size units = bits / MP_LOG_SCALE;
+  mp_size rem = bits % MP_LOG_SCALE;
+  mp_size log2 = s_logv_2[r];
+
+  return convert(size_t, units * log2 + (rem * log2 + (MP_LOG_SCALE - 1)) / MP_LOG_SCALE);
 }