Keith Thompson wrote:
> exp(b*ln(a)) is the obvious and mathematically correct way to
> implement pow(a, b), but I've been told that it's not the best way to
> implement it, probably for reasons having to do with numerical
> stability. Sorry, I don't have any details. (Elsethread,
> P.J. Plauger mentions using this formula with extra precision.)
My portable freestanding version of pow,
which uses my freestanding versions of exp and log,
has trouble matching the accuracy my implementation's pow.
pow(0.0001, 0.25)  10 is 0.000000e+000
fs_pow(0.0001, 0.25)  10 is 3.552714e015
There's a lot of squaring in my fs_exp function,
which I suspect is the main problem.
Portable C code doesn't have access to any types
which are guaranteed to be more precise than double.
double fs_exp(double x)
{
unsigned n, square;
double b, e;
static double x_max;
if (1 > x_max) {
x_max = fs_log(DBL_MAX);
}
if (x_max >= x && x >= x_max) {
for (square = 0; x > 1; x /= 2) {
++square;
}
while (1 > x) {
++square;
x /= 2;
}
e = b = n = 1;
do {
b /= n++;
b *= x;
e += b;
b /= n++;
b *= x;
e += b;
} while (b > DBL_EPSILON / 4);
while (square != 0) {
e *= e;
}
} else {
e = x > 0 ? DBL_MAX : 0;
}
return e;
}

pete
