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