1 /*
  2  *  bignum_div.c
  3  *  BigMath
  4  *
  5  *  Created by Curtis Jones on 2007.08.05.
  6  *  Copyright 2007 __MyCompanyName__. All rights reserved.
  7  *
  8  */
  9 
 10 #include "bignum.h"
 11 
 12 // TODO: make sure that divs->oofm == 1
 13 //
 14 bignum_t *
 15 bignum_div_aux (bignum_t *quot, bignum_t *rmdr, bignum_t *divd, bignum_t *divs)
 16 {
 17   u_int64_t r = 0;
 18   u_int64_t n = divd->oofm;
 19   u_int64_t j = n - 1;
 20   u_int64_t b = divd->base;
 21   
 22   if (divd->sign == BIGNUM_NEGATIVE) {
 23     printf("%s().. divd is negative! ", __PRETTY_FUNCTION__); bignum_print(divd);
 24     return NULL;
 25   }
 26   
 27   if (divs->sign == BIGNUM_NEGATIVE) {
 28     printf("%s().. divs is negative! ", __PRETTY_FUNCTION__); bignum_print(divs);
 29     return NULL;
 30   }
 31   
 32   if (quot != NULL)
 33     __bignum_zero(quot);
 34   
 35   do {
 36     u_int64_t y = r * b + (u_int64_t)divd->numb[j];
 37     
 38     if (quot != NULL) {
 39       __bignum_set(quot, j, y / (u_int64_t)divs->numb[0]);
 40       
 41       if (quot->numb[j] != 0)
 42         quot->oofm = BIGNUM_MAX(quot->oofm, (j+1));
 43     }
 44     
 45     r = y % (u_int64_t)divs->numb[0];
 46   }
 47   while (--j >= 0 && j < n);
 48   
 49   if (rmdr != NULL) {
 50     __bignum_set(rmdr, 0, r);
 51     
 52     if (rmdr->numb[0] != 0)
 53       rmdr->oofm = 1;
 54   }
 55   
 56   return quot;
 57 }
 58 
 59 // Knuth 4.3.1 - Algorithm D
 60 //
 61 // a / b = c
 62 //
 63 // dividend (a) / divisor (b) = quotient (c)
 64 //
 65 bignum_t *
 66 bignum_div (bignum_t *quot, bignum_t *rmdr, bignum_t *divd, bignum_t *divs)
 67 {
 68   u_int64_t i = 0;
 69   u_int64_t m = divd->oofm - divs->oofm;
 70   u_int64_t n = divs->oofm;
 71   u_int64_t b = divd->base;
 72   u_int64_t d = ( b / ((u_int64_t)divs->numb[divs->oofm-1] + 1) );
 73   bignum_t *u = NULL, *v = NULL;
 74   u_int32_t divd_sign, divs_sign, quot_sign;
 75   
 76   // set divd to positive
 77   // set divs to positive
 78   //
 79   // if divd is negative xor divs is negative
 80   //   quot is negative
 81   //
 82   {
 83     if (BIGNUM_NEGATIVE == (divd_sign = divd->sign) ^ BIGNUM_NEGATIVE == (divs_sign = divs->sign))
 84       quot_sign = BIGNUM_NEGATIVE;
 85     else
 86       quot_sign = BIGNUM_POSITIVE;
 87     
 88     divd->sign = BIGNUM_POSITIVE;
 89     divs->sign = BIGNUM_POSITIVE;
 90   }
 91   
 92   if (quot != NULL)
 93     __bignum_zero( quot );
 94   
 95   if (rmdr != NULL)
 96     __bignum_zero( rmdr );
 97   
 98   // if the dividend is less than the divisor then the
 99   // quotient is zero and the remainder is the dividend.
100   //
101   if (bignum_lt(divd, divs)) {
102     if (quot != NULL)
103       bignum_set(quot, 0);
104     
105     if (rmdr != NULL)
106       bignum_set_bn(rmdr, divd);
107     
108     divd->sign = divd_sign;
109     divs->sign = divs_sign;
110     
111     return quot;
112   }
113   
114   // if the dividend is equal to the divisor then the
115   // quotient is one and the remainder is zero.
116   //
117   else if (bignum_eq(divd, divs)) {
118     if (quot != NULL)
119       bignum_set(quot, 1);
120     
121     if (rmdr != NULL)
122       bignum_set(rmdr, 0);
123     
124     divd->sign = divd_sign;
125     divs->sign = divs_sign;
126     
127     return quot;
128   }
129   
130   // TODO: It is possible that divs->oofm is greater than
131   // one, but that the value is still wholly contained in
132   // the first magnitude. we could test the value of divs
133   // relative to the maximum value that can be held in
134   // a single magnitude....
135   //
136   else if (divs->oofm == 1) {
137     bignum_div_aux(quot, rmdr, divd, divs);
138     
139     divd->sign = divd_sign;
140     divs->sign = divs_sign;
141     
142     return quot;
143   }
144   
145   if (NULL == (u = bignum_alloc())) {
146     printf("%s().. [1] bignum_alloc failed\n", __PRETTY_FUNCTION__);
147     return NULL;
148   }
149   
150   if (NULL == (v = bignum_alloc())) {
151     printf("%s().. [2] bignum_alloc failed\n", __PRETTY_FUNCTION__);
152     bignum_free( u );
153     return NULL;
154   }
155   
156   // STEP D1 - Normalize
157   //
158   // Set u = divisor * d
159   // Set v = dividend * d
160   //
161   // If u is still of the same magnitude, then add
162   // an additional zero to the end - the high end.
163   // This does not change its value, but the additional
164   // space is required.
165   //
166   {
167     bignum_t * tmpd = bignum_alloc();
168     
169     bignum_set(tmpd, d);
170     
171     bignum_mul(u, divd, tmpd);
172     bignum_mul(v, divs, tmpd);
173     
174     if (u->oofm == divd->oofm)
175       u->oofm++;
176     
177     bignum_free( tmpd );
178   }
179   
180   // STEP D2 - Initialize j
181   //
182   // Set j to m which is the difference in magnitude between
183   // the dividend and the divisor. We'll be decrementing j
184   // each time through this loop.
185   //
186   {
187     u_int64_t j = m;
188     
189     bignum_t * bn_q    = bignum_alloc();
190     bignum_t * bn_tmp1 = bignum_alloc();
191     bignum_t * bn_tmp2 = bignum_alloc();
192     bignum_t * bn_tmpu = bignum_alloc();
193     
194     // STEP D3 - Part 1 - Calculate q
195     //
196     // q is going to be the j'th digit of the quotient. first
197     // we estimate q and then possibly alter it slightly if 
198     // our estimate is off by one or two.
199     //
200     do {
201       u_int64_t x = (u_int64_t)u->numb[j+n] * b + (u_int64_t)u->numb[j+n-1];
202       u_int64_t q = x / (u_int64_t)v->numb[n-1];
203       u_int64_t r = x - (q * (u_int64_t)v->numb[n-1]);
204       u_int64_t o = 0;
205       
206       // STEP D3 - Part 2 - Verify q
207       //
208       // It is possible that the estimate of q is off. This
209       // code which I'm not about to explain, alters q as
210       // necessary.
211       //
212       for (i = 0; i < 2; ++i) {
213         if (q < b && q * (u_int64_t)v->numb[n-2] <= b * r + u->numb[j+n-2])
214           break;
215         
216         q -= 1;
217         r += (u_int64_t)v->numb[n-1];
218         
219         if (r >= b)
220           break;
221       }
222       
223       // STEP D4 - Multiply and Subtract
224       //
225       // Now that our q correct, multiply "digits" j through
226       // j + n of u by v. Then replace "digits" j through 
227       // j + n of u with the product. Convoluted subtraction.
228       //
229       {
230         int i = 0;
231         
232         bignum_set(bn_q, q);
233         bignum_mul(bn_tmp1, v, bn_q);
234         bignum_set(bn_tmp2, 0);
235         bignum_set(bn_tmpu, 0);
236         
237         for (i = 0; i <= n; ++i)
238           __bignum_set(bn_tmpu, bn_tmpu->oofm++, u->numb[j+i]);
239         
240         if (bignum_lt(bn_tmpu, bn_tmp1)) {
241           bignum_set(bn_q, --q);
242           bignum_mul(bn_tmp1, v, bn_q);
243         }
244         
245         if (bignum_lt(bn_tmpu, bn_tmp1))
246           printf("bignum::bignum_div().. warning: q is still too large.\n");
247         
248         o = bn_tmpu->oofm;
249         bignum_sub(bn_tmp2, bn_tmpu, bn_tmp1);
250         bignum_set_bn(bn_tmpu, bn_tmp2);
251         
252         for (i = 0; i < bn_tmpu->oofm; ++i)
253           __bignum_set(u, i+j, bn_tmpu->numb[i]);
254         
255         for (; i <= n; ++i)
256           __bignum_set(u, i+j, 0);
257         
258         u->oofm -= o - bn_tmpu->oofm;
259       }
260       
261       if (quot != NULL) {
262         __bignum_set(quot, j, q);
263         
264         if (quot->numb[j] != 0)
265           quot->oofm = BIGNUM_MAX(quot->oofm, j+1);
266       }
267     }
268     while (--j >= 0 && j <= m);
269     
270     bignum_free( bn_q    );
271     bignum_free( bn_tmp1 );
272     bignum_free( bn_tmp2 );
273     bignum_free( bn_tmpu );
274   }
275   
276   if (rmdr != NULL) {
277     bignum_t *bn_tmpd = bignum_alloc();
278     bignum_set(bn_tmpd, d);
279     bignum_div_aux(rmdr, NULL, u, bn_tmpd);
280     bignum_free( bn_tmpd );
281   }
282   
283   bignum_free( u );
284   bignum_free( v );
285   
286   divd->sign = divd_sign;
287   divs->sign = divs_sign;
288   
289   if (quot != NULL)
290     quot->sign = quot_sign;
291   
292   return quot;
293 }


syntax highlighted by Code2HTML, v. 0.9.1