1 /*
2 * bignum_root.c
3 * BigMath
4 *
5 * Created by Curtis Jones on 2008.03.19.
6 * Copyright 2008 __MyCompanyName__. All rights reserved.
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