Date: Wed, 11 Mar 2015 06:35:32 +1100 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Steve Kargl <kargl@freebsd.org> Cc: svn-src-head@freebsd.org, svn-src-all@freebsd.org, src-committers@freebsd.org Subject: Re: svn commit: r279856 - in head/lib/msun: man src Message-ID: <20150311051030.D2799@besplex.bde.org> In-Reply-To: <201503101710.t2AHAt0V028926@svn.freebsd.org> References: <201503101710.t2AHAt0V028926@svn.freebsd.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Tue, 10 Mar 2015, Steve Kargl wrote: > Log: > According to POSIX.1-2008, the Bessel functions of second kind > should raise a divide-by-zero floating point exception for x = +-0 > and an invalid floating point exception for x < 0 including x = -Inf. > Update the code to raise the exception and update the documentation > with hopefully better description of the behavior. > > Reviewed by: bde (code only) Er, the functions have always (since Ng wrote them in 1992 or earlier) raised the exceptions, except with broken compilers like clang. > Modified: head/lib/msun/src/e_j0.c > ============================================================================== > --- head/lib/msun/src/e_j0.c Tue Mar 10 17:04:11 2015 (r279855) > +++ head/lib/msun/src/e_j0.c Tue Mar 10 17:10:54 2015 (r279856) > @@ -64,6 +64,8 @@ __FBSDID("$FreeBSD$"); > > static double pzero(double), qzero(double); > > +static const volatile double vone = 1, vzero = 0; > + This volatile hack works around the compiler bug. > static const double > huge = 1e300, > one = 1.0, > @@ -150,10 +152,16 @@ __ieee754_y0(double x) > > EXTRACT_WORDS(hx,lx,x); > 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; -one/zero must be evaluated at compile time, since it generates the divide-by-zero exception if it is evaluated correctly. But clang evaluates it to -Inf at compile time and doesn't generate the exception. > - if(hx<0) return zero/zero; zero/zero must be evaluated at compile time, since it generates the invalid operand exception if it is evaluated correctly. But clang evaluates it to the default quiet NaN at compile time and doesn't generate the exception. > + /* > + * y0(NaN) = NaN. > + * y0(Inf) = 0. > + * y0(-Inf) = NaN and raise invalid exception. > + */ > + if(ix>=0x7ff00000) return vone/(x+x*x); > + /* y0(+-0) = -inf and raise divide-by-zero exception. */ > + if((ix|lx)==0) return -one/vzero; > + /* y0(x<0) = NaN and raise invalid exception. */ > + if(hx<0) return vzero/vzero; The volatile variables are a trick to prevent compile-time evaluation. My review said to not organize the comments like this. They are verbose yet only hint that the return expressions are more subtle than they appear to be. > 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 Volatile hacks are already used a lot to work around this compiler bug. Sometimes even gcc needs them. This is a bug in clang because although translation-time evaluation is permitted unless the state of the FENV_ACCESS pragma is "on", clang doesn't support this pragma and its behaviour is fail-unsafe and the opposite of gcc's. gcc also doesn't support this pragma, but it does document its non-support. Implementations are required to document the default state of this pragma, but clang has almost null installed documentation. I have to look at the gcc documentation to see what arcane flags or pragmas like this do or should do. gcc48.info says: X '-frounding-math' X Disable transformations and optimizations that assume default X floating-point rounding behavior. This is round-to-zero for all X floating point to integer conversions, and round-to-nearest for all X other arithmetic truncations. This option should be specified for X programs that change the FP rounding mode dynamically, or that may X be executed with a non-default rounding mode. This option disables X constant folding of floating-point expressions at compile time X (which may be affected by rounding mode) and arithmetic X transformations that are unsafe in the presence of sign-dependent X rounding modes. X X The default is '-fno-rounding-math'. X X This option is experimental and does not currently guarantee to X disable all GCC optimizations that are affected by rounding mode. X Future versions of GCC may provide finer control of this setting X using C99's 'FENV_ACCESS' pragma. This command-line option will be X used to specify the default state for 'FENV_ACCESS'. X ... X * 'The default state for the 'FENV_ACCESS' pragma (C99 7.6.1).' X X This pragma is not implemented, but the default is to "off" unless X '-frounding-math' is used in which case it is "on". Note that this is mostly about rounding. With FENV_ACCESS on, even 1.0/3.0 cannot be evaluated at compile time. Not evaluating this at compile time would usually just be inefficient. So gcc defaults to FENV_ACCESS off. But gcc still doesn't evaluate at compile time expressions like 1.0/0.0 and 0.0/0.0 that would generate significant exeptions, even if you write these expressions with literal constants. -frounding-math works as advertised in gcc-4.2 and gcc-4.8. It causes the full inefficiencies, with 1.0/3.0 evaluated at run time with it but not without it, while 1.0/0.0 and 0.0/0.0 are evaluated at run time unconditionally. It acts sort of as if it turns FENV_ACCESS on. However, FENV_ACCESS itself is still broken -- it is silently ignored. -frounding-math is not supported in gcc-3.3.3. It is broken (doesn't give runtime evaluation of 1.0/0.0) in gcc-3.3.4. It is broken (unsupported = incompatible) in clang. FENV_ACCESS is broken (silently ignored) in clang. msun is being sloppy by using magic expressions to generate exceptions. It should sprinkle FENV_ACCESS as required. However, FENV_ACCESS requires a C99 compiler, and no C99 compilers are available. So msun has to use some unportable hacks. The ones that worked 20-30 years ago are at least short. Some of the other expessions that should cause exceptions are: (double) (1e300 * 1e300) This remains working in gcc and broken (no overflow exception) in clang. msun actually mostly uses 'return 1e300*1e300'. This causes different (worse) problems with gcc on i386. The expression is evaluated in extra precision/exponent range) and the extras are returned (except C11 breaks this). Then if the caller stores the result 1e600L to a double, an overflow exception occurs after all, but if the caller uses the extra precision then the garbage value 1e600L is used. (double) (1e-300 * 1e-300) This remains broken (no underflow exception) in both gcc and clang. This brokenness goes back to at least gcc-2.95. msun uses lots of volatile hacks to generate overflow and underflow, but it doesn't do this systematically so many of the less important functions remain unfixed. Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20150311051030.D2799>