Subject: Re: lib/6314: Bug in libc with `double' => `long long' conversions
To: Krister Walfridsson <cato@df.lth.se>
From: Torbjorn Granlund <tege@matematik.su.se>
List: netbsd-bugs
Date: 10/18/1998 04:36:14
Krister,

I think your analysis is quite right.

Would it be OK to assume IEEE in these routines?  I think much could be
gained that way.  Here is my shoot at __fixunsdfdi:

---------------------------------------------------------------------------
union real_extract
{
  double d;
  struct
    {
#if defined (__i386__) || defined (__i860__) || defined (__alpha__)
      unsigned int l:32, h:32;
#else
      unsigned int h:32, l:32;
#endif
    } ii;
  unsigned long long ll;
};

unsigned long long
__fixunsdfdi (double df)
{
  union real_extract re;
  unsigned int mh, ml;
  int e;

  re.d = df;
  e = ((signed int) re.ii.h >> 20) - 1022;      /* exponent and sign bit */
  mh = (re.ii.h & 0xfffff) | 0x100000;  /* upper mantissa + implicit bit */
  ml = re.ii.l;                                        /* lower mantissa */
  if (e < 53)
    {
      int downshift = 53 - e;
      if (downshift < 32)
        {
          ml = (mh << 32 - downshift) | (ml >> downshift);
          mh = mh >> downshift;
        }
      else
        {
          if (e < 0)
            return 0;
          return mh >> downshift - 32;
        }
    }
  else
    {
      int upshift;
      if (e > 64)
        return 0xffffffffffffffffLL;    /* overflow */
      upshift = e - 53;
      mh = (mh << upshift) | (ml >> 32 - upshift);
      ml = ml << upshift;
    }

  re.ii.l = ml;
  re.ii.h = mh;
  return re.ll;
}
---------------------------------------------------------------------------

But why doesn't NetBSD use the code from GCC for all the libc/quad routines?
The GCC routines use a license that allows linking to any commercial
program, as long as GCC is used.  Since GCC is the system compiler, and
since these routines will just be used by GCC anyway, that doesn't appear as
a real restriction.

There are some advantages with using the functions from GCC.  1) They work
correctly.  2) They are efficient (or so I think, since I wrote most of them
:-).

Now, I don't want to get into a political argument.  You may have your
reasons not to use the functions provided with.  Fine with me.

Note also that __anddi3, __iordi3, and __xordi3 are no longer called by GCC.
We fixed that before GCC 2.0 was released, if my memory serves me right.
Also, __lshldi3 is not used.  GCC will call __ashldi3 for all 64-bit left
shifts.  negdi2.c and notdi2 can also be deleted.

I took a close look at /usr/src/lib/libc/quad/fixunsdfdi.c.  This code is
written by somebody with poor understanding of computer architecture.
Here is a replacement.  Note that the bit-fiddling version above is
faster, but the one above is more in line with what you already have.

u_quad_t
__fixunsdfdi(x)
	double x;
{
	union uu t;
	u_long tmp;

	if (x < 0)
		return (UQUAD_MAX);	/* ??? should be 0?  ERANGE??? */
#ifdef notdef				/* this falls afoul of a GCC bug */
	if (x >= UQUAD_MAX)
		return (UQUAD_MAX);
#else					/* so we wire in 2^64-1 instead */
	if (x >= 18446744073709551615.0)	/* XXX */
		return (UQUAD_MAX);
#endif
	/*
	 * Now we know that 0 <= x <= 18446744073709549568.  The upper
	 * limit is one ulp less than 18446744073709551615 tested for above.
	 * Dividing this by 2^32 will *not* round irrespective of any
	 * rounding modes (except if the result is an IEEE denorm).
	 * Furthermore, the quotient will fit into a 32-bit integer.
	 */

	tmp = x / ONE;
	t.ul[L] = (u_long) (x - tmp * ONE);
	t.ul[H] = tmp;
	return (t.uq);
}

I tested this code well.  Here is my test routine.  It uses a proven
test method for this type of routines.

___________________________________________________________________________
#include <stdlib.h>

unsigned long long __fixunsdfdi (double);
unsigned long int random_bitstring ();

main (int argc, char **argv)
{
  double df;
  unsigned long long ll, ll2;
  if (argc == 1)
    {
      do
        {
          /* Make `ll' a 64-bit pattern.  */
          ll = ((long long) random_bitstring () << 32) | random_bitstring ();
          /* Some shifts to make `ll' have at most 53 significands... */
          ll >>= 11;  ll <<= random () % 12;
          df = (double) ll;
          ll2 = __fixunsdfdi (df);
        }
      while (ll2 == ll);
      printf ("%.0f\n%qx %qu\n%qx %qu\n", df, ll, ll, ll2, ll2);
      abort ();
    }
  else
    {
      df = strtod (argv[1], 0);
      ll = __fixunsdfdi (df);
      printf ("%.0f\n%qu\n%qx\n", df, ll, ll);
    }
}

unsigned long int
random_bitstring ()
{
  unsigned long int x;
  int ran, n_bits;
  int tot_bits = 0;

  x = 0;
  for (;;)
    {
      ran = random ();
      n_bits = (ran >> 1) % 16;
      tot_bits += n_bits;

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

          if (tot_bits > 8 * sizeof (long) + 6)
            return x;
        }
    }
}
___________________________________________________________________________

Sorry if this got a bit excessive...

Torbjörn