Possible mistake in expf

Joseph Myers joseph@codesourcery.com
Mon Jun 15 10:02:00 GMT 2015

On Thu, 11 Jun 2015, ethereal@ethv.net wrote:

> Hello everyone,
> I believe I may have stumbled across a minor mistake in the expf()
> implementation in newlib's libm.
> newlib/libm/mathfp/sf_pow.c:99 [1] is:
> > x = exp (t);
> However, I don't see why the function exp should be used instead of
> expf, as x itself is a single-precision float. At a quick glance I don't
> see any reason why this would benefit greatly from the increased
> precision the double-precision version of exp would provide, but I'm far
> from an expert on this sort of thing. If someone could confirm or refute
> my suspicions, that would be appreciated!

Since t is float, it won't actually help - but to get good results, t 
ought to be double (note that it's computed using the double version of 
log).  To get results within a few ulp for pow (x, y) as exp (y * log (x)), 
you need y * log (x) to have as many fractional bits as there are bits in 
the mantissa of the result type (that is, for results with large 
exponents, the mantissa precision needed for y * log (x) is about the sum 
of mantissa and exponent bits in the type in question, or about the number 
of bits in the representation for an IEEE type - so 32 bits intermediate 
precision is desirable for powf, rather than just the 24 bits of float).

Joseph S. Myers

More information about the Newlib mailing list