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]

Re: pow() problem



Shaobo Hou wrote:[color=blue]
> 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.[/color]

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(fa bs(x), 1.0/3.0), x)

--
Eric.Sosman@sun .com

[Michael Coyne]

Re: pow() problem

On Tue, 22 Feb 2005 15:12:57 +0000, Shaobo Hou said to the parser:
[color=blue]
> 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.[/color]

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;
}

[Walter Roberson]

Re: pow() problem

In article <cvfi5p$9c$1@wa pping.cs.man.ac .uk>,
Shaobo Hou <hous1@cs.man.a c.uk> wrote:
: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.


--
Those were borogoves and the momerathsoutgra be completely mimsy.

[Shaobo Hou]

Re: pow() problem

Eric Sosman wrote:
[color=blue]
> - If cbrt() isn't available, try something like
> (x >= 0.0) ? pow(x, 1.0/3.0) : -pow(-x, 1.0/3.0)[/color]

This is essentially what I did. thanks

[kyle york]

Re: pow() problem

Shaobo Hou wrote:[color=blue]
> 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?
>[/color]

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.
[color=blue]
> 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.[/color]

C99 has cbrt()

--
Kyle A. York
Sr. Subordinate Grunt

[Tom St Denis]

Re: pow() problem


Walter Roberson wrote:[color=blue]
> In article <cvfi5p$9c$1@wa pping.cs.man.ac .uk>,
> Shaobo Hou <hous1@cs.man.a c.uk> wrote:
> :Can anyone tell me why pow(-8.0, 1.0 / 3.0) (cubic root of -8)[/color]
returns[color=blue]
> :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[/color]
represented[color=blue]
> :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.[/color]

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

[Chris Croughton]

Re: pow() problem

On Tue, 22 Feb 2005 15:12:57 +0000, Shaobo Hou
<hous1@cs.man.a c.uk> wrote:
[color=blue]
> 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.[/color]

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(doub le 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

[Eric Sosman]

Re: pow() problem



Tom St Denis wrote:[color=blue]
> [...]
> That's nonsense. pow is normally implemented as
>
> pow(a, b) = exp(ln(a) * b)[/color]

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

--
Eric.Sosman@sun .com

[AC]

Re: pow() problem


"Shaobo Hou" <hous1@cs.man.a c.uk> wrote in message
news:cvfi5p$9c$ 1@wapping.cs.ma n.ac.uk...[color=blue]
> 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.[/color]

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]

Re: pow() problem

On Tue, 22 Feb 2005 13:38:10 -0500, AC
<test@test.test > wrote:
[color=blue]
> "Shaobo Hou" <hous1@cs.man.a c.uk> wrote in message
> news:cvfi5p$9c$ 1@wapping.cs.ma n.ac.uk...[color=green]
>> 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.[/color]
>
> The first argument has to be positive.[/color]

Or the second argument has to be finite and an integer value. Read
what you quoted:
[color=blue]
> 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.[/color]
^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^ ^^^^^^^^^^
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]

Re: pow() problem


"Chris Croughton" <chris@keristor .net> wrote in message
news:slrnd1n3f0 .bse.chris@ccse rver.keris.net. ..
^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^ ^^^^^^^^^^[color=blue]
> 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[/color]

You are right, I should have mentioned that too.

[Tom St Denis]

Re: pow() problem


Eric Sosman wrote:[color=blue]
> Tom St Denis wrote:[color=green]
> > [...]
> > That's nonsense. pow is normally implemented as
> >
> > pow(a, b) = exp(ln(a) * b)[/color]
>
> (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.[/color]

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]

Re: pow() problem



Tom St Denis wrote:[color=blue]
> Eric Sosman wrote:
>[color=green]
>>Tom St Denis wrote:
>>[color=darkred]
>>>[...]
>>>That's nonsense. pow is normally implemented as
>>>
>>>pow(a, b) = exp(ln(a) * b)[/color]
>>
>> (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
>>understandabl e explanation of the problems and exhibits an
>>implementatio n that mitigates them.[/color]
>
>
> pow(a, b) = exp(ln(a) * b)
>
> is EXACTLY correct.[/color]

My apologies; I'd assumed that because you said
"implemente d as" you intended this as a C expression
(with `log' misspelled) rather than as a mathematical
expression.
[color=blue]
> ... however as implemented [and executed] it may not be ideal.[/color]

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

--
Eric.Sosman@sun .com