Add a new mp_todec_fast...

bn_mp_to_decimal.c

@@ -0,0 +1,23 @@
+#include "tommath_private.h"
+#ifdef BN_MP_TO_DECIMAL_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* stores a bignum as a decimal ASCII string, using Barrett
+ * reduction if available.
+ */
+
+mp_err mp_to_decimal(const mp_int *a, char *str, size_t maxlen)
+{
+   mp_err err;
+
+   if (MP_HAS(S_MP_TO_DECIMAL_FAST) && (a->used > 10)) {
+      err = s_mp_to_decimal_fast(a, str, maxlen);
+   } else {
+      err = mp_to_radix(a, str, maxlen, 10);
+   }
+
+   return err;
+}
+
+#endif

bn_s_mp_to_decimal_fast.c

@@ -0,0 +1,246 @@
+#include "tommath_private.h"
+#include <string.h>
+#ifdef BN_S_MP_TO_DECIMAL_FAST_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* store a bignum as a decimal ASCII string */
+static mp_err s_mp_to_decimal_fast_rec(const mp_int *number, mp_int *nL, mp_int *shiftL, mp_int *mL,
+                                       int precalc_array_index,
+                                       int left,
+                                       char **result,
+                                       size_t *maxlen)
+{
+   mp_int q, nLq, r;
+   mp_err err;
+
+   if (precalc_array_index < 0) {
+      int n = mp_get_i32(number), n_orig = n, digits_to_copy = 0, copy_counter;
+      char *result_str = *result;
+      char sprintf_str[4] = "000";
+
+      while (n) {
+         sprintf_str[2 - digits_to_copy] = mp_s_rmap[n % 10];
+         digits_to_copy++;
+         n /= 10;
+      }
+
+      if (!left && n_orig < 100) {
+         digits_to_copy++;
+         if (n_orig < 10) {
+            digits_to_copy++;
+         }
+         if (n_orig == 0) {
+            digits_to_copy++;
+         }
+      }
+
+      if (*maxlen < ((size_t)digits_to_copy + 1)) {
+         /* no more room */
+         return MP_VAL;
+      }
+
+      for (copy_counter = 0; copy_counter < digits_to_copy; copy_counter++) {
+         result_str[copy_counter] = sprintf_str[3 - digits_to_copy + copy_counter];
+      }
+
+      *maxlen -= (size_t)digits_to_copy;
+      *result += digits_to_copy;
+
+      return MP_OKAY;
+   }
+
+   if ((err = mp_init_multi(&q, &nLq, &r, NULL)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+   if ((err = mp_mul(number, &mL[precalc_array_index], &q)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+   if ((err = mp_div_2d(&q, mp_get_i32(&shiftL[precalc_array_index]), &q, NULL)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   if ((err = mp_mul(&nL[precalc_array_index], &q, &nLq)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   if ((err = mp_sub(number, &nLq, &r)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   if (mp_isneg(&r)) {
+      if ((err = mp_decr(&q)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = mp_add(&r, &nL[precalc_array_index], &r)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   }
+
+   --precalc_array_index;
+   if (left && mp_iszero(&q)) {
+      if ((err = s_mp_to_decimal_fast_rec(&r, nL, shiftL, mL, precalc_array_index, 1, result, maxlen)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   } else {
+      if ((err = s_mp_to_decimal_fast_rec(&q, nL, shiftL, mL, precalc_array_index, left, result, maxlen)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = s_mp_to_decimal_fast_rec(&r, nL, shiftL, mL, precalc_array_index, 0, result, maxlen)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   }
+
+   err = MP_OKAY;
+
+LBL_ERR:
+   mp_clear_multi(&q, &nLq, &r, NULL);
+   return err;
+}
+
+mp_err s_mp_to_decimal_fast(const mp_int *a, char *result, size_t maxlen)
+{
+   mp_int number, n, shift, M, M2, M22, M4, M44;
+   mp_int nL[20], shiftL[20], mL[20];
+   mp_err err;
+   char **result_addr = &result;
+   int precalc_array_index = 1, c;
+
+   /* check range of the maxlen */
+   if (maxlen < 2) {
+      return MP_VAL;
+   }
+
+   MP_ZERO_BUFFER(nL, sizeof(nL));
+   MP_ZERO_BUFFER(shiftL, sizeof(shiftL));
+   MP_ZERO_BUFFER(mL, sizeof(mL));
+
+   if ((err = mp_init_multi(&number, &n, &shift, &M, &M2, &M22, &M4, &M44, &nL[0], &shiftL[0], &mL[0], NULL)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   if ((err = mp_copy(a, &number)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+   if (mp_isneg(&number)) {
+      if ((err = mp_neg(&number, &number)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      result[0] = '-';
+      *result_addr += 1;
+      maxlen -= 1;
+   }
+   mp_set_u32(&n, 1000);
+
+   if ((err = mp_copy(&n, &nL[0])) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   mp_set_u32(&shift, 20);
+
+   if ((err = mp_copy(&shift, &shiftL[0])) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   /* (8 * 2**$shift) / $n rounded up */
+   mp_set_u32(&M, 8389);
+
+   /* $M / 8, rounded up */
+   mp_set_u32(&mL[0], 1049);
+
+   for (precalc_array_index = 1; precalc_array_index < 20; precalc_array_index++) {
+      if ((err = mp_sqr(&n, &n)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if (mp_cmp(&n, &number) == MP_GT) {
+         break;
+      }
+
+      if ((err = mp_mul_2(&shift, &shift)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      /* The following is a Newton-Raphson step, to restore the invariant
+       * that $M is (8 * 2**$shift) / $n, rounded up. */
+      {
+         if ((err = mp_sqr(&M, &M2)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_sqr(&M2, &M4)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+
+         if ((err = mp_mul(&M4, &n, &M4)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_div_2d(&M4, mp_get_i32(&shift) + 6, &M4, NULL)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_mul_2(&M2, &M2)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_sub(&M4, &M2, &M4)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_incr(&M4)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_div_2d(&M4, 3, &M4, NULL)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_decr(&M4)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_neg(&M4, &M)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+      }
+
+      if ((err = mp_init_copy(&nL[precalc_array_index], &n)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = mp_init_copy(&shiftL[precalc_array_index], &shift)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+
+      /* Divide by 8, round up */
+      {
+         if ((err = mp_add_d(&M4, 1, &M4)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_div_2d(&M4, 3, &M4, NULL)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_sub_d(&M4, 1, &M4)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+         if ((err = mp_neg(&M4, &M4)) != MP_OKAY) {
+            goto LBL_ERR;
+         }
+      }
+      if ((err = mp_init_copy(&mL[precalc_array_index], &M4)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+   }
+   if (precalc_array_index >= 20) {
+      err = MP_VAL;
+      goto LBL_ERR;
+   }
+
+   if ((err = s_mp_to_decimal_fast_rec(&number, nL, shiftL, mL, precalc_array_index - 1, 1, result_addr,
+                                       &maxlen)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+   *result_addr[0] = '\0';
+
+   err = MP_OKAY;
+
+LBL_ERR:
+   mp_clear_multi(&number, &n, &shift, &M, &M2, &M22, &M4, &M44, &nL[0], &shiftL[0], &mL[0], NULL);
+   for (c = 1; c < precalc_array_index; c++) {
+      mp_clear_multi(&nL[c], &shiftL[c], &mL[c], NULL);
+   }
+   return err;
+}
+
+#endif

demo/test.c

@@ -2214,6 +2214,54 @@ static int test_s_mp_toom_sqr(void)
    return EXIT_FAILURE;
 }
 
+static int test_mp_to_decimal(void)
+{
+   mp_int a, b;
+   int size, err, strlength;
+   char *str;
+
+   if ((err = mp_init_multi(&a, &b, NULL)) != MP_OKAY) {
+      goto LTM_ERR;
+   }
+   for (size = 1; size <= 1001; size += 10) {
+      int times;
+      printf("Testing mp_to_decimal: %5d digits    \r", size);
+      fflush(stdout);
+      for (times = 0; times < 5; times++) {
+         if ((err = mp_rand(&a, size)) != MP_OKAY) {
+            goto LTM_ERR;
+         }
+         if (times % 2) {
+            /* also test some negative numbers */
+            if ((err = mp_neg(&a, &a)) != MP_OKAY) {
+               goto LTM_ERR;
+            }
+         }
+         if ((err = mp_radix_size(&a, 10, &strlength)) != MP_OKAY) {
+            goto LTM_ERR;
+         }
+         str = (char *)malloc((size_t)strlength);
+         if ((err = mp_to_decimal(&a, str, (size_t)strlength)) != MP_OKAY) {
+            goto LTM_ERR;
+         }
+         if ((err = mp_read_radix(&b, str, 10)) != MP_OKAY) {
+            goto LTM_ERR;
+         }
+         free(str);
+         if (mp_cmp(&a, &b) != MP_EQ) {
+            fprintf(stderr, "mp_to_decimal failed at size %d\n", size);
+            goto LTM_ERR;
+         }
+      }
+   }
+
+   mp_clear_multi(&a, &b, NULL);
+   return EXIT_SUCCESS;
+LTM_ERR:
+   mp_clear_multi(&a, &b, NULL);
+   return EXIT_FAILURE;
+}
+
 int unit_tests(int argc, char **argv)
 {
    static const struct {
@@ -2264,8 +2312,10 @@ int unit_tests(int argc, char **argv)
       T1(s_mp_karatsuba_sqr, S_MP_KARATSUBA_SQR),
       T1(s_mp_toom_mul, S_MP_TOOM_MUL),
       T1(s_mp_toom_sqr, S_MP_TOOM_SQR),
+      T1(mp_to_decimal, S_MP_TO_DECIMAL_FAST)
 #undef T2
 #undef T1
+#undef T0
    };
    unsigned long i, ok, fail, nop;
    uint64_t t;

libtommath_VS2008.vcproj

@@ -816,6 +816,10 @@
 			RelativePath="bn_mp_submod.c"
 			>
 		</File>
+		<File
+			RelativePath="bn_mp_to_decimal.c"
+			>
+		</File>
 		<File
 			RelativePath="bn_mp_to_radix.c"
 			>
@@ -936,6 +940,10 @@
 			RelativePath="bn_s_mp_sub.c"
 			>
 		</File>
+		<File
+			RelativePath="bn_s_mp_to_decimal_fast.c"
+			>
+		</File>
 		<File
 			RelativePath="bn_s_mp_toom_mul.c"
 			>

makefile

@@ -48,13 +48,14 @@ bn_mp_reduce_is_2k.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_root_u32.
 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_bin_size.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_signed_bin.o bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin.o \
-bn_mp_to_unsigned_bin_n.o bn_mp_unsigned_bin_size.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_submod.o bn_mp_to_decimal.o bn_mp_to_radix.o bn_mp_to_signed_bin.o bn_mp_to_signed_bin_n.o \
+bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o bn_mp_unsigned_bin_size.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_to_decimal_fast.o bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o
 
 #END_INS