Addition of faster read_radix method

Author: czurnieden

demo/test.c

@@ -1151,12 +1151,17 @@ static int test_mp_montgomery_reduce(void)
 
 }
 
+#include <time.h>
 static int test_mp_read_radix(void)
 {
    char buf[4096];
-   size_t written;
+   size_t written, maxlen;
    int bignum, i;
 
+   char *buffer, *bcpy;
+
+   clock_t start, stop, t_slow, t_fast;
+
    mp_int a, b;
    DOR(mp_init_multi(&a, &b, NULL));
 
@@ -1187,25 +1192,73 @@ static int test_mp_read_radix(void)
    /* Test the fast method with a slightly larger number */
 
    /* Must be bigger than the cut-off value, of course */
-   bignum = 2* (2 * s_mp_radix_exponent_y[2] * MP_RADIX_BARRETT_START_MULTIPLICATOR);
-   printf("Size of bignum_size = %d\n", bignum);
-   /* Check if "bignum" is small enough for the result to fit into "buf"
-      otherwise lead tester to this function */
-   if (bignum >= 4096) {
-      fprintf(stderr, "Buffer too small, please check function \"test_mp_read_radix\" in \"test.c\"");
+   bignum = (2 * 20 * MP_RADIX_BARRETT_START_MULTIPLICATOR)  * 10;
+   buffer = (char *)malloc((size_t)(bignum + 2));
+   if (buffer == NULL) {
       goto LBL_ERR;
    }
-   /* Produce a random number */
-   bignum /= MP_DIGIT_BIT;
-   DO(mp_rand(&b, bignum));
-   /* Check if it makes the round */
-   printf("Number of limbs in &b = %d, bit_count of &b = %d\n", bignum, mp_count_bits(&b));
+   DO(mp_rand(&a, bignum / MP_DIGIT_BIT));
+   printf("\nNumber of limbs in &b = %d, bit_count of &b = %d\n", bignum / MP_DIGIT_BIT, mp_count_bits(&a));
+   start = clock();
+   for (i = 2; i < 65; i++) {
+      /* printf("FAST radix = %d\n",i); */
+      DO(mp_to_radix(&a, buffer, (size_t)(bignum + 1), &written, i));
+      DO(mp_read_radix(&b, buffer, i));
+      EXPECT(mp_cmp(&a, &b) == MP_EQ);
+   }
+   stop = clock();
+   t_fast = stop - start;
+
+   printf("Same number, slow radix conversions\n");
+   start = clock();
    for (i = 2; i < 65; i++) {
-      DO(mp_to_radix(&b, buf, sizeof(buf), &written, i));
-      DO(mp_read_radix(&a, buf, i));
+      /* printf("SLOW radix = %d\n",i); */
+      maxlen = (size_t)(bignum + 1);
+      bcpy = buffer;
+      DO(s_mp_slower_to_radix(&a, &bcpy, &maxlen, &written, i, false));
+      DO(s_mp_slower_read_radix(&b, bcpy, 0, strlen(bcpy), i));
       EXPECT(mp_cmp(&a, &b) == MP_EQ);
-      /* fprintf(stderr,"radix = %d\n",i); */
    }
+   stop = clock();
+   t_slow = stop - start;
+
+   /* It is "long int" in GLibC but can be bigger and/or even a floating point elsewhere */
+   printf("SLOW: %.10f, FAST: %.10f\n", (double)t_slow/(double)CLOCKS_PER_SEC, (double)t_fast/(double)CLOCKS_PER_SEC);
+
+   /* Check if the branching works. */
+   if (MP_HAS(S_MP_FASTER_READ_RADIX) && MP_HAS(S_MP_FASTER_TO_RADIX)) {
+      if (t_fast > t_slow) {
+         fprintf(stderr, "Timing suspicious in test_mp_read_radix. No fast multiplication? Cut-off too low?\n");
+         goto LBL_ERR;
+      }
+   }
+
+
+   free(buffer);
+
+#if ((MP_DIGIT_BIT <= 16) && (defined MP_CHECK_RADIX_OVF))
+   /* Check a number of size (MP_MAX_DIGIT_COUNT * MP_DIGIT_BIT - 1) at fixed radix "10". */
+   /* Will not work if test is run on platforms with larger int's because
+         #define MP_MAX_DIGIT_COUNT ((INT_MAX - 2) / MP_DIGIT_BIT)
+      So we have to replace the value for INT_MAX with 2^15 - 1 = 32767 to test 16-bit int's. Not
+      very elegant but it works.
+   */
+   bignum = ((32767 - 2) / MP_DIGIT_BIT);
+   bignum = ((bignum - 1) * MP_DIGIT_BIT) + (MP_DIGIT_BIT - 1);
+   /* Manual computation because the automatic methods might not have been included in the build */
+   buffer = (char *)malloc(((bignum + 2)/1000) * 333);
+   if (buffer == NULL) {
+      goto LBL_ERR;
+   }
+   DO(mp_2expt(&a, bignum));
+   DO(mp_decr(&a));
+   printf("Number of limbs in &b = %d, bit_count of &b = %d\n", bignum / MP_DIGIT_BIT, mp_count_bits(&a));
+   DO(mp_to_radix(&a, buffer, ((bignum + 2)/1000) * 333, &written, 10));
+   DO(mp_read_radix(&b, buffer, 10));
+   EXPECT(mp_cmp(&a, &b) == MP_EQ);
+   free(buffer);
+#endif
+
 
 
    while (0) {
@@ -2504,7 +2557,7 @@ static int unit_tests(int argc, char **argv)
       T1(mp_prime_next_prime, MP_PRIME_NEXT_PRIME),
       T1(mp_prime_rand, MP_PRIME_RAND),
       T1(mp_rand, MP_RAND),
-      T1(mp_read_radix, MP_READ_RADIX),
+      T3(mp_read_radix, ONLY_PUBLIC_API, MP_READ_RADIX, MP_TO_RADIX),
       T1(mp_read_write_ubin, MP_TO_UBIN),
       T1(mp_read_write_sbin, MP_TO_SBIN),
       T1(mp_reduce_2k, MP_REDUCE_2K),

mp_read_radix.c

@@ -3,11 +3,23 @@
 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
 /* SPDX-License-Identifier: Unlicense */
 
+static size_t s_mp_strlen(const char *s)
+{
+   const char *p;
+   p = s;
+   while (*p != '\0') {
+      p++;
+   }
+   return (size_t)(p - s);
+}
+
 /* read a string [ASCII] in a given radix */
 mp_err mp_read_radix(mp_int *a, const char *str, int radix)
 {
+
    mp_err   err;
    mp_sign  sign = MP_ZPOS;
+   size_t slen;
 
    /* make sure the radix is ok */
    if ((radix < 2) || (radix > 64)) {
@@ -22,48 +34,24 @@ mp_err mp_read_radix(mp_int *a, const char *str, int radix)
       sign = MP_NEG;
    }
 
-   /* set the integer to the default of zero */
-   mp_zero(a);
-
-   /* process each digit of the string */
-   while (*str != '\0') {
-      /* if the radix <= 36 the conversion is case insensitive
-       * this allows numbers like 1AB and 1ab to represent the same  value
-       * [e.g. in hex]
-       */
-      uint8_t y;
-      char ch = (radix <= 36) ? (char)MP_TOUPPER((int)*str) : *str;
-      unsigned pos = (unsigned)(ch - '+');
-      if (MP_RADIX_MAP_REVERSE_SIZE <= pos) {
-         break;
-      }
-      y = s_mp_radix_map_reverse[pos];
+   slen = s_mp_strlen(str);
 
-      /* if the char was found in the map
-       * and is less than the given radix add it
-       * to the number, otherwise exit the loop.
-       */
-      if (y >= radix) {
-         break;
-      }
-      if ((err = mp_mul_d(a, (mp_digit)radix, a)) != MP_OKAY) {
-         return err;
-      }
-      if ((err = mp_add_d(a, y, a)) != MP_OKAY) {
-         return err;
-      }
-      ++str;
-   }
+   mp_zero(a);
 
-   /* if an illegal character was found, fail. */
-   if ((*str != '\0') && (*str != '\r') && (*str != '\n')) {
-      return MP_VAL;
+   /* Try faster version first */
+   if (MP_HAS(S_MP_FASTER_READ_RADIX)) {
+      if ((err = s_mp_faster_read_radix(a, str, 0, slen, radix)) != MP_OKAY)                             goto LTM_ERR;
+   } else if (MP_HAS(S_MP_SLOWER_READ_RADIX)) {
+      if ((err = s_mp_slower_read_radix(a, str, 0, slen, radix)) != MP_OKAY)                             goto LTM_ERR;
    }
 
    /* set the sign only if a != 0 */
    if (!mp_iszero(a)) {
       a->sign = sign;
    }
-   return MP_OKAY;
+
+LTM_ERR:
+   return err;
 }
+
 #endif

s_mp_faster_read_radix.c

@@ -0,0 +1,69 @@
+#include "tommath_private.h"
+#ifdef S_MP_FASTER_READ_RADIX_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* TODO: It is tunable */
+#define MP_READ_RADIX_CUTOFF_MULTIPLICATOR 3
+
+/* for(n=2,64,t = ceil(log(2^100)/log(n)); printf(t", "); ) */
+static const uint8_t s_read_radix_cutoff[65] = { 0, 0,                            /*  0 .. 1*/
+                                                 100, 64, 50, 44, 39, 36, 34, 32, /*  2 .. 9 */
+                                                 31, 29, 28, 28, 27, 26, 25, 25,  /* 10 .. 17 */
+                                                 24, 24, 24, 23, 23, 23, 22, 22,  /* 18 .. 25 */
+                                                 22, 22, 21, 21, 21, 21, 20, 20,  /* 26 .. 33 */
+                                                 20, 20, 20, 20, 20, 19, 19, 19,  /* 34 .. 41 */
+                                                 19, 19, 19, 19, 19, 19, 18, 18,  /* 42 .. 49 */
+                                                 18, 18, 18, 18, 18, 18, 18, 18,  /* 51 .. 57 */
+                                                 18, 17, 17, 17, 17, 17, 17       /* 58 .. 64 */
+                                               };
+
+/* This is in mp_prime_is_prime.c and can be reused */
+static int s_floor_ilog2(int value)
+{
+   int r = 0;
+   while ((value >>= 1) != 0) {
+      r++;
+   }
+   return r;
+}
+
+mp_err s_mp_faster_read_radix(mp_int *a, const char *str, size_t start, size_t end, int radix)
+{
+   size_t len, mid;
+   mp_int A, B, m;
+   mp_err err = MP_OKAY;
+
+   len = end - start;
+
+   if (len < (size_t)(s_read_radix_cutoff[radix] * MP_READ_RADIX_CUTOFF_MULTIPLICATOR)) {
+      return s_mp_slower_read_radix(a, str, start, end, radix);
+   }
+
+   mid = len / 2u;
+
+   if ((err = mp_init_set(&m, (mp_digit)radix)) != MP_OKAY) {
+      return err;
+   }
+   if ((err = mp_init_multi(&A, &B, NULL)) != MP_OKAY) {
+      mp_clear(&m);
+      return err;
+   }
+
+   if ((err = s_mp_slower_read_radix(&A, str, start, start + mid + 1, radix)) != MP_OKAY)                goto LTM_ERR;
+   if ((err = s_mp_slower_read_radix(&B, str, start + mid +1, end, radix)) != MP_OKAY)                   goto LTM_ERR;
+
+   if (MP_IS_2EXPT((unsigned int)radix)) {
+      if ((err = mp_mul_2d(&A, (int)(((len - mid) - 1u) * (size_t)s_floor_ilog2(radix)), &A)) != MP_OKAY)goto LTM_ERR;
+   } else {
+      if ((err = mp_expt_n(&m, (int)((len - mid) - 1u), &m)) != MP_OKAY)                                 goto LTM_ERR;
+      if ((err = mp_mul(&A, &m, &A)) != MP_OKAY)                                                         goto LTM_ERR;
+   }
+   if ((err = mp_add(&A, &B, a)) != MP_OKAY)                                                             goto LTM_ERR;
+
+LTM_ERR:
+   mp_clear_multi(&A, &B, &m, NULL);
+   return err;
+}
+
+#endif