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