R help - Dec 2010 - NaN not expected

Dear mailing list,

Why does the following code produce numerical results for x.pos.l, but NaNs for
x.neg.l?

x.pos <- function(tau.e, tau.c){
  tau.e + ((-tau.e^3 + 3*tau.e^2*tau.c - 
  3*tau.e*tau.c^2 + tau.c^3)^(1/3))/(2*2^(1/3))
}
(x.pos.l <- x.pos(1, 2))

x.neg <- function(tau.c, tau.i){
  tau.i + ((tau.c^3 - 3*tau.c^2*tau.i + 
  3*tau.c*tau.i^2 - tau.i^3)^(1/3))/(2*2^(1/3))
}
(x.neg.l <- x.neg(2, 3))

Many thanks for help,

Alrik

On 17/12/2010 11:40 AM, Thiem Alrik wrote:
> Dear mailing list,
>
> Why does the following code produce numerical results for x.pos.l, but NaNs
for x.neg.l?
>
> x.pos<- function(tau.e, tau.c){
>    tau.e + ((-tau.e^3 + 3*tau.e^2*tau.c -
>    3*tau.e*tau.c^2 + tau.c^3)^(1/3))/(2*2^(1/3))
> }
> (x.pos.l<- x.pos(1, 2))
>
> x.neg<- function(tau.c, tau.i){
>    tau.i + ((tau.c^3 - 3*tau.c^2*tau.i +
>    3*tau.c*tau.i^2 - tau.i^3)^(1/3))/(2*2^(1/3))
> }
> (x.neg.l<- x.neg(2, 3))

The expression you're taking a cube root of

tau.c^3 - 3*tau.c^2*tau.i +
   3*tau.c*tau.i^2 - tau.i^3

is negative.  You can't raise negative numbers to fractional powers in R. 
(This has been discussed lots of times on the list; use Google if you don't
understand why it has to be that way.)  What you should do is write a cuberoot
function, something like

cuberoot<- function(x) ifelse(x<  0, -(-x)^(1/3), x^(1/3))

and use that instead.

Duncan Murdoch