home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / LIBSRC.ZOO / libsrc / longlong / bdiv.c < prev    next >
C/C++ Source or Header  |  1992-02-22  |  4KB  |  186 lines
  1. #include "longlong.h"
  2.  
  3. static int bshift ();
  4.  
  5. /* Divide a by b, producing quotient q and remainder r.
  6.  
  7.        sizeof a is m
  8.        sizeof b is n
  9.        sizeof q is m - n
  10.        sizeof r is n
  11.  
  12.    The quotient must fit in m - n bytes, i.e., the most significant
  13.    n digits of a must be less than b, and m must be greater than n.  */
  14.  
  15. /* The name of this used to be __div_internal,
  16.    but that is too long for SYSV.  */
  17.  
  18. void 
  19. __bdiv (a, b, q, r, m, n)
  20.      unsigned short *a, *b, *q, *r;
  21.      size_t m, n;
  22. {
  23.   unsigned long qhat, rhat;
  24.   unsigned long acc;
  25.   unsigned short *u = (unsigned short *) alloca (m);
  26.   unsigned short *v = (unsigned short *) alloca (n);
  27.   unsigned short *u0, *u1, *u2;
  28.   unsigned short *v0;
  29.   int d, qn;
  30.   int i, j;
  31.  
  32.   m /= sizeof *a;
  33.   n /= sizeof *b;
  34.   qn = m - n;
  35.  
  36.   /* Remove leading zero digits from divisor, and the same number of
  37.      digits (which must be zero) from dividend.  */
  38.  
  39.   while (b[big_end (n)] == 0)
  40.     {
  41.       r[big_end (n)] = 0;
  42.  
  43.       a += little_end (2);
  44.       b += little_end (2);
  45.       r += little_end (2);
  46.       m--;
  47.       n--;
  48.  
  49.       /* Check for zero divisor.  */
  50.       if (n == 0)
  51.     abort ();
  52.     }
  53.       
  54.   /* If divisor is a single digit, do short division.  */
  55.  
  56.   if (n == 1)
  57.     {
  58.       acc = a[big_end (m)];
  59.       a += little_end (2);
  60.       for (j = big_end (qn); is_not_lsd (j, qn); j = next_lsd (j))
  61.     {
  62.       acc = (acc << 16) | a[j];
  63.       q[j] = acc / *b;
  64.       acc = acc % *b;
  65.     }
  66.       *r = acc;
  67.       return;
  68.     }
  69.  
  70.   /* No such luck, must do long division. Shift divisor and dividend
  71.      left until the high bit of the divisor is 1.  */
  72.  
  73.   for (d = 0; d < 16; d++)
  74.     if (b[big_end (n)] & (1 << (16 - 1 - d)))
  75.       break;
  76.  
  77.   bshift (a, d, u, 0, m);
  78.   bshift (b, d, v, 0, n);
  79.  
  80.   /* Get pointers to the high dividend and divisor digits.  */
  81.  
  82.   u0 = u + big_end (m) - big_end (qn);
  83.   u1 = next_lsd (u0);
  84.   u2 = next_lsd (u1);
  85.   u += little_end (2);
  86.  
  87.   v0 = v + big_end (n);
  88.  
  89.   /* Main loop: find a quotient digit, multiply it by the divisor,
  90.      and subtract that from the dividend, shifted over the right amount. */
  91.  
  92.   for (j = big_end (qn); is_not_lsd (j, qn); j = next_lsd (j))
  93.     {
  94.       /* Quotient digit initial guess: high 2 dividend digits over high
  95.      divisor digit.  */
  96.  
  97.       if (u0[j] == *v0)
  98.     {
  99.       qhat = B - 1;
  100.       rhat = (unsigned long) *v0 + u1[j];
  101.     }
  102.       else
  103.     {
  104.       unsigned long numerator = ((unsigned long) u0[j] << 16) | u1[j];
  105.       qhat = numerator / *v0;
  106.       rhat = numerator % *v0;
  107.     }
  108.  
  109.       /* Now get the quotient right for high 3 dividend digits over
  110.      high 2 divisor digits.  */
  111.  
  112.       while (rhat < B && qhat * *next_lsd (v0) > ((rhat << 16) | u2[j]))
  113.     {
  114.       qhat -= 1;
  115.       rhat += *v0;
  116.     }
  117.         
  118.       /* Multiply quotient by divisor, subtract from dividend.  */
  119.  
  120.       acc = 0;
  121.       for (i = little_end (n); is_not_msd (i, n); i = next_msd (i))
  122.     {
  123.       acc += (unsigned long) (u + j)[i] - v[i] * qhat;
  124.       (u + j)[i] = acc & low16;
  125.       if (acc < B)
  126.         acc = 0;
  127.       else
  128.         acc = (acc >> 16) | -B;
  129.     }
  130.  
  131.       q[j] = qhat;
  132.  
  133.       /* Quotient may have been too high by 1.  If dividend went negative,
  134.      decrement the quotient by 1 and add the divisor back.  */
  135.  
  136.       if ((signed long) (acc + u0[j]) < 0)
  137.     {
  138.       q[j] -= 1;
  139.       acc = 0;
  140.       for (i = little_end (n); is_not_msd (i, n); i = next_msd (i))
  141.         {
  142.           acc += (unsigned long) (u + j)[i] + v[i];
  143.           (u + j)[i] = acc & low16;
  144.           acc = acc >> 16;
  145.         }
  146.     }
  147.     }
  148.  
  149.   /* Now the remainder is what's left of the dividend, shifted right
  150.      by the amount of the normalizing left shift at the top.  */
  151.  
  152.   r[big_end (n)] = bshift (u + 1 + little_end (j - 1),
  153.                16 - d,
  154.                r + little_end (2),
  155.                u[little_end (m - 1)] >> d,
  156.                n - 1);
  157. }
  158.  
  159. /* Left shift U by K giving W; fill the introduced low-order bits with
  160.    CARRY_IN.  Length of U and W is N.  Return carry out.  K must be
  161.    in 0 .. 16.  */
  162.  
  163. static int
  164. bshift (u, k, w, carry_in, n)
  165.      unsigned short *u, *w, carry_in;
  166.      int k, n;
  167. {
  168.   unsigned long acc;
  169.   int i;
  170.  
  171.   if (k == 0)
  172.     {
  173.       bcopy (u, w, n * sizeof *u);
  174.       return 0;
  175.     }
  176.  
  177.   acc = carry_in;
  178.   for (i = little_end (n); is_not_msd (i, n); i = next_msd (i))
  179.     {
  180.       acc |= (unsigned long) u[i] << k;
  181.       w[i] = acc & low16;
  182.       acc = acc >> 16;
  183.     }
  184.   return acc;
  185. }
  186.