pow() problem

[Shaobo Hou]

Can anyone tell me why pow(-8.0, 1.0 / 3.0) (cubic root of -8) returns
nan (in linux) and negative infinity or something (in devcpp in
windows), instead of -2?

The problem seems to be that pow can't handle cubic root of negative
number, I mean, a calculator could do it or am I using the wrong
function here. My current solution is a hack which detects negative
input and handle it differently.

I suspect it is something to do with the fact the power is represented
as a floating point value.

 

[Eric Sosman]

Can anyone tell me why pow(-8.0, 1.0 / 3.0) (cubic root of -8) returns
nan (in linux) and negative infinity or something (in devcpp in
windows), instead of -2?

The problem seems to be that pow can't handle cubic root of negative
number, I mean, a calculator could do it or am I using the wrong
function here. My current solution is a hack which detects negative
input and handle it differently.

I suspect it is something to do with the fact the power is represented
as a floating point value.

In a way, yes: the second argument to pow() is not
exactly one-third, but a nearby value. You aren't actually
calculating "the cube root of -8," but "-8 raised to a
rational exponent close to one-third." That exponent value
is a fraction of the form U/V where U,V are integers and V is
almost certainly a large power of two. Mathematically
speaking,

pow(-8, U/V)
== pow(pow(-8, 1/V), U)
== pow(sqrt(sqrt(...(-8)...)), U)

(since V is a power of two), and this can't be evaluated
in real arithmetic.

Some suggestions:

- The latest "C99" Standard defines a cbrt() function.
Even if you don't have access to a full-blown C99
implementation, you may find that cbrt() is present.

- If cbrt() isn't available, try something like
(x >= 0.0) ? pow(x, 1.0/3.0) : -pow(-x, 1.0/3.0)

- If your system has copysign(), try writing the above as
copysign(pow(fabs(x), 1.0/3.0), x)

 

[Michael Coyne]

Can anyone tell me why pow(-8.0, 1.0 / 3.0) (cubic root of -8) returns nan
(in linux) and negative infinity or something (in devcpp in windows),
instead of -2?

The problem seems to be that pow can't handle cubic root of negative
number, I mean, a calculator could do it or am I using the wrong function
here. My current solution is a hack which detects negative input and
handle it differently.

I suspect it is something to do with the fact the power is represented as
a floating point value.

It has more to do with the -8...

The below demonstrates one method for handling this, which may be what you
did already, when you mention your current solution detects negative input
and handles it differently.


#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(void)
{
float num = -8;
float result;

if (num >= 0)
result = pow(num, 1.0/3.0);
else
result = -pow(-num, 1.0/3.0);

printf("%f\n", result);

return EXIT_SUCCESS;
}

 

[Shaobo Hou]

- If cbrt() isn't available, try something like
(x >= 0.0) ? pow(x, 1.0/3.0) : -pow(-x, 1.0/3.0)

This is essentially what I did. thanks

 

[Chris Croughton]

Can anyone tell me why pow(-8.0, 1.0 / 3.0) (cubic root of -8) returns
nan (in linux) and negative infinity or something (in devcpp in
windows), instead of -2?

The problem seems to be that pow can't handle cubic root of negative
number, I mean, a calculator could do it or am I using the wrong
function here. My current solution is a hack which detects negative
input and handle it differently.

I suspect it is something to do with the fact the power is represented
as a floating point value.

The parameters to pow(x, y) are always floating point. However:

7.12.7.4 The pow functions

2 The pow functions compute x raised to the power y. A domain error
occurs if x is negative and y is finite and not an integer value.

Among other things, 1.0 / 3.0 is not 1/3, it's an approximation, so the
function can't tell that it has a real answer at all (only fractional
powers with an odd integral divisor have real roots, others have complex
roots.

C99 has the function cbrt() which calculates a cube root specifically
(as, I suspect, does your calculator). However, there are few
conforming C99 libraries yet (after all, that specification has only
been out for 6 years!).

In general it is better to give an error if the programmer tries to do
something which has an unrepresentable result.

I would implement my own cube root function, either "from scratch"
(using Newton/Raphson or a better algorithm) or as:

double myCubeRoot(double x)
{
return (x >= 0 ? pow(x, 1.0/3.0) : -pow(-x, 1.0/3.0));
}

(or as a macro).

Note that since pow() uses exp(log(x)*y) (or a faster equivalent) it
will in general be slower than a special-purpose cube root function, as
well as being less accurate (see above about binary approcimation to
fractions).

Chris C

 

[Walter Roberson]

:Can anyone tell me why pow(-8.0, 1.0 / 3.0) (cubic root of -8) returns
:nan (in linux) and negative infinity or something (in devcpp in
:windows), instead of -2?

:I suspect it is something to do with the fact the power is represented
:as a floating point value.

Plausibly. 1/3 is not exactly representable in binary, so whether
1/3 comes out "odd" or "even" (which would allow you to square the
x first) would be dependant on the precision you are working with.

 

[kyle york]

Can anyone tell me why pow(-8.0, 1.0 / 3.0) (cubic root of -8) returns
nan (in linux) and negative infinity or something (in devcpp in
windows), instead of -2?

From C99 section 7.12.7.4, the pow functions:

``A domain error occurs if x is finite and negative and y is finite and
not an integer value.''

pow() is a generic power function and in general a negative number
raised to a non-integral power will result in a complex number. It would
need to special case x**(1/(2**n)) and x**(1/(2**n+1)) which would be
tricky.

The problem seems to be that pow can't handle cubic root of negative
number, I mean, a calculator could do it or am I using the wrong
function here. My current solution is a hack which detects negative
input and handle it differently.

C99 has cbrt()

 

[Tom St Denis]

:Can anyone tell me why pow(-8.0, 1.0 / 3.0) (cubic root of -8) returns
:nan (in linux) and negative infinity or something (in devcpp in
:windows), instead of -2?

:I suspect it is something to do with the fact the power is represented
:as a floating point value.

Plausibly. 1/3 is not exactly representable in binary, so whether
1/3 comes out "odd" or "even" (which would allow you to square the
x first) would be dependant on the precision you are working with.

That's nonsense. pow is normally implemented as

pow(a, b) = exp(ln(a) * b)

It may or may not abs() "a" which is where the problems will occur
because the ln of -a is not defined.

Tom

 

[Eric Sosman]

[...]
That's nonsense. pow is normally implemented as

pow(a, b) = exp(ln(a) * b)

(Slight topic drift): For suitable values of "normally,"
with implications of "quick and dirty." This is not a very
accurate implementation of pow(), and I wouldn't expect to
find it in a high-quality math library. Reference: "The
Standard C Library" by P.J. Plauger gives a brief but quite
understandable explanation of the problems and exhibits an
implementation that mitigates them.

 

[AC]

Can anyone tell me why pow(-8.0, 1.0 / 3.0) (cubic root of -8) returns nan
(in linux) and negative infinity or something (in devcpp in windows),
instead of -2?

The problem seems to be that pow can't handle cubic root of negative
number, I mean, a calculator could do it or am I using the wrong function
here. My current solution is a hack which detects negative input and
handle it differently.

I suspect it is something to do with the fact the power is represented as
a floating point value.

The first argument has to be positive.

7.12.7.4 The pow functions

Synopsis

1 #include <math.h>

double pow(double x, double y);

float powf(float x, float y);

long double powl(long double x, long double y);

Description

The pow functions compute x raised to the power y. *A domain error occurs
if x is finite and negative* and y is finite and not an integer value. A
domain error may occur if x is zero and y is less than or equal to zero. A
range error may occur.

 

[Chris Croughton]

The first argument has to be positive.

Or the second argument has to be finite and an integer value. Read
what you quoted:

7.12.7.4 The pow functions

Synopsis

1 #include <math.h>

double pow(double x, double y);

float powf(float x, float y);

long double powl(long double x, long double y);

Description

The pow functions compute x raised to the power y. *A domain error occurs
if x is finite and negative* and y is finite and not an integer value.

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
There's an 'and' clause in there, pow(-8.0, 3) is valid, for instance
(the 3 is promoted to a double as long as the prototype is in scope).

Chris C

 

[AC]

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

There's an 'and' clause in there, pow(-8.0, 3) is valid, for instance
(the 3 is promoted to a double as long as the prototype is in scope).

Chris C

You are right, I should have mentioned that too.

 

[Tom St Denis]

[...]
That's nonsense. pow is normally implemented as

pow(a, b) = exp(ln(a) * b)

(Slight topic drift): For suitable values of "normally,"
with implications of "quick and dirty." This is not a very
accurate implementation of pow(), and I wouldn't expect to
find it in a high-quality math library. Reference: "The
Standard C Library" by P.J. Plauger gives a brief but quite
understandable explanation of the problems and exhibits an
implementation that mitigates them.

pow(a, b) = exp(ln(a) * b)

is EXACTLY correct.

.... however as implemented [and executed] it may not be ideal. The
"crux" of what the pow function does is probably along the lines of the
equivalence since you can easily find exp and ln from two convergent
series.

Tom

 

[Eric Sosman]

[...]
That's nonsense. pow is normally implemented as

pow(a, b) = exp(ln(a) * b)

(Slight topic drift): For suitable values of "normally,"
with implications of "quick and dirty." This is not a very
accurate implementation of pow(), and I wouldn't expect to
find it in a high-quality math library. Reference: "The
Standard C Library" by P.J. Plauger gives a brief but quite
understandable explanation of the problems and exhibits an
implementation that mitigates them.

pow(a, b) = exp(ln(a) * b)

is EXACTLY correct.

My apologies; I'd assumed that because you said
"implemented as" you intended this as a C expression
(with `log' misspelled) rather than as a mathematical
expression.

... however as implemented [and executed] it may not be ideal.

This was the point of my post: Interpreted as a C
expression, exp(log(a) * b) is a poor way to implement
pow(a,b).