Implementation of expl()

[Bruce Evans]

On Mon, 10 Dec 2007, David Schultz wrote:

> On Mon, Dec 10, 2007, Bruce Evans wrote:

>> 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.

But every different precision requires a significantly different
implementation (mainly for magic number, but table-driven methods
give a huge number of those) since we don't have algorithms or
tools to generate all the magic numbers from the precision parameter.

> 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.

The problems are well understood :-).  I didn't understand them when
you started working on exp2() but do now.  Essentially, they are due
to the strange phenomenon that linear functions are usually harder to
approximate (to within 1 ULP) than transcendental functions expanded
in a Taylor series about 0.  The linear scaling step with a not-exactly
representable factor like log2(e) ends up being the most difficult point
(unless it is done inaccurately).

> 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.

I only agree if "working" includes no new bugs and good accuracy at little
cost.

We still don't support 64-bit precision on i386.  Fixing that has higher
priority.

Bruce