1 2 /*============================================================================ 3 4 This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic 5 Package, Release 3a, by John R. Hauser. 6 7 Copyright 2011, 2012, 2013, 2014 The Regents of the University of California. 8 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 44 float128_t f128_rem( float128_t a, float128_t b ) 45 { 46 union ui128_f128 uA; 47 uint_fast64_t uiA64, uiA0; 48 bool signA; 49 int_fast32_t expA; 50 struct uint128 sigA; 51 union ui128_f128 uB; 52 uint_fast64_t uiB64, uiB0; 53 bool signB; 54 int_fast32_t expB; 55 struct uint128 sigB; 56 struct exp32_sig128 normExpSig; 57 struct uint128 rem; 58 int_fast32_t expDiff; 59 uint_fast32_t q, recip32; 60 uint_fast64_t q64; 61 struct uint128 term, altRem, meanRem; 62 bool signRem; 63 struct uint128 uiZ; 64 union ui128_f128 uZ; 65 66 /*------------------------------------------------------------------------ 67 *------------------------------------------------------------------------*/ 68 uA.f = a; 69 uiA64 = uA.ui.v64; 70 uiA0 = uA.ui.v0; 71 signA = signF128UI64( uiA64 ); 72 expA = expF128UI64( uiA64 ); 73 sigA.v64 = fracF128UI64( uiA64 ); 74 sigA.v0 = uiA0; 75 uB.f = b; 76 uiB64 = uB.ui.v64; 77 uiB0 = uB.ui.v0; 78 signB = signF128UI64( uiB64 ); 79 expB = expF128UI64( uiB64 ); 80 sigB.v64 = fracF128UI64( uiB64 ); 81 sigB.v0 = uiB0; 82 /*------------------------------------------------------------------------ 83 *------------------------------------------------------------------------*/ 84 if ( expA == 0x7FFF ) { 85 if ( 86 (sigA.v64 | sigA.v0) || ((expB == 0x7FFF) && (sigB.v64 | sigB.v0)) 87 ) { 88 goto propagateNaN; 89 } 90 goto invalid; 91 } 92 if ( expB == 0x7FFF ) { 93 if ( sigB.v64 | sigB.v0 ) goto propagateNaN; 94 return a; 95 } 96 /*------------------------------------------------------------------------ 97 *------------------------------------------------------------------------*/ 98 if ( ! expB ) { 99 if ( ! (sigB.v64 | sigB.v0) ) goto invalid; 100 normExpSig = softfloat_normSubnormalF128Sig( sigB.v64, sigB.v0 ); 101 expB = normExpSig.exp; 102 sigB = normExpSig.sig; 103 } 104 if ( ! expA ) { 105 if ( ! (sigA.v64 | sigA.v0) ) return a; 106 normExpSig = softfloat_normSubnormalF128Sig( sigA.v64, sigA.v0 ); 107 expA = normExpSig.exp; 108 sigA = normExpSig.sig; 109 } 110 /*------------------------------------------------------------------------ 111 *------------------------------------------------------------------------*/ 112 sigA.v64 |= UINT64_C( 0x0001000000000000 ); 113 sigB.v64 |= UINT64_C( 0x0001000000000000 ); 114 rem = sigA; 115 expDiff = expA - expB; 116 if ( expDiff < 1 ) { 117 if ( expDiff < -1 ) return a; 118 if ( expDiff ) { 119 --expB; 120 sigB = softfloat_add128( sigB.v64, sigB.v0, sigB.v64, sigB.v0 ); 121 q = 0; 122 } else { 123 q = softfloat_le128( sigB.v64, sigB.v0, rem.v64, rem.v0 ); 124 if ( q ) { 125 rem = softfloat_sub128( rem.v64, rem.v0, sigB.v64, sigB.v0 ); 126 } 127 } 128 } else { 129 recip32 = softfloat_approxRecip32_1( sigB.v64>>17 ); 130 expDiff -= 30; 131 for (;;) { 132 q64 = (uint_fast64_t) (uint32_t) (rem.v64>>19) * recip32; 133 if ( expDiff < 0 ) break; 134 q = (q64 + 0x80000000)>>32; 135 rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, 29 ); 136 term = softfloat_mul128By32( sigB.v64, sigB.v0, q ); 137 rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 ); 138 if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) { 139 rem = softfloat_add128( rem.v64, rem.v0, sigB.v64, sigB.v0 ); 140 } 141 expDiff -= 29; 142 } 143 /*-------------------------------------------------------------------- 144 | (`expDiff' cannot be less than -29 here.) 145 *--------------------------------------------------------------------*/ 146 q = (uint32_t) (q64>>32)>>(~expDiff & 31); 147 rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, expDiff + 30 ); 148 term = softfloat_mul128By32( sigB.v64, sigB.v0, q ); 149 rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 ); 150 if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) { 151 altRem = softfloat_add128( rem.v64, rem.v0, sigB.v64, sigB.v0 ); 152 goto selectRem; 153 } 154 } 155 /*------------------------------------------------------------------------ 156 *------------------------------------------------------------------------*/ 157 do { 158 altRem = rem; 159 ++q; 160 rem = softfloat_sub128( rem.v64, rem.v0, sigB.v64, sigB.v0 ); 161 } while ( ! (rem.v64 & UINT64_C( 0x8000000000000000 )) ); 162 selectRem: 163 meanRem = softfloat_add128( rem.v64, rem.v0, altRem.v64, altRem.v0 ); 164 if ( 165 (meanRem.v64 & UINT64_C( 0x8000000000000000 )) 166 || (! (meanRem.v64 | meanRem.v0) && (q & 1)) 167 ) { 168 rem = altRem; 169 } 170 signRem = signA; 171 if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) { 172 signRem = ! signRem; 173 rem = softfloat_sub128( 0, 0, rem.v64, rem.v0 ); 174 } 175 return softfloat_normRoundPackToF128( signRem, expB - 1, rem.v64, rem.v0 ); 176 /*------------------------------------------------------------------------ 177 *------------------------------------------------------------------------*/ 178 propagateNaN: 179 uiZ = softfloat_propagateNaNF128UI( uiA64, uiA0, uiB64, uiB0 ); 180 goto uiZ; 181 /*------------------------------------------------------------------------ 182 *------------------------------------------------------------------------*/ 183 invalid: 184 softfloat_raiseFlags( softfloat_flag_invalid ); 185 uiZ.v64 = defaultNaNF128UI64; 186 uiZ.v0 = defaultNaNF128UI0; 187 uiZ: 188 uZ.ui = uiZ; 189 return uZ.f; 190 191 } 192 193