userapps/opensource/sshd/libtommath/bn_mp_div.c

   1 /* LibTomMath, multiple-precision integer library -- Tom St Denis
   2  *
   3  * LibTomMath is library that provides for multiple-precision
   4  * integer arithmetic as well as number theoretic functionality.
   5  *
   6  * The library is designed directly after the MPI library by
   7  * Michael Fromberger but has been written from scratch with
   8  * additional optimizations in place.
   9  *
  10  * The library is free for all purposes without any express
  11  * guarantee it works.
  12  *
  13  * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
  14  */
  15 #include <tommath.h>
  16
  17 /* integer signed division.
  18  * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
  19  * HAC pp.598 Algorithm 14.20
  20  *
  21  * Note that the description in HAC is horribly
  22  * incomplete.  For example, it doesn't consider
  23  * the case where digits are removed from 'x' in
  24  * the inner loop.  It also doesn't consider the
  25  * case that y has fewer than three digits, etc..
  26  *
  27  * The overall algorithm is as described as
  28  * 14.20 from HAC but fixed to treat these cases.
  29 */
  30 int
  31 mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
  32 {
  33   mp_int  q, x, y, t1, t2;
  34   int     res, n, t, i, norm, neg;
  35
  36   /* is divisor zero ? */
  37   if (mp_iszero (b) == 1) {
  38     return MP_VAL;
  39   }
  40
  41   /* if a < b then q=0, r = a */
  42   if (mp_cmp_mag (a, b) == MP_LT) {
  43     if (d != NULL) {
  44       res = mp_copy (a, d);
  45     } else {
  46       res = MP_OKAY;
  47     }
  48     if (c != NULL) {
  49       mp_zero (c);
  50     }
  51     return res;
  52   }
  53
  54   if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) {
  55     return res;
  56   }
  57   q.used = a->used + 2;
  58
  59   if ((res = mp_init (&t1)) != MP_OKAY) {
  60     goto __Q;
  61   }
  62
  63   if ((res = mp_init (&t2)) != MP_OKAY) {
  64     goto __T1;
  65   }
  66
  67   if ((res = mp_init_copy (&x, a)) != MP_OKAY) {
  68     goto __T2;
  69   }
  70
  71   if ((res = mp_init_copy (&y, b)) != MP_OKAY) {
  72     goto __X;
  73   }
  74
  75   /* fix the sign */
  76   neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
  77   x.sign = y.sign = MP_ZPOS;
  78
  79   /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
  80   norm = mp_count_bits(&y) % DIGIT_BIT;
  81   if (norm < (int)(DIGIT_BIT-1)) {
  82      norm = (DIGIT_BIT-1) - norm;
  83      if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) {
  84        goto __Y;
  85      }
  86      if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) {
  87        goto __Y;
  88      }
  89   } else {
  90      norm = 0;
  91   }
  92
  93   /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
  94   n = x.used - 1;
  95   t = y.used - 1;
  96
  97   /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
  98   if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
  99     goto __Y;
 100   }
 101
 102   while (mp_cmp (&x, &y) != MP_LT) {
 103     ++(q.dp[n - t]);
 104     if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) {
 105       goto __Y;
 106     }
 107   }
 108
 109   /* reset y by shifting it back down */
 110   mp_rshd (&y, n - t);
 111
 112   /* step 3. for i from n down to (t + 1) */
 113   for (i = n; i >= (t + 1); i--) {
 114     if (i > x.used)
 115       continue;
 116
 117     /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
 118      * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
 119     if (x.dp[i] == y.dp[t]) {
 120       q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
 121     } else {
 122       mp_word tmp;
 123       tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
 124       tmp |= ((mp_word) x.dp[i - 1]);
 125       tmp /= ((mp_word) y.dp[t]);
 126       if (tmp > (mp_word) MP_MASK)
 127         tmp = MP_MASK;
 128       q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
 129     }
 130
 131     /* while (q{i-t-1} * (yt * b + y{t-1})) >
 132              xi * b**2 + xi-1 * b + xi-2
 133
 134        do q{i-t-1} -= 1;
 135     */
 136     q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
 137     do {
 138       q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;
 139
 140       /* find left hand */
 141       mp_zero (&t1);
 142       t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
 143       t1.dp[1] = y.dp[t];
 144       t1.used = 2;
 145       if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {
 146         goto __Y;
 147       }
 148
 149       /* find right hand */
 150       t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
 151       t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
 152       t2.dp[2] = x.dp[i];
 153       t2.used = 3;
 154     } while (mp_cmp_mag(&t1, &t2) == MP_GT);
 155
 156     /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
 157     if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
 158       goto __Y;
 159     }
 160
 161     if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
 162       goto __Y;
 163     }
 164
 165     if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) {
 166       goto __Y;
 167     }
 168
 169     /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
 170     if (x.sign == MP_NEG) {
 171       if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
 172         goto __Y;
 173       }
 174       if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
 175         goto __Y;
 176       }
 177       if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
 178         goto __Y;
 179       }
 180
 181       q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK;
 182     }
 183   }
 184
 185   /* now q is the quotient and x is the remainder
 186    * [which we have to normalize]
 187    */
 188
 189   /* get sign before writing to c */
 190   x.sign = a->sign;
 191
 192   if (c != NULL) {
 193     mp_clamp (&q);
 194     mp_exch (&q, c);
 195     c->sign = neg;
 196   }
 197
 198   if (d != NULL) {
 199     mp_div_2d (&x, norm, &x, NULL);
 200     mp_exch (&x, d);
 201   }
 202
 203   res = MP_OKAY;
 204
 205 __Y:mp_clear (&y);
 206 __X:mp_clear (&x);
 207 __T2:mp_clear (&t2);
 208 __T1:mp_clear (&t1);
 209 __Q:mp_clear (&q);
 210   return res;
 211 }