Date:Mon, 3 Jun 2013 06:50:12 +1000 (EST)From:Bruce Evans <brde@optusnet.com.au>To:numerics@freebsd.orgSubject:interesting bug in cexpf() caused by extra exponent rangeMessage-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: <http://docs.FreeBSD.org/cgi/mid.cgi?20130603061718.M7113>