|
282
|
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, [email protected], http://math.libtomcrypt.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 */ |
|
|
28 int 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 |
|
|
123 LBL_T3:mp_clear (&t3); |
|
|
124 LBL_T2:mp_clear (&t2); |
|
|
125 LBL_T1:mp_clear (&t1); |
|
|
126 return res; |
|
|
127 } |
|
|
128 #endif |
