[PATCH] Fix ppc{32,64} lrintl/lroundl

[Jakub Jelinek <jakub@redhat.com>]

Hi!

lrintl and lroundl are broken on both ppc32 and ppc64.
For ppc64 the problem is that GCC optimizes llroundl call to lroundl
and llrintl to lrintl when long int is as wide as long long int,
but lrintl resp. lroundl are on ppc implemented as calls to
llrintl resp. llroundl.  That means endless recursion.
For ppc32 the problem was that for x such that llrintl (x) != (long) llrintl (x)
(and similarly l*roundl) FE_INVALID wasn't raised and for values where
FE_INVALID was raised even for llrintl/llroundl lrintl/lroundl would
return unexpected return values (as (long) __LONG_LONG_MAX__ is not
__LONG_MAX__ and ((long) (-__LONG_LONG_MAX__ - 1)) is not
(-__LONG_MAX__ - 1).

2006-09-14  Jakub Jelinek  <jakub@redhat.com>

	* sysdeps/ieee754/ldbl-128ibm/s_llrintl.c: Comment fixes.
	* sysdeps/ieee754/ldbl-128ibm/s_llroundl.c: Comment fixes.
	* sysdeps/ieee754/ldbl-128ibm/s_lrintl.c: Rewritten.
	* sysdeps/ieee754/ldbl-128ibm/s_lroundl.c: Rewritten.

--- libc/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c.jj	2006-03-22 23:16:58.000000000 +0100
+++ libc/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c	2006-09-14 12:29:37.000000000 +0200
@@ -1,4 +1,4 @@
-/* Round to int long double floating-point values.
+/* Round to long long int long double floating-point values.
    IBM extended format long double version.
    Copyright (C) 2006 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
@@ -18,9 +18,6 @@
    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
    02111-1307 USA.  */
 
-/* This has been coded in assembler because GCC makes such a mess of it
-   when it's coded in C.  */
-
 #include <math.h>
 #include <fenv_libc.h>
 #include <math_ldbl_opt.h>
--- libc/sysdeps/ieee754/ldbl-128ibm/s_llroundl.c.jj	2006-03-22 23:16:58.000000000 +0100
+++ libc/sysdeps/ieee754/ldbl-128ibm/s_llroundl.c	2006-09-14 12:29:37.000000000 +0200
@@ -1,4 +1,4 @@
-/* Round to long long, long double floating-point values.
+/* Round to long long int long double floating-point values.
    IBM extended format long double version.
    Copyright (C) 2006 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
@@ -18,9 +18,6 @@
    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
    02111-1307 USA.  */
 
-/* This has been coded in assembler because GCC makes such a mess of it
-   when it's coded in C.  */
-
 #include <math.h>
 #include <fenv_libc.h>
 #include <math_ldbl_opt.h>
--- libc/sysdeps/ieee754/ldbl-128ibm/s_lrintl.c.jj	2006-01-28 01:07:25.000000000 +0100
+++ libc/sysdeps/ieee754/ldbl-128ibm/s_lrintl.c	2006-09-14 12:29:37.000000000 +0200
@@ -1,9 +1,162 @@
-/* FIXME */
+/* Round to long int long double floating-point values.
+   IBM extended format long double version.
+   Copyright (C) 2006 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
 #include <math.h>
+#include <fenv_libc.h>
 #include <math_ldbl_opt.h>
+#include <float.h>
+#include <ieee754.h>
+
 
-long int __lrintl (long double d)
+#ifdef __STDC__
+long
+__lrintl (long double x)
+#else
+long
+__lrintl (x)
+     long double x;
+#endif
 {
-  return llrintl (d);
+  double xh, xl;
+  long res, hi, lo;
+  int save_round;
+
+  ldbl_unpack (x, &xh, &xl);
+
+  /* Limit the range of values handled by the conversion to long.
+     We do this because we aren't sure whether that conversion properly
+     raises FE_INVALID.  */
+  if (
+#if __LONG_MAX__ == 2147483647
+      __builtin_expect
+      ((__builtin_fabs (xh) <= (double) __LONG_MAX__ + 2), 1)
+#else
+      __builtin_expect
+      ((__builtin_fabs (xh) <= -(double) (-__LONG_MAX__ - 1)), 1)
+#endif
+#if !defined (FE_INVALID)
+      || 1
+#endif
+    )
+    {
+      save_round = fegetround ();
+
+#if __LONG_MAX__ == 2147483647
+      long long llhi = (long long) xh;
+      if (llhi != (long) llhi)
+	hi = llhi < 0 ? -__LONG_MAX__ - 1 : __LONG_MAX__;
+      else
+	hi = llhi;
+      xh -= hi;
+#else
+      if (__builtin_expect ((xh == -(double) (-__LONG_MAX__ - 1)), 0))
+	{
+	  /* When XH is 9223372036854775808.0, converting to long long will
+	     overflow, resulting in an invalid operation.  However, XL might
+	     be negative and of sufficient magnitude that the overall long
+	     double is in fact in range.  Avoid raising an exception.  In any
+	     case we need to convert this value specially, because
+	     the converted value is not exactly represented as a double
+	     thus subtracting HI from XH suffers rounding error.  */
+	  hi = __LONG_MAX__;
+	  xh = 1.0;
+	}
+      else
+	{
+	  hi = (long) xh;
+	  xh -= hi;
+	}
+#endif
+      ldbl_canonicalize (&xh, &xl);
+
+      lo = (long) xh;
+
+      /* Peg at max/min values, assuming that the above conversions do so.
+         Strictly speaking, we can return anything for values that overflow,
+         but this is more useful.  */
+      res = hi + lo;
+
+      /* This is just sign(hi) == sign(lo) && sign(res) != sign(hi).  */
+      if (__builtin_expect (((~(hi ^ lo) & (res ^ hi)) < 0), 0))
+	goto overflow;
+
+      xh -= lo;
+      ldbl_canonicalize (&xh, &xl);
+
+      hi = res;
+      switch (save_round)
+	{
+	case FE_TONEAREST:
+	  if (fabs (xh) < 0.5
+	      || (fabs (xh) == 0.5
+		  && ((xh > 0.0 && xl < 0.0)
+		      || (xh < 0.0 && xl > 0.0)
+		      || (xl == 0.0 && (res & 1) == 0))))
+	    return res;
+
+	  if (xh < 0.0)
+	    res -= 1;
+	  else
+	    res += 1;
+	  break;
+
+	case FE_TOWARDZERO:
+	  if (res > 0 && (xh < 0.0 || (xh == 0.0 && xl < 0.0)))
+	    res -= 1;
+	  else if (res < 0 && (xh > 0.0 || (xh == 0.0 && xl > 0.0)))
+	    res += 1;
+	  return res;
+	  break;
+
+	case FE_UPWARD:
+	  if (xh > 0.0 || (xh == 0.0 && xl > 0.0))
+	    res += 1;
+	  break;
+
+	case FE_DOWNWARD:
+	  if (xh < 0.0 || (xh == 0.0 && xl < 0.0))
+	    res -= 1;
+	  break;
+	}
+
+      if (__builtin_expect (((~(hi ^ (res - hi)) & (res ^ hi)) < 0), 0))
+	goto overflow;
+
+      return res;
+    }
+  else
+    {
+      if (xh > 0.0)
+	hi = __LONG_MAX__;
+      else if (xh < 0.0)
+	hi = -__LONG_MAX__ - 1;
+      else
+	/* Nan */
+	hi = 0;
+    }
+
+overflow:
+#ifdef FE_INVALID
+  feraiseexcept (FE_INVALID);
+#endif
+  return hi;
 }
+
 long_double_symbol (libm, __lrintl, lrintl);
--- libc/sysdeps/ieee754/ldbl-128ibm/s_lroundl.c.jj	2006-01-28 01:07:25.000000000 +0100
+++ libc/sysdeps/ieee754/ldbl-128ibm/s_lroundl.c	2006-09-14 12:29:37.000000000 +0200
@@ -1,9 +1,142 @@
-/* FIXME */
+/* Round to long int long double floating-point values.
+   IBM extended format long double version.
+   Copyright (C) 2006 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
 #include <math.h>
+#include <fenv_libc.h>
 #include <math_ldbl_opt.h>
+#include <float.h>
+#include <ieee754.h>
 
-long int __lroundl (long double d)
+#ifdef __STDC__
+long
+__lroundl (long double x)
+#else
+long
+__lroundl (x)
+     long double x;
+#endif
 {
-  return llroundl (d);
+  double xh, xl;
+  long res, hi, lo;
+
+  ldbl_unpack (x, &xh, &xl);
+
+  /* Limit the range of values handled by the conversion to long.
+     We do this because we aren't sure whether that conversion properly
+     raises FE_INVALID.  */
+  if (
+#if __LONG_MAX__ == 2147483647
+      __builtin_expect
+      ((__builtin_fabs (xh) <= (double) __LONG_MAX__ + 2), 1)
+#else
+      __builtin_expect
+      ((__builtin_fabs (xh) <= -(double) (-__LONG_MAX__ - 1)), 1)
+#endif
+#if !defined (FE_INVALID)
+      || 1
+#endif
+    )
+    {
+#if __LONG_MAX__ == 2147483647
+      long long llhi = (long long) xh;
+      if (llhi != (long) llhi)
+	hi = llhi < 0 ? -__LONG_MAX__ - 1 : __LONG_MAX__;
+      else
+	hi = llhi;
+      xh -= hi;
+#else
+      if (__builtin_expect ((xh == -(double) (-__LONG_MAX__ - 1)), 0))
+	{
+	  /* When XH is 9223372036854775808.0, converting to long long will
+	     overflow, resulting in an invalid operation.  However, XL might
+	     be negative and of sufficient magnitude that the overall long
+	     double is in fact in range.  Avoid raising an exception.  In any
+	     case we need to convert this value specially, because
+	     the converted value is not exactly represented as a double
+	     thus subtracting HI from XH suffers rounding error.  */
+	  hi = __LONG_MAX__;
+	  xh = 1.0;
+	}
+      else
+	{
+	  hi = (long) xh;
+	  xh -= hi;
+	}
+#endif
+      ldbl_canonicalize (&xh, &xl);
+
+      lo = (long) xh;
+
+      /* Peg at max/min values, assuming that the above conversions do so.
+         Strictly speaking, we can return anything for values that overflow,
+         but this is more useful.  */
+      res = hi + lo;
+
+      /* This is just sign(hi) == sign(lo) && sign(res) != sign(hi).  */
+      if (__builtin_expect (((~(hi ^ lo) & (res ^ hi)) < 0), 0))
+	goto overflow;
+
+      xh -= lo;
+      ldbl_canonicalize (&xh, &xl);
+
+      hi = res;
+      if (xh > 0.5)
+	{
+	  res += 1;
+	}
+      else if (xh == 0.5)
+	{
+	  if (xl > 0.0 || (xl == 0.0 && res >= 0))
+	    res += 1;
+	}
+      else if (-xh > 0.5)
+	{
+	  res -= 1;
+	}
+      else if (-xh == 0.5)
+	{
+	  if (xl < 0.0 || (xl == 0.0 && res <= 0))
+	    res -= 1;
+	}
+
+      if (__builtin_expect (((~(hi ^ (res - hi)) & (res ^ hi)) < 0), 0))
+	goto overflow;
+
+      return res;
+    }
+  else
+    {
+      if (xh > 0.0)
+	hi = __LONG_MAX__;
+      else if (xh < 0.0)
+	hi = -__LONG_MAX__ - 1;
+      else
+	/* Nan */
+	hi = 0;
+    }
+
+overflow:
+#ifdef FE_INVALID
+  feraiseexcept (FE_INVALID);
+#endif
+  return hi;
 }
+
 long_double_symbol (libm, __lroundl, lroundl);

	Jakub