softfloat/source/extF80M_rem.c

   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 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 #ifdef SOFTFLOAT_FAST_INT64
  45
  46 void
  47  extF80M_rem(
  48      const extFloat80_t *aPtr, const extFloat80_t *bPtr, extFloat80_t *zPtr )
  49 {
  50
  51     *zPtr = extF80_rem( *aPtr, *bPtr );
  52
  53 }
  54
  55 #else
  56
  57 void
  58  extF80M_rem(
  59      const extFloat80_t *aPtr, const extFloat80_t *bPtr, extFloat80_t *zPtr )
  60 {
  61     const struct extFloat80M *aSPtr, *bSPtr;
  62     struct extFloat80M *zSPtr;
  63     uint_fast16_t uiA64;
  64     int32_t expA, expB;
  65     uint64_t x64;
  66     bool signRem;
  67     uint64_t sigA;
  68     int32_t expDiff;
  69     uint32_t rem[3], x[3], sig32B, q, recip32, rem2[3], *remPtr, *altRemPtr;
  70     uint32_t *newRemPtr, wordMeanRem;
  71
  72     /*------------------------------------------------------------------------
  73     *------------------------------------------------------------------------*/
  74     aSPtr = (const struct extFloat80M *) aPtr;
  75     bSPtr = (const struct extFloat80M *) bPtr;
  76     zSPtr = (struct extFloat80M *) zPtr;
  77     /*------------------------------------------------------------------------
  78     *------------------------------------------------------------------------*/
  79     uiA64 = aSPtr->signExp;
  80     expA = expExtF80UI64( uiA64 );
  81     expB = expExtF80UI64( bSPtr->signExp );
  82     /*------------------------------------------------------------------------
  83     *------------------------------------------------------------------------*/
  84     if ( (expA == 0x7FFF) || (expB == 0x7FFF) ) {
  85         if ( softfloat_tryPropagateNaNExtF80M( aSPtr, bSPtr, zSPtr ) ) return;
  86         if ( expA == 0x7FFF ) goto invalid;
  87         /*--------------------------------------------------------------------
  88         | If we get here, then argument b is an infinity and `expB' is 0x7FFF;
  89         | Doubling `expB' is an easy way to ensure that `expDiff' later is
  90         | less than -1, which will result in returning a canonicalized version
  91         | of argument a.
  92         *--------------------------------------------------------------------*/
  93         expB += expB;
  94     }
  95     /*------------------------------------------------------------------------
  96     *------------------------------------------------------------------------*/
  97     if ( ! expB ) expB = 1;
  98     x64 = bSPtr->signif;
  99     if ( ! (x64 & UINT64_C( 0x8000000000000000 )) ) {
 100         if ( ! x64 ) goto invalid;
 101         expB += softfloat_normExtF80SigM( &x64 );
 102     }
 103     signRem = signExtF80UI64( uiA64 );
 104     if ( ! expA ) expA = 1;
 105     sigA = aSPtr->signif;
 106     if ( ! (sigA & UINT64_C( 0x8000000000000000 )) ) {
 107         if ( ! sigA ) {
 108             expA = 0;
 109             goto copyA;
 110         }
 111         expA += softfloat_normExtF80SigM( &sigA );
 112     }
 113     /*------------------------------------------------------------------------
 114     *------------------------------------------------------------------------*/
 115     expDiff = expA - expB;
 116     if ( expDiff < -1 ) goto copyA;
 117     rem[indexWord( 3, 2 )] = sigA>>34;
 118     rem[indexWord( 3, 1 )] = sigA>>2;
 119     rem[indexWord( 3, 0 )] = (uint32_t) sigA<<30;
 120     x[indexWord( 3, 0 )] = (uint32_t) x64<<30;
 121     sig32B = x64>>32;
 122     x64 >>= 2;
 123     x[indexWord( 3, 2 )] = x64>>32;
 124     x[indexWord( 3, 1 )] = x64;
 125     if ( expDiff < 1 ) {
 126         if ( expDiff ) {
 127             --expB;
 128             softfloat_add96M( x, x, x );
 129             q = 0;
 130         } else {
 131             q = (softfloat_compare96M( x, rem ) <= 0);
 132             if ( q ) softfloat_sub96M( rem, x, rem );
 133         }
 134     } else {
 135         recip32 = softfloat_approxRecip32_1( sig32B );
 136         expDiff -= 30;
 137         for (;;) {
 138             x64 = (uint64_t) rem[indexWordHi( 3 )] * recip32;
 139             if ( expDiff < 0 ) break;
 140             q = (x64 + 0x80000000)>>32;
 141             softfloat_remStep96MBy32( rem, 29, x, q, rem );
 142             if ( rem[indexWordHi( 3 )] & 0x80000000 ) {
 143                 softfloat_add96M( rem, x, rem );
 144             }
 145             expDiff -= 29;
 146         }
 147         /*--------------------------------------------------------------------
 148         | (`expDiff' cannot be less than -29 here.)
 149         *--------------------------------------------------------------------*/
 150         q = (uint32_t) (x64>>32)>>(~expDiff & 31);
 151         softfloat_remStep96MBy32( rem, expDiff + 30, x, q, rem );
 152         if ( rem[indexWordHi( 3 )] & 0x80000000 ) {
 153             remPtr = rem;
 154             altRemPtr = rem2;
 155             softfloat_add96M( remPtr, x, altRemPtr );
 156             goto selectRem;
 157         }
 158     }
 159     /*------------------------------------------------------------------------
 160     *------------------------------------------------------------------------*/
 161     remPtr = rem;
 162     altRemPtr = rem2;
 163     do {
 164         ++q;
 165         newRemPtr = altRemPtr;
 166         softfloat_sub96M( remPtr, x, newRemPtr );
 167         altRemPtr = remPtr;
 168         remPtr = newRemPtr;
 169     } while ( ! (remPtr[indexWordHi( 3 )] & 0x80000000) );
 170  selectRem:
 171     softfloat_add96M( remPtr, altRemPtr, x );
 172     wordMeanRem = x[indexWordHi( 3 )];
 173     if (
 174         (wordMeanRem & 0x80000000)
 175             || (! wordMeanRem && (q & 1) && ! x[indexWord( 3, 0 )]
 176                     && ! x[indexWord( 3, 1 )])
 177     ) {
 178         remPtr = altRemPtr;
 179     }
 180     if ( remPtr[indexWordHi( 3 )] & 0x80000000 ) {
 181         signRem = ! signRem;
 182         softfloat_negX96M( remPtr );
 183     }
 184     softfloat_normRoundPackMToExtF80M( signRem, expB + 2, remPtr, 80, zSPtr );
 185     return;
 186     /*------------------------------------------------------------------------
 187     *------------------------------------------------------------------------*/
 188  invalid:
 189     softfloat_invalidExtF80M( zSPtr );
 190     return;
 191     /*------------------------------------------------------------------------
 192     *------------------------------------------------------------------------*/
 193  copyA:
 194     if ( expA < 1 ) {
 195         sigA >>= 1 - expA;
 196         expA = 0;
 197     }
 198     zSPtr->signExp = packToExtF80UI64( signRem, expA );
 199     zSPtr->signif = sigA;
 200
 201 }
 202
 203 #endif
 204