aboutsummaryrefslogtreecommitdiff
path: root/softfloat/f64_rem.c
diff options
context:
space:
mode:
Diffstat (limited to 'softfloat/f64_rem.c')
-rw-r--r--[-rwxr-xr-x]softfloat/f64_rem.c160
1 files changed, 117 insertions, 43 deletions
diff --git a/softfloat/f64_rem.c b/softfloat/f64_rem.c
index f738e5c..ca5350c 100755..100644
--- a/softfloat/f64_rem.c
+++ b/softfloat/f64_rem.c
@@ -1,10 +1,42 @@
-// See LICENSE.SoftFloat for license details.
+/*============================================================================
+
+This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic
+Package, Release 3e, by John R. Hauser.
+
+Copyright 2011, 2012, 2013, 2014 The Regents of the University of California.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+ 1. Redistributions of source code must retain the above copyright notice,
+ this list of conditions, and the following disclaimer.
+
+ 2. Redistributions in binary form must reproduce the above copyright notice,
+ this list of conditions, and the following disclaimer in the documentation
+ and/or other materials provided with the distribution.
+
+ 3. Neither the name of the University nor the names of its contributors may
+ be used to endorse or promote products derived from this software without
+ specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
+DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+=============================================================================*/
#include <stdbool.h>
#include <stdint.h>
#include "platform.h"
-#include "primitives.h"
#include "internals.h"
#include "specialize.h"
#include "softfloat.h"
@@ -18,35 +50,44 @@ float64_t f64_rem( float64_t a, float64_t b )
uint_fast64_t sigA;
union ui64_f64 uB;
uint_fast64_t uiB;
- bool signB;
int_fast16_t expB;
uint_fast64_t sigB;
struct exp16_sig64 normExpSig;
+ uint64_t rem;
int_fast16_t expDiff;
- uint_fast64_t q, alternateSigA;
- uint64_t sigMean;
- bool signZ;
+ uint32_t q, recip32;
+ uint_fast64_t q64;
+ uint64_t altRem, meanRem;
+ bool signRem;
uint_fast64_t uiZ;
union ui64_f64 uZ;
+ /*------------------------------------------------------------------------
+ *------------------------------------------------------------------------*/
uA.f = a;
uiA = uA.ui;
signA = signF64UI( uiA );
- expA = expF64UI( uiA );
- sigA = fracF64UI( uiA );
+ expA = expF64UI( uiA );
+ sigA = fracF64UI( uiA );
uB.f = b;
uiB = uB.ui;
- signB = signF64UI( uiB );
expB = expF64UI( uiB );
sigB = fracF64UI( uiB );
+ /*------------------------------------------------------------------------
+ *------------------------------------------------------------------------*/
if ( expA == 0x7FF ) {
- if ( sigA || ( ( expB == 0x7FF ) && sigB ) ) goto propagateNaN;
+ if ( sigA || ((expB == 0x7FF) && sigB) ) goto propagateNaN;
goto invalid;
}
if ( expB == 0x7FF ) {
if ( sigB ) goto propagateNaN;
return a;
}
+ /*------------------------------------------------------------------------
+ *------------------------------------------------------------------------*/
+ if ( expA < expB - 1 ) return a;
+ /*------------------------------------------------------------------------
+ *------------------------------------------------------------------------*/
if ( ! expB ) {
if ( ! sigB ) goto invalid;
normExpSig = softfloat_normSubnormalF64Sig( sigB );
@@ -59,48 +100,81 @@ float64_t f64_rem( float64_t a, float64_t b )
expA = normExpSig.exp;
sigA = normExpSig.sig;
}
+ /*------------------------------------------------------------------------
+ *------------------------------------------------------------------------*/
+ rem = sigA | UINT64_C( 0x0010000000000000 );
+ sigB |= UINT64_C( 0x0010000000000000 );
expDiff = expA - expB;
- sigA = ( sigA | UINT64_C( 0x0010000000000000 ) )<<11;
- sigB = ( sigB | UINT64_C( 0x0010000000000000 ) )<<11;
- if ( expDiff < 0 ) {
+ if ( expDiff < 1 ) {
if ( expDiff < -1 ) return a;
- sigA >>= 1;
- }
- q = ( sigB <= sigA );
- if ( q ) sigA -= sigB;
- expDiff -= 64;
- while ( 0 < expDiff ) {
- q = softfloat_estimateDiv128To64( sigA, 0, sigB );
- q = ( 2 < q ) ? q - 2 : 0;
- sigA = - ( ( sigB>>2 ) * q );
- expDiff -= 62;
- }
- expDiff += 64;
- if ( 0 < expDiff ) {
- q = softfloat_estimateDiv128To64( sigA, 0, sigB );
- q = ( 2 < q ) ? q - 2 : 0;
- q >>= 64 - expDiff;
- sigB >>= 2;
- sigA = ( ( sigA>>1 )<<( expDiff - 1 ) ) - sigB * q;
+ sigB <<= 9;
+ if ( expDiff ) {
+ rem <<= 8;
+ q = 0;
+ } else {
+ rem <<= 9;
+ q = (sigB <= rem);
+ if ( q ) rem -= sigB;
+ }
} else {
- sigA >>= 2;
- sigB >>= 2;
+ recip32 = softfloat_approxRecip32_1( sigB>>21 );
+ /*--------------------------------------------------------------------
+ | Changing the shift of `rem' here requires also changing the initial
+ | subtraction from `expDiff'.
+ *--------------------------------------------------------------------*/
+ rem <<= 9;
+ expDiff -= 30;
+ /*--------------------------------------------------------------------
+ | The scale of `sigB' affects how many bits are obtained during each
+ | cycle of the loop. Currently this is 29 bits per loop iteration,
+ | the maximum possible.
+ *--------------------------------------------------------------------*/
+ sigB <<= 9;
+ for (;;) {
+ q64 = (uint32_t) (rem>>32) * (uint_fast64_t) recip32;
+ if ( expDiff < 0 ) break;
+ q = (q64 + 0x80000000)>>32;
+#ifdef SOFTFLOAT_FAST_INT64
+ rem <<= 29;
+#else
+ rem = (uint_fast64_t) (uint32_t) (rem>>3)<<32;
+#endif
+ rem -= q * (uint64_t) sigB;
+ if ( rem & UINT64_C( 0x8000000000000000 ) ) rem += sigB;
+ expDiff -= 29;
+ }
+ /*--------------------------------------------------------------------
+ | (`expDiff' cannot be less than -29 here.)
+ *--------------------------------------------------------------------*/
+ q = (uint32_t) (q64>>32)>>(~expDiff & 31);
+ rem = (rem<<(expDiff + 30)) - q * (uint64_t) sigB;
+ if ( rem & UINT64_C( 0x8000000000000000 ) ) {
+ altRem = rem + sigB;
+ goto selectRem;
+ }
}
+ /*------------------------------------------------------------------------
+ *------------------------------------------------------------------------*/
do {
- alternateSigA = sigA;
+ altRem = rem;
++q;
- sigA -= sigB;
- } while ( sigA < UINT64_C( 0x8000000000000000 ) );
- sigMean = sigA + alternateSigA;
+ rem -= sigB;
+ } while ( ! (rem & UINT64_C( 0x8000000000000000 )) );
+ selectRem:
+ meanRem = rem + altRem;
if (
- ( UINT64_C( 0x8000000000000000 ) <= sigMean )
- || ( ! sigMean && ( q & 1 ) )
+ (meanRem & UINT64_C( 0x8000000000000000 )) || (! meanRem && (q & 1))
) {
- sigA = alternateSigA;
+ rem = altRem;
+ }
+ signRem = signA;
+ if ( rem & UINT64_C( 0x8000000000000000 ) ) {
+ signRem = ! signRem;
+ rem = -rem;
}
- signZ = ( UINT64_C( 0x8000000000000000 ) <= sigA );
- if ( signZ ) sigA = - sigA;
- return softfloat_normRoundPackToF64( signA ^ signZ, expB, sigA );
+ return softfloat_normRoundPackToF64( signRem, expB, rem );
+ /*------------------------------------------------------------------------
+ *------------------------------------------------------------------------*/
propagateNaN:
uiZ = softfloat_propagateNaNF64UI( uiA, uiB );
goto uiZ;
generated by cgit v1.1 at 2024-10-10 01:25:14 +0000