Upgrade to SoftFloat 3e

1 files changed, 96 insertions, 54 deletions

diff --git a/softfloat/f32_rem.c b/softfloat/f32_rem.c
index 598ee8a..2d74c8c 100755..100644
--- 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 <stdbool.h>
 #include <stdint.h>
 #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;