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]

Re: y0/y1/yn and exceptions - GCC 3.4 changes


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...


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


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'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.


Roger
--


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]