[BZ,#18875] Avoid excess precision in check against float zero
Commit Message
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
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>.)
new file mode 100644
@@ -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)
new file mode 100644
@@ -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)
new file mode 100644
@@ -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)