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) + isin(atan2(y,x)/2)) √(x+yi) = (x²+y²)^(1/4) * (√((√(x²+y²)+x)/2) / (x²+y²)^0.25 + isgn(y) * √((√(x²+y²)-x)/2) / (x²+y²)^0.25) √(x+yi) = ( √((√(x²+y²)+x)/2) + isgn(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 that 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 x is a large negative number