arch/mips/math-emu/ieee754dp.c

   1 /* IEEE754 floating point arithmetic
   2  * double precision: common utilities
   3  */
   4 /*
   5  * MIPS floating point support
   6  * Copyright (C) 1994-2000 Algorithmics Ltd.
   7  *
   8  * ########################################################################
   9  *
  10  *  This program is free software; you can distribute it and/or modify it
  11  *  under the terms of the GNU General Public License (Version 2) as
  12  *  published by the Free Software Foundation.
  13  *
  14  *  This program is distributed in the hope it will be useful, but WITHOUT
  15  *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  16  *  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  17  *  for more details.
  18  *
  19  *  You should have received a copy of the GNU General Public License along
  20  *  with this program; if not, write to the Free Software Foundation, Inc.,
  21  *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
  22  *
  23  * ########################################################################
  24  */
  25
  26 #include <stdarg.h>
  27 #include <linux/compiler.h>
  28
  29 #include "ieee754dp.h"
  30
  31 int ieee754dp_class(union ieee754dp x)
  32 {
  33         COMPXDP;
  34         EXPLODEXDP;
  35         return xc;
  36 }
  37
  38 int ieee754dp_isnan(union ieee754dp x)
  39 {
  40         return ieee754dp_class(x) >= IEEE754_CLASS_SNAN;
  41 }
  42
  43 static inline int ieee754dp_issnan(union ieee754dp x)
  44 {
  45         assert(ieee754dp_isnan(x));
  46         return ((DPMANT(x) & DP_MBIT(DP_FBITS-1)) == DP_MBIT(DP_FBITS-1));
  47 }
  48
  49
  50 union ieee754dp __cold ieee754dp_xcpt(union ieee754dp r, const char *op, ...)
  51 {
  52         struct ieee754xctx ax;
  53         if (!ieee754_tstx())
  54                 return r;
  55
  56         ax.op = op;
  57         ax.rt = IEEE754_RT_DP;
  58         ax.rv.dp = r;
  59         va_start(ax.ap, op);
  60         ieee754_xcpt(&ax);
  61         va_end(ax.ap);
  62         return ax.rv.dp;
  63 }
  64
  65 union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r, const char *op, ...)
  66 {
  67         struct ieee754xctx ax;
  68
  69         assert(ieee754dp_isnan(r));
  70
  71         if (!ieee754dp_issnan(r))       /* QNAN does not cause invalid op !! */
  72                 return r;
  73
  74         if (!ieee754_setandtestcx(IEEE754_INVALID_OPERATION)) {
  75                 /* not enabled convert to a quiet NaN */
  76                 DPMANT(r) &= (~DP_MBIT(DP_FBITS-1));
  77                 if (ieee754dp_isnan(r))
  78                         return r;
  79                 else
  80                         return ieee754dp_indef();
  81         }
  82
  83         ax.op = op;
  84         ax.rt = 0;
  85         ax.rv.dp = r;
  86         va_start(ax.ap, op);
  87         ieee754_xcpt(&ax);
  88         va_end(ax.ap);
  89         return ax.rv.dp;
  90 }
  91
  92 static u64 get_rounding(int sn, u64 xm)
  93 {
  94         /* inexact must round of 3 bits
  95          */
  96         if (xm & (DP_MBIT(3) - 1)) {
  97                 switch (ieee754_csr.rm) {
  98                 case IEEE754_RZ:
  99                         break;
 100                 case IEEE754_RN:
 101                         xm += 0x3 + ((xm >> 3) & 1);
 102                         /* xm += (xm&0x8)?0x4:0x3 */
 103                         break;
 104                 case IEEE754_RU:        /* toward +Infinity */
 105                         if (!sn)        /* ?? */
 106                                 xm += 0x8;
 107                         break;
 108                 case IEEE754_RD:        /* toward -Infinity */
 109                         if (sn) /* ?? */
 110                                 xm += 0x8;
 111                         break;
 112                 }
 113         }
 114         return xm;
 115 }
 116
 117
 118 /* generate a normal/denormal number with over,under handling
 119  * sn is sign
 120  * xe is an unbiased exponent
 121  * xm is 3bit extended precision value.
 122  */
 123 union ieee754dp ieee754dp_format(int sn, int xe, u64 xm)
 124 {
 125         assert(xm);             /* we don't gen exact zeros (probably should) */
 126
 127         assert((xm >> (DP_FBITS + 1 + 3)) == 0);        /* no execess */
 128         assert(xm & (DP_HIDDEN_BIT << 3));
 129
 130         if (xe < DP_EMIN) {
 131                 /* strip lower bits */
 132                 int es = DP_EMIN - xe;
 133
 134                 if (ieee754_csr.nod) {
 135                         ieee754_setcx(IEEE754_UNDERFLOW);
 136                         ieee754_setcx(IEEE754_INEXACT);
 137
 138                         switch(ieee754_csr.rm) {
 139                         case IEEE754_RN:
 140                         case IEEE754_RZ:
 141                                 return ieee754dp_zero(sn);
 142                         case IEEE754_RU:    /* toward +Infinity */
 143                                 if (sn == 0)
 144                                         return ieee754dp_min(0);
 145                                 else
 146                                         return ieee754dp_zero(1);
 147                         case IEEE754_RD:    /* toward -Infinity */
 148                                 if (sn == 0)
 149                                         return ieee754dp_zero(0);
 150                                 else
 151                                         return ieee754dp_min(1);
 152                         }
 153                 }
 154
 155                 if (xe == DP_EMIN - 1
 156                                 && get_rounding(sn, xm) >> (DP_FBITS + 1 + 3))
 157                 {
 158                         /* Not tiny after rounding */
 159                         ieee754_setcx(IEEE754_INEXACT);
 160                         xm = get_rounding(sn, xm);
 161                         xm >>= 1;
 162                         /* Clear grs bits */
 163                         xm &= ~(DP_MBIT(3) - 1);
 164                         xe++;
 165                 }
 166                 else {
 167                         /* sticky right shift es bits
 168                          */
 169                         xm = XDPSRS(xm, es);
 170                         xe += es;
 171                         assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
 172                         assert(xe == DP_EMIN);
 173                 }
 174         }
 175         if (xm & (DP_MBIT(3) - 1)) {
 176                 ieee754_setcx(IEEE754_INEXACT);
 177                 if ((xm & (DP_HIDDEN_BIT << 3)) == 0) {
 178                         ieee754_setcx(IEEE754_UNDERFLOW);
 179                 }
 180
 181                 /* inexact must round of 3 bits
 182                  */
 183                 xm = get_rounding(sn, xm);
 184                 /* adjust exponent for rounding add overflowing
 185                  */
 186                 if (xm >> (DP_FBITS + 3 + 1)) {
 187                         /* add causes mantissa overflow */
 188                         xm >>= 1;
 189                         xe++;
 190                 }
 191         }
 192         /* strip grs bits */
 193         xm >>= 3;
 194
 195         assert((xm >> (DP_FBITS + 1)) == 0);    /* no execess */
 196         assert(xe >= DP_EMIN);
 197
 198         if (xe > DP_EMAX) {
 199                 ieee754_setcx(IEEE754_OVERFLOW);
 200                 ieee754_setcx(IEEE754_INEXACT);
 201                 /* -O can be table indexed by (rm,sn) */
 202                 switch (ieee754_csr.rm) {
 203                 case IEEE754_RN:
 204                         return ieee754dp_inf(sn);
 205                 case IEEE754_RZ:
 206                         return ieee754dp_max(sn);
 207                 case IEEE754_RU:        /* toward +Infinity */
 208                         if (sn == 0)
 209                                 return ieee754dp_inf(0);
 210                         else
 211                                 return ieee754dp_max(1);
 212                 case IEEE754_RD:        /* toward -Infinity */
 213                         if (sn == 0)
 214                                 return ieee754dp_max(0);
 215                         else
 216                                 return ieee754dp_inf(1);
 217                 }
 218         }
 219         /* gen norm/denorm/zero */
 220
 221         if ((xm & DP_HIDDEN_BIT) == 0) {
 222                 /* we underflow (tiny/zero) */
 223                 assert(xe == DP_EMIN);
 224                 if (ieee754_csr.mx & IEEE754_UNDERFLOW)
 225                         ieee754_setcx(IEEE754_UNDERFLOW);
 226                 return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
 227         } else {
 228                 assert((xm >> (DP_FBITS + 1)) == 0);    /* no execess */
 229                 assert(xm & DP_HIDDEN_BIT);
 230
 231                 return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
 232         }
 233 }