mp_log.c

   1 #include "tommath_private.h"
   2 #ifdef MP_LOG_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 #define MP_LOG2_EXPT(a,b) ((mp_count_bits((a)) - 1) / mp_cnt_lsb((b)))
   7
   8 static mp_err s_approx_log_d(const mp_int *a, const mp_int *b, int *lb)
   9 {
  10    mp_word La, Lb;
  11    mp_err err;
  12
  13    /* Approximation of the individual logarithms with low precision */
  14    if ((err = s_mp_fp_log_d(a, &La)) != MP_OKAY)                                                     goto LTM_ERR;
  15    if ((err = s_mp_fp_log_d(b, &Lb)) != MP_OKAY)                                                     goto LTM_ERR;
  16
  17    /* Approximation of log_b(a) with low precision. */
  18    *lb = (int)(((La - (Lb + 1)/2) / Lb) + 1);
  19    /* TODO: just floor it instead? Multiplication is cheaper than division. */
  20    /* *lb = (int)(la_word / lb_word); */
  21    err = MP_OKAY;
  22 LTM_ERR:
  23    return err;
  24 }
  25
  26 static mp_err s_approx_log(const mp_int *a, const mp_int *b, int *lb)
  27 {
  28    mp_int La, Lb, t;
  29    mp_err err;
  30
  31    if ((err = mp_init_multi(&La, &Lb, &t, NULL)) != MP_OKAY) {
  32       return err;
  33    }
  34
  35    if ((err = s_mp_fp_log(a, &La)) != MP_OKAY)                                                            goto LTM_ERR;
  36    if ((err = s_mp_fp_log(b, &Lb)) != MP_OKAY)                                                            goto LTM_ERR;
  37
  38    if ((err = mp_add_d(&Lb, 1u, &t)) != MP_OKAY)                                                          goto LTM_ERR;
  39    if ((err = mp_div_2(&t, &t)) != MP_OKAY)                                                               goto LTM_ERR;
  40    if ((err = mp_sub(&La, &t, &t)) != MP_OKAY)                                                            goto LTM_ERR;
  41    if ((err = mp_div(&t, &Lb, &t, NULL)) != MP_OKAY)                                                      goto LTM_ERR;
  42    if ((err = mp_add_d(&t, 1u, &t)) != MP_OKAY)                                                           goto LTM_ERR;
  43
  44    *lb = mp_get_i32(&t);
  45    err = MP_OKAY;
  46 LTM_ERR:
  47    mp_clear_multi(&t, &Lb, &La, NULL);
  48    return err;
  49 }
  50
  51 mp_err mp_log(const mp_int *a, const mp_int *b, int *lb)
  52 {
  53    mp_int bn;
  54    int n, fla, flb;
  55    mp_err err;
  56    mp_ord cmp;
  57
  58    if (mp_isneg(a) || mp_iszero(a) || (mp_cmp_d(b, 2u) == MP_LT)) {
  59       return MP_VAL;
  60    }
  61
  62    if (MP_IS_POWER_OF_TWO(b)) {
  63       *lb = MP_LOG2_EXPT(a, b);
  64       return MP_OKAY;
  65    }
  66
  67    /* floor(log_2(x)) for cut-off */
  68    fla = mp_count_bits(a) - 1;
  69    flb = mp_count_bits(b) - 1;
  70
  71    cmp = mp_cmp(a, b);
  72
  73    /* "a < b -> 0" and "(b == a) -> 1" */
  74    if ((cmp == MP_LT) || (cmp == MP_EQ)) {
  75       *lb = (cmp == MP_EQ);
  76       return MP_OKAY;
  77    }
  78
  79    /* "a < b^2 -> 1" (bit-count is sufficient, doesn't need to be exact) */
  80    if (((2 * flb)-1) > fla) {
  81       *lb = 1;
  82       return MP_OKAY;
  83    }
  84
  85    if (MP_HAS(S_MP_WORD_TOO_SMALL)) {
  86       err = s_approx_log(a, b, &n);
  87    } else {
  88       err = s_approx_log_d(a, b, &n);
  89    }
  90    if (err != MP_OKAY) {
  91       return err;
  92    }
  93
  94    if ((err = mp_init(&bn)) != MP_OKAY) {
  95       return err;
  96    }
  97
  98    /* Check result. Result is wrong by 2(two) at most. */
  99    if ((err = mp_expt_n(b, n, &bn)) != MP_OKAY) {
 100       /* If the approximation overshot we can give it another try */
 101       if (err == MP_OVF) {
 102          n--;
 103          /* But only one */
 104          if ((err = mp_expt_n(b, n, &bn)) != MP_OKAY)                                                     goto LTM_ERR;
 105       } else {
 106          goto LTM_ERR;
 107       }
 108    }
 109
 110    cmp = mp_cmp(&bn, a);
 111
 112    /* The rare case of a perfect power makes a perfect shortcut, too. */
 113    if (cmp == MP_EQ) {
 114       *lb = n;
 115       goto LTM_OUT;
 116    }
 117
 118    /* We have to make at least one multiplication because it could still be a perfect power. */
 119    if (cmp == MP_LT) {
 120       do {
 121          /* Full big-integer operations are to be avoided if possible */
 122          if (b->used == 1) {
 123             if ((err = mp_mul_d(&bn, b->dp[0], &bn)) != MP_OKAY) {
 124                if (err == MP_OVF) {
 125                   goto LTM_OUT;
 126                }
 127                goto LTM_ERR;
 128             }
 129          } else {
 130             if ((err = mp_mul(&bn, b, &bn)) != MP_OKAY) {
 131                if (err == MP_OVF) {
 132                   goto LTM_OUT;
 133                }
 134                goto LTM_ERR;
 135             }
 136          }
 137          n++;
 138       } while ((cmp = mp_cmp(&bn, a)) == MP_LT);
 139       /* Overshot, take it back. */
 140       if (cmp == MP_GT) {
 141          n--;
 142       }
 143       goto LTM_OUT;
 144    }
 145
 146    /* But it can overestimate, too, for example if "a" is closely below some "b^k" */
 147    if (cmp == MP_GT) {
 148       do {
 149          if (b->used == 1) {
 150             /* These are cheaper exact divisions, but that function is not available in LTM */
 151             if ((err = mp_div_d(&bn, b->dp[0], &bn, NULL)) != MP_OKAY)                                    goto LTM_ERR;
 152
 153          } else {
 154             if ((err = mp_div(&bn, b, &bn, NULL)) != MP_OKAY)                                             goto LTM_ERR;
 155          }
 156          n--;
 157       } while ((cmp = mp_cmp(&bn, a)) == MP_GT);
 158    }
 159
 160 LTM_OUT:
 161    *lb = n;
 162    err = MP_OKAY;
 163 LTM_ERR:
 164    mp_clear(&bn);
 165    return err;
 166 }
 167
 168 #endif