softfloat/f32_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 float32_t f32_rem( float32_t a, float32_t b )
  11 {
  12     union ui32_f32 uA;
  13     uint_fast32_t uiA;
  14     bool signA;
  15     int_fast16_t expA;
  16     uint_fast32_t sigA;
  17     union ui32_f32 uB;
  18     uint_fast32_t uiB;
  19     bool signB;
  20     int_fast16_t expB;
  21     uint_fast32_t sigB;
  22     struct exp16_sig32 normExpSig;
  23     int_fast16_t expDiff;
  24     uint_fast32_t q;
  25     uint_fast64_t sigA64, sigB64, q64;
  26     uint_fast32_t alternateSigA;
  27     uint32_t sigMean;
  28     bool signZ;
  29     uint_fast32_t uiZ;
  30     union ui32_f32 uZ;
  31
  32     uA.f = a;
  33     uiA = uA.ui;
  34     signA = signF32UI( uiA );
  35     expA = expF32UI( uiA );
  36     sigA = fracF32UI( uiA );
  37     uB.f = b;
  38     uiB = uB.ui;
  39     signB = signF32UI( uiB );
  40     expB = expF32UI( uiB );
  41     sigB = fracF32UI( uiB );
  42     if ( expA == 0xFF ) {
  43         if ( sigA || ( ( expB == 0xFF ) && sigB ) ) goto propagateNaN;
  44         goto invalid;
  45     }
  46     if ( expB == 0xFF ) {
  47         if ( sigB ) goto propagateNaN;
  48         return a;
  49     }
  50     if ( ! expB ) {
  51         if ( ! sigB ) goto invalid;
  52         normExpSig = softfloat_normSubnormalF32Sig( sigB );
  53         expB = normExpSig.exp;
  54         sigB = normExpSig.sig;
  55     }
  56     if ( ! expA ) {
  57         if ( ! sigA ) return a;
  58         normExpSig = softfloat_normSubnormalF32Sig( sigA );
  59         expA = normExpSig.exp;
  60         sigA = normExpSig.sig;
  61     }
  62     expDiff = expA - expB;
  63     sigA |= 0x00800000;
  64     sigB |= 0x00800000;
  65     if ( expDiff < 32 ) {
  66         sigA <<= 8;
  67         sigB <<= 8;
  68         if ( expDiff < 0 ) {
  69             if ( expDiff < -1 ) return a;
  70             sigA >>= 1;
  71         }
  72         q = ( sigB <= sigA );
  73         if ( q ) sigA -= sigB;
  74         if ( 0 < expDiff ) {
  75             q = ( (uint_fast64_t) sigA<<32 ) / sigB;
  76             q >>= 32 - expDiff;
  77             sigB >>= 2;
  78             sigA = ( ( sigA>>1 )<<( expDiff - 1 ) ) - sigB * q;
  79         } else {
  80             sigA >>= 2;
  81             sigB >>= 2;
  82         }
  83     } else {
  84         if ( sigB <= sigA ) sigA -= sigB;
  85         sigA64 = (uint_fast64_t) sigA<<40;
  86         sigB64 = (uint_fast64_t) sigB<<40;
  87         expDiff -= 64;
  88         while ( 0 < expDiff ) {
  89             q64 = softfloat_estimateDiv128To64( sigA64, 0, sigB64 );
  90             q64 = ( 2 < q64 ) ? q64 - 2 : 0;
  91             sigA64 = - ( ( sigB * q64 )<<38 );
  92             expDiff -= 62;
  93         }
  94         expDiff += 64;
  95         q64 = softfloat_estimateDiv128To64( sigA64, 0, sigB64 );
  96         q64 = ( 2 < q64 ) ? q64 - 2 : 0;
  97         q = q64>>( 64 - expDiff );
  98         sigB <<= 6;
  99         sigA = ( ( sigA64>>33 )<<( expDiff - 1 ) ) - sigB * q;
 100     }
 101     do {
 102         alternateSigA = sigA;
 103         ++q;
 104         sigA -= sigB;
 105     } while ( sigA < 0x80000000 );
 106     sigMean = sigA + alternateSigA;
 107     if ( ( 0x80000000 <= sigMean ) || ( ! sigMean && ( q & 1 ) ) ) {
 108         sigA = alternateSigA;
 109     }
 110     signZ = ( 0x80000000 <= sigA );
 111     if ( signZ ) sigA = - sigA;
 112     return softfloat_normRoundPackToF32( signA ^ signZ, expB, sigA );
 113  propagateNaN:
 114     uiZ = softfloat_propagateNaNF32UI( uiA, uiB );
 115     goto uiZ;
 116  invalid:
 117     softfloat_raiseFlags( softfloat_flag_invalid );
 118     uiZ = defaultNaNF32UI;
 119  uiZ:
 120     uZ.ui = uiZ;
 121     return uZ.f;
 122
 123 }
 124