release/src/router/dropbear/libtommath/bn_mp_div.c

   1 #include <tommath.h>
   2 #ifdef BN_MP_DIV_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
   4  *
   5  * LibTomMath is a library that provides multiple-precision
   6  * integer arithmetic as well as number theoretic functionality.
   7  *
   8  * The library was designed directly after the MPI library by
   9  * Michael Fromberger but has been written from scratch with
  10  * additional optimizations in place.
  11  *
  12  * The library is free for all purposes without any express
  13  * guarantee it works.
  14  *
  15  * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
  16  */
  17
  18 #ifdef BN_MP_DIV_SMALL
  19
  20 /* slower bit-bang division... also smaller */
  21 int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d)
  22 {
  23    mp_int ta, tb, tq, q;
  24    int    res, n, n2;
  25
  26   /* is divisor zero ? */
  27   if (mp_iszero (b) == 1) {
  28     return MP_VAL;
  29   }
  30
  31   /* if a < b then q=0, r = a */
  32   if (mp_cmp_mag (a, b) == MP_LT) {
  33     if (d != NULL) {
  34       res = mp_copy (a, d);
  35     } else {
  36       res = MP_OKAY;
  37     }
  38     if (c != NULL) {
  39       mp_zero (c);
  40     }
  41     return res;
  42   }
  43
  44   /* init our temps */
  45   if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) {
  46      return res;
  47   }
  48
  49
  50   mp_set(&tq, 1);
  51   n = mp_count_bits(a) - mp_count_bits(b);
  52   if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
  53       ((res = mp_abs(b, &tb)) != MP_OKAY) ||
  54       ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
  55       ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
  56       goto LBL_ERR;
  57   }
  58
  59   while (n-- >= 0) {
  60      if (mp_cmp(&tb, &ta) != MP_GT) {
  61         if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
  62             ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) {
  63            goto LBL_ERR;
  64         }
  65      }
  66      if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
  67          ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
  68            goto LBL_ERR;
  69      }
  70   }
  71
  72   /* now q == quotient and ta == remainder */
  73   n  = a->sign;
  74   n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG);
  75   if (c != NULL) {
  76      mp_exch(c, &q);
  77      c->sign  = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
  78   }
  79   if (d != NULL) {
  80      mp_exch(d, &ta);
  81      d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
  82   }
  83 LBL_ERR:
  84    mp_clear_multi(&ta, &tb, &tq, &q, NULL);
  85    return res;
  86 }
  87
  88 #else
  89
  90 /* integer signed division.
  91  * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
  92  * HAC pp.598 Algorithm 14.20
  93  *
  94  * Note that the description in HAC is horribly
  95  * incomplete.  For example, it doesn't consider
  96  * the case where digits are removed from 'x' in
  97  * the inner loop.  It also doesn't consider the
  98  * case that y has fewer than three digits, etc..
  99  *
 100  * The overall algorithm is as described as
 101  * 14.20 from HAC but fixed to treat these cases.
 102 */
 103 int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
 104 {
 105   mp_int  q, x, y, t1, t2;
 106   int     res, n, t, i, norm, neg;
 107
 108   /* is divisor zero ? */
 109   if (mp_iszero (b) == 1) {
 110     return MP_VAL;
 111   }
 112
 113   /* if a < b then q=0, r = a */
 114   if (mp_cmp_mag (a, b) == MP_LT) {
 115     if (d != NULL) {
 116       res = mp_copy (a, d);
 117     } else {
 118       res = MP_OKAY;
 119     }
 120     if (c != NULL) {
 121       mp_zero (c);
 122     }
 123     return res;
 124   }
 125
 126   if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) {
 127     return res;
 128   }
 129   q.used = a->used + 2;
 130
 131   if ((res = mp_init (&t1)) != MP_OKAY) {
 132     goto LBL_Q;
 133   }
 134
 135   if ((res = mp_init (&t2)) != MP_OKAY) {
 136     goto LBL_T1;
 137   }
 138
 139   if ((res = mp_init_copy (&x, a)) != MP_OKAY) {
 140     goto LBL_T2;
 141   }
 142
 143   if ((res = mp_init_copy (&y, b)) != MP_OKAY) {
 144     goto LBL_X;
 145   }
 146
 147   /* fix the sign */
 148   neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
 149   x.sign = y.sign = MP_ZPOS;
 150
 151   /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
 152   norm = mp_count_bits(&y) % DIGIT_BIT;
 153   if (norm < (int)(DIGIT_BIT-1)) {
 154      norm = (DIGIT_BIT-1) - norm;
 155      if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) {
 156        goto LBL_Y;
 157      }
 158      if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) {
 159        goto LBL_Y;
 160      }
 161   } else {
 162      norm = 0;
 163   }
 164
 165   /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
 166   n = x.used - 1;
 167   t = y.used - 1;
 168
 169   /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
 170   if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
 171     goto LBL_Y;
 172   }
 173
 174   while (mp_cmp (&x, &y) != MP_LT) {
 175     ++(q.dp[n - t]);
 176     if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) {
 177       goto LBL_Y;
 178     }
 179   }
 180
 181   /* reset y by shifting it back down */
 182   mp_rshd (&y, n - t);
 183
 184   /* step 3. for i from n down to (t + 1) */
 185   for (i = n; i >= (t + 1); i--) {
 186     if (i > x.used) {
 187       continue;
 188     }
 189
 190     /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
 191      * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
 192     if (x.dp[i] == y.dp[t]) {
 193       q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
 194     } else {
 195       mp_word tmp;
 196       tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
 197       tmp |= ((mp_word) x.dp[i - 1]);
 198       tmp /= ((mp_word) y.dp[t]);
 199       if (tmp > (mp_word) MP_MASK)
 200         tmp = MP_MASK;
 201       q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
 202     }
 203
 204     /* while (q{i-t-1} * (yt * b + y{t-1})) >
 205              xi * b**2 + xi-1 * b + xi-2
 206
 207        do q{i-t-1} -= 1;
 208     */
 209     q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
 210     do {
 211       q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;
 212
 213       /* find left hand */
 214       mp_zero (&t1);
 215       t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
 216       t1.dp[1] = y.dp[t];
 217       t1.used = 2;
 218       if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {
 219         goto LBL_Y;
 220       }
 221
 222       /* find right hand */
 223       t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
 224       t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
 225       t2.dp[2] = x.dp[i];
 226       t2.used = 3;
 227     } while (mp_cmp_mag(&t1, &t2) == MP_GT);
 228
 229     /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
 230     if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
 231       goto LBL_Y;
 232     }
 233
 234     if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
 235       goto LBL_Y;
 236     }
 237
 238     if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) {
 239       goto LBL_Y;
 240     }
 241
 242     /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
 243     if (x.sign == MP_NEG) {
 244       if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
 245         goto LBL_Y;
 246       }
 247       if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
 248         goto LBL_Y;
 249       }
 250       if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
 251         goto LBL_Y;
 252       }
 253
 254       q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK;
 255     }
 256   }
 257
 258   /* now q is the quotient and x is the remainder
 259    * [which we have to normalize]
 260    */
 261
 262   /* get sign before writing to c */
 263   x.sign = x.used == 0 ? MP_ZPOS : a->sign;
 264
 265   if (c != NULL) {
 266     mp_clamp (&q);
 267     mp_exch (&q, c);
 268     c->sign = neg;
 269   }
 270
 271   if (d != NULL) {
 272     if ((res = mp_div_2d (&x, norm, &x, NULL)) != MP_OKAY) {
 273                 goto LBL_Y;
 274         }
 275     mp_exch (&x, d);
 276   }
 277
 278   res = MP_OKAY;
 279
 280 LBL_Y:mp_clear (&y);
 281 LBL_X:mp_clear (&x);
 282 LBL_T2:mp_clear (&t2);
 283 LBL_T1:mp_clear (&t1);
 284 LBL_Q:mp_clear (&q);
 285   return res;
 286 }
 287
 288 #endif
 289
 290 #endif
 291
 292 /* $Source: /cvs/libtom/libtommath/bn_mp_div.c,v $ */
 293 /* $Revision: 1.3 $ */
 294 /* $Date: 2006/03/31 14:18:44 $ */