pd.if.org Git - zpackage/blob - tomsfastmath/src/divide/fp_div.c

   1 /* TomsFastMath, a fast ISO C bignum library.
   2  *
   3  * This project is meant to fill in where LibTomMath
   4  * falls short.  That is speed ;-)
   5  *
   6  * This project is public domain and free for all purposes.
   7  *
   8  * Tom St Denis, tomstdenis@gmail.com
   9  */
  10 #include <tfm_private.h>
  11
  12 /* a/b => cb + d == a */
  13 int fp_div(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
  14 {
  15   fp_int  q, x, y, t1, t2;
  16   int     n, t, i, norm, neg;
  17
  18   /* is divisor zero ? */
  19   if (fp_iszero (b) == 1) {
  20     return FP_VAL;
  21   }
  22
  23   /* if a < b then q=0, r = a */
  24   if (fp_cmp_mag (a, b) == FP_LT) {
  25     if (d != NULL) {
  26       fp_copy (a, d);
  27     }
  28     if (c != NULL) {
  29       fp_zero (c);
  30     }
  31     return FP_OKAY;
  32   }
  33
  34   fp_init(&q);
  35   q.used = a->used + 2;
  36
  37   fp_init(&t1);
  38   fp_init(&t2);
  39   fp_init_copy(&x, a);
  40   fp_init_copy(&y, b);
  41
  42   /* fix the sign */
  43   neg = (a->sign == b->sign) ? FP_ZPOS : FP_NEG;
  44   x.sign = y.sign = FP_ZPOS;
  45
  46   /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
  47   norm = fp_count_bits(&y) % DIGIT_BIT;
  48   if (norm < (int)(DIGIT_BIT-1)) {
  49      norm = (DIGIT_BIT-1) - norm;
  50      fp_mul_2d (&x, norm, &x);
  51      fp_mul_2d (&y, norm, &y);
  52   } else {
  53      norm = 0;
  54   }
  55
  56   /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
  57   n = x.used - 1;
  58   t = y.used - 1;
  59
  60   /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
  61   fp_lshd (&y, n - t);                                             /* y = y*b**{n-t} */
  62
  63   while (fp_cmp (&x, &y) != FP_LT) {
  64     ++(q.dp[n - t]);
  65     fp_sub (&x, &y, &x);
  66   }
  67
  68   /* reset y by shifting it back down */
  69   fp_rshd (&y, n - t);
  70
  71   /* step 3. for i from n down to (t + 1) */
  72   for (i = n; i >= (t + 1); i--) {
  73     if (i > x.used) {
  74       continue;
  75     }
  76
  77     /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
  78      * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
  79     if (x.dp[i] == y.dp[t]) {
  80       q.dp[i - t - 1] = ((((fp_word)1) << DIGIT_BIT) - 1);
  81     } else {
  82       fp_word tmp;
  83       tmp = ((fp_word) x.dp[i]) << ((fp_word) DIGIT_BIT);
  84       tmp |= ((fp_word) x.dp[i - 1]);
  85       tmp /= ((fp_word) y.dp[t]);
  86       q.dp[i - t - 1] = (fp_digit) (tmp);
  87     }
  88
  89     /* while (q{i-t-1} * (yt * b + y{t-1})) >
  90              xi * b**2 + xi-1 * b + xi-2
  91
  92        do q{i-t-1} -= 1;
  93     */
  94     q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
  95     do {
  96       q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
  97
  98       /* find left hand */
  99       fp_zero (&t1);
 100       t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
 101       t1.dp[1] = y.dp[t];
 102       t1.used = 2;
 103       fp_mul_d (&t1, q.dp[i - t - 1], &t1);
 104
 105       /* find right hand */
 106       t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
 107       t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
 108       t2.dp[2] = x.dp[i];
 109       t2.used = 3;
 110     } while (fp_cmp_mag(&t1, &t2) == FP_GT);
 111
 112     /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
 113     fp_mul_d (&y, q.dp[i - t - 1], &t1);
 114     fp_lshd  (&t1, i - t - 1);
 115     fp_sub   (&x, &t1, &x);
 116
 117     /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
 118     if (x.sign == FP_NEG) {
 119       fp_copy (&y, &t1);
 120       fp_lshd (&t1, i - t - 1);
 121       fp_add (&x, &t1, &x);
 122       q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
 123     }
 124   }
 125
 126   /* now q is the quotient and x is the remainder
 127    * [which we have to normalize]
 128    */
 129
 130   /* get sign before writing to c */
 131   x.sign = x.used == 0 ? FP_ZPOS : a->sign;
 132
 133   if (c != NULL) {
 134     fp_clamp (&q);
 135     fp_copy (&q, c);
 136     c->sign = neg;
 137   }
 138
 139   if (d != NULL) {
 140     fp_div_2d (&x, norm, &x, NULL);
 141
 142 /* the following is a kludge, essentially we were seeing the right remainder but
 143    with excess digits that should have been zero
 144  */
 145     for (i = b->used; i < x.used; i++) {
 146         x.dp[i] = 0;
 147     }
 148     fp_clamp(&x);
 149     fp_copy (&x, d);
 150   }
 151
 152   return FP_OKAY;
 153 }
 154
 155 /* $Source$ */
 156 /* $Revision$ */
 157 /* $Date$ */