Date: Sun, 2 Jun 2013 18:11:30 -0700 From: David Schultz <das@freebsd.org> To: Bruce Evans <brde@optusnet.com.au> Cc: numerics@freebsd.org Subject: Re: interesting bug in cexpf() caused by extra exponent range Message-ID: <20130603011130.GA84261@zim.MIT.EDU> In-Reply-To: <20130603061718.M7113@besplex.bde.org> References: <20130603061718.M7113@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Mon, Jun 03, 2013, Bruce Evans wrote: > 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: I think you have identified two potential bugs, but neither are in cexpf(): one is in the compiler, and the other is in expf(). First the compiler bug: cexpf() is careful to assign the result of expf() to a float variable, and the C standard requires that the compiler not retain any extra precision. Second, the expf() bug: If the compiler decides to return the result in a register that has more than float precision, expf() needs to be prepared to deal with that. It can return INFINITY, or it can return the correct answer (which will overflow to infinity when the extra precision is dropped). Returning 1e60 is bogus. I think just using a macro, e.g., RETURN(huge * huge), to ensure that the extra precision gets dropped, is the most straightforward solution. > 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. Yes, I think this is why it needs to be fixed in expf(), not in callers of expf().
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130603011130.GA84261>