computing t = pow(-11.5, .333)

[Charles Richmond]

If you want real roots (which may not exist for a quartic, but must
exist for a cubic) consider a Newton-Raphson solution. Off to the
books with you. This is now no longer a C question, but
algorithmic, and better suited to comp.programming.

OT: I would recommend Bairstow's method. Google it or see Wikipedia.

 

[pete]

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.552714e-015

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]

Unless you write a portable set of extra precision functions...

Thank you.
That's like the next level of a hard video game for me.
But I only write these things for my own amusement anyway,
so I guess it's time to go to the next level.