softfloat/f64_rem.c

   1
   2 #include <stdbool.h>
   3 #include <stdint.h>
   4 #include "platform.h"
   5 #include "primitives.h"
   6 #include "internals.h"
   7 #include "specialize.h"
   8 #include "softfloat.h"
   9
  10 float64_t f64_rem( float64_t a, float64_t b )
  11 {
  12     union ui64_f64 uA;
  13     uint_fast64_t uiA;
  14     bool signA;
  15     int_fast16_t expA;
  16     uint_fast64_t sigA;
  17     union ui64_f64 uB;
  18     uint_fast64_t uiB;
  19     bool signB;
  20     int_fast16_t expB;
  21     uint_fast64_t sigB;
  22     struct exp16_sig64 normExpSig;
  23     int_fast16_t expDiff;
  24     uint_fast64_t q, alternateSigA;
  25     uint64_t sigMean;
  26     bool signZ;
  27     uint_fast64_t uiZ;
  28     union ui64_f64 uZ;
  29
  30     uA.f = a;
  31     uiA = uA.ui;
  32     signA = signF64UI( uiA );
  33     expA = expF64UI( uiA );
  34     sigA = fracF64UI( uiA );
  35     uB.f = b;
  36     uiB = uB.ui;
  37     signB = signF64UI( uiB );
  38     expB = expF64UI( uiB );
  39     sigB = fracF64UI( uiB );
  40     if ( expA == 0x7FF ) {
  41         if ( sigA || ( ( expB == 0x7FF ) && sigB ) ) goto propagateNaN;
  42         goto invalid;
  43     }
  44     if ( expB == 0x7FF ) {
  45         if ( sigB ) goto propagateNaN;
  46         return a;
  47     }
  48     if ( ! expB ) {
  49         if ( ! sigB ) goto invalid;
  50         normExpSig = softfloat_normSubnormalF64Sig( sigB );
  51         expB = normExpSig.exp;
  52         sigB = normExpSig.sig;
  53     }
  54     if ( ! expA ) {
  55         if ( ! sigA ) return a;
  56         normExpSig = softfloat_normSubnormalF64Sig( sigA );
  57         expA = normExpSig.exp;
  58         sigA = normExpSig.sig;
  59     }
  60     expDiff = expA - expB;
  61     sigA = ( sigA | UINT64_C( 0x0010000000000000 ) )<<11;
  62     sigB = ( sigB | UINT64_C( 0x0010000000000000 ) )<<11;
  63     if ( expDiff < 0 ) {
  64         if ( expDiff < -1 ) return a;
  65         sigA >>= 1;
  66     }
  67     q = ( sigB <= sigA );
  68     if ( q ) sigA -= sigB;
  69     expDiff -= 64;
  70     while ( 0 < expDiff ) {
  71         q = softfloat_estimateDiv128To64( sigA, 0, sigB );
  72         q = ( 2 < q ) ? q - 2 : 0;
  73         sigA = - ( ( sigB>>2 ) * q );
  74         expDiff -= 62;
  75     }
  76     expDiff += 64;
  77     if ( 0 < expDiff ) {
  78         q = softfloat_estimateDiv128To64( sigA, 0, sigB );
  79         q = ( 2 < q ) ? q - 2 : 0;
  80         q >>= 64 - expDiff;
  81         sigB >>= 2;
  82         sigA = ( ( sigA>>1 )<<( expDiff - 1 ) ) - sigB * q;
  83     } else {
  84         sigA >>= 2;
  85         sigB >>= 2;
  86     }
  87     do {
  88         alternateSigA = sigA;
  89         ++q;
  90         sigA -= sigB;
  91     } while ( sigA < UINT64_C( 0x8000000000000000 ) );
  92     sigMean = sigA + alternateSigA;
  93     if (
  94         ( UINT64_C( 0x8000000000000000 ) <= sigMean )
  95             || ( ! sigMean && ( q & 1 ) )
  96     ) {
  97         sigA = alternateSigA;
  98     }
  99     signZ = ( UINT64_C( 0x8000000000000000 ) <= sigA );
 100     if ( signZ ) sigA = - sigA;
 101     return softfloat_normRoundPackToF64( signA ^ signZ, expB, sigA );
 102  propagateNaN:
 103     uiZ = softfloat_propagateNaNF64UI( uiA, uiB );
 104     goto uiZ;
 105  invalid:
 106     softfloat_raiseFlags( softfloat_flag_invalid );
 107     uiZ = defaultNaNF64UI;
 108  uiZ:
 109     uZ.ui = uiZ;
 110     return uZ.f;
 111
 112 }
 113