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