home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / LINUX / MATH_EMU.ZIP / MATH_EMU / FPU_DIV.C < prev    next >
Encoding:
C/C++ Source or Header  |  1979-12-31  |  9.6 KB  |  251 lines
  1. /*        $NetBSD: fpu_div.c,v 1.2 1994/11/20 20:52:38 deraadt Exp $ */
  2.  
  3. /*
  4.  * Copyright (c) 1992, 1993
  5.  *        The Regents of the University of California.  All rights reserved.
  6.  *
  7.  * This software was developed by the Computer Systems Engineering group
  8.  * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
  9.  * contributed to Berkeley.
  10.  *
  11.  * All advertising materials mentioning features or use of this software
  12.  * must display the following acknowledgement:
  13.  *        This product includes software developed by the University of
  14.  *        California, Lawrence Berkeley Laboratory.
  15.  *
  16.  * Redistribution and use in source and binary forms, with or without
  17.  * modification, are permitted provided that the following conditions
  18.  * are met:
  19.  * 1. Redistributions of source code must retain the above copyright
  20.  *    notice, this list of conditions and the following disclaimer.
  21.  * 2. Redistributions in binary form must reproduce the above copyright
  22.  *    notice, this list of conditions and the following disclaimer in the
  23.  *    documentation and/or other materials provided with the distribution.
  24.  * 3. All advertising materials mentioning features or use of this software
  25.  *    must display the following acknowledgement:
  26.  *        This product includes software developed by the University of
  27.  *        California, Berkeley and its contributors.
  28.  * 4. Neither the name of the University nor the names of its contributors
  29.  *    may be used to endorse or promote products derived from this software
  30.  *    without specific prior written permission.
  31.  *
  32.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  33.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  34.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  35.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  36.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  37.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  38.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  39.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  40.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  41.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  42.  * SUCH DAMAGE.
  43.  *
  44.  *        @(#)fpu_div.c       8.1 (Berkeley) 6/11/93
  45.  */
  46.  
  47. /*
  48.  * Perform an FPU divide (return x / y).
  49.  */
  50.  
  51. #include "types.h"
  52.  
  53. #include "reg.h"
  54.  
  55. #include "fpu_arit.h"
  56. #include "fpu_emul.h"
  57.  
  58. /*
  59.  * Division of normal numbers is done as follows:
  60.  *
  61.  * x and y are floating point numbers, i.e., in the form 1.bbbb * 2^e.
  62.  * If X and Y are the mantissas (1.bbbb's), the quotient is then:
  63.  *
  64.  *        q = (X / Y) * 2^((x exponent) - (y exponent))
  65.  *
  66.  * Since X and Y are both in [1.0,2.0), the quotient's mantissa (X / Y)
  67.  * will be in [0.5,2.0).  Moreover, it will be less than 1.0 if and only
  68.  * if X < Y.  In that case, it will have to be shifted left one bit to
  69.  * become a normal number, and the exponent decremented.  Thus, the
  70.  * desired exponent is:
  71.  *
  72.  *        left_shift = x->fp_mant < y->fp_mant;
  73.  *        result_exp = x->fp_exp - y->fp_exp - left_shift;
  74.  *
  75.  * The quotient mantissa X/Y can then be computed one bit at a time
  76.  * using the following algorithm:
  77.  *
  78.  *        Q = 0;                        -- Initial quotient.
  79.  *        R = X;                        -- Initial remainder,
  80.  *        if (left_shift)               --   but fixed up in advance.
  81.  *                  R *= 2;
  82.  *        for (bit = FP_NMANT; --bit >= 0; R *= 2) {
  83.  *                  if (R >= Y) {
  84.  *                            Q |= 1 << bit;
  85.  *                            R -= Y;
  86.  *                  }
  87.  *        }
  88.  *
  89.  * The subtraction R -= Y always removes the uppermost bit from R (and
  90.  * can sometimes remove additional lower-order 1 bits); this proof is
  91.  * left to the reader.
  92.  *
  93.  * This loop correctly calculates the guard and round bits since they are
  94.  * included in the expanded internal representation.  The sticky bit
  95.  * is to be set if and only if any other bits beyond guard and round
  96.  * would be set.  From the above it is obvious that this is true if and
  97.  * only if the remainder R is nonzero when the loop terminates.
  98.  *
  99.  * Examining the loop above, we can see that the quotient Q is built
  100.  * one bit at a time ``from the top down''.  This means that we can
  101.  * dispense with the multi-word arithmetic and just build it one word
  102.  * at a time, writing each result word when it is done.
  103.  *
  104.  * Furthermore, since X and Y are both in [1.0,2.0), we know that,
  105.  * initially, R >= Y.  (Recall that, if X < Y, R is set to X * 2 and
  106.  * is therefore at in [2.0,4.0).)  Thus Q is sure to have bit FP_NMANT-1
  107.  * set, and R can be set initially to either X - Y (when X >= Y) or
  108.  * 2X - Y (when X < Y).  In addition, comparing R and Y is difficult,
  109.  * so we will simply calculate R - Y and see if that underflows.
  110.  * This leads to the following revised version of the algorithm:
  111.  *
  112.  *        R = X;
  113.  *        bit = FP_1;
  114.  *        D = R - Y;
  115.  *        if (D >= 0) {
  116.  *                  result_exp = x->fp_exp - y->fp_exp;
  117.  *                  R = D;
  118.  *                  q = bit;
  119.  *                  bit >>= 1;
  120.  *        } else {
  121.  *                  result_exp = x->fp_exp - y->fp_exp - 1;
  122.  *                  q = 0;
  123.  *        }
  124.  *        R <<= 1;
  125.  *        do  {
  126.  *                  D = R - Y;
  127.  *                  if (D >= 0) {
  128.  *                            q |= bit;
  129.  *                            R = D;
  130.  *                  }
  131.  *                  R <<= 1;
  132.  *        } while ((bit >>= 1) != 0);
  133.  *        Q[0] = q;
  134.  *        for (i = 1; i < 4; i++) {
  135.  *                  q = 0, bit = 1 << 31;
  136.  *                  do {
  137.  *                            D = R - Y;
  138.  *                            if (D >= 0) {
  139.  *                                      q |= bit;
  140.  *                                      R = D;
  141.  *                            }
  142.  *                            R <<= 1;
  143.  *                  } while ((bit >>= 1) != 0);
  144.  *                  Q[i] = q;
  145.  *        }
  146.  *
  147.  * This can be refined just a bit further by moving the `R <<= 1'
  148.  * calculations to the front of the do-loops and eliding the first one.
  149.  * The process can be terminated immediately whenever R becomes 0, but
  150.  * this is relatively rare, and we do not bother.
  151.  */
  152.  
  153. struct fpn *
  154. fpu_div(fe)
  155.           register struct fpemu *fe;
  156. {
  157.           register struct fpn *x = &fe->fe_f1, *y = &fe->fe_f2;
  158.           register u_int q, bit;
  159.           register u_int r0, r1, r2, r3, d0, d1, d2, d3, y0, y1, y2, y3;
  160.           FPU_DECL_CARRY
  161.  
  162.           fe->fe_fpsr &= ~FPSR_EXCP; /* clear all exceptions */
  163.  
  164.           /*
  165.            * Since divide is not commutative, we cannot just use ORDER.
  166.            * Check either operand for NaN first; if there is at least one,
  167.            * order the signalling one (if only one) onto the right, then
  168.            * return it.  Otherwise we have the following cases:
  169.            *
  170.            *        Inf / Inf = NaN, plus NV exception
  171.            *        Inf / num = Inf [i.e., return x]
  172.            *        Inf / 0   = Inf [i.e., return x]
  173.            *        0 / Inf = 0 [i.e., return x]
  174.            *        0 / num = 0 [i.e., return x]
  175.            *        0 / 0   = NaN, plus NV exception
  176.            *        num / Inf = 0
  177.            *        num / num = num (do the divide)
  178.            *        num / 0   = Inf, plus DZ exception
  179.            */
  180.           if (ISNAN(x) || ISNAN(y)) {
  181.                     ORDER(x, y);
  182.                     return (y);
  183.           }
  184.           if (ISINF(x) || ISZERO(x)) {
  185.                     if (x->fp_class == y->fp_class)
  186.                               return (fpu_newnan(fe));
  187.                     return (x);
  188.           }
  189.  
  190.           /* all results at this point use XOR of operand signs */
  191.           x->fp_sign ^= y->fp_sign;
  192.           if (ISINF(y)) {
  193.                     x->fp_class = FPC_ZERO;
  194.                     return (x);
  195.           }
  196.           if (ISZERO(y)) {
  197.                     fe->fe_fpsr |= FPSR_DZ;
  198.                     x->fp_class = FPC_INF;
  199.                     return (x);
  200.           }
  201.  
  202.           /*
  203.            * Macros for the divide.  See comments at top for algorithm.
  204.            * Note that we expand R, D, and Y here.
  205.            */
  206.  
  207. #define   SUBTRACT /* D = R - Y */ FPU_SUBS(d3, r3, y3); FPU_SUBCS(d2, r2, y2); FPU_SUBCS(d1, r1, y1); FPU_SUBC(d0, r0, y0)
  208.  
  209. #define   NONNEGATIVE /* D >= 0 */ ((int)d0 >= 0)
  210.  
  211. #ifdef FPU_SHL1_BY_ADD
  212. #define   SHL1 /* R <<= 1 */ FPU_ADDS(r3, r3, r3); FPU_ADDCS(r2, r2, r2); FPU_ADDCS(r1, r1, r1); FPU_ADDC(r0, r0, r0)
  213. #else
  214. #define SHL1 r0 = (r0 << 1) | (r1 >> 31), r1 = (r1 << 1) | (r2 >> 31), r2 = (r2 << 1) | (r3 >> 31), r3 <<= 1
  215. #endif
  216.  
  217. #define   LOOP /* do ... while (bit >>= 1) */ do { SHL1; SUBTRACT; if (NONNEGATIVE) { q |= bit; r0 = d0, r1 = d1, r2 = d2, r3 = d3; } } while ((bit >>= 1) != 0)
  218.  
  219. #define   WORD(r, i) /* calculate r->fp_mant[i] */ q = 0; bit = 1 << 31; LOOP; (x)->fp_mant[i] = q
  220.  
  221.           /* Setup.  Note that we put our result in x. */
  222.           r0 = x->fp_mant[0];
  223.           r1 = x->fp_mant[1];
  224.           r2 = x->fp_mant[2];
  225.           r3 = x->fp_mant[3];
  226.           y0 = y->fp_mant[0];
  227.           y1 = y->fp_mant[1];
  228.           y2 = y->fp_mant[2];
  229.           y3 = y->fp_mant[3];
  230.  
  231.           bit = FP_1;
  232.           SUBTRACT;
  233.           if (NONNEGATIVE) {
  234.                     x->fp_exp -= y->fp_exp;
  235.                     r0 = d0, r1 = d1, r2 = d2, r3 = d3;
  236.                     q = bit;
  237.                     bit >>= 1;
  238.           } else {
  239.                     x->fp_exp -= y->fp_exp + 1;
  240.                     q = 0;
  241.           }
  242.           LOOP;
  243.           x->fp_mant[0] = q;
  244.           WORD(x, 1);
  245.           WORD(x, 2);
  246.           WORD(x, 3);
  247.           x->fp_sticky = r0 | r1 | r2 | r3;
  248.  
  249.           return (x);
  250. }
  251.