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_IMPL.C < prev    next >
Encoding:
C/C++ Source or Header  |  1979-12-31  |  17.5 KB  |  487 lines
  1. /*        $NetBSD: fpu_implode.c,v 1.2 1994/11/20 20:52:42 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_implode.c   8.1 (Berkeley) 6/11/93
  45.  */
  46.  
  47. /*
  48.  * FPU subroutines: `implode' internal format numbers into the machine's
  49.  * `packed binary' format.
  50.  */
  51.  
  52. #include "types.h"
  53.  
  54. #include "ieee.h"
  55. #include "reg.h"
  56.  
  57. #include "fpu_emul.h"
  58. #include "fpu_arit.h"
  59.  
  60. /* Conversion from internal format -- note asymmetry. */
  61. static u_int        fpu_ftoi __P((struct fpemu *fe, struct fpn *fp));
  62. static u_int        fpu_ftos __P((struct fpemu *fe, struct fpn *fp));
  63. static u_int        fpu_ftod __P((struct fpemu *fe, struct fpn *fp, u_int *));
  64. static u_int        fpu_ftox __P((struct fpemu *fe, struct fpn *fp, u_int *));
  65.  
  66. /*
  67.  * Round a number (algorithm from Motorola MC68882 manual, modified for
  68.  * our internal format).  Set inexact exception if rounding is required.
  69.  * Return true iff we rounded up.
  70.  *
  71.  * After rounding, we discard the guard and round bits by shifting right
  72.  * 2 bits (a la fpu_shr(), but we do not bother with fp->fp_sticky).
  73.  * This saves effort later.
  74.  *
  75.  * Note that we may leave the value 2.0 in fp->fp_mant; it is the caller's
  76.  * responsibility to fix this if necessary.
  77.  */
  78. int
  79. round(register struct fpemu *fe, register struct fpn *fp)
  80. {
  81.           register u_int m0, m1, m2, m3;
  82.           register int gr, s, ret;
  83.  
  84.           m0 = fp->fp_mant[0];
  85.           m1 = fp->fp_mant[1];
  86.           m2 = fp->fp_mant[2];
  87.           m3 = fp->fp_mant[3];
  88.           gr = m3 & 3;
  89.           s = fp->fp_sticky;
  90.  
  91.           /* mant >>= FP_NG */
  92.           m3 = (m3 >> FP_NG) | (m2 << (32 - FP_NG));
  93.           m2 = (m2 >> FP_NG) | (m1 << (32 - FP_NG));
  94.           m1 = (m1 >> FP_NG) | (m0 << (32 - FP_NG));
  95.           m0 >>= FP_NG;
  96.  
  97.           if ((gr | s) == 0)  /* result is exact: no rounding needed */
  98.                     goto rounddown;
  99.  
  100.           fe->fe_fpsr |= FPSR_INEX2;    /* inexact */
  101.  
  102.           /* Go to rounddown to round down; break to round up. */
  103.           switch (fe->fe_fpcr & FPCR_ROUND) {
  104.  
  105.           case FPCR_NEAR:
  106.           default:
  107.                     /*
  108.                      * Round only if guard is set (gr & 2).  If guard is set,
  109.                      * but round & sticky both clear, then we want to round
  110.                      * but have a tie, so round to even, i.e., add 1 iff odd.
  111.                      */
  112.                     if ((gr & 2) == 0)
  113.                               goto rounddown;
  114.                     if ((gr & 1) || fp->fp_sticky || (m3 & 1))
  115.                               break;
  116.                     goto rounddown;
  117.  
  118.           case FPCR_ZERO:
  119.                     /* Round towards zero, i.e., down. */
  120.                     goto rounddown;
  121.  
  122.           case FPCR_MINF:
  123.                     /* Round towards -Inf: up if negative, down if positive. */
  124.                     if (fp->fp_sign)
  125.                               break;
  126.                     goto rounddown;
  127.  
  128.           case FPCR_PINF:
  129.                     /* Round towards +Inf: up if positive, down otherwise. */
  130.                     if (!fp->fp_sign)
  131.                               break;
  132.                     goto rounddown;
  133.           }
  134.  
  135.           /* Bump low bit of mantissa, with carry. */
  136. #ifdef sparc /* ``cheating'' (left out FPU_DECL_CARRY; know this is faster) */
  137.           FPU_ADDS(m3, m3, 1);
  138.           FPU_ADDCS(m2, m2, 0);
  139.           FPU_ADDCS(m1, m1, 0);
  140.           FPU_ADDC(m0, m0, 0);
  141. #else
  142.           if (++m3 == 0 && ++m2 == 0 && ++m1 == 0)
  143.                     m0++;
  144. #endif
  145.           fp->fp_mant[0] = m0;
  146.           fp->fp_mant[1] = m1;
  147.           fp->fp_mant[2] = m2;
  148.           fp->fp_mant[3] = m3;
  149.           return (1);
  150.  
  151. rounddown:
  152.           fp->fp_mant[0] = m0;
  153.           fp->fp_mant[1] = m1;
  154.           fp->fp_mant[2] = m2;
  155.           fp->fp_mant[3] = m3;
  156.           return (0);
  157. }
  158.  
  159. /*
  160.  * For overflow: return true if overflow is to go to +/-Inf, according
  161.  * to the sign of the overflowing result.  If false, overflow is to go
  162.  * to the largest magnitude value instead.
  163.  */
  164. static int
  165. toinf(struct fpemu *fe, int sign)
  166. {
  167.           int inf;
  168.  
  169.           /* look at rounding direction */
  170.           switch (fe->fe_fpcr & FPCR_ROUND) {
  171.  
  172.           default:
  173.           case FPCR_NEAR:               /* the nearest value is always Inf */
  174.                     inf = 1;
  175.                     break;
  176.  
  177.           case FPCR_ZERO:               /* toward 0 => never towards Inf */
  178.                     inf = 0;
  179.                     break;
  180.  
  181.           case FPCR_PINF:               /* toward +Inf iff positive */
  182.                     inf = (sign == 0);
  183.                     break;
  184.  
  185.           case FPCR_MINF:               /* toward -Inf iff negative */
  186.                     inf = sign;
  187.                     break;
  188.           }
  189.           return (inf);
  190. }
  191.  
  192. /*
  193.  * fpn -> int (int value returned as return value).
  194.  *
  195.  * N.B.: this conversion always rounds towards zero (this is a peculiarity
  196.  * of the SPARC instruction set).
  197.  */
  198. static u_int
  199. fpu_ftoi(fe, fp)
  200.           struct fpemu *fe;
  201.           register struct fpn *fp;
  202. {
  203.           register u_int i;
  204.           register int sign, exp;
  205.  
  206.           sign = fp->fp_sign;
  207.           switch (fp->fp_class) {
  208.  
  209.           case FPC_ZERO:
  210.                     return (0);
  211.  
  212.           case FPC_NUM:
  213.                     /*
  214.                      * If exp >= 2^32, overflow.  Otherwise shift value right
  215.                      * into last mantissa word (this will not exceed 0xffffffff),
  216.                      * shifting any guard and round bits out into the sticky
  217.                      * bit.  Then ``round'' towards zero, i.e., just set an
  218.                      * inexact exception if sticky is set (see round()).
  219.                      * If the result is > 0x80000000, or is positive and equals
  220.                      * 0x80000000, overflow; otherwise the last fraction word
  221.                      * is the result.
  222.                      */
  223.                     if ((exp = fp->fp_exp) >= 32)
  224.                               break;
  225.                     /* NB: the following includes exp < 0 cases */
  226.                     if (fpu_shr(fp, FP_NMANT - 1 - FP_NG - exp) != 0)
  227.                               /* m68881/2 do not underflow when
  228.                                  converting to integer */;
  229.                     round(fe, fp);
  230.                     i = fp->fp_mant[3];
  231.                     if (i >= ((u_int)0x80000000 + sign))
  232.                               break;
  233.                     return (sign ? -i : i);
  234.  
  235.           default:            /* Inf, qNaN, sNaN */
  236.                     break;
  237.           }
  238.           /* overflow: replace any inexact exception with invalid */
  239.           fe->fe_fpsr = (fe->fe_fpsr & ~FPSR_INEX2) | FPSR_OPERR;
  240.           return (0x7fffffff + sign);
  241. }
  242.  
  243. /*
  244.  * fpn -> single (32 bit single returned as return value).
  245.  * We assume <= 29 bits in a single-precision fraction (1.f part).
  246.  */
  247. static u_int
  248. fpu_ftos(fe, fp)
  249.           struct fpemu *fe;
  250.           register struct fpn *fp;
  251. {
  252.           register u_int sign = fp->fp_sign << 31;
  253.           register int exp;
  254.  
  255. #define   SNG_EXP(e)          ((e) << SNG_FRACBITS)         /* makes e an exponent */
  256. #define   SNG_MASK  (SNG_EXP(1) - 1)    /* mask for fraction */
  257.  
  258.           /* Take care of non-numbers first. */
  259.           if (ISNAN(fp)) {
  260.                     /*
  261.                      * Preserve upper bits of NaN, per SPARC V8 appendix N.
  262.                      * Note that fp->fp_mant[0] has the quiet bit set,
  263.                      * even if it is classified as a signalling NaN.
  264.                      */
  265.                     (void) fpu_shr(fp, FP_NMANT - 1 - SNG_FRACBITS);
  266.                     exp = SNG_EXP_INFNAN;
  267.                     goto done;
  268.           }
  269.           if (ISINF(fp))
  270.                     return (sign | SNG_EXP(SNG_EXP_INFNAN));
  271.           if (ISZERO(fp))
  272.                     return (sign);
  273.  
  274.           /*
  275.            * Normals (including subnormals).  Drop all the fraction bits
  276.            * (including the explicit ``implied'' 1 bit) down into the
  277.            * single-precision range.  If the number is subnormal, move
  278.            * the ``implied'' 1 into the explicit range as well, and shift
  279.            * right to introduce leading zeroes.  Rounding then acts
  280.            * differently for normals and subnormals: the largest subnormal
  281.            * may round to the smallest normal (1.0 x 2^minexp), or may
  282.            * remain subnormal.  In the latter case, signal an underflow
  283.            * if the result was inexact or if underflow traps are enabled.
  284.            *
  285.            * Rounding a normal, on the other hand, always produces another
  286.            * normal (although either way the result might be too big for
  287.            * single precision, and cause an overflow).  If rounding a
  288.            * normal produces 2.0 in the fraction, we need not adjust that
  289.            * fraction at all, since both 1.0 and 2.0 are zero under the
  290.            * fraction mask.
  291.            *
  292.            * Note that the guard and round bits vanish from the number after
  293.            * rounding.
  294.            */
  295.           if ((exp = fp->fp_exp + SNG_EXP_BIAS) <= 0) {     /* subnormal */
  296.                     /* -NG for g,r; -SNG_FRACBITS-exp for fraction */
  297.                     (void) fpu_shr(fp, FP_NMANT - FP_NG - SNG_FRACBITS - exp);
  298.                     if (round(fe, fp) && fp->fp_mant[3] == SNG_EXP(1))
  299.                               return (sign | SNG_EXP(1) | 0);
  300.                     if (fe->fe_fpsr & FPSR_INEX2)
  301.                               /* mc68881/2 don't underflow when converting */;
  302.                     return (sign | SNG_EXP(0) | fp->fp_mant[3]);
  303.           }
  304.           /* -FP_NG for g,r; -1 for implied 1; -SNG_FRACBITS for fraction */
  305.           (void) fpu_shr(fp, FP_NMANT - FP_NG - 1 - SNG_FRACBITS);
  306. #ifdef DIAGNOSTIC
  307.           if ((fp->fp_mant[3] & SNG_EXP(1 << FP_NG)) == 0)
  308.                   /*  panic("fpu_ftos"); */
  309. #endif
  310.           if (round(fe, fp) && fp->fp_mant[3] == SNG_EXP(2))
  311.                     exp++;
  312.           if (exp >= SNG_EXP_INFNAN) {
  313.                     /* overflow to inf or to max single */
  314.                     fe->fe_fpsr |= FPSR_OPERR | FPSR_INEX2;
  315.                     if (toinf(fe, sign))
  316.                               return (sign | SNG_EXP(SNG_EXP_INFNAN));
  317.                     return (sign | SNG_EXP(SNG_EXP_INFNAN - 1) | SNG_MASK);
  318.           }
  319. done:
  320.           /* phew, made it */
  321.           return (sign | SNG_EXP(exp) | (fp->fp_mant[3] & SNG_MASK));
  322. }
  323.  
  324. /*
  325.  * fpn -> double (32 bit high-order result returned; 32-bit low order result
  326.  * left in res[1]).  Assumes <= 61 bits in double precision fraction.
  327.  *
  328.  * This code mimics fpu_ftos; see it for comments.
  329.  */
  330. static u_int
  331. fpu_ftod(fe, fp, res)
  332.           struct fpemu *fe;
  333.           register struct fpn *fp;
  334.           u_int *res;
  335. {
  336.           register u_int sign = fp->fp_sign << 31;
  337.           register int exp;
  338.  
  339. #define   DBL_EXP(e)          ((e) << (DBL_FRACBITS & 31))
  340. #define   DBL_MASK  (DBL_EXP(1) - 1)
  341.  
  342.           if (ISNAN(fp)) {
  343.                     (void) fpu_shr(fp, FP_NMANT - 1 - DBL_FRACBITS);
  344.                     exp = DBL_EXP_INFNAN;
  345.                     goto done;
  346.           }
  347.           if (ISINF(fp)) {
  348.                     sign |= DBL_EXP(DBL_EXP_INFNAN);
  349.                     res[1] = 0;
  350.                     return (sign);
  351.           }
  352.           if (ISZERO(fp)) {
  353.                     res[1] = 0;
  354.                     return (sign);
  355.           }
  356.  
  357.           if ((exp = fp->fp_exp + DBL_EXP_BIAS) <= 0) {
  358.                     (void) fpu_shr(fp, FP_NMANT - FP_NG - DBL_FRACBITS - exp);
  359.                     if (round(fe, fp) && fp->fp_mant[2] == DBL_EXP(1)) {
  360.                               res[1] = 0;
  361.                               return (sign | DBL_EXP(1) | 0);
  362.                     }
  363.                     if (fe->fe_fpsr & FPSR_INEX2)
  364.                               /* mc68881/2 don't underflow when converting */;
  365.                     exp = 0;
  366.                     goto done;
  367.           }
  368.           (void) fpu_shr(fp, FP_NMANT - FP_NG - 1 - DBL_FRACBITS);
  369.           if (round(fe, fp) && fp->fp_mant[2] == DBL_EXP(2))
  370.                     exp++;
  371.           if (exp >= DBL_EXP_INFNAN) {
  372.                     fe->fe_fpsr |= FPSR_OPERR | FPSR_INEX2;
  373.                     if (toinf(fe, sign)) {
  374.                               res[1] = 0;
  375.                               return (sign | DBL_EXP(DBL_EXP_INFNAN) | 0);
  376.                     }
  377.                     res[1] = ~0;
  378.                     return (sign | DBL_EXP(DBL_EXP_INFNAN) | DBL_MASK);
  379.           }
  380. done:
  381.           res[1] = fp->fp_mant[3];
  382.           return (sign | DBL_EXP(exp) | (fp->fp_mant[2] & DBL_MASK));
  383. }
  384.  
  385. /*
  386.  * fpn -> 68k extended (32 bit high-order result returned; two 32-bit low
  387.  * order result left in res[1] & res[2]).  Assumes == 64 bits in extended
  388.  * precision fraction.
  389.  *
  390.  * This code mimics fpu_ftos; see it for comments.
  391.  */
  392. static u_int
  393. fpu_ftox(fe, fp, res)
  394.           struct fpemu *fe;
  395.           register struct fpn *fp;
  396.           u_int *res;
  397. {
  398.           register u_int sign = fp->fp_sign << 31;
  399.           register int exp;
  400.  
  401. #define   EXT_EXP(e)          ((e) << 16)
  402. #define   EXT_MASK  (EXT_EXP(1) - 1)
  403.  
  404.           if (ISNAN(fp)) {
  405.                     (void) fpu_shr(fp, FP_NMANT - 1 - EXT_FRACBITS);
  406.                     exp = EXT_EXP_INFNAN;
  407.                     goto done;
  408.           }
  409.           if (ISINF(fp)) {
  410.                     sign |= EXT_EXP(EXT_EXP_INFNAN);
  411.                     res[1] = res[2] = 0;
  412.                     return (sign);
  413.           }
  414.           if (ISZERO(fp)) {
  415.                     res[1] = res[2] = 0;
  416.                     return (sign);
  417.           }
  418.  
  419.           if ((exp = fp->fp_exp + EXT_EXP_BIAS) <= 0) {
  420.                     /* I'm not sure about this <=... exp==0 doesn't mean
  421.                        it's a denormal in extended format */
  422.                     (void) fpu_shr(fp, FP_NMANT - FP_NG - EXT_FRACBITS - exp);
  423.                     if (round(fe, fp) && fp->fp_mant[2] == EXT_EXP(1)) {
  424.                               res[1] = res[2] = 0;
  425.                               return (sign | EXT_EXP(1) | 0);
  426.                     }
  427.                     if (fe->fe_fpsr & FPSR_INEX2)
  428.                               /* mc68881/2 don't underflow */;
  429.                     exp = 0;
  430.                     goto done;
  431.           }
  432.           (void) fpu_shr(fp, FP_NMANT - FP_NG - EXT_FRACBITS);
  433.           if (round(fe, fp) && fp->fp_mant[2] == EXT_EXP(2))
  434.                     exp++;
  435.           if (exp >= EXT_EXP_INFNAN) {
  436.                     fe->fe_fpsr |= FPSR_OPERR | FPSR_INEX2;
  437.                     if (toinf(fe, sign)) {
  438.                               res[1] = res[2] = 0;
  439.                               return (sign | EXT_EXP(EXT_EXP_INFNAN) | 0);
  440.                     }
  441.                     res[1] = res[2] = ~0;
  442.                     return (sign | EXT_EXP(EXT_EXP_INFNAN) | EXT_MASK);
  443.           }
  444. done:
  445.           res[1] = fp->fp_mant[2];
  446.           res[2] = fp->fp_mant[3];
  447.           return (sign | EXT_EXP(exp));
  448. }
  449.  
  450. /*
  451.  * Implode an fpn, writing the result into the given space.
  452.  */
  453. void
  454. fpu_implode(fe, fp, type, space)
  455.           struct fpemu *fe;
  456.           register struct fpn *fp;
  457.           int type;
  458.           register u_int *space;
  459. {
  460.           fe->fe_fpsr &= ~FPSR_EXCP;
  461.  
  462.           switch (type) {
  463.           case FTYPE_LNG:
  464.                     space[0] = fpu_ftoi(fe, fp);
  465.                     break;
  466.  
  467.           case FTYPE_SNG:
  468.                     space[0] = fpu_ftos(fe, fp);
  469.                     break;
  470.  
  471.           case FTYPE_DBL:
  472.                     space[0] = fpu_ftod(fe, fp, space);
  473.                     break;
  474.  
  475.           case FTYPE_EXT:
  476.                     /* funky rounding precision options ?? */
  477.                     space[0] = fpu_ftox(fe, fp, space);
  478.                     break;
  479.  
  480.           default:
  481. #ifdef DIAGNOSTIC
  482.                     panic("fpu_implode"); 
  483. #endif
  484.  
  485.           }
  486. }
  487.