diff mbox

[BZ,#18875] Avoid excess precision in check against float zero

Message ID 20150826133745.GA17104@gmail.com
State Dropped
Headers show

Commit Message

H.J. Lu Aug. 26, 2015, 1:37 p.m. UTC
On Tue, Jul 07, 2015 at 03:55:16PM -0400, Carlos O'Donell wrote:
> Joseph,
> 
> I'm seeing these failures on i686 in test-float, have you seen these?
> If not, I'll investigate the failures to see if they are related to
> my compiler or linker.
> 
> glibc-i686/math/test-float.out:
> 
> testing float (without inline functions)
> Failure: exp10 (-0x1.31p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp10 (-0x1.343792p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp10 (-0x1.343794p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp10 (-0x1.344p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp10_downward (-0x1.31p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp10_downward (-0x1.343792p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp10_downward (-0x1.343794p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp10_downward (-0x1.344p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp10_towardzero (-0x1.31p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp10_towardzero (-0x1.343792p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp10_towardzero (-0x1.343794p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp10_towardzero (-0x1.344p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp2 (-0x1.3045fep+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2 (-0x3.fb8p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2 (-0x3.fc8p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2 (-0x3.fd8p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2 (-0x4.01p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp2 (-0x4.32p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_downward (-0x1.3045fep+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_downward (-0x3.fb8p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_downward (-0x3.fc8p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_downward (-0x3.fd8p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_downward (-0x4.01p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_downward (-0x4.32p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_towardzero (-0x1.3045fep+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_towardzero (-0x3.fb8p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_towardzero (-0x3.fc8p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_towardzero (-0x3.fd8p+8): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_towardzero (-0x4.01p+12): errno set to 0, expected 34 (ERANGE)
> Failure: exp2_towardzero (-0x4.32p+8): errno set to 0, expected 34 (ERANGE)
> Failure: pow10 (-0x1.31p+8): errno set to 0, expected 34 (ERANGE)
> Failure: pow10 (-0x1.343792p+12): errno set to 0, expected 34 (ERANGE)
> Failure: pow10 (-0x1.343794p+12): errno set to 0, expected 34 (ERANGE)
> Failure: pow10 (-0x1.344p+12): errno set to 0, expected 34 (ERANGE)
> Failure: pow10_downward (-0x1.31p+8): errno set to 0, expected 34 (ERANGE)
> Failure: pow10_downward (-0x1.343792p+12): errno set to 0, expected 34 (ERANGE)
> Failure: pow10_downward (-0x1.343794p+12): errno set to 0, expected 34 (ERANGE)
> Failure: pow10_downward (-0x1.344p+12): errno set to 0, expected 34 (ERANGE)
> Failure: pow10_towardzero (-0x1.31p+8): errno set to 0, expected 34 (ERANGE)
> Failure: pow10_towardzero (-0x1.343792p+12): errno set to 0, expected 34 (ERANGE)
> Failure: pow10_towardzero (-0x1.343794p+12): errno set to 0, expected 34 (ERANGE)
> Failure: pow10_towardzero (-0x1.344p+12): errno set to 0, expected 34 (ERANGE)
> 
> Test suite completed:
>   38654 test cases plus 36610 tests for exception flags and
>     36610 tests for errno executed.
>   42 errors occurred.
> 

For check against float zero:

  float z;
  z = ...
  if (z == 0)

on x86 with excess precision, we may get

(gdb) info float
=>R7: Valid   0x3d30e0d9b4d78a806800 +6.369668676344082156e-217
  R6: Empty   0x40029000000000000000

and z (+6.369668676344082156e-217) == 0 is false.  We can avoid excess
precision by checking binary representations of float zero directly.

OK for master?

Thanks.

H.J.
---
	[BZ #18875]
	* sysdeps/i386/fpu/w_exp10f.c: New file.
	* sysdeps/i386/fpu/w_exp2f.c: Likwise.
	* sysdeps/i386/fpu/w_expf.c: Likwise.
---
 sysdeps/i386/fpu/w_exp10f.c | 41 +++++++++++++++++++++++++++++++++++++++++
 sysdeps/i386/fpu/w_exp2f.c  | 38 ++++++++++++++++++++++++++++++++++++++
 sysdeps/i386/fpu/w_expf.c   | 37 +++++++++++++++++++++++++++++++++++++
 3 files changed, 116 insertions(+)
 create mode 100644 sysdeps/i386/fpu/w_exp10f.c
 create mode 100644 sysdeps/i386/fpu/w_exp2f.c
 create mode 100644 sysdeps/i386/fpu/w_expf.c

Comments

Joseph Myers Aug. 26, 2015, 2:24 p.m. UTC | #1
On Wed, 26 Aug 2015, H.J. Lu wrote:

> 	[BZ #18875]
> 	* sysdeps/i386/fpu/w_exp10f.c: New file.
> 	* sysdeps/i386/fpu/w_exp2f.c: Likwise.
> 	* sysdeps/i386/fpu/w_expf.c: Likwise.

I don't really like having extra copies of the wrappers, which will mean 
more places to clean up when we eventually obsolete _LIB_VERSION and 
__kernel_standard and make the wrappers set errno directly.

The problem is: the functions return results with excess range and 
precision.  This seems a dubious thing to do anyway in C standard terms 
(not allowed for return statements under Annex F, but there's no actual 
statement that library functions return as if by a C return statement).  
Your fix works around this for the wrappers, so as to set errno and avoid 
test failures.  But if the underflow in the wrappers isn't all the way to 
zero, the wrappers might well return the original value with excess range 
and precision, and confuse any callers in much the same way the wrappers 
were confused.

I think it's better to do what I suggested in the bug and copy the code 
from atan2f to eliminate excess range and precision, and force underflow 
exceptions, for underflowing results, before the functions return to the 
wrappers at all.  You're still storing to the stack in order to get at a 
value without excess range and precision, just before returning to the 
wrapper rather than in the wrapper.  The same thing should be done for the 
double versions of these functions, which could have the issue just as 
much as the float versions depending on what code the compiler chooses to 
generate.  That way you avoid adding any architecture-specific versions of 
files where we don't already have them, just patch existing 
architecture-specific files.

(In either case, the general question of avoiding excess range and 
precision in return values from all libm functions is deferred; the patch 
would just deal with float and double versions of these three functions, 
for underflowing results.  The general issue is already mentioned at 
<https://sourceware.org/glibc/wiki/Development_Todo/Master#libm_itself>.)
diff mbox

Patch

diff --git a/sysdeps/i386/fpu/w_exp10f.c b/sysdeps/i386/fpu/w_exp10f.c
new file mode 100644
index 0000000..bee27b0
--- /dev/null
+++ b/sysdeps/i386/fpu/w_exp10f.c
@@ -0,0 +1,41 @@ 
+/* Copyright (C) 2015 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, see
+   <http://www.gnu.org/licenses/>.  */
+
+
+/*
+ * wrapper exp10f(x)
+ */
+
+#include <math.h>
+#include <math_private.h>
+
+float
+__exp10f (float x)
+{
+  float z = __ieee754_exp10f (x);
+  int32_t iz;
+  GET_FLOAT_WORD (iz, z);
+  if (__builtin_expect (!isfinite (z) || (iz & 0x7fffffff) == 0, 0)
+      && isfinite (x) && _LIB_VERSION != _IEEE_)
+    /* exp10f overflow (146) if x > 0, underflow (147) if x < 0.  */
+    return __kernel_standard_f (x, x, 146 + !!signbit (x));
+
+  return z;
+}
+weak_alias (__exp10f, exp10f)
+strong_alias (__exp10f, __pow10f)
+weak_alias (__pow10f, pow10f)
diff --git a/sysdeps/i386/fpu/w_exp2f.c b/sysdeps/i386/fpu/w_exp2f.c
new file mode 100644
index 0000000..bf86212
--- /dev/null
+++ b/sysdeps/i386/fpu/w_exp2f.c
@@ -0,0 +1,38 @@ 
+/* Copyright (C) 2015 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, see
+   <http://www.gnu.org/licenses/>.  */
+
+/*
+ * wrapper exp2f(x)
+ */
+
+#include <math.h>
+#include <math_private.h>
+
+float
+__exp2f (float x)
+{
+  float z = __ieee754_exp2f (x);
+  int32_t iz;
+  GET_FLOAT_WORD (iz, z);
+  if (__builtin_expect (!isfinite (z) || (iz & 0x7fffffff) == 0, 0)
+      && isfinite (x) && _LIB_VERSION != _IEEE_)
+    /* exp2 overflow: 144, exp2 underflow: 145 */
+    return __kernel_standard_f (x, x, 144 + !!signbit (x));
+
+  return z;
+}
+weak_alias (__exp2f, exp2f)
diff --git a/sysdeps/i386/fpu/w_expf.c b/sysdeps/i386/fpu/w_expf.c
new file mode 100644
index 0000000..5eb7f56
--- /dev/null
+++ b/sysdeps/i386/fpu/w_expf.c
@@ -0,0 +1,37 @@ 
+/* Copyright (C) 2015 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, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include <math_private.h>
+
+/* wrapper expf */
+float
+__expf (float x)
+{
+  float z = __ieee754_expf (x);
+  int32_t iz;
+
+  GET_FLOAT_WORD (iz, z);
+
+  if (__builtin_expect (!isfinite (z) || (iz & 0x7fffffff) == 0, 0)
+      && isfinite (x) && _LIB_VERSION != _IEEE_)
+    return __kernel_standard_f (x, x, 106 + !!signbit (x));
+
+  return z;
+}
+hidden_def (__expf)
+weak_alias (__expf, expf)