Addition of mp_radix_size_overestimate (behaviour like the function in GMP)

bn_mp_fwrite.c

@@ -8,12 +8,13 @@ mp_err mp_fwrite(const mp_int *a, int radix, FILE *stream)
 {
    char *buf;
    mp_err err;
-   int len;
+   int len = 0;
    size_t written;
 
-   /* TODO: this function is not in this PR */
    if (MP_HAS(MP_RADIX_SIZE_OVERESTIMATE)) {
-      /* if ((err = mp_radix_size_overestimate(&t, base, &len)) != MP_OKAY)      goto LBL_ERR; */
+      if ((err = mp_radix_size_overestimate(a, radix, &len)) != MP_OKAY) {
+         return err;
+      }
    } else {
       if ((err = mp_radix_size(a, radix, &len)) != MP_OKAY) {
          return err;
@@ -33,8 +34,8 @@ mp_err mp_fwrite(const mp_int *a, int radix, FILE *stream)
       err = MP_ERR;
       goto LBL_ERR;
    }
-   err = MP_OKAY;
 
+   err = MP_OKAY;
 
 LBL_ERR:
    MP_FREE_BUFFER(buf, (size_t)len);

bn_mp_radix_size.c

@@ -6,12 +6,14 @@
 /* returns size of ASCII representation */
 mp_err mp_radix_size(const mp_int *a, int radix, int *size)
 {
-   mp_err  err;
-   int digs;
-   mp_int   t;
-   mp_digit d;
+   mp_err err;
+   mp_int a_, b;
 
-   *size = 0;
+   if (size == NULL) {
+      return MP_VAL;
+   } else {
+      *size = 0;
+   }
 
    /* make sure the radix is in range */
    if ((radix < 2) || (radix > 64)) {
@@ -23,42 +25,38 @@ mp_err mp_radix_size(const mp_int *a, int radix, int *size)
       return MP_OKAY;
    }
 
-   /* special case for binary */
-   if (radix == 2) {
-      *size = (mp_count_bits(a) + ((a->sign == MP_NEG) ? 1 : 0) + 1);
+   /* A small shortcut for powers of two. */
+   if (!(radix&(radix-1))) {
+      unsigned int x = (unsigned int)radix;
+      int y, bit_count, rem;
+      for (y=0; (y < 7) && !(x & 1u); y++) {
+         x >>= 1;
+      }
+      bit_count = mp_count_bits(a);
+      *size = bit_count/y;
+      rem = bit_count - ((*size) * y);
+      /* Add 1 for the remainder if any and 1 for "\0". */
+      *size += ((rem == 0) ? 1 : 2) + (a->sign == MP_NEG);
       return MP_OKAY;
    }
 
-   /* digs is the digit count */
-   digs = 0;
-
-   /* if it's negative add one for the sign */
-   if (a->sign == MP_NEG) {
-      ++digs;
+   if ((err = mp_init(&b)) != MP_OKAY) {
+      goto LBL_ERR;
    }
 
-   /* init a copy of the input */
-   if ((err = mp_init_copy(&t, a)) != MP_OKAY) {
-      return err;
+   a_ = *a;
+   a_.sign = MP_ZPOS;
+   if ((err = mp_ilogb(&a_, (uint32_t)radix, &b)) != MP_OKAY) {
+      goto LBL_ERR;
    }
 
-   /* force temp to positive */
-   t.sign = MP_ZPOS;
-
-   /* fetch out all of the digits */
-   while (!MP_IS_ZERO(&t)) {
-      if ((err = mp_div_d(&t, (mp_digit)radix, &t, &d)) != MP_OKAY) {
-         goto LBL_ERR;
-      }
-      ++digs;
-   }
+   *size = (int)mp_get_l(&b);
 
-   /* return digs + 1, the 1 is for the NULL byte that would be required. */
-   *size = digs + 1;
-   err = MP_OKAY;
+   /* mp_ilogb truncates to zero, hence we need one extra put on top and one for `\0`. */
+   *size += 2 + (a->sign == MP_NEG);
 
 LBL_ERR:
-   mp_clear(&t);
+   mp_clear(&b);
    return err;
 }
 

bn_mp_radix_size_overestimate.c

@@ -0,0 +1,70 @@
+#include "tommath_private.h"
+#ifdef BN_MP_RADIX_SIZE_OVERESTIMATE_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+
+/*
+   Table of {0, log_2([1..64])} times 2^p where p is the scale
+   factor defined in LTM_RADIX_SIZE_SCALE.
+ */
+/* *INDENT-OFF* */
+#define LTM_RADIX_SIZE_SCALE 13
+static const uint16_t logbases[65] = {
+       0u,     0u,  8192u, 12984u, 16384u,
+   19021u, 21176u, 22997u, 24576u, 25968u,
+   27213u, 28339u, 29368u, 30314u, 31189u,
+   32005u, 32768u, 33484u, 34160u, 34799u,
+   35405u, 35981u, 36531u, 37057u, 37560u,
+   38042u, 38506u, 38952u, 39381u, 39796u,
+   40197u, 40584u, 40960u, 41323u, 41676u,
+   42019u, 42352u, 42675u, 42991u, 43298u,
+   43597u, 43889u, 44173u, 44451u, 44723u,
+   44989u, 45249u, 45503u, 45752u, 45995u,
+   46234u, 46468u, 46698u, 46923u, 47144u,
+   47360u, 47573u, 47783u, 47988u, 48190u,
+   48389u, 48584u, 48776u, 48965u, 49152u
+};
+/* *INDENT-ON* */
+mp_err mp_radix_size_overestimate(const mp_int *a, const int radix, int *size)
+{
+   mp_int bi_bit_count, bi_k;
+   int bit_count;
+   mp_err err = MP_OKAY;
+
+   *size = 0;
+
+   if ((radix < 2) || (radix > 64)) {
+      return MP_VAL;
+   }
+
+   bit_count = mp_count_bits(a) + 1;
+
+   if (bit_count == 0) {
+      *size = 2;
+      return MP_OKAY;
+   }
+
+   if ((err = mp_init_multi(&bi_bit_count, &bi_k, NULL)) != MP_OKAY) {
+      return err;
+   }
+
+   mp_set_l(&bi_bit_count, bit_count);
+   mp_set_u32(&bi_k, logbases[radix]);
+   if ((err = mp_mul_2d(&bi_bit_count, LTM_RADIX_SIZE_SCALE, &bi_bit_count)) != MP_OKAY)  goto LTM_ERR;
+   if ((err = mp_div(&bi_bit_count, &bi_k, &bi_bit_count, NULL)) != MP_OKAY)              goto LTM_ERR;
+
+   *size = (int)mp_get_l(&bi_bit_count) + 4;
+#if ( (defined MP_8BIT) && (INT_MAX > 0xFFFF))
+   *size += 3;
+#endif
+
+LTM_ERR:
+   mp_clear_multi(&bi_bit_count, &bi_k, NULL);
+   return err;
+}
+
+
+
+
+#endif

demo/test.c

@@ -2296,6 +2296,114 @@ static int test_s_mp_toom_sqr(void)
    mp_clear_multi(&a, &b, &c, NULL);
    return EXIT_FAILURE;
 }
+static int test_mp_radix_size(void)
+{
+   mp_err err;
+   mp_int a;
+   int radix, size;
+/* *INDENT-OFF* */
+   int results[65] = {
+       0, 0, 1627, 1027, 814, 702, 630, 581, 543,
+       514, 491, 471, 455, 441, 428, 418, 408, 399,
+       391, 384, 378, 372, 366, 361, 356, 352, 347,
+       343, 340, 336, 333, 330, 327, 324, 321, 318,
+       316, 314, 311, 309, 307, 305, 303, 301, 299,
+       298, 296, 294, 293, 291, 290, 288, 287, 285,
+       284, 283, 281, 280, 279, 278, 277, 276, 275,
+       273, 272
+   };
+/* *INDENT-ON* */
+
+   mp_init(&a);
+
+   /* number to result in a different size for every base: 67^(4 * 67) */
+   mp_set(&a, 67);
+   if ((err = mp_expt_u32(&a, 268u, &a)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+
+   for (radix = 2; radix < 65; radix++) {
+      if ((err = mp_radix_size(&a, radix, &size)) != MP_OKAY) {
+         goto LTM_ERR;
+      }
+      if (size != results[radix]) {
+         fprintf(stderr, "mp_radix_size: result for base %d was %d instead of %d\n",
+                 radix, size, results[radix]);
+         goto LTM_ERR;
+      }
+      a.sign = MP_NEG;
+      if ((err = mp_radix_size(&a, radix, &size)) != MP_OKAY) {
+         goto LTM_ERR;
+      }
+      if (size != (results[radix] + 1)) {
+         fprintf(stderr, "mp_radix_size: result for base %d was %d instead of %d\n",
+                 radix, size, results[radix]);
+         goto LTM_ERR;
+      }
+      a.sign = MP_ZPOS;
+   }
+
+   mp_clear(&a);
+   return EXIT_SUCCESS;
+LTM_ERR:
+   mp_clear(&a);
+   return EXIT_FAILURE;
+}
+
+static int test_mp_radix_size_overestimate(void)
+{
+
+   mp_err err;
+   mp_int a;
+   int radix, size;
+/* *INDENT-OFF* */
+   int results[65] = {
+       0, 0, 1627, 1027, 814, 702, 630, 581, 543,
+       514, 491, 471, 455, 441, 428, 418, 408, 399,
+       391, 384, 378, 372, 366, 361, 356, 352, 347,
+       343, 340, 336, 333, 330, 327, 324, 321, 318,
+       316, 314, 311, 309, 307, 305, 303, 301, 299,
+       298, 296, 294, 293, 291, 290, 288, 287, 285,
+       284, 283, 281, 280, 279, 278, 277, 276, 275,
+       273, 272
+   };
+/* *INDENT-ON* */
+
+   mp_init(&a);
+
+   /* number to result in a different size for every base: 67^(4 * 67) */
+   mp_set(&a, 67);
+   if ((err = mp_expt_u32(&a, 268u, &a)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+
+   for (radix = 2; radix < 65; radix++) {
+      if ((err = mp_radix_size_overestimate(&a, radix, &size)) != MP_OKAY) {
+         goto LTM_ERR;
+      }
+      if (size < results[radix]) {
+         fprintf(stderr, "mp_radix_size_overestimate: result for base %d was %d instead of %d\n",
+                 radix, size, results[radix]);
+         goto LTM_ERR;
+      }
+      a.sign = MP_NEG;
+      if ((err = mp_radix_size_overestimate(&a, radix, &size)) != MP_OKAY) {
+         goto LTM_ERR;
+      }
+      if (size < results[radix]) {
+         fprintf(stderr, "mp_radix_size_overestimate: result for base %d was %d instead of %d\n",
+                 radix, size, results[radix]);
+         goto LTM_ERR;
+      }
+      a.sign = MP_ZPOS;
+   }
+
+   mp_clear(&a);
+   return EXIT_SUCCESS;
+LTM_ERR:
+   mp_clear(&a);
+   return EXIT_FAILURE;
+}
 
 static int test_mp_read_write_ubin(void)
 {
@@ -2426,6 +2534,7 @@ static int unit_tests(int argc, char **argv)
       const char *name;
       int (*fn)(void);
    } test[] = {
+
 #define T0(n)           { #n, test_##n }
 #define T1(n, o)        { #n, MP_HAS(o) ? test_##n : NULL }
 #define T2(n, o1, o2)   { #n, MP_HAS(o1) && MP_HAS(o2) ? test_##n : NULL }
@@ -2461,6 +2570,8 @@ static int unit_tests(int argc, char **argv)
       T1(mp_read_write_sbin, MP_TO_SBIN),
       T1(mp_reduce_2k, MP_REDUCE_2K),
       T1(mp_reduce_2k_l, MP_REDUCE_2K_L),
+      T1(mp_radix_size, MP_RADIX_SIZE),
+      T1(mp_radix_size_overestimate, MP_RADIX_SIZE_OVERESTIMATE),
 #if defined(__STDC_IEC_559__) || defined(__GCC_IEC_559)
       T1(mp_set_double, MP_SET_DOUBLE),
 #endif

doc/bn.tex

@@ -2118,8 +2118,16 @@ \subsection{To ASCII}
 \begin{alltt}
 int mp_radix_size (mp_int * a, int radix, int *size)
 \end{alltt}
-This stores in ``size'' the number of characters (including space for the NUL terminator) required.  Upon error this
-function returns an error code and ``size'' will be zero.
+This stores in ``size'' the exact number of characters (including space for the NUL terminator) required.  Upon error this function returns an error code and ``size'' will be zero. It uses \texttt{mp\_ilogb} to compute the size in O(log n) with quite a large overhead. This function is deprecated in favor of
+\texttt{mp\_radix\_sizeinbase} described below. If the exact size is still required after the removal of this function use \texttt{mp\_ilogb} directly.
+
+\index{mp\_radix\_overestimate}
+\begin{alltt}
+int mp_radix_size_overestimate (mp_int * a, int radix, int *size)
+\end{alltt}
+This stores in ``size'' the number of characters (including space for the NUL terminator) required.  Upon error this function returns an error code and ``size'' will be zero. This function can overshoot by one (two in case of \texttt{MP\_8BIT} and also for \texttt{MP\_16BIT} if the type \texttt{int} has more than 16 bits). The results for bases that are powers-of-two are exact. It uses tables to compute the result in O(1) with quite a small overhead.
+
+The behavior is similar to the behavior of the function \texttt{mpn\_sizeinbase} in GMP for \texttt{MP\_32BIT}, and \texttt{MP\_64BIT}.
 
 If \texttt{LTM\_NO\_FILE} is not defined a function to write to a file is also available.
 \index{mp\_fwrite}

libtommath_VS2008.vcproj

@@ -684,6 +684,10 @@
 			RelativePath="bn_mp_radix_size.c"
 			>
 		</File>
+		<File
+			RelativePath="bn_mp_radix_size_overestimate.c"
+			>
+		</File>
 		<File
 			RelativePath="bn_mp_radix_smap.c"
 			>

makefile

@@ -42,19 +42,19 @@ bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o
 bn_mp_or.o bn_mp_pack.o bn_mp_pack_count.o bn_mp_prime_fermat.o bn_mp_prime_frobenius_underwood.o \
 bn_mp_prime_is_prime.o bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o \
 bn_mp_prime_rabin_miller_trials.o bn_mp_prime_rand.o bn_mp_prime_strong_lucas_selfridge.o \
-bn_mp_radix_size.o bn_mp_radix_smap.o bn_mp_rand.o bn_mp_read_radix.o bn_mp_reduce.o bn_mp_reduce_2k.o \
-bn_mp_reduce_2k_l.o bn_mp_reduce_2k_setup.o bn_mp_reduce_2k_setup_l.o bn_mp_reduce_is_2k.o \
-bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_root_u32.o bn_mp_rshd.o bn_mp_sbin_size.o bn_mp_set.o \
-bn_mp_set_double.o bn_mp_set_i32.o bn_mp_set_i64.o bn_mp_set_l.o bn_mp_set_ll.o bn_mp_set_u32.o \
-bn_mp_set_u64.o bn_mp_set_ul.o bn_mp_set_ull.o bn_mp_shrink.o bn_mp_signed_rsh.o bn_mp_sqr.o \
-bn_mp_sqrmod.o bn_mp_sqrt.o bn_mp_sqrtmod_prime.o bn_mp_sub.o bn_mp_sub_d.o bn_mp_submod.o \
-bn_mp_to_radix.o bn_mp_to_sbin.o bn_mp_to_ubin.o bn_mp_ubin_size.o bn_mp_unpack.o bn_mp_xor.o bn_mp_zero.o \
-bn_prime_tab.o bn_s_mp_add.o bn_s_mp_balance_mul.o bn_s_mp_exptmod.o bn_s_mp_exptmod_fast.o \
-bn_s_mp_get_bit.o bn_s_mp_invmod_fast.o bn_s_mp_invmod_slow.o bn_s_mp_karatsuba_mul.o \
-bn_s_mp_karatsuba_sqr.o bn_s_mp_montgomery_reduce_fast.o bn_s_mp_mul_digs.o bn_s_mp_mul_digs_fast.o \
-bn_s_mp_mul_high_digs.o bn_s_mp_mul_high_digs_fast.o bn_s_mp_prime_is_divisible.o \
-bn_s_mp_rand_jenkins.o bn_s_mp_rand_platform.o bn_s_mp_reverse.o bn_s_mp_sqr.o bn_s_mp_sqr_fast.o \
-bn_s_mp_sub.o bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o
+bn_mp_radix_size.o bn_mp_radix_size_overestimate.o bn_mp_radix_smap.o bn_mp_rand.o bn_mp_read_radix.o \
+bn_mp_reduce.o bn_mp_reduce_2k.o bn_mp_reduce_2k_l.o bn_mp_reduce_2k_setup.o bn_mp_reduce_2k_setup_l.o \
+bn_mp_reduce_is_2k.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_root_u32.o bn_mp_rshd.o \
+bn_mp_sbin_size.o bn_mp_set.o bn_mp_set_double.o bn_mp_set_i32.o bn_mp_set_i64.o bn_mp_set_l.o \
+bn_mp_set_ll.o bn_mp_set_u32.o bn_mp_set_u64.o bn_mp_set_ul.o bn_mp_set_ull.o bn_mp_shrink.o \
+bn_mp_signed_rsh.o bn_mp_sqr.o bn_mp_sqrmod.o bn_mp_sqrt.o bn_mp_sqrtmod_prime.o bn_mp_sub.o bn_mp_sub_d.o \
+bn_mp_submod.o bn_mp_to_radix.o bn_mp_to_sbin.o bn_mp_to_ubin.o bn_mp_ubin_size.o bn_mp_unpack.o \
+bn_mp_xor.o bn_mp_zero.o bn_prime_tab.o bn_s_mp_add.o bn_s_mp_balance_mul.o bn_s_mp_exptmod.o \
+bn_s_mp_exptmod_fast.o bn_s_mp_get_bit.o bn_s_mp_invmod_fast.o bn_s_mp_invmod_slow.o \
+bn_s_mp_karatsuba_mul.o bn_s_mp_karatsuba_sqr.o bn_s_mp_montgomery_reduce_fast.o bn_s_mp_mul_digs.o \
+bn_s_mp_mul_digs_fast.o bn_s_mp_mul_high_digs.o bn_s_mp_mul_high_digs_fast.o \
+bn_s_mp_prime_is_divisible.o bn_s_mp_rand_jenkins.o bn_s_mp_rand_platform.o bn_s_mp_reverse.o \
+bn_s_mp_sqr.o bn_s_mp_sqr_fast.o bn_s_mp_sub.o bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o
 
 #END_INS