Line data Source code
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 : */
28 0 : 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 0 : if ((b & 1) == 0 && a->sign == MP_NEG) {
35 0 : return MP_VAL;
36 : }
37 :
38 0 : if ((res = mp_init (&t1)) != MP_OKAY) {
39 0 : return res;
40 : }
41 :
42 0 : if ((res = mp_init (&t2)) != MP_OKAY) {
43 0 : goto LBL_T1;
44 : }
45 :
46 0 : if ((res = mp_init (&t3)) != MP_OKAY) {
47 0 : goto LBL_T2;
48 : }
49 :
50 : /* if a is negative fudge the sign but keep track */
51 0 : neg = a->sign;
52 0 : a->sign = MP_ZPOS;
53 :
54 : /* t2 = 2 */
55 0 : mp_set (&t2, 2);
56 :
57 : do {
58 : /* t1 = t2 */
59 0 : if ((res = mp_copy (&t2, &t1)) != MP_OKAY) {
60 0 : goto LBL_T3;
61 : }
62 :
63 : /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
64 :
65 : /* t3 = t1**(b-1) */
66 0 : if ((res = mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) {
67 0 : goto LBL_T3;
68 : }
69 :
70 : /* numerator */
71 : /* t2 = t1**b */
72 0 : if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) {
73 0 : goto LBL_T3;
74 : }
75 :
76 : /* t2 = t1**b - a */
77 0 : if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) {
78 0 : goto LBL_T3;
79 : }
80 :
81 : /* denominator */
82 : /* t3 = t1**(b-1) * b */
83 0 : if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) {
84 0 : goto LBL_T3;
85 : }
86 :
87 : /* t3 = (t1**b - a)/(b * t1**(b-1)) */
88 0 : if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) {
89 0 : goto LBL_T3;
90 : }
91 :
92 0 : if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) {
93 0 : goto LBL_T3;
94 : }
95 0 : } while (mp_cmp (&t1, &t2) != MP_EQ);
96 :
97 : /* result can be off by a few so check */
98 : for (;;) {
99 0 : if ((res = mp_expt_d (&t1, b, &t2)) != MP_OKAY) {
100 0 : goto LBL_T3;
101 : }
102 :
103 0 : if (mp_cmp (&t2, a) == MP_GT) {
104 0 : if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) {
105 0 : goto LBL_T3;
106 : }
107 : } else {
108 0 : break;
109 : }
110 : }
111 :
112 : /* reset the sign of a first */
113 0 : a->sign = neg;
114 :
115 : /* set the result */
116 0 : mp_exch (&t1, c);
117 :
118 : /* set the sign of the result */
119 0 : c->sign = neg;
120 :
121 0 : res = MP_OKAY;
122 :
123 0 : LBL_T3:mp_clear (&t3);
124 0 : LBL_T2:mp_clear (&t2);
125 0 : LBL_T1:mp_clear (&t1);
126 0 : 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 $ */
|