diff options
Diffstat (limited to 'libtommath/bn_mp_n_root.c')
| -rw-r--r-- | libtommath/bn_mp_n_root.c | 118 |
1 files changed, 8 insertions, 110 deletions
diff --git a/libtommath/bn_mp_n_root.c b/libtommath/bn_mp_n_root.c index e6bf725..a14ee67 100644 --- a/libtommath/bn_mp_n_root.c +++ b/libtommath/bn_mp_n_root.c @@ -1,4 +1,4 @@ -#include <tommath.h> +#include <tommath_private.h> #ifdef BN_MP_N_ROOT_C /* LibTomMath, multiple-precision integer library -- Tom St Denis * @@ -12,121 +12,19 @@ * The library is free for all purposes without any express * guarantee it works. * - * Tom St Denis, tomstdenis@gmail.com, http://libtom.org + * Tom St Denis, tstdenis82@gmail.com, http://libtom.org */ -/* find the n'th root of an integer - * - * Result found such that (c)**b <= a and (c+1)**b > a - * - * This algorithm uses Newton's approximation - * x[i+1] = x[i] - f(x[i])/f'(x[i]) - * which will find the root in log(N) time where - * each step involves a fair bit. This is not meant to - * find huge roots [square and cube, etc]. +/* wrapper function for mp_n_root_ex() + * computes c = (a)**(1/b) such that (c)**b <= a and (c+1)**b > a */ int mp_n_root (mp_int * a, mp_digit b, mp_int * c) { - mp_int t1, t2, t3; - int res, neg; - - /* input must be positive if b is even */ - if ((b & 1) == 0 && a->sign == MP_NEG) { - return MP_VAL; - } - - if ((res = mp_init (&t1)) != MP_OKAY) { - return res; - } - - if ((res = mp_init (&t2)) != MP_OKAY) { - goto LBL_T1; - } - - if ((res = mp_init (&t3)) != MP_OKAY) { - goto LBL_T2; - } - - /* if a is negative fudge the sign but keep track */ - neg = a->sign; - a->sign = MP_ZPOS; - - /* t2 = 2 */ - mp_set (&t2, 2); - - do { - /* t1 = t2 */ - if ((res = mp_copy (&t2, &t1)) != MP_OKAY) { - goto LBL_T3; - } - - /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */ - - /* t3 = t1**(b-1) */ - if ((res = mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) { - goto LBL_T3; - } - - /* numerator */ - /* t2 = t1**b */ - if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) { - goto LBL_T3; - } - - /* t2 = t1**b - a */ - if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) { - goto LBL_T3; - } - - /* denominator */ - /* t3 = t1**(b-1) * b */ - if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) { - goto LBL_T3; - } - - /* t3 = (t1**b - a)/(b * t1**(b-1)) */ - if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) { - goto LBL_T3; - } - - if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) { - goto LBL_T3; - } - } while (mp_cmp (&t1, &t2) != MP_EQ); - - /* result can be off by a few so check */ - for (;;) { - if ((res = mp_expt_d (&t1, b, &t2)) != MP_OKAY) { - goto LBL_T3; - } - - if (mp_cmp (&t2, a) == MP_GT) { - if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) { - goto LBL_T3; - } - } else { - break; - } - } - - /* reset the sign of a first */ - a->sign = neg; - - /* set the result */ - mp_exch (&t1, c); - - /* set the sign of the result */ - c->sign = neg; - - res = MP_OKAY; - -LBL_T3:mp_clear (&t3); -LBL_T2:mp_clear (&t2); -LBL_T1:mp_clear (&t1); - return res; + return mp_n_root_ex(a, b, c, 0); } + #endif /* $Source$ */ -/* $Revision: 0.41 $ */ -/* $Date: 2007-04-18 09:58:18 +0000 $ */ +/* $Revision$ */ +/* $Date$ */ |