-
Notifications
You must be signed in to change notification settings - Fork 0
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
All complex numbers can be written in rectangular form:
1+i
Or in polar notation:
√(2) ∠ π/4
If you put any number
r ∠ θ
To the power of a real exponent k, your result will be:
r^k ∠ k*θ
Thus, to take the square root of x+yi, first convert to polar, raise that to the power of 1/2, then convert that back to rectangular.
Convert to polar:
x+yi = √(x²+y²) ∠ atan2(y,x)
Raise to the power of 1/2:
√(x+yi) = ∜(x²+y²) ∠ atan2(y,x)/2
Convert back to rectangular:
√(x+yi) = ∜(x²+y²) * (cos(atan2(y,x)/2) + i*sin(atan2(y,x)/2))
According to half-angle formulas,
cos(θ/2) = √((1+cosθ)/2)
sin(θ/2) = csgn(θ)*√((1-cosθ)/2)
(Assuming -π < θ ≤ π, with csgn(θ) being 1 if θ ≥ 0 and -1 if θ < 0)
And since
cos(atan2(y,x)) = x/√(x²+y²)
sin(atan2(y,x)) = y/√(x²+y²)
that means:
cos(atan2(y,x)/2) = √(( 1 + x/√(x²+y²) )/2) = √((√(x²+y²)+x)/2) / ∜(x²+y²)
sin(atan2(y,x)/2) = csgn(atan2(y,x))*√(( 1 - x/√(x²+y²) )/2) = csgn(y)*√((√(x²+y²)-x)/2) / ∜(x²+y²)
(Feel free to ponder on the csgn part to make sure it makes sense. Keep in mind, -π ≤ atan2(y,x) < π)
√(x+yi) = ∜(x²+y²) * (cos(atan2(y,x)/2) + i*sin(atan2(y,x)/2))
√(x+yi) = ∜(x²+y²) * (√((√(x²+y²)+x)/2) / ∜(x²+y²) + i * csgn(y) * √((√(x²+y²)-x)/2) / ∜(x²+y²))
√(x+yi) = ( √((√(x²+y²)+x)/2) + i * csgn(y) * √((√(x²+y²)-x)/2) )
√(x²+y²) = |x+yi|
√(x+yi) = √((|x+yi|+x)/2) + csgn(y)i√((|x+yi|-x)/2)
That's most of the algorithm. In case you're not quite convinced this works, a proof is shown at the bottom of this section (right before "What about overflow?") that this expression, when squared, always gives you x+yi.
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, then 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: round-off.
Suppose you want to compute √(-1E17+0.01i)
Using the given algorithm, we'd get an absolute value of slightly less than 1.0000000000000000000...0005E17, according to Wolfram Alpha. More concisely, we'd get a value of about 1E17+5E-22. Under the double floating point system, however, 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. That's because u is calculated by adding the absolute value to the real part (then doing some other stuff). Our issue arises because we added a large positive number to a large negative number, resulting in underflow. The same issue cannot happen if we add two positive numbers, and since the absolute value is always positive, this means this issue can only happen if the real part is negative.
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 find v, then solve for u with division?". Watch how that affects our previous example:
v = csgn(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 round-off 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 round-off 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.
As stated, the square root can be calculated as follows:
√(x+yi) = √((|x+yi|+x)/2) + csgn(y)i√((|x+yi|-x)/2)
If we square both sides, we get:
x+yi = (√((|x+yi|+x)/2) + csgn(y)i√((|x+yi|-x)/2))²
Thus, if we simplify the right hand side, we should get x+yi. But first, a small lemma:
(a+b)² = a² + 2ab + b²
Using this lemma, we can evaluate the expression squared:
(√((|x+yi|+x)/2) + csgn(y)i√((|x+yi|-x)/2))² =
(√((|x+yi|+x)/2))² + 2 (√((|x+yi|+x)/2)) (csgn(y)i√((|x+yi|-x)/2)) + (csgn(y)i√((|x+yi|-x)/2))² =
(|x+yi|+x)/2 + 2csgn(y)i * √((|x+yi|+x)/2) * √((|x+yi|-x)/2) + csgn(y)²i²(|x+yi|-x)/2 =
(|x+yi|+x)/2 + 2csgn(y)i * √((|x+yi|+x)/2 * (|x+yi|-x)/2) + 1*-1*(|x+yi|-x)/2 =
(|x+yi|+x)/2 + 2csgn(y)i * √((|x+yi|+x)(|x+yi|-x)/4) - (|x+yi|-x)/2 =
(|x+yi|+x)/2 - (|x+yi|-x)/2 + 2csgn(y)i * √((|x+yi|+x)(|x+yi|-x))/2 =
(|x+yi|+x-|x+yi|+x)/2 + csgn(y)i * √(|x+yi|²-x²) =
(2x)/2 + csgn(y)i * √(x²+y²-x²) =
x + csgn(y)i * √(y²) =
x + csgn(y)i * |y|
If y ≥ 0, this evaluates to x + 1i*y = x+yi, which is the correct answer.
If y < 0, this evaluates to x + (-1)i*(-y) = x+yi, which is the correct answer.
Furthermore, just to quickly prove that this is always the principal square root:
The real part is the square root of a non-negative real number (since |x+yi| ≥ |x|). The square root of a non-negative real number is always positive or zero. If it's positive, this is the principal square root, because the real part is positive. If it's zero, however, that means |x+yi|+x is 0. That means x = -|x+yi|, or x = -√(x²+y²). That means x is negative or zero. If we square both sides, x²=x²+y², thus y=0. Thus this case only occurs when x is non-positive and y is 0. In that case, csgn(y) = 1, and the imaginary part is 1 times the square root of a non-negative real, which in turn is a non-negative real. Thus, either the real part is positive, or it's 0 and the imaginary part is non-negative.
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)=csgn(y)*∞i
√(x+∞i)=∞+∞i
√(x-∞i)=∞-∞i