Suppose you’re given a complex number

*z* = *x* + *i**y*

and you want to find a complex number

*w* = *u* + *iv*

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

ℓ = √(*x*² + *y*²)

*u* = √((ℓ + *x*)/2)

*v* = sign(*y*) √((ℓ − *x*)/2).

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

(3 + 2*i*)² = 5 + 12*i*.

(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*² + *y*² 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*² + *y*²) 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^{154}, 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.