1/* $NetBSD: bn_mp_n_root.c,v 1.2 2017/01/28 21:31:47 christos Exp $ */ 2 3#include <tommath.h> 4#ifdef BN_MP_N_ROOT_C 5/* LibTomMath, multiple-precision integer library -- Tom St Denis 6 * 7 * LibTomMath is a library that provides multiple-precision 8 * integer arithmetic as well as number theoretic functionality. 9 * 10 * The library was designed directly after the MPI library by 11 * Michael Fromberger but has been written from scratch with 12 * additional optimizations in place. 13 * 14 * The library is free for all purposes without any express 15 * guarantee it works. 16 * 17 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org 18 */ 19 20/* find the n'th root of an integer 21 * 22 * Result found such that (c)**b <= a and (c+1)**b > a 23 * 24 * This algorithm uses Newton's approximation 25 * x[i+1] = x[i] - f(x[i])/f'(x[i]) 26 * which will find the root in log(N) time where 27 * each step involves a fair bit. This is not meant to 28 * find huge roots [square and cube, etc]. 29 */ 30int mp_n_root (mp_int * a, mp_digit b, mp_int * c) 31{ 32 mp_int t1, t2, t3; 33 int res, neg; 34 35 /* input must be positive if b is even */ 36 if ((b & 1) == 0 && a->sign == MP_NEG) { 37 return MP_VAL; 38 } 39 40 if ((res = mp_init (&t1)) != MP_OKAY) { 41 return res; 42 } 43 44 if ((res = mp_init (&t2)) != MP_OKAY) { 45 goto LBL_T1; 46 } 47 48 if ((res = mp_init (&t3)) != MP_OKAY) { 49 goto LBL_T2; 50 } 51 52 /* if a is negative fudge the sign but keep track */ 53 neg = a->sign; 54 a->sign = MP_ZPOS; 55 56 /* t2 = 2 */ 57 mp_set (&t2, 2); 58 59 do { 60 /* t1 = t2 */ 61 if ((res = mp_copy (&t2, &t1)) != MP_OKAY) { 62 goto LBL_T3; 63 } 64 65 /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */ 66 67 /* t3 = t1**(b-1) */ 68 if ((res = mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) { 69 goto LBL_T3; 70 } 71 72 /* numerator */ 73 /* t2 = t1**b */ 74 if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) { 75 goto LBL_T3; 76 } 77 78 /* t2 = t1**b - a */ 79 if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) { 80 goto LBL_T3; 81 } 82 83 /* denominator */ 84 /* t3 = t1**(b-1) * b */ 85 if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) { 86 goto LBL_T3; 87 } 88 89 /* t3 = (t1**b - a)/(b * t1**(b-1)) */ 90 if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) { 91 goto LBL_T3; 92 } 93 94 if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) { 95 goto LBL_T3; 96 } 97 } while (mp_cmp (&t1, &t2) != MP_EQ); 98 99 /* result can be off by a few so check */ 100 for (;;) { 101 if ((res = mp_expt_d (&t1, b, &t2)) != MP_OKAY) { 102 goto LBL_T3; 103 } 104 105 if (mp_cmp (&t2, a) == MP_GT) { 106 if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) { 107 goto LBL_T3; 108 } 109 } else { 110 break; 111 } 112 } 113 114 /* reset the sign of a first */ 115 a->sign = neg; 116 117 /* set the result */ 118 mp_exch (&t1, c); 119 120 /* set the sign of the result */ 121 c->sign = neg; 122 123 res = MP_OKAY; 124 125LBL_T3:mp_clear (&t3); 126LBL_T2:mp_clear (&t2); 127LBL_T1:mp_clear (&t1); 128 return res; 129} 130#endif 131 132/* Source: /cvs/libtom/libtommath/bn_mp_n_root.c,v */ 133/* Revision: 1.4 */ 134/* Date: 2006/12/28 01:25:13 */ 135