This is the mail archive of the libc-alpha@sources.redhat.com mailing list for the glibc project.
| Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
|---|---|---|
| Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
| Other format: | [Raw text] | |
Roger Sayle <roger@eyesopen.com> writes:
> On Sat, 20 Dec 2003, Andreas Jaeger wrote:
>> Roger Sayle <roger@eyesopen.com> writes:
>> > On Fri, 19 Dec 2003, Andreas Jaeger wrote:
>> >> I analyzed the y0 (0.0) failure for double:
>> >>
>> >> In e_j0.c we check for y0(0.0) with:
>> >> if((ix|lx)==0) return -one/zero;
>> >>
>> >> Previously GCC evaluated this at compiletime and generated a NaN but
>> >> since this was done at compile-time, no exception was set at run-time
>> >>
>> >> Now GCC does the calculation at run-time and sets the exception
>> >> (Roger, is this correct? I remember you doing work in this area).
>> >>
>> >> What do you think? Is GCC correct and we should change glibc? In
>> >> that case I'll make the changes to the math testsuite.
>> >
>> > Hi Andreas,
>> >
>> > Indeed the change was made to GCC to not evaluate floating point
>> > division by zero in in-line code at compile-time, precisely to allow
>> > the observable floating point exception to be raised if that division
>> > ever gets executed.
>> >
>> > The solution, given that this code is supposed to return a constant,
>> > is to change the glibc source code such that the division is performed
>> > in the initializer of a global or static variable. In these cases, the
>> > division by zero, and conversion into +-Inf or NaN is performed at
>> > compile-time. This is the idiom used by many of the testcases in
>> > gcc.c-torture/execute/ieee that needs such non-finite values.
>>
>> I'm not sure what's the correct thing is for glibc, the standards do
>> not mention the behaviour of y0/y1/yn AFAIK. So, how shall we update
>> glibc?
>
>
> My apologies in advance, if the above question was directed at the
> libc-alpha list...
No problem at all - I wasn't clear...
>
> The first step is to get a consensus from the glibc folks on what the
> correct behaviour should be. Whether y0(0.0) is expected to return
> -Inf with signal, or -Inf without signal. Versions of GCC prior
> to 3.4 have a bug which makes it difficult to ensure that the correct
> exception is generated at run-time, but this has now been fixed.
>
> A good example of a well defined interface is atanh which includes:
>
> * Special cases:
> * atanh(x) is NaN if |x| > 1 with signal;
> * atanh(NaN) is that NaN with no signal;
> * atanh(+-1) is +-INF with signal.
>
>
>
> If the consensus is that y0(0.0) shouldn't signal, then the solution
> is to use a patch such as the one below:
>
>
> Index: e_j0.c
> ===================================================================
> RCS file: /cvs/glibc/libc/sysdeps/ieee754/dbl-64/e_j0.c,v
> retrieving revision 1.4
> diff -c -3 -p -r1.4 e_j0.c
> *** e_j0.c 13 Feb 2001 01:23:45 -0000 1.4
> --- e_j0.c 20 Dec 2003 17:09:36 -0000
> *************** S[] = {0.0, 1.56191029464890010492e-02
> *** 92,101 ****
> --- 92,106 ----
>
> #ifdef __STDC__
> static const double zero = 0.0;
> + static const double neginf = -1.0 / 0.0;
> + static const double nan = 0.0 / 0.0;
> #else
> static double zero = 0.0;
> + static double neginf = -1.0 / 0.0;
> + static double nan = 0.0 / 0.0;
> #endif
>
> +
> #ifdef __STDC__
> double __ieee754_j0(double x)
> #else
> *************** V[] = {1.27304834834123699328e-02, /*
> *** 187,194 ****
> ix = 0x7fffffff&hx;
> /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
> if(ix>=0x7ff00000) return one/(x+x*x);
> ! if((ix|lx)==0) return -one/zero;
> ! if(hx<0) return zero/zero;
> if(ix >= 0x40000000) { /* |x| >= 2.0 */
> /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
> * where x0 = x-pi/4
> --- 192,199 ----
> ix = 0x7fffffff&hx;
> /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
> if(ix>=0x7ff00000) return one/(x+x*x);
> ! if((ix|lx)==0) return neginf; /* Return -Inf without exception */
> ! if(hx<0) return nan; /* Return NaN without exception */
> if(ix >= 0x40000000) { /* |x| >= 2.0 */
> /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
> * where x0 = x-pi/4
See below for comments.
>
>
> Of course, this issue applies to y0, y1 and yn, and the above changes
> will need to be made consistently to flt-32, dbl-64, ldbl-96 and ldbl-128.
I'll take care of this once we agree on the course of action. Thanks
Roger!
> I'm happy to prepare this large patch myself, once Uli and the other
> glibc folks have decided what the intended behaviour is expected to be.
> Clearly, any code that depends upon y0 trapping or not trapping without
> an applicable standard, is not going to be portable. But I also
> appreciate the argument for keeping y0 non-signaling to preserve the
> backwards compatible behaviour of previous glibc/gcc combinations.
Reading through the Posix Standard, I noticed that:
- y0 (-1.0) has to set the "Invalid operation" exception
- y0 (0.0) has to to set the "Overflow exception" (and not the divide
by zero one!)
From the Posix documentation on y0/y1/yn:
RETURN VALUE
Upon successful completion, these functions shall return the relevant Bessel value of x of the
second kind.
If x is NaN, NaN shall be returned.
If the x argument to these functions is negative, -HUGE_VAL or NaN shall be returned, and a
domain error may occur.
If x is 0.0, -HUGE_VAL shall be returned and a range error may occur.
If the correct result would cause underflow, 0.0 shall be returned and a range error may occur.
If the correct result would cause overflow, -HUGE_VAL or 0.0 shall be returned and a range
error may occur.
ERRORS
These functions may fail if:
Domain Error The value of x is negative.
If the integer expression (math_errhandling & MATH_ERRNO) is non-zero,
then errno shall be set to [EDOM]. If the integer expression (math_errhandling
& MATH_ERREXCEPT) is non-zero, then the invalid floating-point exception
shall be raised.
Range Error The value of x is 0.0, or the correct result would cause overflow.
If the integer expression (math_errhandling & MATH_ERRNO) is non-zero,
then errno shall be set to [ERANGE]. If the integer expression
(math_errhandling & MATH_ERREXCEPT) is non-zero, then the overflow
floating-point exception shall be raised.
Range Error The value of x is too large in magnitude, or the correct result would cause
underflow.
If the integer expression (math_errhandling & MATH_ERRNO) is non-zero,
then errno shall be set to [ERANGE]. If the integer expression
(math_errhandling & MATH_ERREXCEPT) is non-zero, then the underflow
floating-point exception shall be raised.
I'm appending a patch for the testsuite and for the dbl-64 vrsion of
y0. If this is ok, I'll change all other functions accordingly. The
patch below gives correct answers now for y0 (double) with both GCC
3.4 and 3.3.
Uli, what do you suggest? Ok to commit and change the other files
accordingly?
Andreas
2003-12-28 Andreas Jaeger <aj@suse.de>
* math/libm-test.inc (yn_test): Expect invalid exception for
negative arguments.
(y0_test): Likewise.
(y1_test): Likewise.
* sysdeps/ieee754/dbl-64/e_j0.c (y0): Raise only overflow for y0
(0). Raise Invalid exception for negative args.
============================================================
Index: math/libm-test.inc
--- math/libm-test.inc 7 Dec 2003 03:24:31 -0000 1.55
+++ math/libm-test.inc 28 Dec 2003 10:40:38 -0000
@@ -4145,7 +4145,7 @@ y0_test (void)
/* y0 is the Bessel function of the second kind of order 0 */
START (y0);
- TEST_f_f (y0, -1.0, minus_infty);
+ TEST_f_f (y0, -1.0, minus_infty, INVALID_EXCEPTION);
TEST_f_f (y0, 0.0, minus_infty);
TEST_f_f (y0, nan_value, nan_value);
TEST_f_f (y0, plus_infty, 0);
@@ -4179,7 +4179,7 @@ y1_test (void)
/* y1 is the Bessel function of the second kind of order 1 */
START (y1);
- TEST_f_f (y1, -1.0, minus_infty);
+ TEST_f_f (y1, -1.0, minus_infty, INVALID_EXCEPTION);
TEST_f_f (y1, 0.0, minus_infty);
TEST_f_f (y1, plus_infty, 0);
TEST_f_f (y1, nan_value, nan_value);
@@ -4214,7 +4214,7 @@ yn_test (void)
START (yn);
/* yn (0, x) == y0 (x) */
- TEST_ff_f (yn, 0, -1.0, minus_infty);
+ TEST_ff_f (yn, 0, -1.0, minus_infty, INVALID_EXCEPTION);
TEST_ff_f (yn, 0, 0.0, minus_infty);
TEST_ff_f (yn, 0, nan_value, nan_value);
TEST_ff_f (yn, 0, plus_infty, 0);
@@ -4228,7 +4228,7 @@ yn_test (void)
TEST_ff_f (yn, 0, 10.0, 0.0556711672835993914244598774101900481L);
/* yn (1, x) == y1 (x) */
- TEST_ff_f (yn, 1, -1.0, minus_infty);
+ TEST_ff_f (yn, 1, -1.0, minus_infty, INVALID_EXCEPTION);
TEST_ff_f (yn, 1, 0.0, minus_infty);
TEST_ff_f (yn, 1, plus_infty, 0);
TEST_ff_f (yn, 1, nan_value, nan_value);
============================================================
Index: sysdeps/ieee754/dbl-64/e_j0.c
--- sysdeps/ieee754/dbl-64/e_j0.c 13 Feb 2001 01:23:45 -0000 1.4
+++ sysdeps/ieee754/dbl-64/e_j0.c 28 Dec 2003 11:00:23 -0000
@@ -185,9 +185,9 @@ V[] = {1.27304834834123699328e-02, /*
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
- /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
+ /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */
if(ix>=0x7ff00000) return one/(x+x*x);
- if((ix|lx)==0) return -one/zero;
+ if((ix|lx)==0) return -HUGE_VALF+x; /* -inf and overflow exception. */
if(hx<0) return zero/zero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
--
Andreas Jaeger, aj@suse.de, http://www.suse.de/~aj
SuSE Linux AG, Maxfeldstr. 5, 90409 Nürnberg, Germany
GPG fingerprint = 93A3 365E CE47 B889 DF7F FED1 389A 563C C272 A126
Attachment:
pgp00000.pgp
Description: PGP signature
| Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
|---|---|---|
| Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |