sqrt

Preface

sqrt is short for "square root".

All complex numbers have two square roots (except 0, which just has 0).

√(1) = 1, -1

√(4) = 2, -2

√(-1) = i, -i

√(2i) = 1+i, -1-i

√(3+4i) = 2+i, -2-i

Both square roots will always be arithmetic opposites. That is, one square root will always be the other times -1.

For the sake of simplicity, mathematicians have defined something called the "principal square root", which only yields one of the two square roots. When dealing strictly with positive reals, the principal square root is always whichever of the two roots is positive.

√(1) = 1

√(4) = 2

When dealing with complex numbers, there technically isn't one preferred convention, but the most common one to go by (as well as the one this library is designed under) is to choose whichever of the two has a positive real part.

√(2i) = 1+i

√(3+4i) = 2+i

In the special case that our input is negative, both square roots will be imaginary and have 0 as their real part. In this circumstance, it is most common to pick whichever one has a positive imaginary part and assign that as the principal square root.

√(-1) = i

√(-4) = 2i

Algorithm

Thanks to Euler's identity, Complex Numbers can be written in rectangular form:

1+i

Or in polar notation:

√(2) * e^(i*π/4)

If you put any number:

r * e^(Θi)

To the power of a real exponent k, your result will be:

r^k * e^(k*Θi)

Thus, if we wanted to take the square root of x+yi, we would first convert to polar, raise that to the power of 0.5, and convert back to rectangular

x+yi = √(x²+y²) * e^(atan2(y,x)*i)

√(x+yi) = (x²+y²)^(1/4) * e^(atan2(y,x)/2 * i)

Convert back to rectangular:

(x²+y²)^(1/4) * (cos(atan2(y,x)/2) + i*sin(atan2(y,x)/2))

Thanks to our half angle formulas, we get

cos(atan2(y,x)/2) = cos(sgn(y)*acos(x/√(x²+y²))/2) = √(( x/√(x²+y²) + 1)/2) = √((√(x²+y²+x))/2) / (x²+y²)^0.25

sin(atan2(y,x)/2) = sin(sgn(y)*acos(x/√(x²+y²))/2) = sgn(y) * √(( 1 - x/√(x²+y²) )/2) = sgn(y) * √((√(x²+y²)-x)/2) / (x²+y²)^0.25

(Where sgn returns 1 if the input is positive, -1 if the input is negative. While sgn(0) is technically 0, for the sake of this calculation, as well as the convention that √(-1) = +i, it's best to assume here that sgn(0)=1)

√(x+yi) = (x²+y²)^(1/4) * (cos(atan2(y,x)/2) + i*sin(atan2(y,x)/2))

√(x+yi) = (x²+y²)^(1/4) * (√((√(x²+y²)+x)/2) / (x²+y²)^0.25 + i*sgn(y) * √((√(x²+y²)-x)/2) / (x²+y²)^0.25)

√(x+yi) = ( √((√(x²+y²)+x)/2) + i*sgn(y)* √((√(x²+y²)-x)/2) )

√(x²+y²) = |x+yi|

√(x+yi) = √((|x+yi|+x)/2) + isgn(y)√((|x+yi|-x)/2)

That's most of the algorithm. However, there are a few slight technicalities and optimizations. For one thing, if you multiply together the real and imaginary parts of √(x+yi), you are guaranteed to get y/2.

PROOF: Let's say u and v are the real and imaginary parts of √(x+yi) respectively

√(x+yi)=u+vi

√(x+yi)²=(u+vi)²=u²+2uvi-v²

x+yi = (u²-v²) + (2uv)i

y=2uv

y/2=uv

Since division is much less computationally expensive than square roots, the Complex Number Library takes advantage of the above fact by only computing u, the solving v by dividing y/(2*u).

This speeds up calculations considerably. Whereas before, you'd be performing 3 square roots (one of which is used to find the absolute value), this way, you only need 2 square roots, albeit at the cost of one additional division.

However, there is another technicality that must be considered: roundoff.

Suppose you want to compute √(-1E17+0.01i)

Using the given algorithm, we'd get an absolute value of about 1.0000000000000000000...0005E17. More precisely, we'd get a value of about 1E17+5E-22. Under the double floating point system, that will only register as 1E17.

Using the algorithm, we in theory get a real part of √((1E17+5E-22+(-1E17))/2) = √(5E-22/2) = 1.58113883E-11. Using the division formula, we get an imaginary part of 1E17/(2*1.58113883E-11) = 3.16227766E27. And, thus, we'd get our correct solution of 1.581133883E-11+3.16227766E27*i.

Of course, that's not what the computer actually does. Like I said, the absolute value will only be stored as 1E17. The real part will be computed as √((1E17+(-1E17))/2) = 0, and the imaginary part will be computed as 1E17/(2*0) = ∞. And we get an incorrect answer of 0+∞i.

How do we solve this? Well, for starters, we have to address that this can only ever be an issue if the real part is negative. I'd explain why, but I honestly think it'd be more effective to look back at the calculations and analyze what would happen if the real part was positive. Remember, the absolute value is always positive!

With that said, there's a simple fix. Whenever the real part is negative, we ask the computer, "hey, instead of finding u, then solving for v with division, could you compute v, then solve for u with division?". Watch how that effects our previous example:

v = sgn(0.01)*√((|-1E17+0.01i|-(-1E17))/2) = 1*√((1E17+1E17)/2) = √(1E17) = 3.16227766E27

u = y/(2*v) = 0.01/(2*3.16227766E27) = 1.58113883E-11.

Thus, we get the correct solution of 1.58113883E-11+3.16227766E27i

Now, you might ask "Even if this is faster, wouldn't it still be safer to apply the original algorithm and avoid roundoff altogether?" And to that I must say "...no".

If we were to apply the original formula, where we solve both u and v independently of each other, we'd get

√(-1E17+0.01i) = 0 + 3.16227766E27i

The imaginary component is correct, but the real component isn't. Thus, our optimization not only improves efficiency, but actually removes roundoff error. It makes sure that, of the two components, we pick whichever one is guaranteed never to cancel out and create and underflow, then solve for the other component as a function of that. It's the best of both worlds.

What about overflow?

The square root doesn't have much room for overflow error. There are very few cases where our input isn't infinite but our output still overflows. But, it can still happen. If |x+yi|+|x| surpasses the overflow limit, even though x+yi doesn't, then that's a special case. Our formula will fail, since even though (|x+yi|±x)/2 won't overflow, |x+yi|±x will. It's that crucial divide by two step that trips it up. How do we solve this? Well...lazily. We just divide our number by 4, take the square root, then multiply by 2. Simple...but effective.

In the case that the input is an overflow, the output can also be expected to overflow. Admittedly, this library has very, very poor support for special cases involving infinity. However, for square roots, all the cases are covered.

√(∞+yi)=∞

√(-∞+yi)=sgn(y)*∞i

√(x+∞i)=∞+∞i

√(x-∞i)=∞-∞i