util/src/TclTk/tcl8.6.12/libtommath/bn_mp_log_u32.c

   1 #include "tommath_private.h"
   2 #ifdef BN_MP_LOG_U32_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 /* Compute log_{base}(a) */
   7 static mp_word s_pow(mp_word base, mp_word exponent)
   8 {
   9    mp_word result = 1;
  10    while (exponent != 0u) {
  11       if ((exponent & 1u) == 1u) {
  12          result *= base;
  13       }
  14       exponent >>= 1;
  15       base *= base;
  16    }
  17
  18    return result;
  19 }
  20
  21 static mp_digit s_digit_ilogb(mp_digit base, mp_digit n)
  22 {
  23    mp_word bracket_low = 1, bracket_mid, bracket_high, N;
  24    mp_digit ret, high = 1uL, low = 0uL, mid;
  25
  26    if (n < base) {
  27       return 0uL;
  28    }
  29    if (n == base) {
  30       return 1uL;
  31    }
  32
  33    bracket_high = (mp_word) base ;
  34    N = (mp_word) n;
  35
  36    while (bracket_high < N) {
  37       low = high;
  38       bracket_low = bracket_high;
  39       high <<= 1;
  40       bracket_high *= bracket_high;
  41    }
  42
  43    while (((mp_digit)(high - low)) > 1uL) {
  44       mid = (low + high) >> 1;
  45       bracket_mid = bracket_low * s_pow(base, (mp_word)(mid - low));
  46
  47       if (N < bracket_mid) {
  48          high = mid ;
  49          bracket_high = bracket_mid ;
  50       }
  51       if (N > bracket_mid) {
  52          low = mid ;
  53          bracket_low = bracket_mid ;
  54       }
  55       if (N == bracket_mid) {
  56          return (mp_digit) mid;
  57       }
  58    }
  59
  60    if (bracket_high == N) {
  61       ret = high;
  62    } else {
  63       ret = low;
  64    }
  65
  66    return ret;
  67 }
  68
  69 /* TODO: output could be "int" because the output of mp_radix_size is int, too,
  70          as is the output of mp_bitcount.
  71          With the same problem: max size is INT_MAX * MP_DIGIT not INT_MAX only!
  72 */
  73 mp_err mp_log_u32(const mp_int *a, unsigned int base, unsigned int *c)
  74 {
  75    mp_err err;
  76    mp_ord cmp;
  77    unsigned int high, low, mid;
  78    mp_int bracket_low, bracket_high, bracket_mid, t, bi_base;
  79
  80    err = MP_OKAY;
  81
  82    if (a->sign == MP_NEG) {
  83       return MP_VAL;
  84    }
  85
  86    if (MP_IS_ZERO(a)) {
  87       return MP_VAL;
  88    }
  89
  90    if (base < 2u) {
  91       return MP_VAL;
  92    }
  93
  94    /* A small shortcut for bases that are powers of two. */
  95    if ((base & (base - 1u)) == 0u) {
  96       int y, bit_count;
  97       for (y=0; (y < 7) && ((base & 1u) == 0u); y++) {
  98          base >>= 1;
  99       }
 100       bit_count = mp_count_bits(a) - 1;
 101       *c = (unsigned int)(bit_count/y);
 102       return MP_OKAY;
 103    }
 104
 105    if (a->used == 1) {
 106       *c = (unsigned int)s_digit_ilogb(base, a->dp[0]);
 107       return err;
 108    }
 109
 110    cmp = mp_cmp_d(a, base);
 111    if ((cmp == MP_LT) || (cmp == MP_EQ)) {
 112       *c = cmp == MP_EQ;
 113       return err;
 114    }
 115
 116    if ((err =
 117            mp_init_multi(&bracket_low, &bracket_high,
 118                          &bracket_mid, &t, &bi_base, NULL)) != MP_OKAY) {
 119       return err;
 120    }
 121
 122    low = 0u;
 123    mp_set(&bracket_low, 1uL);
 124    high = 1u;
 125
 126    mp_set(&bracket_high, base);
 127
 128    /*
 129        A kind of Giant-step/baby-step algorithm.
 130        Idea shamelessly stolen from https://programmingpraxis.com/2010/05/07/integer-logarithms/2/
 131        The effect is asymptotic, hence needs benchmarks to test if the Giant-step should be skipped
 132        for small n.
 133     */
 134    while (mp_cmp(&bracket_high, a) == MP_LT) {
 135       low = high;
 136       if ((err = mp_copy(&bracket_high, &bracket_low)) != MP_OKAY) {
 137          goto LBL_ERR;
 138       }
 139       high <<= 1;
 140       if ((err = mp_sqr(&bracket_high, &bracket_high)) != MP_OKAY) {
 141          goto LBL_ERR;
 142       }
 143    }
 144    mp_set(&bi_base, base);
 145
 146    while ((high - low) > 1u) {
 147       mid = (high + low) >> 1;
 148
 149       if ((err = mp_expt_u32(&bi_base, mid - low, &t)) != MP_OKAY) {
 150          goto LBL_ERR;
 151       }
 152       if ((err = mp_mul(&bracket_low, &t, &bracket_mid)) != MP_OKAY) {
 153          goto LBL_ERR;
 154       }
 155       cmp = mp_cmp(a, &bracket_mid);
 156       if (cmp == MP_LT) {
 157          high = mid;
 158          mp_exch(&bracket_mid, &bracket_high);
 159       }
 160       if (cmp == MP_GT) {
 161          low = mid;
 162          mp_exch(&bracket_mid, &bracket_low);
 163       }
 164       if (cmp == MP_EQ) {
 165          *c = mid;
 166          goto LBL_END;
 167       }
 168    }
 169
 170    *c = (mp_cmp(&bracket_high, a) == MP_EQ) ? high : low;
 171
 172 LBL_END:
 173 LBL_ERR:
 174    mp_clear_multi(&bracket_low, &bracket_high, &bracket_mid,
 175                   &t, &bi_base, NULL);
 176    return err;
 177 }
 178
 179
 180 #endif