84 lines
3.2 KiB
C
84 lines
3.2 KiB
C
#include "tommath_private.h"
|
|
#ifdef S_MP_FP_LOG_C
|
|
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
|
|
/* SPDX-License-Identifier: Unlicense */
|
|
|
|
/* Fixed point bitwise logarithm base two of "a" with precision "p", Big-integer version. */
|
|
static mp_err s_mp_fp_log_fraction(const mp_int *a, int p, mp_int *c)
|
|
{
|
|
mp_int b, L_out, twoep, a_bar;
|
|
int L, pmL;
|
|
int i;
|
|
mp_err err;
|
|
|
|
if ((err = mp_init_multi(&b, &L_out, &twoep, &a_bar, NULL)) != MP_OKAY) {
|
|
return err;
|
|
}
|
|
|
|
L = mp_count_bits(a) - 1;
|
|
pmL = (p < L) ? L - p: p - L;
|
|
if ((err = mp_mul_2d(a, pmL, &a_bar)) != MP_OKAY) goto LTM_ERR;
|
|
if ((err = mp_2expt(&b, p - 1)) != MP_OKAY) goto LTM_ERR;
|
|
mp_set_i32(&L_out, L);
|
|
if ((err = mp_mul_2d(&L_out, p, &L_out)) != MP_OKAY) goto LTM_ERR;
|
|
|
|
if ((err = mp_2expt(&twoep, p + 1)) != MP_OKAY) goto LTM_ERR;
|
|
|
|
for (i = 0; i < p; i++) {
|
|
if ((err = mp_sqr(&a_bar, &a_bar)) != MP_OKAY) goto LTM_ERR;
|
|
if ((err = mp_div_2d(&a_bar, p, &a_bar, NULL)) != MP_OKAY) goto LTM_ERR;
|
|
if (mp_cmp(&a_bar, &twoep) != MP_LT) {
|
|
if ((err = mp_div_2(&a_bar, &a_bar)) != MP_OKAY) goto LTM_ERR;
|
|
if ((err = mp_add(&L_out, &b, &L_out)) != MP_OKAY) goto LTM_ERR;
|
|
}
|
|
if ((err = mp_div_2(&b, &b)) != MP_OKAY) goto LTM_ERR;
|
|
}
|
|
|
|
mp_exch(c, &L_out);
|
|
|
|
LTM_ERR:
|
|
mp_clear_multi(&b, &L_out, &twoep, &a_bar, NULL);
|
|
return err;
|
|
}
|
|
|
|
mp_err s_mp_fp_log(const mp_int *a, mp_int *c)
|
|
{
|
|
mp_int La, t;
|
|
mp_err err;
|
|
int fla;
|
|
/*
|
|
We have arbitrary precision here and could adapt "prec" to actual precision,
|
|
there is no need for lowering precision for small numbers. Some quick runtime
|
|
tests revealed a gain that can only be called marginally at best.
|
|
But: YMMV, as always.
|
|
*/
|
|
int prec = MP_PRECISION_FIXED_LOG;
|
|
|
|
fla = mp_count_bits(a) - 1;
|
|
|
|
if ((err = mp_init_multi(&La, &t, NULL)) != MP_OKAY) {
|
|
return err;
|
|
}
|
|
|
|
if (fla > prec) {
|
|
if ((err = mp_div_2d(a, fla - prec, &t, NULL)) != MP_OKAY) goto LTM_ERR;
|
|
if ((err = s_mp_fp_log_fraction(&t, prec,
|
|
&La)) != MP_OKAY) goto LTM_ERR;
|
|
mp_set_i32(&t,fla - prec);
|
|
if ((err = mp_mul_2d(&t,prec, &t)) != MP_OKAY) goto LTM_ERR;
|
|
if ((err = mp_add(&La, &t, &La)) != MP_OKAY) goto LTM_ERR;
|
|
} else {
|
|
if ((err = s_mp_fp_log_fraction(a, prec,
|
|
&La)) != MP_OKAY) goto LTM_ERR;
|
|
}
|
|
|
|
mp_exch(&La, c);
|
|
|
|
LTM_ERR:
|
|
mp_clear_multi(&La, &t, NULL);
|
|
return err;
|
|
}
|
|
|
|
|
|
#endif
|