# How can I more efficiently define a square root with the branch cut along the positive imaginary axis?

I have a system of equations I need to solve involving square roots. I don't always have a good guess what the solutions are, and so I am fortunate that NSolve can handle these equations. I want to be able to try moving the branch cut of the square root function to the positive imaginary axis since I want continuous functions as I cross the positive and negative real axis. I can easily accomplish this as follows:

Define my new square root function as

MySqrt[z_] := If[Re[z] < 0 && Im[z] >= 0, -Sqrt[z], +Sqrt[z]]


Now for some equations I pull out of a hat:

EQNS = {x + MySqrt[y] - 2/MySqrt[z] + z - 1,
y - x + MySqrt[x] + 1/MySqrt[z] - 2,
MySqrt[x] + 1/MySqrt[y] - 2 y + z};
NSolve[EQNS == 0, {x, y, z}] // Timing


Yields:

  Out={0.92901,
{{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I},
{x -> 2.72138, y -> 1.42833, z -> 0.37027}}}


Now to check this against the standard definition:

EQNS = {x + Sqrt[y] - 2/Sqrt[z] + z - 1,
y - x + Sqrt[x] + 1/Sqrt[z] - 2,
Sqrt[x] + 1/Sqrt[y] - 2 y + z};
NSolve[EQNS == 0, {x, y, z}] // Timing


Yields:

Out={0.279246, {
{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I},
{x -> 0.698929 - 3.41798 I, y -> 1.03321 - 0.455899 I, z -> -0.301342 + 0.075144 I},
{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}}}


Which brings me to my question: Why does my definition for the square root make it so that it takes NSolve so much longer? It is about a factor of 4 here, which isn't a huge deal, but I am solving even more complicated equations with more variables and I am finding the increase to be a factor of 10 or greater. Can anyone suggest a better way to define a new (faster?) square root function? Or make any other suggestions to make this process faster?

Multiplying a number by $e^{\pi i/2} = i$ is equivalent to rotating the phase by $90^\circ$. So, if you pre-multiply your number by $i$ and take the square root, you will be using the same branch cut as the default. Then, to compensate, you will need to multiply by $e^{-\pi i/4}$ to rotate the phase backwards by $45^\circ$. Here is a function that does this:

isqrt[d_] := -Sqrt[I d](-1)^(3/4)


First, let's compare isqrt with MySqrt:

sample = RandomComplex[{-5 - 5 I, 5 + 5 I}, 100];
MinMax @ Abs[isqrt[sample] - MySqrt /@ sample]


{0., 6.28037*10^-16}

So, basically the same (I check the norm of the differences to account for differences in numerics). Next, let's check the timing:

EQNS={x+MySqrt[y]-2/MySqrt[z]+z-1,y-x+MySqrt[x]+1/MySqrt[z]-2,MySqrt[x]+1/MySqrt[y]-2 y+z};
NSolve[EQNS==0,{x,y,z}] //Timing

EQNS={x+isqrt[y]-2/isqrt[z]+z-1,y-x+isqrt[x]+1/isqrt[z]-2,isqrt[x]+1/isqrt[y]-2 y+z};
NSolve[EQNS==0,{x,y,z}] //Timing


{0.783682, {{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}, {x -> 2.72138, y -> 1.42833, z -> 0.37027}}}

{0.228611, {{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}, {x -> 2.72138, y -> 1.42833, z -> 0.37027}}}

The reason that MySqrt is so slow is because MySqrt[y] contains an "unevaluated" If object, while isqrt[y] is a simple arithmetic object:

MySqrt[y]
isqrt[y]


If[Re[y] < 0 && Im[y] >= 0, -Sqrt[y], +Sqrt[y]]

-(-1)^(3/4) Sqrt[I y]

NSolve probably has to use more complicated algebraic methods when presented with these If objects.