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