target/m68k/softfloat.c

   1 /*
   2  * Ported from a work by Andreas Grabher for Previous, NeXT Computer Emulator,
   3  * derived from NetBSD M68040 FPSP functions,
   4  * derived from release 2a of the SoftFloat IEC/IEEE Floating-point Arithmetic
   5  * Package. Those parts of the code (and some later contributions) are
   6  * provided under that license, as detailed below.
   7  * It has subsequently been modified by contributors to the QEMU Project,
   8  * so some portions are provided under:
   9  *  the SoftFloat-2a license
  10  *  the BSD license
  11  *  GPL-v2-or-later
  12  *
  13  * Any future contributions to this file will be taken to be licensed under
  14  * the Softfloat-2a license unless specifically indicated otherwise.
  15  */
  16
  17 /* Portions of this work are licensed under the terms of the GNU GPL,
  18  * version 2 or later. See the COPYING file in the top-level directory.
  19  */
  20
  21 #include "qemu/osdep.h"
  22 #include "softfloat.h"
  23 #include "fpu/softfloat-macros.h"
  24
  25 static floatx80 propagateFloatx80NaNOneArg(floatx80 a, float_status *status)
  26 {
  27     if (floatx80_is_signaling_nan(a, status)) {
  28         float_raise(float_flag_invalid, status);
  29     }
  30
  31     if (status->default_nan_mode) {
  32         return floatx80_default_nan(status);
  33     }
  34
  35     return floatx80_maybe_silence_nan(a, status);
  36 }
  37
  38 /*----------------------------------------------------------------------------
  39  | Returns the modulo remainder of the extended double-precision floating-point
  40  | value `a' with respect to the corresponding value `b'.
  41  *----------------------------------------------------------------------------*/
  42
  43 floatx80 floatx80_mod(floatx80 a, floatx80 b, float_status *status)
  44 {
  45     flag aSign, zSign;
  46     int32_t aExp, bExp, expDiff;
  47     uint64_t aSig0, aSig1, bSig;
  48     uint64_t qTemp, term0, term1;
  49
  50     aSig0 = extractFloatx80Frac(a);
  51     aExp = extractFloatx80Exp(a);
  52     aSign = extractFloatx80Sign(a);
  53     bSig = extractFloatx80Frac(b);
  54     bExp = extractFloatx80Exp(b);
  55
  56     if (aExp == 0x7FFF) {
  57         if ((uint64_t) (aSig0 << 1)
  58             || ((bExp == 0x7FFF) && (uint64_t) (bSig << 1))) {
  59             return propagateFloatx80NaN(a, b, status);
  60         }
  61         goto invalid;
  62     }
  63     if (bExp == 0x7FFF) {
  64         if ((uint64_t) (bSig << 1)) {
  65             return propagateFloatx80NaN(a, b, status);
  66         }
  67         return a;
  68     }
  69     if (bExp == 0) {
  70         if (bSig == 0) {
  71         invalid:
  72             float_raise(float_flag_invalid, status);
  73             return floatx80_default_nan(status);
  74         }
  75         normalizeFloatx80Subnormal(bSig, &bExp, &bSig);
  76     }
  77     if (aExp == 0) {
  78         if ((uint64_t) (aSig0 << 1) == 0) {
  79             return a;
  80         }
  81         normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
  82     }
  83     bSig |= LIT64(0x8000000000000000);
  84     zSign = aSign;
  85     expDiff = aExp - bExp;
  86     aSig1 = 0;
  87     if (expDiff < 0) {
  88         return a;
  89     }
  90     qTemp = (bSig <= aSig0);
  91     if (qTemp) {
  92         aSig0 -= bSig;
  93     }
  94     expDiff -= 64;
  95     while (0 < expDiff) {
  96         qTemp = estimateDiv128To64(aSig0, aSig1, bSig);
  97         qTemp = (2 < qTemp) ? qTemp - 2 : 0;
  98         mul64To128(bSig, qTemp, &term0, &term1);
  99         sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1);
 100         shortShift128Left(aSig0, aSig1, 62, &aSig0, &aSig1);
 101     }
 102     expDiff += 64;
 103     if (0 < expDiff) {
 104         qTemp = estimateDiv128To64(aSig0, aSig1, bSig);
 105         qTemp = (2 < qTemp) ? qTemp - 2 : 0;
 106         qTemp >>= 64 - expDiff;
 107         mul64To128(bSig, qTemp << (64 - expDiff), &term0, &term1);
 108         sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1);
 109         shortShift128Left(0, bSig, 64 - expDiff, &term0, &term1);
 110         while (le128(term0, term1, aSig0, aSig1)) {
 111             ++qTemp;
 112             sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1);
 113         }
 114     }
 115     return
 116         normalizeRoundAndPackFloatx80(
 117             80, zSign, bExp + expDiff, aSig0, aSig1, status);
 118 }
 119
 120 /*----------------------------------------------------------------------------
 121  | Returns the mantissa of the extended double-precision floating-point
 122  | value `a'.
 123  *----------------------------------------------------------------------------*/
 124
 125 floatx80 floatx80_getman(floatx80 a, float_status *status)
 126 {
 127     flag aSign;
 128     int32_t aExp;
 129     uint64_t aSig;
 130
 131     aSig = extractFloatx80Frac(a);
 132     aExp = extractFloatx80Exp(a);
 133     aSign = extractFloatx80Sign(a);
 134
 135     if (aExp == 0x7FFF) {
 136         if ((uint64_t) (aSig << 1)) {
 137             return propagateFloatx80NaNOneArg(a , status);
 138         }
 139         float_raise(float_flag_invalid , status);
 140         return floatx80_default_nan(status);
 141     }
 142
 143     if (aExp == 0) {
 144         if (aSig == 0) {
 145             return packFloatx80(aSign, 0, 0);
 146         }
 147         normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
 148     }
 149
 150     return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign,
 151                                 0x3FFF, aSig, 0, status);
 152 }
 153
 154 /*----------------------------------------------------------------------------
 155  | Returns the exponent of the extended double-precision floating-point
 156  | value `a' as an extended double-precision value.
 157  *----------------------------------------------------------------------------*/
 158
 159 floatx80 floatx80_getexp(floatx80 a, float_status *status)
 160 {
 161     flag aSign;
 162     int32_t aExp;
 163     uint64_t aSig;
 164
 165     aSig = extractFloatx80Frac(a);
 166     aExp = extractFloatx80Exp(a);
 167     aSign = extractFloatx80Sign(a);
 168
 169     if (aExp == 0x7FFF) {
 170         if ((uint64_t) (aSig << 1)) {
 171             return propagateFloatx80NaNOneArg(a , status);
 172         }
 173         float_raise(float_flag_invalid , status);
 174         return floatx80_default_nan(status);
 175     }
 176
 177     if (aExp == 0) {
 178         if (aSig == 0) {
 179             return packFloatx80(aSign, 0, 0);
 180         }
 181         normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
 182     }
 183
 184     return int32_to_floatx80(aExp - 0x3FFF, status);
 185 }
 186
 187 /*----------------------------------------------------------------------------
 188  | Scales extended double-precision floating-point value in operand `a' by
 189  | value `b'. The function truncates the value in the second operand 'b' to
 190  | an integral value and adds that value to the exponent of the operand 'a'.
 191  | The operation performed according to the IEC/IEEE Standard for Binary
 192  | Floating-Point Arithmetic.
 193  *----------------------------------------------------------------------------*/
 194
 195 floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status)
 196 {
 197     flag aSign, bSign;
 198     int32_t aExp, bExp, shiftCount;
 199     uint64_t aSig, bSig;
 200
 201     aSig = extractFloatx80Frac(a);
 202     aExp = extractFloatx80Exp(a);
 203     aSign = extractFloatx80Sign(a);
 204     bSig = extractFloatx80Frac(b);
 205     bExp = extractFloatx80Exp(b);
 206     bSign = extractFloatx80Sign(b);
 207
 208     if (bExp == 0x7FFF) {
 209         if ((uint64_t) (bSig << 1) ||
 210             ((aExp == 0x7FFF) && (uint64_t) (aSig << 1))) {
 211             return propagateFloatx80NaN(a, b, status);
 212         }
 213         float_raise(float_flag_invalid , status);
 214         return floatx80_default_nan(status);
 215     }
 216     if (aExp == 0x7FFF) {
 217         if ((uint64_t) (aSig << 1)) {
 218             return propagateFloatx80NaN(a, b, status);
 219         }
 220         return packFloatx80(aSign, floatx80_infinity.high,
 221                             floatx80_infinity.low);
 222     }
 223     if (aExp == 0) {
 224         if (aSig == 0) {
 225             return packFloatx80(aSign, 0, 0);
 226         }
 227         if (bExp < 0x3FFF) {
 228             return a;
 229         }
 230         normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
 231     }
 232
 233     if (bExp < 0x3FFF) {
 234         return a;
 235     }
 236
 237     if (0x400F < bExp) {
 238         aExp = bSign ? -0x6001 : 0xE000;
 239         return roundAndPackFloatx80(status->floatx80_rounding_precision,
 240                                     aSign, aExp, aSig, 0, status);
 241     }
 242
 243     shiftCount = 0x403E - bExp;
 244     bSig >>= shiftCount;
 245     aExp = bSign ? (aExp - bSig) : (aExp + bSig);
 246
 247     return roundAndPackFloatx80(status->floatx80_rounding_precision,
 248                                 aSign, aExp, aSig, 0, status);
 249 }