Implementation of expl()

[David Schultz]

On Mon, Dec 10, 2007, Bruce Evans wrote:
> 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.

I'm not worried about maintenance, since I don't expect God to add
any major new features to e^x any time soon. Writing exp2() took a
lot of reading papers and tinkering, and it's a pain in the neck
to generate the constants and figure out the resulting error for
each interval. I seem to recall trying the same tricks for expl(),
but there were problems with rescaling in base e instead of base 2
without losing accuracy. In any case, I don't have the kind of time
needed to fix all of that stuff now.

I don't really care how expl() is implemented; anything that works
is better than anything that doesn't. If someone comes up with a
better, cleverer scheme later, that's great, but we've been
talking about a lot of this stuff forever and there's still
nothing in the tree.