TAOCP 4.2.1 Exercise 16
We wish to compute the real and imaginary parts of $\frac{a + b i}{c + d i}$ without forming $c^2 + d^2$, which may overflow even if $|c|$ or $|d|$ is near the largest representable floating point num...
Section 4.2.1: Single-Precision Calculations
Exercise 16. [HM21] (Robert L. Smith.) Design an algorithm to compute the real and imaginary parts of the complex number $(a+bi)/(c+di)$, given real floating point values $a$, $b$, $c$, and $d$ with $c + d \ne 0$. Avoid the computation of $c^2 + d^2$, since it would cause floating point overflow even when $|c|$ or $|d|$ is approximately the square root of the maximum allowable floating point value.
Verified: yes
Solve time: 2m02s
Solution
We wish to compute the real and imaginary parts of
$\frac{a + b i}{c + d i}$
without forming $c^2 + d^2$, which may overflow even if $|c|$ or $|d|$ is near the largest representable floating point number. Let $a$, $b$, $c$, $d$ be normalized floating point numbers, with $c + d \ne 0$.
Begin by writing the usual formula in terms of multiplication and division:
$\frac{a + b i}{c + d i} = \frac{(a + b i)(c - d i)}{c^2 + d^2} = \frac{ac + bd}{c^2 + d^2} + i \frac{bc - ad}{c^2 + d^2}.$
Direct evaluation of $c^2 + d^2$ can overflow. To avoid this, factor the larger of $|c|$ and $|d|$ from numerator and denominator. Assume, without loss of generality, that $|c| \ge |d|$. Set
$r = \frac{d}{c}.$
Then $|r| \le 1$, and
$c^2 + d^2 = c^2 (1 + r^2),$
so that
$\frac{1}{c^2 + d^2} = \frac{1}{c^2 (1 + r^2)} = \frac{1}{c^2} \cdot \frac{1}{1 + r^2}.$
The real and imaginary parts can then be expressed as
\begin{align*}
\operatorname{Re}\frac{a + b i}{c + d i} &= \frac{ac + bd}{c^2 + d^2} = \frac{a + b r}{c (1 + r^2)},\
\operatorname{Im}\frac{a + b i}{c + d i} &= \frac{bc - ad}{c^2 + d^2} = \frac{b - a r}{c (1 + r^2)}.
\end{align*}
If instead $|d| > |c|$, set
$r = \frac{c}{d},$
so that $|r| \le 1$, and
$c^2 + d^2 = d^2 (1 + r^2),$
yielding
\begin{align*}
\operatorname{Re}\frac{a + b i}{c + d i} &= \frac{ac + bd}{c^2 + d^2} = \frac{b + a r}{d (1 + r^2)},\
\operatorname{Im}\frac{a + b i}{c + d i} &= \frac{bc - ad}{c^2 + d^2} = \frac{-a + b r}{d (1 + r^2)}.
\end{align*}
These formulas avoid any intermediate quantity exceeding the magnitude of $|c|$ or $|d|$, since $|r| \le 1$ and $|1 + r^2| \le 2$, thus preventing overflow.
The algorithm can now be stated in stepwise form:
Step 1. Test whether $|c| \ge |d|$.
Step 2. If $|c| \ge |d|$, compute $r \leftarrow d/c$, then
$\operatorname{Re} \leftarrow \frac{a + b r}{c (1 + r^2)}, \quad \operatorname{Im} \leftarrow \frac{b - a r}{c (1 + r^2)}.$
Step 3. If $|d| > |c|$, compute $r \leftarrow c/d$, then
$\operatorname{Re} \leftarrow \frac{b + a r}{d (1 + r^2)}, \quad \operatorname{Im} \leftarrow \frac{-a + b r}{d (1 + r^2)}.$
All divisions are by numbers of magnitude at least $\max(|c|,|d|)$, which is nonzero by hypothesis, and all multiplications and additions involve numbers of magnitude at most $2 \max(|a|, |b|)$ times a factor $\le 2$, ensuring that no overflow occurs unless the final result itself cannot be represented. Each step preserves the correct ratio of numerator to denominator.
This completes the proof. ∎