| 1 |
|
|---|
| 2 | /*============================================================================
|
|---|
| 3 |
|
|---|
| 4 | This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic
|
|---|
| 5 | Package, Release 3e, by John R. Hauser.
|
|---|
| 6 |
|
|---|
| 7 | Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017 The Regents of the
|
|---|
| 8 | University of California. All rights reserved.
|
|---|
| 9 |
|
|---|
| 10 | Redistribution and use in source and binary forms, with or without
|
|---|
| 11 | modification, are permitted provided that the following conditions are met:
|
|---|
| 12 |
|
|---|
| 13 | 1. Redistributions of source code must retain the above copyright notice,
|
|---|
| 14 | this list of conditions, and the following disclaimer.
|
|---|
| 15 |
|
|---|
| 16 | 2. Redistributions in binary form must reproduce the above copyright notice,
|
|---|
| 17 | this list of conditions, and the following disclaimer in the documentation
|
|---|
| 18 | and/or other materials provided with the distribution.
|
|---|
| 19 |
|
|---|
| 20 | 3. Neither the name of the University nor the names of its contributors may
|
|---|
| 21 | be used to endorse or promote products derived from this software without
|
|---|
| 22 | specific prior written permission.
|
|---|
| 23 |
|
|---|
| 24 | THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
|
|---|
| 25 | EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
|
|---|
| 26 | WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
|
|---|
| 27 | DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
|
|---|
| 28 | DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
|
|---|
| 29 | (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
|
|---|
| 30 | LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
|
|---|
| 31 | ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
|---|
| 32 | (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
|
|---|
| 33 | SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|---|
| 34 |
|
|---|
| 35 | =============================================================================*/
|
|---|
| 36 |
|
|---|
| 37 | #include <stdbool.h>
|
|---|
| 38 | #include <stdint.h>
|
|---|
| 39 | #include "platform.h"
|
|---|
| 40 | #include "internals.h"
|
|---|
| 41 | #include "specialize.h"
|
|---|
| 42 | #include "softfloat.h"
|
|---|
| 43 | #include <iprt/cdefs.h> /* RT_BIT_64 */
|
|---|
| 44 | #include <iprt/x86.h> /* X86_FSW_C? */
|
|---|
| 45 | #include <iprt/assert.h>
|
|---|
| 46 |
|
|---|
| 47 | /** VBox: Copy of extF80_rem modified to fit FPREM and FPREM1. */
|
|---|
| 48 | extFloat80_t extF80_partialRem( extFloat80_t a, extFloat80_t b, uint8_t roundingMode,
|
|---|
| 49 | uint16_t *pfCxFlags, softfloat_state_t *pState )
|
|---|
| 50 | {
|
|---|
| 51 | union { struct extFloat80M s; extFloat80_t f; } uA;
|
|---|
| 52 | uint_fast16_t uiA64;
|
|---|
| 53 | uint_fast64_t uiA0;
|
|---|
| 54 | bool signA;
|
|---|
| 55 | int_fast32_t expA;
|
|---|
| 56 | uint_fast64_t sigA;
|
|---|
| 57 | union { struct extFloat80M s; extFloat80_t f; } uB;
|
|---|
| 58 | uint_fast16_t uiB64;
|
|---|
| 59 | uint_fast64_t uiB0;
|
|---|
| 60 | int_fast32_t expB;
|
|---|
| 61 | uint_fast64_t sigB;
|
|---|
| 62 | struct exp32_sig64 normExpSig;
|
|---|
| 63 | int_fast32_t expDiff;
|
|---|
| 64 | struct uint128 rem, shiftedSigB;
|
|---|
| 65 | uint_fast32_t q, recip32;
|
|---|
| 66 | uint_fast64_t q64;
|
|---|
| 67 | struct uint128 term, altRem, meanRem;
|
|---|
| 68 | bool signRem;
|
|---|
| 69 | struct uint128 uiZ;
|
|---|
| 70 | uint_fast16_t uiZ64;
|
|---|
| 71 | uint_fast64_t uiZ0;
|
|---|
| 72 | union { struct extFloat80M s; extFloat80_t f; } uZ;
|
|---|
| 73 |
|
|---|
| 74 | *pfCxFlags = 0; /* C2=0 - complete */ /* VBox */
|
|---|
| 75 |
|
|---|
| 76 | /*------------------------------------------------------------------------
|
|---|
| 77 | *------------------------------------------------------------------------*/
|
|---|
| 78 | uA.f = a;
|
|---|
| 79 | uiA64 = uA.s.signExp;
|
|---|
| 80 | uiA0 = uA.s.signif;
|
|---|
| 81 | signA = signExtF80UI64( uiA64 );
|
|---|
| 82 | expA = expExtF80UI64( uiA64 );
|
|---|
| 83 | sigA = uiA0;
|
|---|
| 84 | uB.f = b;
|
|---|
| 85 | uiB64 = uB.s.signExp;
|
|---|
| 86 | uiB0 = uB.s.signif;
|
|---|
| 87 | expB = expExtF80UI64( uiB64 );
|
|---|
| 88 | sigB = uiB0;
|
|---|
| 89 | /*------------------------------------------------------------------------
|
|---|
| 90 | *------------------------------------------------------------------------*/
|
|---|
| 91 | if ( expA == 0x7FFF ) {
|
|---|
| 92 | if (
|
|---|
| 93 | (sigA & UINT64_C( 0x7FFFFFFFFFFFFFFF )) /* VBox: NaN or Indefinite */
|
|---|
| 94 | || ((expB == 0x7FFF) && (sigB & UINT64_C( 0x7FFFFFFFFFFFFFFF ))) /* VBox: NaN or Indefinite */
|
|---|
| 95 | ) {
|
|---|
| 96 | goto propagateNaN;
|
|---|
| 97 | }
|
|---|
| 98 | goto invalid; /* VBox: Infinity */
|
|---|
| 99 | }
|
|---|
| 100 | if ( expB == 0x7FFF ) {
|
|---|
| 101 | if ( sigB & UINT64_C( 0x7FFFFFFFFFFFFFFF ) ) goto propagateNaN; /* VBox: NaN or Indefinite */
|
|---|
| 102 | /*--------------------------------------------------------------------
|
|---|
| 103 | | Argument b is an infinity. Doubling `expB' is an easy way to ensure
|
|---|
| 104 | | that `expDiff' later is less than -1, which will result in returning
|
|---|
| 105 | | a canonicalized version of argument a.
|
|---|
| 106 | *--------------------------------------------------------------------*/
|
|---|
| 107 | expB += expB;
|
|---|
| 108 | }
|
|---|
| 109 | /*------------------------------------------------------------------------
|
|---|
| 110 | *------------------------------------------------------------------------*/
|
|---|
| 111 | if ( ! expB ) expB = 1;
|
|---|
| 112 | if ( ! (sigB & UINT64_C( 0x8000000000000000 )) ) {
|
|---|
| 113 | if ( ! sigB ) goto invalid; /* VBox: Zero -> /0 -> invalid. */
|
|---|
| 114 | normExpSig = softfloat_normSubnormalExtF80Sig( sigB );
|
|---|
| 115 | expB += normExpSig.exp;
|
|---|
| 116 | sigB = normExpSig.sig;
|
|---|
| 117 | }
|
|---|
| 118 | if ( ! expA ) expA = 1;
|
|---|
| 119 | if ( ! (sigA & UINT64_C( 0x8000000000000000 )) ) {
|
|---|
| 120 | if ( ! sigA ) {
|
|---|
| 121 | #if 0 /* A is zero. VBox: Don't mix denormals and zero returns! */ /* VBox */
|
|---|
| 122 | expA = 0;
|
|---|
| 123 | goto copyA;
|
|---|
| 124 | #else /* VBox */
|
|---|
| 125 | uiZ64 = packToExtF80UI64( signA, 0); /* VBox */
|
|---|
| 126 | uiZ0 = 0; /* VBox */
|
|---|
| 127 | goto uiZ; /* VBox */
|
|---|
| 128 | #endif /* VBox */
|
|---|
| 129 | }
|
|---|
| 130 | normExpSig = softfloat_normSubnormalExtF80Sig( sigA );
|
|---|
| 131 | expA += normExpSig.exp;
|
|---|
| 132 | sigA = normExpSig.sig;
|
|---|
| 133 | }
|
|---|
| 134 | /*------------------------------------------------------------------------
|
|---|
| 135 | *------------------------------------------------------------------------*/
|
|---|
| 136 | expDiff = expA - expB;
|
|---|
| 137 |
|
|---|
| 138 | /*------------------------------------------------------------------------ VBox
|
|---|
| 139 | | Do at most 63 rounds. If exponent difference is 64 or higher, return VBox
|
|---|
| 140 | | a partial remainder. VBox
|
|---|
| 141 | *------------------------------------------------------------------------*/ /* VBox */
|
|---|
| 142 | bool const fPartial = expDiff >= 64; /* VBox */
|
|---|
| 143 | if ( fPartial ) { /* VBox */ /* VBox */
|
|---|
| 144 | unsigned N = 32 + ((unsigned)expDiff % 32); /* (Amount of work documented by AMD.) */ /* VBox */
|
|---|
| 145 | expB = expA - N; /* VBox */
|
|---|
| 146 | expDiff = N; /* VBox */
|
|---|
| 147 | roundingMode = softfloat_round_minMag; /* VBox */
|
|---|
| 148 | } /* VBox */
|
|---|
| 149 |
|
|---|
| 150 | /*------------------------------------------------------------------------
|
|---|
| 151 | | The final rounds.
|
|---|
| 152 | *------------------------------------------------------------------------*/
|
|---|
| 153 | if ( expDiff < -1 ) goto copyA;
|
|---|
| 154 | rem = softfloat_shortShiftLeft128( 0, sigA, 32 );
|
|---|
| 155 | shiftedSigB = softfloat_shortShiftLeft128( 0, sigB, 32 );
|
|---|
| 156 | uint64_t quotient = 0; /* VBox */
|
|---|
| 157 | if ( expDiff < 1 ) {
|
|---|
| 158 | if ( expDiff ) {
|
|---|
| 159 | --expB;
|
|---|
| 160 | shiftedSigB = softfloat_shortShiftLeft128( 0, sigB, 33 );
|
|---|
| 161 | q = 0;
|
|---|
| 162 | } else {
|
|---|
| 163 | q = (sigB <= sigA);
|
|---|
| 164 | quotient = q; /* VBox */
|
|---|
| 165 | if ( q ) {
|
|---|
| 166 | rem =
|
|---|
| 167 | softfloat_sub128(
|
|---|
| 168 | rem.v64, rem.v0, shiftedSigB.v64, shiftedSigB.v0 );
|
|---|
| 169 | }
|
|---|
| 170 | }
|
|---|
| 171 | } else {
|
|---|
| 172 | recip32 = softfloat_approxRecip32_1( sigB>>32 );
|
|---|
| 173 | expDiff -= 30;
|
|---|
| 174 | for (;;) {
|
|---|
| 175 | q64 = (uint_fast64_t) (uint32_t) (rem.v64>>2) * recip32;
|
|---|
| 176 | if ( expDiff < 0 ) break;
|
|---|
| 177 | q = (q64 + 0x80000000)>>32;
|
|---|
| 178 | quotient = (quotient << 29) + q; /* VBox */
|
|---|
| 179 | rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, 29 );
|
|---|
| 180 | term = softfloat_mul64ByShifted32To128( sigB, q );
|
|---|
| 181 | rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
|
|---|
| 182 | if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
|
|---|
| 183 | rem =
|
|---|
| 184 | softfloat_add128(
|
|---|
| 185 | rem.v64, rem.v0, shiftedSigB.v64, shiftedSigB.v0 );
|
|---|
| 186 | quotient--; /* VBox */
|
|---|
| 187 | }
|
|---|
| 188 | expDiff -= 29;
|
|---|
| 189 | }
|
|---|
| 190 | /*--------------------------------------------------------------------
|
|---|
| 191 | | (`expDiff' cannot be less than -29 here.)
|
|---|
| 192 | *--------------------------------------------------------------------*/
|
|---|
| 193 | q = (uint32_t) (q64>>32)>>(~expDiff & 31);
|
|---|
| 194 | quotient = (quotient << (expDiff + 30)) + q; /* VBox */
|
|---|
| 195 | rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, expDiff + 30 );
|
|---|
| 196 | term = softfloat_mul64ByShifted32To128( sigB, q );
|
|---|
| 197 | rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
|
|---|
| 198 | if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
|
|---|
| 199 | altRem =
|
|---|
| 200 | softfloat_add128(
|
|---|
| 201 | rem.v64, rem.v0, shiftedSigB.v64, shiftedSigB.v0 );
|
|---|
| 202 | goto selectRem;
|
|---|
| 203 | }
|
|---|
| 204 | }
|
|---|
| 205 | /*------------------------------------------------------------------------
|
|---|
| 206 | *------------------------------------------------------------------------*/
|
|---|
| 207 | do {
|
|---|
| 208 | altRem = rem;
|
|---|
| 209 | ++q;
|
|---|
| 210 | quotient++; /* VBox */
|
|---|
| 211 | rem =
|
|---|
| 212 | softfloat_sub128(
|
|---|
| 213 | rem.v64, rem.v0, shiftedSigB.v64, shiftedSigB.v0 );
|
|---|
| 214 | } while ( ! (rem.v64 & UINT64_C( 0x8000000000000000 )) );
|
|---|
| 215 | selectRem:
|
|---|
| 216 | if (roundingMode == softfloat_round_near_even) { /* VBox */
|
|---|
| 217 | meanRem = softfloat_add128( rem.v64, rem.v0, altRem.v64, altRem.v0 );
|
|---|
| 218 | if (
|
|---|
| 219 | (meanRem.v64 & UINT64_C( 0x8000000000000000 ))
|
|---|
| 220 | || (! (meanRem.v64 | meanRem.v0) && (q & 1))
|
|---|
| 221 | ) {
|
|---|
| 222 | rem = altRem;
|
|---|
| 223 | quotient--; /* VBox */
|
|---|
| 224 | }
|
|---|
| 225 | } /* VBox */
|
|---|
| 226 | signRem = signA;
|
|---|
| 227 | if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
|
|---|
| 228 | if (roundingMode != softfloat_round_near_even) { /* VBox */
|
|---|
| 229 | rem = altRem; /* VBox */
|
|---|
| 230 | quotient--; /* VBox */
|
|---|
| 231 | } else { /* VBox */
|
|---|
| 232 | signRem = ! signRem;
|
|---|
| 233 | rem = softfloat_sub128( 0, 0, rem.v64, rem.v0 );
|
|---|
| 234 | Assert(!fPartial); /* VBox */
|
|---|
| 235 | } /* VBox */
|
|---|
| 236 | } else Assert(roundingMode == softfloat_round_near_even); /* VBox */
|
|---|
| 237 |
|
|---|
| 238 | /* VBox: Set pfCxFlags */ /* VBox */
|
|---|
| 239 | if ( fPartial ) { /* VBox */
|
|---|
| 240 | *pfCxFlags = X86_FSW_C2; /* C2 = 1 - incomplete */ /* VBox */
|
|---|
| 241 | } else { /* VBox */
|
|---|
| 242 | *pfCxFlags = X86_FSW_CX_FROM_QUOTIENT( quotient ); /* C2 = 0 - complete */ /* VBox */
|
|---|
| 243 | } /* VBox */
|
|---|
| 244 |
|
|---|
| 245 | return
|
|---|
| 246 | softfloat_normRoundPackToExtF80(
|
|---|
| 247 | signRem, rem.v64 | rem.v0 ? expB + 32 : 0, rem.v64, rem.v0, 80 SOFTFLOAT_STATE_ARG_COMMA );
|
|---|
| 248 | /*------------------------------------------------------------------------
|
|---|
| 249 | *------------------------------------------------------------------------*/
|
|---|
| 250 | propagateNaN:
|
|---|
| 251 | uiZ = softfloat_propagateNaNExtF80UI( uiA64, uiA0, uiB64, uiB0 SOFTFLOAT_STATE_ARG_COMMA );
|
|---|
| 252 | uiZ64 = uiZ.v64;
|
|---|
| 253 | uiZ0 = uiZ.v0;
|
|---|
| 254 | goto uiZ;
|
|---|
| 255 | /*------------------------------------------------------------------------
|
|---|
| 256 | *------------------------------------------------------------------------*/
|
|---|
| 257 | invalid:
|
|---|
| 258 | softfloat_raiseFlags( softfloat_flag_invalid SOFTFLOAT_STATE_ARG_COMMA );
|
|---|
| 259 | uiZ64 = defaultNaNExtF80UI64;
|
|---|
| 260 | uiZ0 = defaultNaNExtF80UI0;
|
|---|
| 261 | goto uiZ;
|
|---|
| 262 | /*------------------------------------------------------------------------
|
|---|
| 263 | *------------------------------------------------------------------------*/
|
|---|
| 264 | copyA:
|
|---|
| 265 | if ( expA < 1 ) {
|
|---|
| 266 | #if 0 /* VBox */
|
|---|
| 267 | sigA >>= 1 - expA;
|
|---|
| 268 | expA = 0;
|
|---|
| 269 | #else /* VBox */
|
|---|
| 270 | Assert(sigA != 0); /* We don't get here for zero values, only denormals. */ /* VBox */
|
|---|
| 271 | /* Apply the bias adjust if underflows exceptions aren't masked, unless VBox
|
|---|
| 272 | the divisor is +/-Infinity. VBox
|
|---|
| 273 | Note! extB has been tweaked, so don't use if for Inf classification. */ /* VBox */
|
|---|
| 274 | if ( (pState->exceptionMask & softfloat_flag_underflow) /* VBox */
|
|---|
| 275 | || (expExtF80UI64( b.signExp ) == 0x7fff && !(sigB & (RT_BIT_64( 63 ) - 1))) ) { /* VBox */
|
|---|
| 276 | sigA >>= 1 - expA; /* VBox */
|
|---|
| 277 | expA = 0; /* VBox */
|
|---|
| 278 | } else { /* VBox */
|
|---|
| 279 | softfloat_raiseFlags( softfloat_flag_underflow SOFTFLOAT_STATE_ARG_COMMA ); /* VBox */
|
|---|
| 280 | expA = (expA + RTFLOAT80U_EXP_BIAS_ADJUST) & RTFLOAT80U_EXP_MAX; /* VBox */
|
|---|
| 281 | } /* VBox */
|
|---|
| 282 | #endif /* VBox */
|
|---|
| 283 | }
|
|---|
| 284 | uiZ64 = packToExtF80UI64( signA, expA );
|
|---|
| 285 | uiZ0 = sigA;
|
|---|
| 286 | uiZ:
|
|---|
| 287 | uZ.s.signExp = uiZ64;
|
|---|
| 288 | uZ.s.signif = uiZ0;
|
|---|
| 289 | return uZ.f;
|
|---|
| 290 |
|
|---|
| 291 | }
|
|---|
| 292 |
|
|---|
