LXR linux/arch/mips/math-emu/ieee754sp.c

   1/* IEEE754 floating point arithmetic
   2 * single precision
   3 */
   4/*
   5 * MIPS floating point support
   6 * Copyright (C) 1994-2000 Algorithmics Ltd.
   7 *
   8 *  This program is free software; you can distribute it and/or modify it
   9 *  under the terms of the GNU General Public License (Version 2) as
  10 *  published by the Free Software Foundation.
  11 *
  12 *  This program is distributed in the hope it will be useful, but WITHOUT
  13 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14 *  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15 *  for more details.
  16 *
  17 *  You should have received a copy of the GNU General Public License along
  18 *  with this program; if not, write to the Free Software Foundation, Inc.,
  19 *  51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA.
  20 */
  21
  22#include <linux/compiler.h>
  23
  24#include "ieee754sp.h"
  25
  26int ieee754sp_class(union ieee754sp x)
  27{
  28        COMPXSP;
  29        EXPLODEXSP;
  30        return xc;
  31}
  32
  33static inline int ieee754sp_isnan(union ieee754sp x)
  34{
  35        return ieee754_class_nan(ieee754sp_class(x));
  36}
  37
  38static inline int ieee754sp_issnan(union ieee754sp x)
  39{
  40        int qbit;
  41
  42        assert(ieee754sp_isnan(x));
  43        qbit = (SPMANT(x) & SP_MBIT(SP_FBITS - 1)) == SP_MBIT(SP_FBITS - 1);
  44        return ieee754_csr.nan2008 ^ qbit;
  45}
  46
  47
  48/*
  49 * Raise the Invalid Operation IEEE 754 exception
  50 * and convert the signaling NaN supplied to a quiet NaN.
  51 */
  52union ieee754sp __cold ieee754sp_nanxcpt(union ieee754sp r)
  53{
  54        assert(ieee754sp_issnan(r));
  55
  56        ieee754_setcx(IEEE754_INVALID_OPERATION);
  57        if (ieee754_csr.nan2008) {
  58                SPMANT(r) |= SP_MBIT(SP_FBITS - 1);
  59        } else {
  60                SPMANT(r) &= ~SP_MBIT(SP_FBITS - 1);
  61                if (!ieee754sp_isnan(r))
  62                        SPMANT(r) |= SP_MBIT(SP_FBITS - 2);
  63        }
  64
  65        return r;
  66}
  67
  68static unsigned ieee754sp_get_rounding(int sn, unsigned xm)
  69{
  70        /* inexact must round of 3 bits
  71         */
  72        if (xm & (SP_MBIT(3) - 1)) {
  73                switch (ieee754_csr.rm) {
  74                case FPU_CSR_RZ:
  75                        break;
  76                case FPU_CSR_RN:
  77                        xm += 0x3 + ((xm >> 3) & 1);
  78                        /* xm += (xm&0x8)?0x4:0x3 */
  79                        break;
  80                case FPU_CSR_RU:        /* toward +Infinity */
  81                        if (!sn)        /* ?? */
  82                                xm += 0x8;
  83                        break;
  84                case FPU_CSR_RD:        /* toward -Infinity */
  85                        if (sn) /* ?? */
  86                                xm += 0x8;
  87                        break;
  88                }
  89        }
  90        return xm;
  91}
  92
  93
  94/* generate a normal/denormal number with over,under handling
  95 * sn is sign
  96 * xe is an unbiased exponent
  97 * xm is 3bit extended precision value.
  98 */
  99union ieee754sp ieee754sp_format(int sn, int xe, unsigned xm)
 100{
 101        assert(xm);             /* we don't gen exact zeros (probably should) */
 102
 103        assert((xm >> (SP_FBITS + 1 + 3)) == 0);        /* no excess */
 104        assert(xm & (SP_HIDDEN_BIT << 3));
 105
 106        if (xe < SP_EMIN) {
 107                /* strip lower bits */
 108                int es = SP_EMIN - xe;
 109
 110                if (ieee754_csr.nod) {
 111                        ieee754_setcx(IEEE754_UNDERFLOW);
 112                        ieee754_setcx(IEEE754_INEXACT);
 113
 114                        switch(ieee754_csr.rm) {
 115                        case FPU_CSR_RN:
 116                        case FPU_CSR_RZ:
 117                                return ieee754sp_zero(sn);
 118                        case FPU_CSR_RU:      /* toward +Infinity */
 119                                if (sn == 0)
 120                                        return ieee754sp_min(0);
 121                                else
 122                                        return ieee754sp_zero(1);
 123                        case FPU_CSR_RD:      /* toward -Infinity */
 124                                if (sn == 0)
 125                                        return ieee754sp_zero(0);
 126                                else
 127                                        return ieee754sp_min(1);
 128                        }
 129                }
 130
 131                if (xe == SP_EMIN - 1 &&
 132                    ieee754sp_get_rounding(sn, xm) >> (SP_FBITS + 1 + 3))
 133                {
 134                        /* Not tiny after rounding */
 135                        ieee754_setcx(IEEE754_INEXACT);
 136                        xm = ieee754sp_get_rounding(sn, xm);
 137                        xm >>= 1;
 138                        /* Clear grs bits */
 139                        xm &= ~(SP_MBIT(3) - 1);
 140                        xe++;
 141                } else {
 142                        /* sticky right shift es bits
 143                         */
 144                        xm = XSPSRS(xm, es);
 145                        xe += es;
 146                        assert((xm & (SP_HIDDEN_BIT << 3)) == 0);
 147                        assert(xe == SP_EMIN);
 148                }
 149        }
 150        if (xm & (SP_MBIT(3) - 1)) {
 151                ieee754_setcx(IEEE754_INEXACT);
 152                if ((xm & (SP_HIDDEN_BIT << 3)) == 0) {
 153                        ieee754_setcx(IEEE754_UNDERFLOW);
 154                }
 155
 156                /* inexact must round of 3 bits
 157                 */
 158                xm = ieee754sp_get_rounding(sn, xm);
 159                /* adjust exponent for rounding add overflowing
 160                 */
 161                if (xm >> (SP_FBITS + 1 + 3)) {
 162                        /* add causes mantissa overflow */
 163                        xm >>= 1;
 164                        xe++;
 165                }
 166        }
 167        /* strip grs bits */
 168        xm >>= 3;
 169
 170        assert((xm >> (SP_FBITS + 1)) == 0);    /* no excess */
 171        assert(xe >= SP_EMIN);
 172
 173        if (xe > SP_EMAX) {
 174                ieee754_setcx(IEEE754_OVERFLOW);
 175                ieee754_setcx(IEEE754_INEXACT);
 176                /* -O can be table indexed by (rm,sn) */
 177                switch (ieee754_csr.rm) {
 178                case FPU_CSR_RN:
 179                        return ieee754sp_inf(sn);
 180                case FPU_CSR_RZ:
 181                        return ieee754sp_max(sn);
 182                case FPU_CSR_RU:        /* toward +Infinity */
 183                        if (sn == 0)
 184                                return ieee754sp_inf(0);
 185                        else
 186                                return ieee754sp_max(1);
 187                case FPU_CSR_RD:        /* toward -Infinity */
 188                        if (sn == 0)
 189                                return ieee754sp_max(0);
 190                        else
 191                                return ieee754sp_inf(1);
 192                }
 193        }
 194        /* gen norm/denorm/zero */
 195
 196        if ((xm & SP_HIDDEN_BIT) == 0) {
 197                /* we underflow (tiny/zero) */
 198                assert(xe == SP_EMIN);
 199                if (ieee754_csr.mx & IEEE754_UNDERFLOW)
 200                        ieee754_setcx(IEEE754_UNDERFLOW);
 201                return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm);
 202        } else {
 203                assert((xm >> (SP_FBITS + 1)) == 0);    /* no excess */
 204                assert(xm & SP_HIDDEN_BIT);
 205
 206                return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
 207        }
 208}
 209