fp_div.c

00001 /* TomsFastMath, a fast ISO C bignum library.
00002  * 
00003  * This project is meant to fill in where LibTomMath
00004  * falls short.  That is speed ;-)
00005  *
00006  * This project is public domain and free for all purposes.
00007  * 
00008  * Tom St Denis, tomstdenis@gmail.com
00009  */
00010 #include <tfm.h>
00011 
00012 /* a/b => cb + d == a */
00013 int fp_div(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
00014 {
00015   fp_int  q, x, y, t1, t2;
00016   int     n, t, i, norm, neg;
00017 
00018   /* is divisor zero ? */
00019   if (fp_iszero (b) == 1) {
00020     return FP_VAL;
00021   }
00022 
00023   /* if a < b then q=0, r = a */
00024   if (fp_cmp_mag (a, b) == FP_LT) {
00025     if (d != NULL) {
00026       fp_copy (a, d);
00027     } 
00028     if (c != NULL) {
00029       fp_zero (c);
00030     }
00031     return FP_OKAY;
00032   }
00033 
00034   fp_init(&q);
00035   q.used = a->used + 2;
00036 
00037   fp_init(&t1);
00038   fp_init(&t2);
00039   fp_init_copy(&x, a);
00040   fp_init_copy(&y, b);
00041 
00042   /* fix the sign */
00043   neg = (a->sign == b->sign) ? FP_ZPOS : FP_NEG;
00044   x.sign = y.sign = FP_ZPOS;
00045 
00046   /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
00047   norm = fp_count_bits(&y) % DIGIT_BIT;
00048   if (norm < (int)(DIGIT_BIT-1)) {
00049      norm = (DIGIT_BIT-1) - norm;
00050      fp_mul_2d (&x, norm, &x);
00051      fp_mul_2d (&y, norm, &y);
00052   } else {
00053      norm = 0;
00054   }
00055 
00056   /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
00057   n = x.used - 1;
00058   t = y.used - 1;
00059 
00060   /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
00061   fp_lshd (&y, n - t);                                             /* y = y*b**{n-t} */
00062 
00063   while (fp_cmp (&x, &y) != FP_LT) {
00064     ++(q.dp[n - t]);
00065     fp_sub (&x, &y, &x);
00066   }
00067 
00068   /* reset y by shifting it back down */
00069   fp_rshd (&y, n - t);
00070 
00071   /* step 3. for i from n down to (t + 1) */
00072   for (i = n; i >= (t + 1); i--) {
00073     if (i > x.used) {
00074       continue;
00075     }
00076 
00077     /* step 3.1 if xi == yt then set q{i-t-1} to b-1, 
00078      * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
00079     if (x.dp[i] == y.dp[t]) {
00080       q.dp[i - t - 1] = ((((fp_word)1) << DIGIT_BIT) - 1);
00081     } else {
00082       fp_word tmp;
00083       tmp = ((fp_word) x.dp[i]) << ((fp_word) DIGIT_BIT);
00084       tmp |= ((fp_word) x.dp[i - 1]);
00085       tmp /= ((fp_word) y.dp[t]);
00086       q.dp[i - t - 1] = (fp_digit) (tmp);
00087     }
00088 
00089     /* while (q{i-t-1} * (yt * b + y{t-1})) > 
00090              xi * b**2 + xi-1 * b + xi-2 
00091      
00092        do q{i-t-1} -= 1; 
00093     */
00094     q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
00095     do {
00096       q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
00097 
00098       /* find left hand */
00099       fp_zero (&t1);
00100       t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
00101       t1.dp[1] = y.dp[t];
00102       t1.used = 2;
00103       fp_mul_d (&t1, q.dp[i - t - 1], &t1);
00104 
00105       /* find right hand */
00106       t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
00107       t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
00108       t2.dp[2] = x.dp[i];
00109       t2.used = 3;
00110     } while (fp_cmp_mag(&t1, &t2) == FP_GT);
00111 
00112     /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
00113     fp_mul_d (&y, q.dp[i - t - 1], &t1);
00114     fp_lshd  (&t1, i - t - 1);
00115     fp_sub   (&x, &t1, &x);
00116 
00117     /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
00118     if (x.sign == FP_NEG) {
00119       fp_copy (&y, &t1);
00120       fp_lshd (&t1, i - t - 1);
00121       fp_add (&x, &t1, &x);
00122       q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
00123     }
00124   }
00125 
00126   /* now q is the quotient and x is the remainder 
00127    * [which we have to normalize] 
00128    */
00129   
00130   /* get sign before writing to c */
00131   x.sign = x.used == 0 ? FP_ZPOS : a->sign;
00132 
00133   if (c != NULL) {
00134     fp_clamp (&q);
00135     fp_copy (&q, c);
00136     c->sign = neg;
00137   }
00138 
00139   if (d != NULL) {
00140     fp_div_2d (&x, norm, &x, NULL);
00141 
00142 /* the following is a kludge, essentially we were seeing the right remainder but 
00143    with excess digits that should have been zero
00144  */
00145     for (i = b->used; i < x.used; i++) {
00146         x.dp[i] = 0;
00147     }
00148     fp_clamp(&x);
00149     fp_copy (&x, d);
00150   }
00151 
00152   return FP_OKAY;
00153 }
00154 
00155 /* $Source: /cvs/libtom/tomsfastmath/src/divide/fp_div.c,v $ */
00156 /* $Revision: 1.1 $ */
00157 /* $Date: 2006/12/31 21:25:53 $ */