#define _GNU_SOURCE
#include <stdlib.h>
#include <mpfr.h>
#include <ieee754.h>
#include <stdio.h>
#include <math.h>

#ifdef DO_FLT
#define MANT_DIG __FLT_MANT_DIG__
#define MAX_EXP __FLT_MAX_EXP__
#define union_type union ieee754_float
#define suff(x) x##d
#define fma_suff fmaf
#define BIAS 0x7f
#define FMT "%a"
#define fltfield f
#elif defined DO_DBL
#define MANT_DIG __DBL_MANT_DIG__
#define MAX_EXP __DBL_MAX_EXP__
#define union_type union ieee754_double
#define suff(x) x##d
#define fma_suff fma
#define BIAS 0x3ff
#define FMT "%a"
#define fltfield d
#elif defined DO_LDBL
#define MANT_DIG __LDBL_MANT_DIG__
#define MAX_EXP __LDBL_MAX_EXP__
#define union_type union ieee854_long_double
#define suff(x) x##ld
#define fma_suff fmal
#define BIAS 0x3fff
#if MANT_DIG == 64
#define FMT "%La"
#else
#define FMT "%.40La"
#endif
#define fltfield d
#else
#error Define DO_FLT, DO_DBL or DO_LDBL
#endif

unsigned int
randbits (int count)
{
  static unsigned int val, bits;
  unsigned int ret = 0;

  if (count > bits)
    {
      if (count == 32 && bits == 0)
	{
	  ret = ((unsigned int) random ()) << 1;
	  count -= 31;
	}
      else
	{
	  ret = val << (count - bits);
	  count -= bits;
	}
      val = (random () & 0x7fffffffU);
      bits = 31;
    }
  ret |= val & ((1U << count) - 1);
  val >>= count;
  bits -= count;
  return ret;
}

void
randflt (int mode, union_type *p, unsigned int exp)
{
  if (mode == 1 || (mode >= 2 && randbits (4) == 0))
    {
      unsigned int mantissa[4] = { 0, 0, 0, 0 };
      unsigned int bits = 0, cur = randbits (1);
      do
	{
	  unsigned int curbits = randbits (8) % (MANT_DIG - bits + 6);
	  if (curbits > MANT_DIG - bits)
	    curbits = MANT_DIG - bits;
	  if (cur && curbits)
	    {
	      do
		{
		  unsigned int thisbits = 32 - (bits & 31);
		  if (curbits < thisbits)
		    thisbits = curbits;
		  mantissa[bits / 32] |= (~0U << (32 - thisbits))
					 >> (32 - thisbits - (bits & 31));
		  bits += thisbits;
		  curbits -= thisbits;
		}
	      while (curbits);
	    }
	  else
	    bits += curbits;
	  cur ^= 1;
	}
      while (bits < MANT_DIG);
#ifdef DO_FLT
      p->ieee.mantissa = mantissa[0];
#elif defined DO_DBL
      p->ieee.mantissa1 = mantissa[0];
      p->ieee.mantissa0 = mantissa[1];
#elif defined DO_LDBL
#if MANT_DIG == 64
      p->ieee.mantissa1 = mantissa[0];
      p->ieee.mantissa0 = mantissa[1] & 0x7fffffffU;
#elif MANT_DIG == 113
      p->ieee.mantissa3 = mantissa[0];
      p->ieee.mantissa2 = mantissa[1];
      p->ieee.mantissa1 = mantissa[2];
      p->ieee.mantissa0 = mantissa[3];
#endif
#endif
    }
  else
    {
#ifdef DO_FLT
      p->ieee.mantissa = randbits (23);
#elif defined DO_DBL
      p->ieee.mantissa1 = randbits (32);
      p->ieee.mantissa0 = randbits (20);
#elif defined DO_LDBL
#if MANT_DIG == 64
      p->ieee.mantissa1 = randbits (32);
      p->ieee.mantissa0 = randbits (31);
#elif MANT_DIG == 113
      p->ieee.mantissa3 = randbits (32);
      p->ieee.mantissa2 = randbits (32);
      p->ieee.mantissa1 = randbits (32);
      p->ieee.mantissa0 = randbits (16);
#endif
#endif
    }
  p->ieee.negative = randbits (1);
  p->ieee.exponent = exp;
#if defined (DO_LDBL) && MANT_DIG == 64
  p->ieee.empty = 0;
  if (exp != 0)
    p->ieee.mantissa0 |= 0x80000000;
#endif
}

int
main (int argc, const char **argv)
{
  long int i, count = 100;

  if (argc >= 2)
    count = strtoul (argv[1], NULL, 0);
  mpfr_set_default_prec (MANT_DIG);
  mpfr_set_emin (-MAX_EXP-MANT_DIG+4);
  mpfr_set_emax (MAX_EXP);
  for (i = 0; i < count; i++)
    {
      union_type u, v, w, x, y;
      int mode = randbits (2);
      randflt (mode, &u, randbits (15));
      randflt (mode, &v, randbits (15));
      if (randbits (4))
	randflt (mode, &w, ((int) u.ieee.exponent + (int) v.ieee.exponent)
			   - BIAS - MANT_DIG + randbits (18) % (2 * MANT_DIG));
      else
	randflt (mode, &w, randbits (15));
      mpfr_t xx, xy, xz;
      suff (mpfr_init_set_) (xx, u.fltfield, GMP_RNDN);
      suff (mpfr_init_set_) (xy, v.fltfield, GMP_RNDN);
      suff (mpfr_init_set_) (xz, w.fltfield, GMP_RNDN);
      int i = mpfr_fma (xx, xx, xy, xz, GMP_RNDN);
      mpfr_subnormalize (xx, i, GMP_RNDN);
      x.fltfield = suff (mpfr_get_) (xx, GMP_RNDN);
      mpfr_clear (xx); mpfr_clear (xy); mpfr_clear (xz);
      y.fltfield = fma_suff (u.fltfield, v.fltfield, w.fltfield);
      if (x.fltfield != y.fltfield
	  && !(isnan (x.fltfield) && isnan (y.fltfield)))
	printf ("fma ("FMT", "FMT", "FMT") = "FMT" mpfr_fma "FMT"\n",
		u.fltfield, v.fltfield, w.fltfield, y.fltfield, x.fltfield);
    }
  return 0;
}