RFH: testers for ieee test

[Joern Rennecke]

I'm trying to write a test to check I get all the subnormal
corner cases for ieee-754 double precision floating point division right.
The usual thing to do with a new test is to try it on a known good
platform.  Unfortunately, the x86 processor of my workstation is actually
known to get these calculations wrong.  Could someone with a fully
ieee compliant FPU try this test?

I've currently made it output lots of debug information and to continue
in the faceof errors; to get only a single result via the exit code,
just disable the DEBUG define.

It appears my machine gets about one in 3.7 of the fraction approximation
tests wrong :-(

FWIW, the random stuff is cribbed from arith-rand-ll.c.
-------------- next part --------------
/* Test proper rounding for division when generating subnormal numbers.  */

/* The debug code is only meant to run on linux with long double support.  */
#define DEBUG

long long
simple_rand ()
{
  static unsigned long long seed = 47114711;
  unsigned long long this = seed * 1103515245 + 12345;
  seed = this;
  return this >> 8;
}

unsigned long long int
random_bitstring (int *lenp)
{
  unsigned long long int x;
  int n_bits;
  int tot_bits = 0;
  int len = 0;
  int limit = *lenp;

  x = 0;
  for (;;)
    {
      long long ran = simple_rand ();
      int ones_p = ran & 1;
      n_bits = (ran >> 1) % 16;
      if (x || ones_p)
	{
	  if (n_bits > limit - len)
	    n_bits = limit - len;
	  *lenp = len += n_bits;
	}
      tot_bits += n_bits;

      if (n_bits == 0 && len)
	return x;
      else if (n_bits)
	{
	  x <<= n_bits;
	  if (ones_p)
	    x |= (1 << n_bits) - 1;

	  if (len == limit
	      || (len && tot_bits > 8 * sizeof (long long) + 6))
	    return x;
	}
    }
}

main ()
{
  int i;
  double e0, e1, ie1, i2e1;
  /* Hack to get x86 math to work.  */
  volatile double eval_tmp;
#define EVAL(x) (eval_tmp = (x))

  for (i = 0, e1 = 0.5; EVAL (1. + e1) > 1.; e1 *= 0.5, i++);
  e1 *= 2.;
  ie1 = 1./e1;
  i2e1 = 2./e1;
  if (i < 52 || (sizeof (double) == 8 && i > 52))
    abort ();
  for (i = 1, e0 = 0.5; EVAL (e0 * e0); e0 *= e0, i <<= 1);
  for (;EVAL (e0 * 0.5); e0 *= 0.5, i++);
  if (i < 0x3ff + 51 || (sizeof (double) == 8 && i > 0x3ff + 51))
    abort ();
  
  /* First, check that a quotient that can be computed exactly is properly
     rounded, and try variantions on the fraction to do some simple
     round-to-nearest checks for inexact results.  */
  for (i = 0; i < 10000; i++)
    {
      unsigned long long x, y;
      int xlen, ylen;
      long long ran;
      double xd, xr, yd, yd2, e2y, pd, ep;

      xlen = 53;
      x = random_bitstring (&xlen);
      ylen = 54-xlen;
      y = random_bitstring (&ylen);
      y |= 1;
      if (x > 2 &&
	  (double) (x|3) * y >= (double) (1LL << 53))
	{
	  x >>= 1;
	  xlen--;
	}
      x |= 1;
      ran = simple_rand ();
      x ^= ran & 2;
      xd = (double)x * e0;
      yd = (double)y * ie1;
      yd2 = yd * 2.;
      e2y = (double) (1LL << ylen);
      if (EVAL (yd2 + e2y) == yd2)
	abort ();
      if ((yd2 + e2y) / yd2 > (1.+e1)/1.)
	abort ();
      pd = xd * yd;
      ep = e0 * (1LL << xlen-1) * (1LL << ylen-1);
      if (EVAL (pd + ep) == pd)
	ep += ep;
      if (pd + ep == pd)
	abort ();
      if (EVAL((pd + ep) / pd) > EVAL (1 + e1))
	abort ();
      if (EVAL (pd / yd) != x * e0)
	abort ();
      /* Round to even.  */
      xr = ((x & 2) + x) >> 1;
      if (EVAL (pd / yd2) != xr * e0)
	abort ();
      /* Round to nearest - upwards.  */
      xr = x+1 >> 1;
      if (EVAL (pd / (yd2-e2y)) < xr * e0)
	abort ();
      if (EVAL ((pd + ep) / yd2) < xr * e0)
	abort ();
      /* Round to nearest - downwards.  */
      xr = x >> 1;
      if (EVAL (pd / (yd2+e2y)) > xr * e0)
	abort ();
      if (EVAL ((pd - ep) / yd2) > xr * e0)
	abort ();
    }
  /* Now generate a set of 53 bit random numbers, calculate a fractional
     approximation which is likely to be hard to distinguish from the
     exact result, and check for proper rounding.  */
  for (i = 0; i < 10000; i++)
    {
      unsigned long long x, y, x0, y0, x1, y1, x2, y2, x3, tmp;
      int rest_sign;
      long long ran;
      int xlen;
      long long a[10];
      int j, k;

      do
	{
	  xlen = 53;
	  x = random_bitstring (&xlen);
	}
      while (xlen < 10);
      x |= 1;
      x0 = x;
      /* Look for a close, but inexact approximation that fits in 53 bits
	 numerator / denominator.  */
      y = 1LL << xlen - 1;
      y0 = y;
      rest_sign = 0;
      for (j = 0; j < 20; j++)
	{
	  a[j] = x/y;
	  x1 = a[j], y1 = 1;
	  for (k = j - 1; k >= 0; k--)
	    {
	      tmp = a[k] * x1 + y1;
	      if (tmp >= 1LL << 53 || tmp/a[k] < x1)
		goto end_approx;
	      y1 = x1;
	      x1 = tmp;
	    }
	  tmp = x - a[j] * y;
	  if (!tmp)
	    break;
	  rest_sign = j & 1 ? -1 : 1;
	  x2 = x1, y2 = y1;
	  x = y;
	  y = tmp;
	}
      end_approx:
      if (!rest_sign)
	continue;
#ifdef DEBUG
      printf ("%d %d %f %f %e\n", j, rest_sign,
	      (double)x0/y0, (double)x2/y2,
	      (double)((long double)x0/y0-((long double)x2/y2)));
#endif
      x3 = EVAL (x2*e0*y0/(2.*y2)) / e0 * 2.;
#ifdef  DEBUG
      printf ("%llx %llx\n", x0, x3);
#endif
      if (rest_sign > 0 ? x3 >= x0 : x3 <= x0)
#ifdef DEBUG
	printf("ERROR!\n");
#else
	abort ();
#endif
    }

  exit (0);
}