svn commit: r279856 - in head/lib/msun: man src

[Bruce Evans]

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