softfloat/f64_div.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_div( 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     bool signZ;
  23     struct exp16_sig64 normExpSig;
  24     int_fast16_t expZ;
  25     uint_fast64_t sigZ;
  26     struct uint128 term, rem;
  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     signZ = signA ^ signB;
  41     if ( expA == 0x7FF ) {
  42         if ( sigA ) goto propagateNaN;
  43         if ( expB == 0x7FF ) {
  44             if ( sigB ) goto propagateNaN;
  45             goto invalid;
  46         }
  47         goto infinity;
  48     }
  49     if ( expB == 0x7FF ) {
  50         if ( sigB ) goto propagateNaN;
  51         goto zero;
  52     }
  53     if ( ! expB ) {
  54         if ( ! sigB ) {
  55             if ( ! ( expA | sigA ) ) goto invalid;
  56             softfloat_raiseFlags( softfloat_flag_infinity );
  57             goto infinity;
  58         }
  59         normExpSig = softfloat_normSubnormalF64Sig( sigB );
  60         expB = normExpSig.exp;
  61         sigB = normExpSig.sig;
  62     }
  63     if ( ! expA ) {
  64         if ( ! sigA ) goto zero;
  65         normExpSig = softfloat_normSubnormalF64Sig( sigA );
  66         expA = normExpSig.exp;
  67         sigA = normExpSig.sig;
  68     }
  69     expZ = expA - expB + 0x3FD;
  70     sigA = ( sigA | UINT64_C( 0x0010000000000000 ) )<<10;
  71     sigB = ( sigB | UINT64_C( 0x0010000000000000 ) )<<11;
  72     if ( sigB <= ( sigA + sigA ) ) {
  73         ++expZ;
  74         sigA >>= 1;
  75     }
  76     sigZ = softfloat_estimateDiv128To64( sigA, 0, sigB );
  77     if ( ( sigZ & 0x1FF ) <= 2 ) {
  78         term = softfloat_mul64To128( sigB, sigZ );
  79         rem = softfloat_sub128( sigA, 0, term.v64, term.v0 );
  80         while ( UINT64_C( 0x8000000000000000 ) <= rem.v64 ) {
  81             --sigZ;
  82             rem = softfloat_add128( rem.v64, rem.v0, 0, sigB );
  83         }
  84         sigZ |= ( rem.v0 != 0 );
  85     }
  86     return softfloat_roundPackToF64( signZ, expZ, sigZ );
  87  propagateNaN:
  88     uiZ = softfloat_propagateNaNF64UI( uiA, uiB );
  89     goto uiZ;
  90  invalid:
  91     softfloat_raiseFlags( softfloat_flag_invalid );
  92     uiZ = defaultNaNF64UI;
  93     goto uiZ;
  94  infinity:
  95     uiZ = packToF64UI( signZ, 0x7FF, 0 );
  96     goto uiZ;
  97  zero:
  98     uiZ = packToF64UI( signZ, 0, 0 );
  99  uiZ:
 100     uZ.ui = uiZ;
 101     return uZ.f;
 102
 103 }
 104