modules/src/flt_arith/flt_div.c

   1 /*
   2   (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3   See the copyright notice in the ACK home directory, in the file "Copyright".
   4 */
   5
   6 /* $Id: flt_div.c,v 1.13 1997/03/13 18:38:24 ceriel Exp $ */
   7
   8 #include "flt_misc.h"
   9
  10 void
  11 flt_div(e1,e2,e3)
  12         register flt_arith *e1,*e2,*e3;
  13 {
  14         long    result[2];
  15         register long   *lp;
  16         unsigned short u[9], v[5];
  17         register int j;
  18         register unsigned short *u_p = u;
  19         int maxv = 4;
  20         flt_arith cpe1 = *e1;
  21
  22         flt_status = 0;
  23         e3->flt_sign = e1->flt_sign ^ e2->flt_sign;
  24         if ((e2->m1 | e2->m2) == 0) {
  25                 flt_status = FLT_DIV0;
  26                 e3->flt_exp = EXT_MAX;
  27                 e3->m1 = 0x80000000;
  28                 e3->m2 = 0L;
  29                 return;
  30         }
  31         if ((e1->m1 | e1->m2) == 0) {   /* 0 / anything == 0 */
  32                 e3->flt_exp = 0;        /* make sure */
  33                 e3->m1 = e3->m2 = 0L;
  34                 e3->flt_sign = 0;
  35                 return;
  36         }
  37         e3->flt_exp = e1->flt_exp - e2->flt_exp;
  38
  39         u[4] = (e1->m2 & 1) << 15;
  40         flt_b64_sft(&(cpe1.flt_mantissa),1);
  41         flt_split(&cpe1, u);
  42         u[5] = 0; u[6] = 0; u[7] = 0;
  43         flt_split(e2, &v[1]);
  44         while (! v[maxv]) maxv--;
  45         result[0] = 0;
  46         result[1] = 0;
  47         lp = result;
  48
  49         /*
  50          * Use an algorithm of Knuth (The art of programming, Seminumerical
  51          * algorithms), to divide u by v. u and v are both seen as numbers
  52          * with base 65536.
  53          */
  54         for (j = 0; j <= 3; j++, u_p++) {
  55                 long q_est, temp;
  56                 long v1 = v[1];
  57
  58                 if (j == 2) lp++;
  59                 if (u_p[0] == 0 && u_p[1] < v[1]) continue;
  60                 temp = ((long)u_p[0] << 16) + u_p[1];
  61                 if (u_p[0] >= v[1]) {
  62                         q_est = 0x0000FFFFL;
  63                 }
  64                 else if (v[1] == 1) {
  65                         q_est = temp;
  66                 }
  67                 else if (temp >= 0) {
  68                         q_est = temp / v1;
  69                 }
  70                 else {
  71                         long rem;
  72                         q_est = (0x7FFFFFFF/v1)+((temp&0x7FFFFFFF)/v1);
  73                         rem = (0x7FFFFFFF%v1)+((temp&0x7FFFFFFF)%v1)+1;
  74                         while (rem >= v1) {
  75                                 q_est++;
  76                                 rem -= v1;
  77                         }
  78                 }
  79                 temp -= q_est * v1;
  80                 while (!(temp&0xFFFF0000) &&
  81                        ucmp((long)v[2]*q_est,(temp<<16)+(long)u_p[2]) > 0) {
  82                         q_est--;
  83                         temp += v1;
  84                 }
  85                 /*      Now, according to Knuth, we have an estimate of the
  86                         quotient, that is either correct or one too big, but
  87                         almost always correct.
  88                 */
  89                 if (q_est != 0)  {
  90                         int i;
  91                         long k = 0;
  92                         int borrow = 0;
  93
  94                         for (i = maxv; i > 0; i--) {
  95                                 long tmp = q_est * (long)v[i] + k + borrow;
  96                                 unsigned short md = tmp & 0xFFFF;
  97
  98                                 borrow = (md > u_p[i]);
  99                                 u_p[i] -= md;
 100                                 k = (tmp >> 16) & 0xFFFF;
 101                         }
 102                         k += borrow;
 103                         borrow = (long)u_p[0] < k;
 104                         u_p[0] -= k;
 105
 106                         if (borrow) {
 107                                 /* So, this does not happen often; the estimate
 108                                    was one too big; correct this
 109                                 */
 110                                 q_est--;
 111                                 borrow = 0;
 112                                 for (i = maxv; i > 0; i--) {
 113                                         long tmp
 114                                             = v[i]+(long)u_p[i]+borrow;
 115
 116                                         u_p[i] = tmp & 0xFFFF;
 117                                         borrow = (tmp >> 16) & 0xFFFF;
 118                                 }
 119                                 u_p[0] += borrow;
 120                         }
 121                         *lp |= (j & 1) ? q_est : (q_est<<16);
 122                 }
 123         }
 124         e3->m1 = result[0];
 125         e3->m2 = result[1];
 126
 127         flt_nrm(e3);
 128         flt_chk(e3);
 129 }