Repository navigation
sqrt
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, some authors go by different conventions. However, by far 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, the most common convention is to pick whichever one has a positive imaginary part and assign that as the principal square root.
√(-1) = i
√(-4) = 2i
Thanks to Euler's formula, 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))
According to tangent half-angle formulas,
cos(θ/2) = √((1+cos(θ))/2)
sin(θ/2) = sgn(sin(θ))*√((1-cos(θ))/2)
And since
cos(atan2(y,x)) = x/√(x²+y²)
sin(atan2(y,x)) = y/√(x²+y²)
that means:
cos(atan2(y,x)/2) = √(( x/√(x²+y²) + 1)/2) = √((√(x²+y²+x))/2) / (x²+y²)^0.25
sin(atan2(y,x)/2) = sgn(y/√(x²+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. This should hopefully make sense, since last time we tried to find the real part with square roots, we got the wrong answer of 0. 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 an underflow, then solve for the other component as a function of that. It's the best of both worlds.
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 isn't an overflow, |x+yi|±x is. 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 back 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, most of the cases are covered pretty well.
√(∞+yi)=∞
√(-∞+yi)=sgn(y)*∞i
√(x+∞i)=∞+∞i
√(x-∞i)=∞-∞i