```  1 /*
2  *  bignum_root.c
3  *  BigMath
4  *
5  *  Created by Curtis Jones on 2008.03.19.
7  *
8  */
9
10 #include "bignum.h"
11
12 /**
13  * Sets root equal to the nth root of num. Returns root.
14  *
15  * http://en.wikipedia.org/wiki/Nth_root_algorithm
16  *
17  * 1 / n *
18  *   (n - 1) * x_sub_k +
19  *   num / (x_sub_k ^ (n - 1))
20  *
21  * n /
22  *   (n - 1) * x_sub_k +
23  *   num / (x_sub_k ^ (n - 1))
24  *
25  * TODO - how do we pick a good first x_sub_k? it seems that using a multi-bit
26  *        shift is the way to go for that. we might be fine with a hard-coded
27  *        initial x_sub_k for now. adding floating point support might change
28  *        that.
29  *
30  * TODO - special cases:
31  *   if num is zero...
32  *   if num is negative...
33  *   if n is zero...
34  *   if n is one... root = num
35  *   if n is two...
36  *   if n is negative...
37  *
38  */
39 bignum_t *
40 bignum_root (bignum_t *root, bignum_t *n, bignum_t *num)
41 {
42   int first_pass = 0;
43
44   bignum_t * result  = NULL;
45   bignum_t * x_sub_k = NULL;
46   bignum_t * n_min_1 = NULL;
47   bignum_t * part_1  = NULL;
48   bignum_t * part_2  = NULL;
49   bignum_t * part_3  = NULL;
50   bignum_t * part_4  = NULL;
51   bignum_t * part_5  = NULL;
52   bignum_t * part_6  = NULL;
53
54   if (NULL == (x_sub_k = bignum_alloc()))
55     goto finish;
56
57   if (NULL == (n_min_1 = bignum_alloc()))
58     goto finish;
59
60   if (NULL == (part_1 = bignum_alloc()))
61     goto finish;
62
63   if (NULL == (part_2 = bignum_alloc()))
64     goto finish;
65
66   if (NULL == (part_3 = bignum_alloc()))
67     goto finish;
68
69   if (NULL == (part_4 = bignum_alloc()))
70     goto finish;
71
72   if (NULL == (part_5 = bignum_alloc()))
73     goto finish;
74
75   if (NULL == (part_6 = bignum_alloc()))
76     goto finish;
77
78   bignum_set(x_sub_k, 2);
79
80   // n __-__ 1
81   bignum_sub(n_min_1, n, bignum_one());
82
83   int i = 0;
84
85   // we loop forever, and hope that we'll detect that we've found the best
86   // possible answer and break out of the loop. otherwise, there's a short-
87   // circuit that'll stop this from looping infinitely.
88   //
89   while (1)
90   {
91     // short circuit in case something goes wrong
92     if (i++ > 1000)
93       break;
94
95     // (n - 1) __*__ x_sub_k
96     bignum_mul(part_1, n_min_1, x_sub_k);
97
98     // x_sub_k __^__ (n - 1)
99     bignum_pow(part_2, x_sub_k, n_min_1);
100
101     // num __/__ (x_sub_k ^ (n - 1))
102     bignum_div(part_3, NULL, num, part_2);
103
104     // ((n - 1) * x_sub_k) __+__ (num / (x_sub_k ^ (n - 1)))
105     bignum_add(part_4, part_1, part_3);
106
107     // [...] __/__ n
108     bignum_div(part_5, part_6, part_4, n);
109
110     // if the result is zero then something has gone very, very wrong.
111     if (bignum_eq(part_5, bignum_zero()))
112       break;
113
114     // if we're getting a repeating answer, then we have arrived.
115     else if (bignum_eq(x_sub_k, part_5))
116       break;
117
118     // this process should always result in a lower and lower answer, slowly
119     // approaching the correct root value. however, if the result is higher
120     // than it was last time, then we have used a value for x_sub_k that is
121     // less than the correct root value.
122     //
123     // if this is the first time through this process, or this also happened
124     // last time, then multiply x_sub_k by one hundred and try again.
125     //
126     // if this is not the first time through and this did not also happen
127     // last time, then we have gone as far as we can go. the current x_sub_k is
128     // the answer we're sticking with.
129     //
130     else if (bignum_lt(x_sub_k, part_5)) {
131       if (first_pass || i == 1) {
132         bignum_mul(part_6, x_sub_k, bignum_one_hundred());
133         bignum_set_bn(x_sub_k, part_6);
134         first_pass++;
135       }
136       else {
137         bignum_t * tmp = part_5;
138         part_5 = x_sub_k;
139         x_sub_k = tmp;
140         break;
141       }
142     }
143     else {
144       first_pass = 0;
145       bignum_set_bn(x_sub_k, part_5);
146     }
147   }
148
149   bignum_set_bn(root, part_5);
150
151   result = root;
152
153 finish:
154   if (x_sub_k != NULL)
155     bignum_free(x_sub_k);
156
157   if (n_min_1 != NULL)
158     bignum_free(n_min_1);
159
160   if (part_1 != NULL)
161     bignum_free(part_1);
162
163   if (part_2 != NULL)
164     bignum_free(part_2);
165
166   if (part_3 != NULL)
167     bignum_free(part_3);
168
169   if (part_4 != NULL)
170     bignum_free(part_4);
171
172   if (part_5 != NULL)
173     bignum_free(part_5);
174
175   if (part_6 != NULL)
176     bignum_free(part_6);
177
178   return result;
179 }

```

syntax highlighted by Code2HTML, v. 0.9.1