Wikipedia:Reference desk/Archives/Mathematics/2023 August 18

Welcome to the Wikipedia Mathematics Reference Desk Archives
The page you are currently viewing is a transcluded archive page. While you can leave answers for any questions shown below, please ask new questions on one of the current reference desk pages.

August 18

Code for complex square root

I'm trying to roll my own complex square root formula in a spreadsheet. (Actually, Calc has a built-in function IMSQRT which does this, but it takes strings as input and output and I'd rather use pairs of cells instead of trying to assemble and parse strings.) Our square root article gives a formula, but it has accuracy issues due to rounding when the real part is large compared to the imaginary part. For example when I use it to compute the square root of -10000000000+2i it gives 100000i rather than the more accurate value of .00001+100000i. Stack exchange gives some different formulas but they seem to have the rounding error issue as well. I was wondering if anyone knows how this is implemented in actual function libraries such as Python cmath. My own version isn't quite done, but it's turning out to be more complicated than I though it would be. --RDBury (talk) 06:36, 18 August 2023 (UTC)

I suspect this is the same problem as loss of accuracy when solving the quadratic equation using the familiar quadratic formula; see Quadratic formula § Muller's method. A better method is to use

$x_{1}={\frac {-b-\operatorname {sgn}(b){\sqrt {b^{2}-4ac}}}{2a}},$

$x_{2}={\frac {c}{ax_{1}}}.$

Solving $(x+iy)^{2}=(u+iv)$ for $x$ and $y$ gives rise, after expansion and elimination of $y$, to a quadratic equation in $x^{2}$:

$(x^{2})^{2}-u(x^{2})-{\frac {v^{2}}{4}}=0.$

Use the positive root of this equation to find $x={\sqrt {x^{2}}}$ and compute $y={\frac {v}{2x}}.$ (Disclaimer: I have not actually tried this out.)  --Lambiam 11:13, 18 August 2023 (UTC)

That's more or less the approach I'm taking. First compute $r={\sqrt {x^{2}+y^{2}}},$ then ${\sqrt {\frac {r+|x|}{2}}}.$ At this point u and v are either ±r, ±x/(2r) or 0. It's then a matter or working out which expression to use for which variable when, which sign to use, and when to use a special case to avoid division by 0. It's very easy to get the conjugate of a square root or a result that's not the principal branch. From a mathematician's point of view, you can reduce to the case where z is in the first quadrant by using ${\sqrt {-z}}=i{\sqrt {z}}$ and ${\sqrt {\bar {z}}}={\overline {\sqrt {z}}}$ when applicable. --RDBury (talk) 12:15, 18 August 2023 (UTC)

Here is Python-like pseudo code; test before use.

Input: (u, v)

Returns: (x, y) such that x + iy is the principal square root of u + iv

if v = 0:

if u >= 0:

return (sqrt(u), 0)

return (0, sqrt(−u))

if u = 0:

set x = sqrt(abs(v)/2)

return (x, sign(v) × x)

if u > 0:

set x = sqrt((u + sqrt(u² + v²)) / 2)

return (x, v / (2 × x))

if u < 0:

set y = sqrt((−u + sqrt(u² + v²)) / 2)

set x = v / (2 × y)

if x < 0:

return (−x, −y)

return (x, y)

 --Lambiam 14:35, 18 August 2023 (UTC)

In the code above, the cases u = 0 and u > 0 can be merged in an obvious way into one case u >= 0. A separate case for v = 0 would not be necessary except for the possibility that (u, v) = (0, 0), which would otherwise lead to the undefined division of 0/0. But if the case u = 0 remains, the separate case for v = 0 is not needed. (However, if in practice the supplied argument is often real, keeping it may improve the efficiency.)  --Lambiam 22:40, 18 August 2023 (UTC)

Lots of spreadsheets have an atan2 function. Presumably you could do r_out = sqrt(sqrt(r_in^2 + i_in^2)) * cos(atan2(i_in,r_in)/2) and i_out = sqrt(sqrt(r_in^2 + i_in^2)) * sin(atan2(i_in,r_in)/2). --Amble (talk) 17:09, 18 August 2023 (UTC)

Throwing trigonometry into basic vector operations is a recipe for significant slowdown and a pile of tricky edge cases. –jacobolus (t) 17:13, 18 August 2023 (UTC)

Since this is all going to live in a spreadsheet, I imagine that raw speed is not so important, but that it's handy to be able to write a formula for each output cell as a one-liner with no intermediate variables or conditionals. --Amble (talk) 17:23, 18 August 2023 (UTC)

It is possible to invoke user-defined functions written in Python from an Excel worksheet.^[1][2]  --Lambiam 22:18, 18 August 2023 (UTC)

Actually the original question was about how it's done in polished code libraries. I was able to do some "quick and dirty" formulas myself. --RDBury (talk) 13:08, 19 August 2023 (UTC)

The code implementing the Python function cmath.sqrt can be seen here, lines 739–805.  --Lambiam 14:23, 19 August 2023 (UTC)

Thanks. Except for a few functions I'm not familiar with (e.g. hypot) is looka a lot like your version. RDBury (talk) 16:27, 19 August 2023 (UTC)

Wikipedia has an article about it hypot. 😀 NadVolum (talk) 18:31, 19 August 2023 (UTC)