1 /* $NetBSD: bn_mp_n_root.c,v 1.1.1.2 2014/04/24 12:45:31 pettai 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 */
mp_n_root(mp_int * a,mp_digit b,mp_int * c)30 int 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
125 LBL_T3:mp_clear (&t3);
126 LBL_T2:mp_clear (&t2);
127 LBL_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