math: Set errno to ERANGE for logb (+-0) [BZ #6793]

logb (+-0) is a pole error: it returns -Inf and raises the
divide-by-zero exception, but it never set errno, even though glibc
defines math_errhandling to include MATH_ERRNO.

Set errno in the zero branch that already exists in every logb
implementation, instead of adding a w_logb wrapper; the
USE_LOGB*_BUILTIN paths have no such branch, so add one there.  The
double and float versions use __math_divzero and __math_divzerof.
There is no long double equivalent, so those keep the explicit
division and use math_opt_barrier to stop the compiler from folding
it away.

The i386 fxtract implementations of logb and logbf cannot set errno,
and adding the error handling to the assembly is not worthwhile, so
they are removed in favour of the generic C ones.  s_logbl.c moves to
sysdeps/x86/fpu, replacing the x86_64 copy that only included it.

The manual described logb (0) as returning +Inf without signalling,
which was wrong in both respects.

Tested on x86_64-linux-gnu.

Signed-off-by: Shamil Abdulaev <ashamil435@gmail.com>
Reviewed-by: Adhemerval Zanella  <adhemerval.zanella@linaro.org>

manual/math.texi

@@ -669,7 +669,8 @@ due to intermediate rounding.
 If @var{x} is de-normalized, @code{logb} returns the exponent @var{x}
 would have if it were normalized.  If @var{x} is infinity (positive or
 negative), @code{logb} returns @math{@infinity{}}.  If @var{x} is zero,
-@code{logb} returns @math{@infinity{}}.  It does not signal.
+a pole error occurs: @code{logb} raises the divide-by-zero exception,
+sets @code{errno} to @code{ERANGE} and returns @math{-@infinity{}}.
 @end deftypefun
 
 @deftypefun int ilogb (double @var{x})

math/libm-test-logb.inc

@@ -23,10 +23,8 @@ static const struct test_f_f_data logb_test_data[] =
     TEST_f_f (logb, plus_infty, plus_infty, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (logb, minus_infty, plus_infty, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
 
-    /* Bug 6793: errno setting may be missing.  */
+    TEST_f_f (logb, 0, minus_infty, NO_INEXACT_EXCEPTION|DIVIDE_BY_ZERO_EXCEPTION|ERRNO_ERANGE),
-    TEST_f_f (logb, 0, minus_infty, NO_INEXACT_EXCEPTION|DIVIDE_BY_ZERO_EXCEPTION),
+    TEST_f_f (logb, minus_zero, minus_infty, NO_INEXACT_EXCEPTION|DIVIDE_BY_ZERO_EXCEPTION|ERRNO_ERANGE),
-
-    TEST_f_f (logb, minus_zero, minus_infty, NO_INEXACT_EXCEPTION|DIVIDE_BY_ZERO_EXCEPTION),
     TEST_f_f (logb, qnan_value, qnan_value, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (logb, -qnan_value, qnan_value, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED),
     TEST_f_f (logb, snan_value, qnan_value, NO_INEXACT_EXCEPTION|INVALID_EXCEPTION),

sysdeps/i386/fpu/s_logb.S

@@ -1,16 +0,0 @@
-/*
- * Public domain.
- */
-
-#include <machine/asm.h>
-#include <libm-alias-double.h>
-
-RCSID("$NetBSD: s_logb.S,v 1.4 1995/05/09 00:14:30 jtc Exp $")
-
-ENTRY(__logb)
-	fldl	4(%esp)
-	fxtract
-	fstp	%st
-	ret
-END (__logb)
-libm_alias_double (__logb, logb)

sysdeps/i386/fpu/s_logbf.S

@@ -1,16 +0,0 @@
-/*
- * Public domain.
- */
-
-#include <machine/asm.h>
-#include <libm-alias-float.h>
-
-RCSID("$NetBSD: s_logbf.S,v 1.3 1995/05/09 00:15:12 jtc Exp $")
-
-ENTRY(__logbf)
-	flds	4(%esp)
-	fxtract
-	fstp	%st
-	ret
-END (__logbf)
-libm_alias_float (__logb, logb)

sysdeps/ieee754/dbl-64/s_logb.c

@@ -20,11 +20,15 @@
 #include <math_private.h>
 #include <libm-alias-double.h>
 #include <fix-int-fp-convert-zero.h>
+#include "math_config.h"
 
 double
 __logb (double x)
 {
 #if USE_LOGB_BUILTIN
+  if (__glibc_unlikely (x == 0))
+    /* Pole error: logb (+-0).  */
+    return __math_divzero (1);
   return __builtin_logb (x);
 #else
   int64_t ix, ex;
@@ -32,7 +36,8 @@ __logb (double x)
   EXTRACT_WORDS64 (ix, x);
   ix &= UINT64_C(0x7fffffffffffffff);
   if (ix == 0)
-    return -1.0 / fabs (x);
+    /* Pole error: logb (+-0).  */
+    return __math_divzero (1);
   ex = ix >> 52;
   if (ex == 0x7ff)
     return x * x;

sysdeps/ieee754/flt-32/s_logbf.c

@@ -16,11 +16,15 @@
 #include <math_private.h>
 #include <libm-alias-float.h>
 #include <fix-int-fp-convert-zero.h>
+#include "math_config.h"
 
 float
 __logbf (float x)
 {
 #if USE_LOGBF_BUILTIN
+  if (__glibc_unlikely (x == 0))
+    /* Pole error: logbf (+-0).  */
+    return __math_divzerof (1);
   return __builtin_logbf (x);
 #else
   int32_t ix, rix;
@@ -28,7 +32,8 @@ __logbf (float x)
   GET_FLOAT_WORD (ix, x);
   ix &= 0x7fffffff;		/* high |x| */
   if (ix == 0)
-    return (float) -1.0 / fabsf (x);
+    /* Pole error: logbf (+-0).  */
+    return __math_divzerof (1);
   if (ix >= 0x7f800000)
     return x * x;
   if (__glibc_unlikely ((rix = ix >> 23) == 0))

sysdeps/ieee754/ldbl-128/s_logbl.c

@@ -23,6 +23,8 @@ static char rcsid[] = "$NetBSD: $";
  */
 
 #include <math.h>
+#include <errno.h>
+#include <math-barriers.h>
 #include <math_private.h>
 #include <libm-alias-ldouble.h>
 
@@ -30,6 +32,12 @@ _Float128
 __logbl (_Float128 x)
 {
 #if USE_LOGBL_BUILTIN
+  if (__glibc_unlikely (x == 0))
+    {
+      /* Pole error: logbl (+-0).  */
+      __set_errno (ERANGE);
+      return math_opt_barrier (-1.0) / fabsl (x);
+    }
   return __builtin_logbl (x);
 #else
   /* Use generic implementation.  */
@@ -38,7 +46,11 @@ __logbl (_Float128 x)
   GET_LDOUBLE_WORDS64 (hx, lx, x);
   hx &= 0x7fffffffffffffffLL;	/* high |x| */
   if ((hx | lx) == 0)
-    return -1.0 / fabsl (x);
+    {
+      /* Pole error: logbl (+-0).  */
+      __set_errno (ERANGE);
+      return math_opt_barrier (-1.0) / fabsl (x);
+    }
   if (hx >= 0x7fff000000000000LL)
     return x * x;
   if ((ex = hx >> 48) == 0)	/* IEEE 754 logb */

sysdeps/ieee754/ldbl-128ibm/s_logbl.c

@@ -19,6 +19,8 @@
  */
 
 #include <math.h>
+#include <errno.h>
+#include <math-barriers.h>
 #include <math_private.h>
 #include <math_ldbl_opt.h>
 #include <fix-int-fp-convert-zero.h>
@@ -34,7 +36,11 @@ __logbl (long double x)
   hxs = hx;
   hx &= 0x7fffffffffffffffLL;	/* high |x| */
   if (hx == 0)
-    return -1.0 / fabs (x);
+    {
+      /* Pole error: logbl (+-0).  */
+      __set_errno (ERANGE);
+      return math_opt_barrier (-1.0) / fabs (x);
+    }
   if (hx >= 0x7ff0000000000000LL)
     return x * x;
   if (__glibc_unlikely ((rhx = hx >> 52) == 0))

sysdeps/m68k/m680x0/fpu/s_logbl.c

@@ -19,6 +19,8 @@
  */
 
 #include <math.h>
+#include <errno.h>
+#include <math-barriers.h>
 #include <math_private.h>
 
 long double
@@ -29,7 +31,11 @@ __logbl (long double x)
   GET_LDOUBLE_WORDS (es, ix, lx, x);
   es &= 0x7fff;			/* exponent */
   if ((es | ix | lx) == 0)
-    return -1.0 / fabsl (x);
+    {
+      /* Pole error: logbl (+-0).  */
+      __set_errno (ERANGE);
+      return math_opt_barrier (-1.0) / fabsl (x);
+    }
   if (es == 0x7fff)
     return x * x;
   if (es == 0)			/* IEEE 754 logb */

sysdeps/powerpc/fpu/s_logb.c

@@ -22,6 +22,8 @@
 # include <sysdeps/ieee754/dbl-64/s_logb.c>
 #else
 # include <math.h>
+# include <errno.h>
+# include <math-barriers.h>
 # include <math_private.h>
 # include <math_ldbl_opt.h>
 # include <libm-alias-double.h>
@@ -34,8 +36,11 @@ __logb (double x)
   double ret;
 
   if (__glibc_unlikely (x == 0.0))
-    /* Raise FE_DIVBYZERO and return -HUGE_VAL[LF].  */
+    {
-    return -1.0 / fabs (x);
+      /* Pole error: raise FE_DIVBYZERO, set errno and return -HUGE_VAL.  */
+      __set_errno (ERANGE);
+      return math_opt_barrier (-1.0) / fabs (x);
+    }
 
   /* Mask to extract the exponent.  */
   asm ("xxland %x0,%x1,%x2\n"

sysdeps/powerpc/fpu/s_logbf.c

@@ -22,6 +22,8 @@
 # include <sysdeps/ieee754/flt-32/s_logbf.c>
 #else
 # include <math.h>
+# include <errno.h>
+# include <math-barriers.h>
 # include <libm-alias-float.h>
 /* This implementation avoids FP to INT conversions by using VSX
    bitwise instructions over FP values.  */
@@ -32,8 +34,11 @@ __logbf (float x)
   double ret;
 
   if (__glibc_unlikely (x == 0.0))
-    /* Raise FE_DIVBYZERO and return -HUGE_VAL[LF].  */
+    {
-    return -1.0 / fabs (x);
+      /* Pole error: raise FE_DIVBYZERO, set errno and return -HUGE_VALF.  */
+      __set_errno (ERANGE);
+      return math_opt_barrier (-1.0) / fabs (x);
+    }
 
   /* mask to extract the exponent.  */
   asm ("xxland %x0,%x1,%x2\n"

sysdeps/powerpc/fpu/s_logbl.c

@@ -22,6 +22,8 @@
 # include <./sysdeps/ieee754/ldbl-128ibm/s_logbl.c>
 #else
 # include <math.h>
+# include <errno.h>
+# include <math-barriers.h>
 # include <math_private.h>
 # include <math_ldbl_opt.h>
 
@@ -35,8 +37,11 @@ __logbl (long double x)
   int64_t hx;
 
   if (__glibc_unlikely (x == 0.0))
-    /* Raise FE_DIVBYZERO and return -HUGE_VAL[LF].  */
+    {
-    return -1.0L / __builtin_fabsl (x);
+      /* Pole error: raise FE_DIVBYZERO, set errno and return -HUGE_VALL.  */
+      __set_errno (ERANGE);
+      return math_opt_barrier (-1.0L) / __builtin_fabsl (x);
+    }
 
   ldbl_unpack (x, &xh, &xl);
   EXTRACT_WORDS64 (hx, xh);

sysdeps/x86/fpu/s_logbl.c

@@ -2,6 +2,9 @@
  * Public domain.
  */
 
+#include <math.h>
+#include <errno.h>
+#include <math-barriers.h>
 #include <libm-alias-ldouble.h>
 
 long double
@@ -9,6 +12,13 @@ __logbl (long double x)
 {
   long double res;
 
+  if (__glibc_unlikely (x == 0))
+    {
+      /* Pole error: logbl (+-0).  */
+      __set_errno (ERANGE);
+      return math_opt_barrier (-1.0L) / 0.0L;
+    }
+
   asm ("fxtract\n"
        "fstp	%%st" : "=t" (res) : "0" (x));
   return res;

sysdeps/x86_64/fpu/s_logbl.c

@@ -1 +0,0 @@
-#include <sysdeps/i386/fpu/s_logbl.c>