why does similar code generate complex numbers?

In the following code, the array 'f1' contains real numbers while the array 'f2' contains complex numbers. From what I understand and am expecting, both arrays should contain identical data, consisting only of real numbers. But this is not what I get. Instead, the array 'f2', generated using exponential notation, contains complex numbers and the question is why is this happening?

x = [-27 -8 -1 1 8 27];
f1   = nthroot(x, 3);
f2   = x .^ (1/3);

The array 'f1' contains the expected results:

-3 -2 -1 1 2 3

For the array 'f2', here are the numbers generated:

1.5000 + 2.5981i 1.0000 + 1.7321i 0.5000 + 0.8660i 1.0000 + 0.0000i 2.0000 + 0.0000i 3.0000 + 0.0000i

The array 'f2' unexpectedly contains complex numbers. The last 3 numbers give the correct/expected result but the first 3 numbers are not the correct/expected result. That is, the first 3 numbers in 'f2' are NOT the cubed root of the first 3 numbers in 'x'. Why is this happening?

In trying to debug/understand the problem, when I manually enter the negative numbers of 'x' in the command window then I get the expected results:

-27 ^ (1/3) ---> ans = -3

-8 ^ (1/3) ---> ans = -2

-1 ^ (1/3) ---> ans = -1

However, when I use array syntax in the command window and select the same negative numbers from 'x' then I get unexpected and incorrect results:

x(1) ^ (1/3) ---> ans = 1.5000 + 2.5981i

x(2) ^ (1/3) ---> ans = 1.0000 + 1.7321i

x(3) ^ (1/3) ---> ans = 0.5000 + 0.8660i

I have no experience working with complex numbers. The question I am asking is why do I get incorrect results (i.e. in complex number format) when using exponential notation in my script?

Thanks to all who respond

Answer 1

[Walter Roberson]

When x is negative, x^n is defined by MATLAB as exp(n*log(x)) . log(x) when x is negative, is log(-x) + pi * 1i . n times that is n*log(-x) + n*pi*1i . exp() of that is going to be complex except when n is an integer (in which case n*pi becomes an exact multiple of pi)

Thus for x negative, f2 = x .^ (1/3); is going to be exp(log(-x)/3 + pi/3*1i) which is going to be complex valued.

nthroot() is defined in terms of real roots, with it being an error for the value to be negative and the power to not be an odd integer.

[Walter Roberson]

 -27 ^ (1/3) 

The - is formally an operator, not a part of the syntax for building numbers. Formally speaking, -27 is not "an indivisible unit with value negative 27": it is "unary minus applied to the expression to the right of it".

And the expression to the right of the - in 27 ^ (1/3) is all of 27 ^ (1/3) -- which is positive 3. - is then applied to that, giving negative 3.

Another way of writing this is to say that ^ and .^ have higher binding priority than - has .

This is not the situation in all programming languages: in some programming languages, - followed by a numeric constant has higher binding precedence than ^ has. But those languages differ on whether the syntax -X ^ Y is otherwise to be interpreted as (-X) ^ Y or as -(X^Y) if X is not a numeric constant.

You need to use

 (-27) ^ (1/3)

[IB00]

Operator precedence! Good catch. My mistake. Thanks Walter.

I just saw your first post. You wrote, "When x is negative, x^n is defined by MATLAB as exp(n*log(x))".

1) Does Matlab define x^n the same way when x is positive? I assume yes but, given I'm in a bit over my head, I am asking just to confirm.

2) Would you be able to refer me to the section in the Matlab documentation where I can read more about how "x^n is defined when x is negative"?

Thanks again Walter

[Walter Roberson]

If I recall correctly, MATLAB actually uses log() and exp() when the base is negative, even for integer exponents . At some point I saw someone demonstrate this based on precision, a case where log() / exp() exactly matched the observed result but repeated multiplication did not exactly match.. but that would have been a number of years ago now.

MATLAB does not definitely use log() and exp() when the base is positive. I believe that it does so when the exponent is not an integer, but for the case where the exponent is a hard-coded integer, it might not. People have done some timing tests, which I have mostly forgotten the details of now. I seem to recall that X.^2 has timing consistent with being implemented as X*X . My memory seems to think that X.^4 had timing consistent with X*X*X*X, and not consistent with temp = X*X; temp*temp. My memory is telling me that the precision results were consistent with repeated multiplication, not with log() / exp() for positive integer expontnets. Someone submitted a FEX contribution to do faster integer exponents by decomposing the power into binary.

I do not know what the algorithm is for nthroot; I have seen cases where nthroot(x, N) did not exactly equal x.^(1/N), and my memory is telling me that x.^(1/N) in that case was exactly equal to the log/exp solution.

Answer 2

[John D'Errico]

What, for example, are the THREE cube roots of a negative number? We can find them all as solutions to the equation

syms x
xsol = solve(x^3 == -4,'maxdegree',3)
xsol =
                        -2^(2/3)
  2^(2/3)*((3^(1/2)*1i)/2 + 1/2)
 -2^(2/3)*((3^(1/2)*1i)/2 - 1/2)

vpa(xsol)
ans =
                          -1.5874010519681994748
 0.79370052598409973738 + 1.3747296369986026264i
 0.79370052598409973738 - 1.3747296369986026264i

So -4 has 3 cube roots. We can plot them in the complex plane.

plot(double(xsol),'o')
xline(0);
yline(0);

Which one has the smallest phase angle in the complex plane?

angle(xsol)
ans =
    pi
  pi/3
 -pi/3

So the second root it found is the one that lives at an angle of pi/3 in the complex plane, or 60 degrees for the degree fanatics. This is the root with the smallest complex phase angle, and is the one chosen by power.

NTHROOT explicitly looks for the real root.

nthroot(-4,3)
ans =
      -1.5874

It tells you that.

 help nthroot
 nthroot Real n-th root of real numbers.
 
    nthroot(X, N) returns the real Nth root of the elements of X.
    Both X and N must be real, and if X is negative, N must be an odd integer.

[IB00]

Thanks for helping out. I may be wrong but I believe you have not understood the question I was asking. I have just edited my original post in order to specfically explain what the problem is. If you have a chance, can you review my edit. Thanks.