FMA based division

Below is a source code for FMA-based division algorithm with integer reciprocal estimation. This implementation gives mostly correctly-rounded results, but not perfect.

#include <math.h>

static double mulsign(double x, double y) {
  union {
    double f;
    uint64_t ul;
  } cnvx, cnvy;

  cnvx.f = x;
  cnvy.f = y;
  cnvy.ul &= 0x8000000000000000ULL;
  cnvx.ul ^= cnvy.ul;
  return cnvx.f;
}

static inline double divide(double x, double y) {
  union {
    double f;
    uint64_t ul;
  } cnvx, cnvy;

  // Adjust the exponents of arguments

  cnvy.f = y;
  cnvy.ul &= 0x7fc0000000000000ULL;

  cnvx.f = x;
  cnvx.ul &= 0x7fc0000000000000ULL;

  cnvx.ul =  0x7fd0000000000000ULL - ((cnvy.ul + cnvx.ul) >> 1);

  y *= cnvx.f;
  x *= cnvx.f;

  //

  cnvy.f = y;
  int o = (cnvy.ul & 0x7fc0000000000000ULL) >= 0x7fc0000000000000ULL;

  // Reciprocal estimation

  int32_t i = 0x7fff & (int32_t)(cnvy.ul >> (52-15));
  cnvy.ul = (uint64_t)0x7fe0000000000000LL - cnvy.ul;

  int d;
  d = (8432 * i - 763081686) >> 14;
  d = (d    * i + 982096282) >> 15;
  d = (d    * i -   2075152);

  cnvy.ul = o ? 0 : (cnvy.ul - ((uint64_t)d << (52 - 23 - 7)));

  // Newton-Raphson

  double t = cnvy.f;
  t = fma(t, fma(t, -y, 1), t);
  t = fma(t, fma(t, -y, 1), t);
  t = fma(t, fma(t, -y, 1), t);

  // Reconstruction

  double u = x * t;
  u = fma(t, fma(u, -y, x), u);

  // Fixup

  if (isnan(u)) u = mulsign(mulsign(INFINITY, x), y);
  if (isinf(y)) u = mulsign(mulsign(0       , x), y);
  if (y == 0 && x == 0) u = NAN;
  if (isnan(x)) u = x;
  if (isnan(y)) u = y;

  //

  return u;
}