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

Improve this question
$\endgroup$
14
  • 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
    Commented Jan 13, 2015 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
    Commented Jan 13, 2015 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
    Commented Jan 13, 2015 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
    Commented Jan 13, 2015 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.
    Commented Jan 20, 2015 at 21:45

4 Answers 4

Reset to default
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 )

Improve this answer
$\endgroup$
3
  • $\begingroup$ Which compiler/flags? $\endgroup$
    – Oleksandr R.
    Commented Jan 20, 2015 at 22:30
  • $\begingroup$ no flags. See the edit though - the results may differ depending on how the compiler optimizes the code. $\endgroup$
    – george2079
    Commented Jan 20, 2015 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.
    Commented Jan 20, 2015 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
Improve this answer
$\endgroup$
4
  • 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.
    Commented Jan 20, 2015 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
    Commented Jan 21, 2015 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.
    Commented Jan 21, 2015 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.
    Commented Jan 21, 2015 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 *)
Improve this answer
$\endgroup$
0
$\begingroup$

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

Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ I get an imaginary part: i.sstatic.net/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
    Commented May 12, 2019 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
    Commented May 12, 2019 at 12:23
  • $\begingroup$ The difference is in the inputs (integer exponent 2 vs real 2.0). $\endgroup$
    – Daniel Lichtblau
    Commented May 12, 2019 at 20:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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