19

I have run across some confusing behaviour with square roots of complex numbers in python. Running this code:

from cmath import sqrt
a = 0.2
b = 0.2 + 0j
print(sqrt(a / (a - 1)))
print(sqrt(b / (b - 1)))

gives the output

0.5j
-0.5j

A similar thing happens with

print(sqrt(-1 * b))
print(sqrt(-b))

It appears these pairs of statements should give the same answer?

CleptoMarcus
  • 193
  • 1
  • 1
  • 5
  • According to [Wolfram](https://www.wolframalpha.com/) you are correct. The first pair ([link](https://www.wolframalpha.com/input/?i=sqrt(0.2+%2F+(0.2+-+1))) and [link](https://www.wolframalpha.com/input/?i=sqrt(+(0.2%2B0i)+%2F+(+(0.2%2B0i)+-1+)))) both should be `0.5i`, and the second pair ([link](https://www.wolframalpha.com/input/?i=sqrt(-1+*+(0.2+%2B+0i))) and [link](https://www.wolframalpha.com/input/?i=sqrt(-1+*+(0.2+%2B+0i)))) should both be `0.447214... i`. The source for `cmath.sqrt()` is [here](https://hg.python.org/cpython/file/tip/Modules/cmathmodule.c#l732)... – Jens Apr 12 '16 at 21:42
  • 7
    Both answers are correct, the question is why it returns different conjugates. – tzaman Apr 12 '16 at 21:45
  • FWIW the behavior appears the same in 2.7 and 3.5. – tzaman Apr 12 '16 at 21:49
  • 2
    There are simply multiple solutions to a complex root. – roadrunner66 Apr 12 '16 at 21:49
  • 3
    @tzaman ...And if it's defined somewhere which one Python should return. If it's not defined, Python has right to choose any. – George Sovetov Apr 12 '16 at 21:51

2 Answers2

12

Both answers (+0.5j and -0.5j) are correct, since they are complex conjugates -- i.e. the real part is identical, and the imaginary part is sign-flipped.

Looking at the code makes the behavior clear - the imaginary part of the result always has the same sign as the imaginary part of the input, as seen in lines 790 and 793:

r.imag = copysign(d, z.imag);

Since a/(a-1) is 0.25 which is implicitly 0.25+0j you get a positive result; b/(b-1) produces 0.25-0j (for some reason; not sure why it doesn't result in 0.25+0j tbh) so your result is similarly negative.

EDIT: This question has some useful discussion on the same issue.

Community
  • 1
  • 1
tzaman
  • 46,925
  • 11
  • 90
  • 115
  • `a / (a - 1)` is `-0.25`. Converting it to complex gives you `-0.25+0j`. With b, you get `-0.25-0j`. – Seth Apr 12 '16 at 22:08
  • 1
    The [documentation](https://docs.python.org/2/library/cmath.html#cmath.sqrt) also specifies that the branch cut is chosen to be along a ray from 0 to -∞, which would give the observed behavior. It's interesting that `b/(b-1)` gives a negative zero for the imaginary component, though. – user2357112 Apr 12 '16 at 22:19
  • @user2357112 yeah, that's pretty weird, I'm not sure why it does that. Simply typing `1-0j` in as a literal into IDLE produces `1+0j` as a value. Interestingly, `-0j` by itself produces `-0-0j` though. – tzaman Apr 12 '16 at 22:23
  • 3
    @tzaman: That's a nasty edge case in the syntax. A complex number with a negative zero imaginary component gets `repr`esented as `0.2-0j`, but the expression `0.2-0j` is actually a subtraction of `0.2` and `0j`. The subtraction produces a *positive* zero imaginary component, because `0j - 0j` is *positive* `0j`. It's [come up](http://bugs.python.org/issue25839) a few times on the bug tracker. There's a related problem where negative zero real components aren't visible in the `repr` output at all. – user2357112 Apr 12 '16 at 22:32
  • @user2357112 that's fascinating, thanks for the link! – tzaman Apr 12 '16 at 22:35
2

I can answer why this is happening, but not why the behavior was chosen.

a/(a - 1)

evaluates to 0.2/-0.8 which is -0.25, which is converted to a complex number by cmath.sqrt, while

b/(b - 1)

evaluates to (0.2+0j)/(-0.8+0j) which is (-0.25-0j), which is converted to a complex number with a negative complex component.

For a simpler example,

cmath.sqrt(0j) == 0j
cmath.sqrt(-0j) == -0j
Matt Jordan
  • 2,133
  • 9
  • 10