s_mp_div_school.c

   1 #include "tommath_private.h"
   2 #ifdef S_MP_DIV_SCHOOL_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 /* integer signed division.
   7  * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
   8  * HAC pp.598 Algorithm 14.20
   9  *
  10  * Note that the description in HAC is horribly
  11  * incomplete.  For example, it doesn't consider
  12  * the case where digits are removed from 'x' in
  13  * the inner loop.  It also doesn't consider the
  14  * case that y has fewer than three digits, etc..
  15  *
  16  * The overall algorithm is as described as
  17  * 14.20 from HAC but fixed to treat these cases.
  18 */
  19 mp_err s_mp_div_school(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
  20 {
  21    mp_int q, x, y, t1, t2;
  22    mp_digit xdpi;
  23    int n, t, i, norm;
  24    bool neg;
  25    mp_err err;
  26
  27    if ((err = mp_init_size(&q, a->used + 2)) != MP_OKAY) {
  28       return err;
  29    }
  30    q.used = a->used + 2;
  31
  32    if ((err = mp_init(&t1)) != MP_OKAY)                           goto LBL_Q;
  33    if ((err = mp_init(&t2)) != MP_OKAY)                           goto LBL_T1;
  34    if ((err = mp_init_copy(&x, a)) != MP_OKAY)                    goto LBL_T2;
  35    if ((err = mp_init_copy(&y, b)) != MP_OKAY)                    goto LBL_X;
  36
  37    /* fix the sign */
  38    neg = (a->sign != b->sign);
  39    x.sign = y.sign = MP_ZPOS;
  40
  41    /* normalize both x and y, ensure that y >= b/2, [b == 2**MP_DIGIT_BIT] */
  42    norm = mp_count_bits(&y) % MP_DIGIT_BIT;
  43    if (norm < (MP_DIGIT_BIT - 1)) {
  44       norm = (MP_DIGIT_BIT - 1) - norm;
  45       if ((err = mp_mul_2d(&x, norm, &x)) != MP_OKAY)             goto LBL_Y;
  46       if ((err = mp_mul_2d(&y, norm, &y)) != MP_OKAY)             goto LBL_Y;
  47    } else {
  48       norm = 0;
  49    }
  50
  51    /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
  52    n = x.used - 1;
  53    t = y.used - 1;
  54
  55    /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
  56    /* y = y*b**{n-t} */
  57    if ((err = mp_lshd(&y, n - t)) != MP_OKAY)                     goto LBL_Y;
  58
  59    while (mp_cmp(&x, &y) != MP_LT) {
  60       ++(q.dp[n - t]);
  61       if ((err = mp_sub(&x, &y, &x)) != MP_OKAY)                  goto LBL_Y;
  62    }
  63
  64    /* reset y by shifting it back down */
  65    mp_rshd(&y, n - t);
  66
  67    /* step 3. for i from n down to (t + 1) */
  68    for (i = n; i >= (t + 1); i--) {
  69       if (i > x.used) {
  70          continue;
  71       }
  72       /* Do not assume that more than enough memory is automatically allocated and set to '0' */
  73       xdpi = (i == x.used) ? 0u : x.dp[i];
  74
  75       /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
  76        * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
  77       if (xdpi == y.dp[t]) {
  78          q.dp[(i - t) - 1] = ((mp_digit)1 << (mp_digit)MP_DIGIT_BIT) - (mp_digit)1;
  79       } else {
  80          mp_word tmp;
  81          tmp = (mp_word)xdpi << (mp_word)MP_DIGIT_BIT;
  82          tmp |= (mp_word)x.dp[i - 1];
  83          tmp /= (mp_word)y.dp[t];
  84          if (tmp > (mp_word)MP_MASK) {
  85             tmp = MP_MASK;
  86          }
  87          q.dp[(i - t) - 1] = (mp_digit)(tmp & (mp_word)MP_MASK);
  88       }
  89
  90       /* while (q{i-t-1} * (yt * b + y{t-1})) >
  91                xi * b**2 + xi-1 * b + xi-2
  92
  93          do q{i-t-1} -= 1;
  94       */
  95       q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1uL) & (mp_digit)MP_MASK;
  96       do {
  97          q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & (mp_digit)MP_MASK;
  98
  99          /* find left hand */
 100          mp_zero(&t1);
 101          t1.dp[0] = ((t - 1) < 0) ? 0u : y.dp[t - 1];
 102          t1.dp[1] = y.dp[t];
 103          t1.used = 2;
 104          if ((err = mp_mul_d(&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY)   goto LBL_Y;
 105
 106          /* find right hand */
 107          t2.dp[0] = ((i - 2) < 0) ? 0u : x.dp[i - 2];
 108          t2.dp[1] = x.dp[i - 1]; /* i >= 1 always holds */
 109          t2.dp[2] = xdpi;
 110          t2.used = 3;
 111       } while (mp_cmp_mag(&t1, &t2) == MP_GT);
 112
 113       /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
 114       if ((err = mp_mul_d(&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY)       goto LBL_Y;
 115       if ((err = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY)                  goto LBL_Y;
 116       if ((err = mp_sub(&x, &t1, &x)) != MP_OKAY)                        goto LBL_Y;
 117
 118       /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
 119       if (mp_isneg(&x)) {
 120          if ((err = mp_copy(&y, &t1)) != MP_OKAY)                        goto LBL_Y;
 121          if ((err = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY)               goto LBL_Y;
 122          if ((err = mp_add(&x, &t1, &x)) != MP_OKAY)                     goto LBL_Y;
 123
 124          q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1uL) & MP_MASK;
 125       }
 126    }
 127
 128    /* now q is the quotient and x is the remainder
 129     * [which we have to normalize]
 130     */
 131
 132    /* get sign before writing to c */
 133    x.sign = mp_iszero(&x) ? MP_ZPOS : a->sign;
 134
 135    if (c != NULL) {
 136       mp_clamp(&q);
 137       mp_exch(&q, c);
 138       c->sign = (neg ? MP_NEG : MP_ZPOS);
 139    }
 140
 141    if (d != NULL) {
 142       if ((err = mp_div_2d(&x, norm, &x, NULL)) != MP_OKAY)       goto LBL_Y;
 143       mp_exch(&x, d);
 144    }
 145
 146 LBL_Y:
 147    mp_clear(&y);
 148 LBL_X:
 149    mp_clear(&x);
 150 LBL_T2:
 151    mp_clear(&t2);
 152 LBL_T1:
 153    mp_clear(&t1);
 154 LBL_Q:
 155    mp_clear(&q);
 156    return err;
 157 }
 158
 159 #endif