17
$\begingroup$

Why (-1.)^2. in Mathematica returns a complex number? It looks like in both C and Fortran it returns 1. Why does Mathematica behave differently than the other systems? Is it because the way Mathematica calculates the power is unique compared to others? If this is true then what are the advantages of this design?

I'm curious because it seems cost a "bug" in my code. Sometimes the code is more than 10X slower than other time, and it took me a while to figure out that the problem comes the following function

f[n_] := Compile[{{x, _Real}}, Cos[x]^n] 

f[2] /@ Range[0, 4 π, π/10]; // AbsoluteTiming
(* {0.000258, Null} *)

f[2.] /@ Range[0, 4 π, π/10]; // AbsoluteTiming
(* {0.140602, Null} *)

Notice that the second evaluation of f is much slower than the first, because cos[x]^2. is treated as a complex function and the uncompiled version is invoked.

enter image description here

Here are some examples from other languages, it looks all of them gives 1 as the answer.

Fortran

enter image description here

C

enter image description here

Python

enter image description here

Matlab/Octave

enter image description here

$\endgroup$
  • 2
    $\begingroup$ The imaginary components which are involved in the internal algorithms don't cancel exactly due to precision issues. Use Chop to remove the imaginary artifacts. $\endgroup$ – Bob Hanlon Jan 13 '15 at 22:52
  • 3
    $\begingroup$ Actually it's nor surprising. (-1)^2.1 and (-1)^1.9 do not have zero imaginary parts. Im[(-1)^x] == Sin[Pi x] crosses zero at x=2, but when using inexact numbers there's always a chance that we'd get a slight deviation from precise zero. It's just numerical error. Why C and Fortran don't do this? If you work with a data type that doesn't support imaginary numbers at all then you can't get imaginary results. pow(-1., 2.01) is nan in C. $\endgroup$ – Szabolcs Jan 13 '15 at 22:53
  • 2
    $\begingroup$ In other words this question is: why is Sin[2. Pi] not exact 0? You don't get exact 0 in C++ either $\endgroup$ – Szabolcs Jan 13 '15 at 22:54
  • 4
    $\begingroup$ A machine double, raised to a machine double power, will in effect compute exp(pow*log(base)). In[2]:= 2.*Log[-1.] Out[2]= 0. + 6.28319 I In[3]:= Exp[%] -16 Out[3]= 1. - 2.44929 10 I $\endgroup$ – Daniel Lichtblau Jan 13 '15 at 23:57
  • 2
    $\begingroup$ @xslittlegrass the other systems follow IEEE754, which requires that they get the answer that they do. Mathematica doesn't. This still doesn't explain why Mathematica feels satisfied with 1.1 ULP error in this result; someone else will have to answer that one for you. $\endgroup$ – Oleksandr R. Jan 20 '15 at 21:45
6
$\begingroup$

lets make a fair fortran test... raise a complex typed number 1 to the power:

   write(*,*)(-1.d0,0.d0)**2.,(-1.d0,0.d0)**2

(1.00000000000,-2.449212707...E-16 ) , ( 1.0000000000,0.000000000 )

also note in fortran the real power can not yield a complex result, that is (-1.)**.5 throws an error or yields NaN depending on the compiler so you see (-1.)**2. works only because (/if) the compiler is smart enough to recognize that 2. should be treated as an integer.

**the first test was with the gfortran compiler, intel fortran gives the expected

( 1.0000000000,0.000000000 ) , ( 1.0000000000,0.000000000 )

$\endgroup$
  • $\begingroup$ Which compiler/flags? $\endgroup$ – Oleksandr R. Jan 20 '15 at 22:30
  • $\begingroup$ no flags. See the edit though - the results may differ depending on how the compiler optimizes the code. $\endgroup$ – george2079 Jan 20 '15 at 22:39
  • $\begingroup$ Interesting. gcc is usually stricter about IEEE compliance than the Intel compiler, so to see this result is quite surprising. On the other hand, Intel has the better math library, so if this is a library call, it's not so surprising after all. $\endgroup$ – Oleksandr R. Jan 20 '15 at 22:56
9
$\begingroup$

My original answer (below) is wrong. Arbitrary precision does not fix this problem, I only fooled myself (and others) into thinking that it does. Now please consider:

Power[-1`5, 2`5]
1.0000 + 0.*10^-5 I

This agrees with the output of the effective computation that Daniel described in a comment:

power[base_, pow_] := Exp[pow*Log[base]]

power[-1`5, 2`5]

1.0000 + 0.*10^-5 I

Following that evaluation we get a complex value from Log when base is negative:

Log[-1`5]

0.*10^-5 + 3.1416 I

Given the method used complex numbers seem inherent in the calculation.


Mistaken assertion

Round-off errors occur in floating point calculations. You can use arbitrary precision to avoid these:

(-1`5)^2
1.0000
$\endgroup$
  • 5
    $\begingroup$ I think the question is really: why is Mathematica's Power function less accurate than those in other languages? Here we have 1.1 ULP of error. Practically speaking we can say that the other languages follow IEEE754 (which requires exact rounding, so <0.5 ULP error in all operations), and Mathematica doesn't. Some justification of this is probably called for so long as Mathematica gives inferior results. But, interestingly, we do not have these problems as soon as you go to arbitrary precision. $\endgroup$ – Oleksandr R. Jan 20 '15 at 21:49
  • $\begingroup$ @Oleksandr Thanks for the comment. By the way I believe I was mistaken regarding arbitrary precision and I revised my answer accordingly. $\endgroup$ – Mr.Wizard Jan 21 '15 at 2:11
  • 1
    $\begingroup$ +1 for Dan's comment explaining it (transcendental functions are not included in the IEEE754 specification, so >0.5 ULP error should be acceptable for them, even if technically the implementation of Power should not be done this way if we are following the standard), but I don't follow why arbitrary precision doesn't resolve the problem? SetPrecision[-1, $MachinePrecision]^SetPrecision[2, $MachinePrecision] gives an exactly rounded result; remove the $ and this is lost. The arbitrary-precision Log and Exp implementations seem more accurate than their machine-precision counterparts. $\endgroup$ – Oleksandr R. Jan 21 '15 at 10:22
  • $\begingroup$ N.B. result = SetPrecision[-1, $MachinePrecision]^SetPrecision[2, $MachinePrecision]; Log[2, 10^$MachinePrecision] - Log[2, 10^Accuracy@Im[result]] gives 3.05. But I don't know what to make of this; I am not sure if it is significance arithmetic telling us that we should have lost this much precision, or if we really have. $\endgroup$ – Oleksandr R. Jan 21 '15 at 10:32
6
$\begingroup$

My guess is that when you write 1. and 2. as opposed to 1 and 2 you are telling Mathematica that the numbers you are using are not integers but that they are rather real numbers whose best decimal representation is, as far as you know or within your accuracy criteria, 1.0 and 2.0. Thus, numerical methods are legitimate and you can expect numerical errors.

It is for a similar reason that1+1. returns 2. (a real) rather than 2 (integer).

Mathematica is aware of the issue as testified by this example:

Simplify[(-1.)^(2 z), Element[z, Integers]]
(* 1. *)
Simplify[(-1)^(2 z), Element[z, Integers]]
(* 1 *)
$\endgroup$
0
$\begingroup$

(-1.)^2 does not return an imaginary number in Version 12.0; it returns 1.

$\endgroup$
  • 2
    $\begingroup$ I get an imaginary part: i.stack.imgur.com/F5UjW.png -- Note that the code to test is (-1.)^2., with a Real number 2., not an integer 2 as you have written. $\endgroup$ – Michael E2 May 12 at 12:21
  • $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is useful for learning how to format your questions and answers. You may also find this meta Q&A helpful $\endgroup$ – Michael E2 May 12 at 12:23
  • $\begingroup$ The difference is in the inputs (integer exponent 2 vs real 2.0). $\endgroup$ – Daniel Lichtblau May 12 at 20:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.