Date: Mon, 3 Jun 2013 06:50:12 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: numerics@freebsd.org Subject: interesting bug in cexpf() caused by extra exponent range Message-ID: <20130603061718.M7113@besplex.bde.org>
next in thread | raw e-mail | index | archive | help
I found that cexpf(192.066 + I*1.176e-38) is broken with clang iff
clang uses the i387 (the default without -march). Only the imaginary
part is detected as being broken by my tests.
This is because the imaginary part is evaluated as expf(192.066) *
sinf(1.176e-38). expf(192.066) is expected to overflow, but it
is evaluated as huge*huge = 1e60 in extra precision and exponent
range. clang but not gcc manages to keep the extras.
sinf(1.176e-38) is tiny and multiplication reduces the value to
below FLT_MAX so that it no longer overflows to Inf when the
extra exponent range is killed.
I think that it is correct and good to evaluate expf(192.066) with
extra precision and exponent range and then use these extras (except
C11 breaks the returning of the extras). However, 1e60 is not the
correct evaluation, which is 2.589e83. With the correct evaluation,
the result of the multiplication (3.047e45) would exceed FLT_MAX
in extra precision/range so it would overflow to Inf in float
precision.
I think the bug is in cexpf(). Its thresholds depend on overflow
occuring, but without a strict assignment there is no gurantee that
overflow occurs:
@ if (hx >= exp_ovfl && hx <= cexp_ovfl) {
@ /*
@ * x is between 88.7 and 192, so we must scale to avoid
@ * overflow in expf(x).
@ */
@ return (__ldexp_cexpf(z, 0));
@ } else {
@ /*
@ * Cases covered here:
@ * - x < exp_ovfl and exp(x) won't overflow (common case)
@ * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0
Actually, it doesn't overflow for all s > 0, because exp(x) is garbage
instead of Inf as expected.
@ * - x = +-Inf (generated by exp())
@ * - x = NaN (spurious inexact exception from y)
@ */
@ exp_x = expf(x);
The garbage would be converted to Inf if the compiler were perfectly pessimized
to C90/99/(11?) spec.
There is the usual bizarreness that if this used exp(x) to intentionally
get extra precision/range, then the assignment would work and and lose the
extras. Here the extra precision would be useful in the non-overflowing
case, but the algorithm depends on overflow occuring. It would have to
do a little more classification to force a value of Inf above cexp_ovfl
if it used exp() and assigned it to a double or long double so as to
keep the extra precision.
@ return (cpackf(exp_x * cosf(y), exp_x * sinf(y)));
@ }
@ }
I think a strict assignment is needed to fix this.
Similarly in cexp().
This is fragile in user code too. Naive code might use expf(190) * 1e-38
(no assignment). Then overflow shouldn't occur, but the result is garbage.
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130603061718.M7113>