How to compute the square root of a complex number

Suppose you're given a complex number

z = x + iy

and you want to find a complex number

w = u + iv

such that w^2 = z. If all goes well, you can compute w as follows:

= (x^2 + y^2)
u = (( + x)/2)
v = sign(y) (( - x)/2).

For example, if z = 5 + 12i, then = 13, u = 3, and v = 2. A quick calculation confirms

(3 + 2i)^2 = 5 + 12i.

(That example worked out very nicely. More on why in the next post.)

Note that >= |x|, and so the square roots above have non-negative arguments.

Numerical issues

I said if all goes well" above because in floating point arithmetic, various things could go wrong. In [1], the authors point out that x^2 + y^2 could overflow. They give a more complicated algorithm that addresses the problem of overflow, and they reference a paper with an even more complicated algorithm to cover all the corner cases - NaNs, infinities, etc.

If you have access to a hypot function, it will let you compute (x^2 + y^2) without overflow. If you don't have access to a hypot function, you can roll your own quickly using the algorithm here.

This may seem unnecessary. Why would you need to worry about overflow? And why would you have access to a real square root function but not a complex square root function?

If you're using standard 64-bit floating point numbers, the hypotenuse calculation is not going to overflow as long as your arguments are less than 10¹⁵⁴, and you can use the function csqrt from complex.h.

But if you're using a 16-bit float or even an 8-bit float, overflow is an ever-present concern, and you might not have the math library support you're accustomed to.

Related posts

• Anatomy of a floating point number
• Anatomy of a posit number
• Math library functions that seem unnecessary

[1] Jean-Michel Muller et al. Handbook of Floating Point Arithmetic.