From 5cd8d2b93fad2011fdf8f49a263d77697af35b32 Mon Sep 17 00:00:00 2001 From: Andrew Waterman Date: Wed, 11 Jul 2018 01:16:28 -0700 Subject: Upgrade to SoftFloat 3e --- softfloat/f32_rem.c | 150 +++++++++++++++++++++++++++++++++------------------- 1 file changed, 96 insertions(+), 54 deletions(-) mode change 100755 => 100644 softfloat/f32_rem.c (limited to 'softfloat/f32_rem.c') diff --git a/softfloat/f32_rem.c b/softfloat/f32_rem.c old mode 100755 new mode 100644 index 598ee8a..2d74c8c --- a/softfloat/f32_rem.c +++ b/softfloat/f32_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 #include #include "platform.h" -#include "primitives.h" #include "internals.h" #include "specialize.h" #include "softfloat.h" @@ -18,37 +50,39 @@ float32_t f32_rem( float32_t a, float32_t b ) uint_fast32_t sigA; union ui32_f32 uB; uint_fast32_t uiB; - bool signB; int_fast16_t expB; uint_fast32_t sigB; struct exp16_sig32 normExpSig; + uint32_t rem; int_fast16_t expDiff; - uint_fast32_t q; - uint_fast64_t sigA64, sigB64, q64; - uint_fast32_t alternateSigA; - uint32_t sigMean; - bool signZ; + uint32_t q, recip32, altRem, meanRem; + bool signRem; uint_fast32_t uiZ; union ui32_f32 uZ; + /*------------------------------------------------------------------------ + *------------------------------------------------------------------------*/ uA.f = a; uiA = uA.ui; signA = signF32UI( uiA ); - expA = expF32UI( uiA ); - sigA = fracF32UI( uiA ); + expA = expF32UI( uiA ); + sigA = fracF32UI( uiA ); uB.f = b; uiB = uB.ui; - signB = signF32UI( uiB ); expB = expF32UI( uiB ); sigB = fracF32UI( uiB ); + /*------------------------------------------------------------------------ + *------------------------------------------------------------------------*/ if ( expA == 0xFF ) { - if ( sigA || ( ( expB == 0xFF ) && sigB ) ) goto propagateNaN; + if ( sigA || ((expB == 0xFF) && sigB) ) goto propagateNaN; goto invalid; } if ( expB == 0xFF ) { if ( sigB ) goto propagateNaN; return a; } + /*------------------------------------------------------------------------ + *------------------------------------------------------------------------*/ if ( ! expB ) { if ( ! sigB ) goto invalid; normExpSig = softfloat_normSubnormalF32Sig( sigB ); @@ -61,57 +95,65 @@ float32_t f32_rem( float32_t a, float32_t b ) expA = normExpSig.exp; sigA = normExpSig.sig; } - expDiff = expA - expB; - sigA |= 0x00800000; + /*------------------------------------------------------------------------ + *------------------------------------------------------------------------*/ + rem = sigA | 0x00800000; sigB |= 0x00800000; - if ( expDiff < 32 ) { - sigA <<= 8; - sigB <<= 8; - if ( expDiff < 0 ) { - if ( expDiff < -1 ) return a; - sigA >>= 1; - } - q = ( sigB <= sigA ); - if ( q ) sigA -= sigB; - if ( 0 < expDiff ) { - q = ( (uint_fast64_t) sigA<<32 ) / sigB; - q >>= 32 - expDiff; - sigB >>= 2; - sigA = ( ( sigA>>1 )<<( expDiff - 1 ) ) - sigB * q; + expDiff = expA - expB; + if ( expDiff < 1 ) { + if ( expDiff < -1 ) return a; + sigB <<= 6; + if ( expDiff ) { + rem <<= 5; + q = 0; } else { - sigA >>= 2; - sigB >>= 2; + rem <<= 6; + q = (sigB <= rem); + if ( q ) rem -= sigB; } } else { - if ( sigB <= sigA ) sigA -= sigB; - sigA64 = (uint_fast64_t) sigA<<40; - sigB64 = (uint_fast64_t) sigB<<40; - expDiff -= 64; - while ( 0 < expDiff ) { - q64 = softfloat_estimateDiv128To64( sigA64, 0, sigB64 ); - q64 = ( 2 < q64 ) ? q64 - 2 : 0; - sigA64 = - ( ( sigB * q64 )<<38 ); - expDiff -= 62; - } - expDiff += 64; - q64 = softfloat_estimateDiv128To64( sigA64, 0, sigB64 ); - q64 = ( 2 < q64 ) ? q64 - 2 : 0; - q = q64>>( 64 - expDiff ); + recip32 = softfloat_approxRecip32_1( sigB<<8 ); + /*-------------------------------------------------------------------- + | Changing the shift of `rem' here requires also changing the initial + | subtraction from `expDiff'. + *--------------------------------------------------------------------*/ + rem <<= 7; + expDiff -= 31; + /*-------------------------------------------------------------------- + | The scale of `sigB' affects how many bits are obtained during each + | cycle of the loop. Currently this is 29 bits per loop iteration, + | which is believed to be the maximum possible. + *--------------------------------------------------------------------*/ sigB <<= 6; - sigA = ( ( sigA64>>33 )<<( expDiff - 1 ) ) - sigB * q; + for (;;) { + q = (rem * (uint_fast64_t) recip32)>>32; + if ( expDiff < 0 ) break; + rem = -(q * (uint32_t) sigB); + expDiff -= 29; + } + /*-------------------------------------------------------------------- + | (`expDiff' cannot be less than -30 here.) + *--------------------------------------------------------------------*/ + q >>= ~expDiff & 31; + rem = (rem<<(expDiff + 30)) - q * (uint32_t) sigB; } + /*------------------------------------------------------------------------ + *------------------------------------------------------------------------*/ do { - alternateSigA = sigA; + altRem = rem; ++q; - sigA -= sigB; - } while ( sigA < 0x80000000 ); - sigMean = sigA + alternateSigA; - if ( ( 0x80000000 <= sigMean ) || ( ! sigMean && ( q & 1 ) ) ) { - sigA = alternateSigA; + rem -= sigB; + } while ( ! (rem & 0x80000000) ); + meanRem = rem + altRem; + if ( (meanRem & 0x80000000) || (! meanRem && (q & 1)) ) rem = altRem; + signRem = signA; + if ( 0x80000000 <= rem ) { + signRem = ! signRem; + rem = -rem; } - signZ = ( 0x80000000 <= sigA ); - if ( signZ ) sigA = - sigA; - return softfloat_normRoundPackToF32( signA ^ signZ, expB, sigA ); + return softfloat_normRoundPackToF32( signRem, expB, rem ); + /*------------------------------------------------------------------------ + *------------------------------------------------------------------------*/ propagateNaN: uiZ = softfloat_propagateNaNF32UI( uiA, uiB ); goto uiZ; -- cgit v1.1