Line data    Source code

       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           0 : static mp_word s_pow(mp_word base, mp_word exponent)
       8             : {
       9           0 :    mp_word result = 1uLL;
      10           0 :    while (exponent != 0u) {
      11           0 :       if ((exponent & 1u) == 1u) {
      12           0 :          result *= base;
      13             :       }
      14           0 :       exponent >>= 1;
      15           0 :       base *= base;
      16             :    }
      17             : 
      18           0 :    return result;
      19             : }
      20             : 
      21           0 : static mp_digit s_digit_ilogb(mp_digit base, mp_digit n)
      22             : {
      23           0 :    mp_word bracket_low = 1uLL, bracket_mid, bracket_high, N;
      24           0 :    mp_digit ret, high = 1uL, low = 0uL, mid;
      25             : 
      26           0 :    if (n < base) {
      27           0 :       return 0uL;
      28             :    }
      29           0 :    if (n == base) {
      30           0 :       return 1uL;
      31             :    }
      32             : 
      33           0 :    bracket_high = (mp_word) base ;
      34           0 :    N = (mp_word) n;
      35             : 
      36           0 :    while (bracket_high < N) {
      37           0 :       low = high;
      38           0 :       bracket_low = bracket_high;
      39           0 :       high <<= 1;
      40           0 :       bracket_high *= bracket_high;
      41             :    }
      42             : 
      43           0 :    while (((mp_digit)(high - low)) > 1uL) {
      44           0 :       mid = (low + high) >> 1;
      45           0 :       bracket_mid = bracket_low * s_pow(base, (mp_word)(mid - low));
      46             : 
      47           0 :       if (N < bracket_mid) {
      48           0 :          high = mid ;
      49           0 :          bracket_high = bracket_mid ;
      50             :       }
      51           0 :       if (N > bracket_mid) {
      52           0 :          low = mid ;
      53           0 :          bracket_low = bracket_mid ;
      54             :       }
      55           0 :       if (N == bracket_mid) {
      56           0 :          return (mp_digit) mid;
      57             :       }
      58             :    }
      59             : 
      60           0 :    if (bracket_high == N) {
      61           0 :       ret = high;
      62             :    } else {
      63           0 :       ret = low;
      64             :    }
      65             : 
      66           0 :    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           0 : mp_err mp_log_u32(const mp_int *a, uint32_t base, uint32_t *c)
      74             : {
      75           0 :    mp_err err;
      76           0 :    mp_ord cmp;
      77           0 :    uint32_t high, low, mid;
      78           0 :    mp_int bracket_low, bracket_high, bracket_mid, t, bi_base;
      79             : 
      80           0 :    err = MP_OKAY;
      81             : 
      82           0 :    if (a->sign == MP_NEG) {
      83           0 :       return MP_VAL;
      84             :    }
      85             : 
      86           0 :    if (MP_IS_ZERO(a)) {
      87           0 :       return MP_VAL;
      88             :    }
      89             : 
      90           0 :    if (base < 2u) {
      91           0 :       return MP_VAL;
      92             :    }
      93             : 
      94             :    /* `base' is at least 2 */
      95             : 
      96             :    /* A small shortcut for bases that are powers of two. */
      97           0 :    if ((base & (base - 1u)) == 0u) {
      98             :       int y, bit_count;
      99             : 
     100           0 :       for (y=0; (y < 7) && ((base & 1u) == 0u); y++) {
     101             :          /* We must go through this loop at least once */
     102           0 :          base >>= 1;
     103             :       }
     104           0 :       bit_count = mp_count_bits(a) - 1;
     105             :       /*
     106             :        * `y' is necessarily at least 1 because `base' is a power of two and
     107             :        * larger than 1, so we must have gone through the loop at least once, so
     108             :        * we can't be dividing by zero.
     109             :        *
     110             :        * scan-build thinks we can be dividing by zero...  WAT.
     111             :        */
     112           0 :       *c = (uint32_t)(bit_count/y);
     113           0 :       return MP_OKAY;
     114             :    }
     115             : 
     116           0 :    if (a->used == 1) {
     117           0 :       *c = (uint32_t)s_digit_ilogb(base, a->dp[0]);
     118           0 :       return err;
     119             :    }
     120             : 
     121           0 :    cmp = mp_cmp_d(a, base);
     122           0 :    if ((cmp == MP_LT) || (cmp == MP_EQ)) {
     123           0 :       *c = cmp == MP_EQ;
     124           0 :       return err;
     125             :    }
     126             : 
     127           0 :    if ((err =
     128           0 :            mp_init_multi(&bracket_low, &bracket_high,
     129             :                          &bracket_mid, &t, &bi_base, NULL)) != MP_OKAY) {
     130           0 :       return err;
     131             :    }
     132             : 
     133           0 :    low = 0u;
     134           0 :    mp_set(&bracket_low, 1uL);
     135           0 :    high = 1u;
     136             : 
     137           0 :    mp_set(&bracket_high, base);
     138             : 
     139             :    /*
     140             :        A kind of Giant-step/baby-step algorithm.
     141             :        Idea shamelessly stolen from https://programmingpraxis.com/2010/05/07/integer-logarithms/2/
     142             :        The effect is asymptotic, hence needs benchmarks to test if the Giant-step should be skipped
     143             :        for small n.
     144             :     */
     145           0 :    while (mp_cmp(&bracket_high, a) == MP_LT) {
     146           0 :       low = high;
     147           0 :       if ((err = mp_copy(&bracket_high, &bracket_low)) != MP_OKAY) {
     148           0 :          goto LBL_ERR;
     149             :       }
     150           0 :       high <<= 1;
     151           0 :       if ((err = mp_sqr(&bracket_high, &bracket_high)) != MP_OKAY) {
     152           0 :          goto LBL_ERR;
     153             :       }
     154             :    }
     155           0 :    mp_set(&bi_base, base);
     156             : 
     157           0 :    while ((high - low) > 1u) {
     158           0 :       mid = (high + low) >> 1;
     159             : 
     160           0 :       if ((err = mp_expt_u32(&bi_base, (uint32_t)(mid - low), &t)) != MP_OKAY) {
     161           0 :          goto LBL_ERR;
     162             :       }
     163           0 :       if ((err = mp_mul(&bracket_low, &t, &bracket_mid)) != MP_OKAY) {
     164           0 :          goto LBL_ERR;
     165             :       }
     166           0 :       cmp = mp_cmp(a, &bracket_mid);
     167           0 :       if (cmp == MP_LT) {
     168           0 :          high = mid;
     169           0 :          mp_exch(&bracket_mid, &bracket_high);
     170             :       }
     171           0 :       if (cmp == MP_GT) {
     172           0 :          low = mid;
     173           0 :          mp_exch(&bracket_mid, &bracket_low);
     174             :       }
     175           0 :       if (cmp == MP_EQ) {
     176           0 :          *c = mid;
     177           0 :          goto LBL_END;
     178             :       }
     179             :    }
     180             : 
     181           0 :    *c = (mp_cmp(&bracket_high, a) == MP_EQ) ? high : low;
     182             : 
     183           0 : LBL_END:
     184           0 : LBL_ERR:
     185           0 :    mp_clear_multi(&bracket_low, &bracket_high, &bracket_mid,
     186             :                   &t, &bi_base, NULL);
     187           0 :    return err;
     188             : }
     189             : 
     190             : 
     191             : #endif