Implementation of expl()
brde at optusnet.com.au
Sun Dec 9 21:18:51 PST 2007
On Sun, 9 Dec 2007, David Schultz wrote:
> Some minor fixes (mostly whitespace) are below.
> Also, don't you lose precision when you do stuff like this?
> z.e = 1 / z.e;
> In any case, if you can get me the appropriate constants for the
> 128-bit format, I'll clean everything up and check it in, which
> will make a lot of ports maintainers happy. That will pave the way
> to other stuff (e.g., MD versions of this) as well. We can worry
> about subnormals later.
Why not convert the fdlibm algorithm for exp() as is done for expf()?
It is much better (*), doesn't need to be debugged (except for the
conversion), and would be easier to maintain.
Better algorithms exist, like someone named das@ used for exp2(), but
would be harder to debug and maintain.
amd64 should probably use hardware for expl(), since long doubles are
much harder to handle efficiently than doubles. Thus the software
algorithm is needed mainly for ia64 and sparc64.
(*) Starting with the transformation of exp(x) to x*(exp(x)+1)/(exp(x)-1).
The power series for the latter converges much faster and has better
numeric properties. Polynomial approximations inherit these properties.
More information about the freebsd-standards