MATHEMATICS OF COMPUTATION Volume 00, Number 0, Pages 000–000 S 0025-5718(XX)0000-0

ERROR BOUNDS ON COMPLEX FLOATING-POINT MULTIPLICATION

RICHARD BRENT, COLIN PERCIVAL, AND PAUL ZIMMERMANN

In memory of Erin Brent (1947–2005)

Abstract. Given floating-point arithmetic with t-digit base-β significands in which all arithmetic operations are performed as if calculated to infinite precision and rounded to a nearest representable value, we prove that the product of complex values z_{0} and z_{1} can be computed with maximum absolute error |z_{0}| |z_{1}| \frac{1}{2} β^{1−t}\sqrt{5}. In particular, this provides relative error bounds of 2^{−24}\sqrt{5} and 2^{−53}\sqrt{5} for IEEE 754 single and double precision arithmetic respectively, provided that overflow, underflow, and denormals do not occur.

We also provide the numerical worst cases for IEEE 754 single and double precision arithmetic.

1. Introduction

In an earlier paper [2], the second author made the claim that the maximum relative error which can occur when computing the product z_{0}z_{1} of two complex values using floating-point arithmetic is ϵ\sqrt{5}, where ϵ is the maximum relative error which can result from rounded floating-point addition, subtraction, or multiplication. While reviewing that paper a few years later, the other two authors noted that the proof given was incorrect, although the result claimed was true.

Since the bound of ϵ\sqrt{8} which is commonly used [1] is suboptimal, we present here a corrected proof of the tighter bound. Interestingly, by explicitly finding worst-case inputs, we can demonstrate that our error bound is effectively optimal.

Throughout this paper, we concern ourselves with floating-point arithmetic with t-digit base-β significands, denote by ulp(x) for x ≠ 0 the (unique) power of β such that β^{t−1} ≤ |x| /ulp(x) < β^{t}, and write ϵ = \frac{1}{2}ulp(1) = \frac{1}{2}β^{1−t}; we also define ulp(0) = 0. The notations x ⊕ y, x ⊖ y, and x ⊗ y represent rounded-to-nearest floating-point addition, subtraction, and multiplication of the values x and y.

2. An Error Bound

Theorem 1. Let z_{0} = a_{0} + b_{0}i and z_{1} = a_{1} + b_{1}i, with a_{0}, b_{0}, a_{1}, b_{1} floatingpoint values with t-digit base-β significands, and let z_{2} = ((a_{0} ⊗ a_{1}) ⊖ (b_{0} ⊗ b_{1})) + ((a_{0} ⊗ b_{1}) ⊕ (b_{0} ⊗ a_{1}))i be computed. Providing that no overflow or underflow occur, no denormal values are produced, arithmetic results are correctly rounded to a nearest representable value, z_{0}z_{1} ≠ 0, and ϵ ≤ 2^{−5}, the relative error

|z_{2}(z_{0}z_{1})^{−1} − 1|

is less than ϵ\sqrt{5} = \frac{1}{2}β^{1−t}\sqrt{5}.

\begin{matrix}c \\ ◯ \end{matrix}1997 American Mathematical Society

Proof. Let a_{0}, b_{0}, a_{1}, and b_{1} be chosen such that the relative error is maximized. By multiplying z_{0} and z_{1} by powers of i and/or taking complex conjugates, we can assume without loss of generality that

(1)  0 ≤ a_{0}, b_{0}, a_{1}, b_{1}
(2)  b_{0}b_{1} ≤ a_{0}a_{1}

and given our assumptions that overflow, underflow, and denormals do not occur, and that rounding is performed to a nearest representable value, we can conclude that for any x occurring in the computation, the error introduced when rounding x is at most \frac{1}{2}ulp(x) and is strictly less than ϵ · x.

We note that the error |ℑ(z_{2} − z_{0}z_{1})| in the imaginary part of z_{2} is bounded as follows:

|ℑ(z_{2} − z_{0}z_{1})| = |((a_{0} ⊗ b_{1}) ⊕ (b_{0} ⊗ a_{1})) − (a_{0}b_{1} + b_{0}a_{1})|
≤ |a_{0} ⊗ b_{1} − a_{0}b_{1}| + |b_{0} ⊗ a_{1} − b_{0}a_{1}|
+ |((a_{0} ⊗ b_{1}) ⊕ (b_{0} ⊗ a_{1})) − (a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1})|

and consider two cases:

Case I1: ulp(a_{0}b_{1} + b_{0}a_{1}) < ulp(a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1})

Using first the definition of ulp and second the assumption above, we must have

a_{0}b_{1} + b_{0}a_{1} < β^{t}ulp(a_{0}b_{1} + b_{0}a_{1}) ≤ a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1}

and therefore

|(a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1}) − β^{t}ulp(a_{0}b_{1} + b_{0}a_{1})| < (a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1}) − (a_{0}b_{1} + b_{0}a_{1})
≤ |a_{0} ⊗ b_{1} − a_{0}b_{1}| + |b_{0} ⊗ a_{1} − b_{0}a_{1}|
≤ ϵ · (a_{0}b_{1} + b_{0}a_{1}).

However, β^{t}ulp(a_{0}b_{1} + b_{0}a_{1}) is a representable floating-point value; so given our assumption that rounding is performed to a nearest representable value, we must now have

|((a_{0} ⊗ b_{1}) ⊕ (b_{0} ⊗ a_{1})) − (a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1})| < ϵ · (a_{0}b_{1} + b_{0}a_{1}).

Case I2: ulp(a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1}) ≤ ulp(a_{0}b_{1} + b_{0}a_{1})

From our assumption that the results of arithmetic operations are correctly rounded, we obtain

|((a_{0} ⊗ b_{1}) ⊕ (b_{0} ⊗ a_{1})) − (a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1})| ≤ \frac{1}{2}ulp(a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1})

≤ \frac{1}{2}ulp(a_{0}b_{1} + b_{0}a_{1})

≤ ϵ · (a_{0}b_{1} + b_{0}a_{1}).

Combining these two cases with the earlier-stated bound, we obtain

|ℑ(z_{2} − z_{0}z_{1})| ≤ |a_{0} ⊗ b_{1} − a_{0}b_{1}| + |b_{0} ⊗ a_{1} − b_{0}a_{1}|
+ |((a_{0} ⊗ b_{1}) ⊕ (b_{0} ⊗ a_{1})) − (a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1})|
< ϵ · (a_{0}b_{1}) + ϵ · (b_{0}a_{1}) + ϵ · (a_{0}b_{1} + b_{0}a_{1})
= ϵ · (2a_{0}b_{1} + 2b_{0}a_{1}).

Now that we have a bound on the imaginary part of the error, we turn our attention to the real part, and consider the following four cases (where the examples given apply to β = 2):

ulp(b_{0}b_{1}) ≤ ulp(a_{0}a_{1}) ≤ ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1})  e.g., z_{0} = z_{1} = 0.8 + 0.1i
ulp(b_{0}b_{1}) < ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1}) < ulp(a_{0}a_{1})  e.g., z_{0} = z_{1} = 0.8 + 0.4i
ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1}) ≤ ulp(b_{0}b_{1}) < ulp(a_{0}a_{1})  e.g., z_{0} = z_{1} = 0.8 + 0.7i
ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1}) < ulp(b_{0}b_{1}) = ulp(a_{0}a_{1})  e.g., z_{0} = z_{1} = 0.8 + 0.8i.

Since we have assumed that b_{0}b_{1} ≤ a_{0}a_{1}, we know that ulp(b_{0}b_{1}) ≤ ulp(a_{0}a_{1}), and thus these four cases cover all possible inputs. Consequently, it suffices to prove the required bound for each of these four cases.

Case R1: ulp(b_{0}b_{1}) ≤ ulp(a_{0}a_{1}) ≤ ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1})

Note that the right inequality can only be strict if a_{0} ⊗ a_{1} rounds up to a power of β and b_{0}b_{1} = 0.

We observe that

a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1} < a_{0}a_{1} − b_{0}b_{1} + ϵ · (a_{0}a_{1} + b_{0}b_{1})

and bound the real part of the complex error as follows:

|ℜ(z_{2} − z_{0}z_{1})| ≤ |a_{0} ⊗ a_{1} − a_{0}a_{1}| + |b_{0} ⊗ b_{1} − b_{0}b_{1}|
+ |((a_{0} ⊗ a_{1}) ⊖ (b_{0} ⊗ b_{1})) − (a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1})|

≤ \frac{1}{2}ulp(a_{0}a_{1}) + \frac{1}{2}ulp(b_{0}b_{1}) + \frac{1}{2}ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1})

≤ \frac{1}{2}ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1}) + \frac{1}{2}ulp(b_{0}b_{1}) + \frac{1}{2}ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1})

< 2ϵ · (a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1}) + ϵ · (b_{0}b_{1})

< ϵ · (2a_{0}a_{1} − b_{0}b_{1}) + ϵ^{2} · (2a_{0}a_{1} + 2b_{0}b_{1}).

Applying the triangle inequality, we now observe that

|z_{2} − z_{0}z_{1}| = \sqrt{ℜ(z_{2} − z_{0}z_{1})^{2} + ℑ(z_{2} − z_{0}z_{1})^{2}}

< ϵ\sqrt{(2a_{0}a_{1} − b_{0}b_{1})^{2} + (2a_{0}b_{1} + 2b_{0}a_{1})^{2}} + ϵ^{2} · (2a_{0}a_{1} + 2b_{0}b_{1})

≤ ϵ\sqrt{\frac{32}{7} |z_{0}z_{1}|^{2} − \frac{4}{7}(a_{0}b_{1} − b_{0}a_{1})^{2} − \frac{1}{7}(2a_{0}a_{1} − 5b_{0}b_{1})^{2}} + 2ϵ^{2} |z_{0}z_{1}|

≤ ϵ (\sqrt{32/7} + 2ϵ) |z_{0}z_{1}| < ϵ\sqrt{5} |z_{0}z_{1}|

as required.

Case R2: ulp(b_{0}b_{1}) < ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1}) < ulp(a_{0}a_{1})

Noting that ulp(x) < ulp(y) implies ulp(x) ≤ β^{−1}ulp(y) ≤ \frac{1}{2}ulp(y), we obtain

|ℜ(z_{2} − z_{0}z_{1})| ≤ \frac{1}{2}ulp(a_{0}a_{1}) + \frac{1}{2}ulp(b_{0}b_{1}) + \frac{1}{2}ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1})

≤ \frac{7}{8}ulp(a_{0}a_{1})

≤ ϵ · (\frac{7}{4}a_{0}a_{1})

and therefore

|z_{2} − z_{0}z_{1}| = \sqrt{ℜ(z_{2} − z_{0}z_{1})^{2} + ℑ(z_{2} − z_{0}z_{1})^{2}}

< ϵ\sqrt{(\frac{7}{4}a_{0}a_{1})^{2} + (2a_{0}b_{1} + 2b_{0}a_{1})^{2}}

= ϵ\sqrt{\frac{1024}{207} |z_{0}z_{1}|^{2} − \frac{196}{207}(a_{0}b_{1} − b_{0}a_{1})^{2} − \frac{1}{3312}(79a_{0}a_{1} − 128b_{0}b_{1})^{2}}

≤ ϵ\sqrt{1024/207} |z_{0}z_{1}| < ϵ\sqrt{5} |z_{0}z_{1}|

as required.

Case R3: ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1}) ≤ ulp(b_{0}b_{1}) < ulp(a_{0}a_{1})

In this case, there is no rounding error introduced in computing the difference between a_{0} ⊗ a_{1} and b_{0} ⊗ b_{1} since ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1}) ≤ ulp(b_{0}b_{1}) ≤ ulp(b_{0} ⊗ b_{1}) and ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1}) < ulp(a_{0}a_{1}) ≤ ulp(a_{0} ⊗ a_{1}). Also,

ulp(b_{0}b_{1}) ≤ \frac{1}{β} ulp(a_{0}a_{1}) ≤ \frac{1}{2}ulp(a_{0}a_{1})

so we have

|ℜ(z_{2} − z_{0}z_{1})| ≤ \frac{1}{2}ulp(a_{0}a_{1}) + \frac{1}{2}ulp(b_{0}b_{1})

≤ \frac{3}{4}ulp(a_{0}a_{1})

≤ ϵ · (\frac{3}{2}a_{0}a_{1})

and consequently

|z_{2} − z_{0}z_{1}| = \sqrt{ℜ(z_{2} − z_{0}z_{1})^{2} + ℑ(z_{2} − z_{0}z_{1})^{2}}

< ϵ\sqrt{(\frac{3}{2}a_{0}a_{1})^{2} + (2a_{0}b_{1} + 2b_{0}a_{1})^{2}}

= ϵ\sqrt{\frac{256}{55} |z_{0}z_{1}|^{2} − \frac{36}{55}(a_{0}b_{1} − b_{0}a_{1})^{2} − \frac{1}{220}(23a_{0}a_{1} − 32b_{0}b_{1})^{2}}

≤ ϵ\sqrt{256/55} |z_{0}z_{1}| < ϵ\sqrt{5} |z_{0}z_{1}|

as required.

Case R4: ulp(a_{0} ⊗ a_{1} − b_{0} ⊗ b_{1}) < ulp(b_{0}b_{1}) = ulp(a_{0}a_{1})

In this case, there is again no rounding error introduced in computing the difference between a_{0} ⊗ a_{1} and b_{0} ⊗ b_{1}, so we obtain

|ℜ(z_{2} − z_{0}z_{1})| ≤ |a_{0} ⊗ a_{1} − a_{0}a_{1}| + |b_{0} ⊗ b_{1} − b_{0}b_{1}|
< ϵ · (a_{0}a_{1} + b_{0}b_{1})

and consequently

|z_{2} − z_{0}z_{1}| = \sqrt{ℜ(z_{2} − z_{0}z_{1})^{2} + ℑ(z_{2} − z_{0}z_{1})^{2}}
< ϵ\sqrt{(a_{0}a_{1} + b_{0}b_{1})^{2} + (2a_{0}b_{1} + 2b_{0}a_{1})^{2}}
= ϵ\sqrt{5 |z_{0}z_{1}|^{2} − (a_{0}b_{1} − b_{0}a_{1})^{2} − 4(a_{0}a_{1} − b_{0}b_{1})^{2}}
≤ ϵ\sqrt{5} |z_{0}z_{1}|

as required. □

3. Worst-case multiplicands for β = 2

Having proved an upper bound on the relative error which can result from floating-point rounding when computing the product of complex values, we now turn to a more number-theoretic problem: finding precise worst-case inputs for β = 2. Starting with the assumption that some inputs produce errors very close to the proven upper bound, we will repeatedly reduce the set of possible inputs until an exhaustive search becomes feasible.

Theorem 2. Let β = 2 and assume that z_{0} = a_{0} + b_{0}i ≠ 0 and z_{1} = a_{1} + b_{1}i ≠ 0, where a_{0}, b_{0}, a_{1}, b_{1} are floating-point values with t-digit base-β significands, and z_{2} = ((a_{0} ⊗ a_{1}) ⊖ (b_{0} ⊗ b_{1})) + ((a_{0} ⊗ b_{1}) ⊕ (b_{0} ⊗ a_{1}))i are such that

(1)  0 ≤ a_{0}, b_{0}, a_{1}, b_{1}
(2)  b_{0}b_{1} ≤ a_{0}a_{1}
(3)  b_{0}a_{1} ≤ a_{0}b_{1}
(4)  1/2 ≤ a_{0}a_{1} < 1

and no overflow, underflow, or denormal values occur during the computation of z_{2}. Assume further that the results of arithmetic operations are correctly rounded to a nearest representable value and that

(5)  \frac{|z_{2} − z_{0}z_{1}|}{|z_{0}z_{1}|} > ϵ\sqrt{5 − nϵ} > ϵ · max (\sqrt{1024/207}, \sqrt{32/7} + 2ϵ)

for some positive integer n. Then

a_{0}a_{1} = 1/2 + (j_{aa} + 1/2)ϵ + k_{aa}ϵ^{2}
a_{0}b_{1} = 1/2 + (j_{ab} + 1/2)ϵ + k_{ab}ϵ^{2}
b_{0}a_{1} = 1/2 + (j_{ba} + 1/2)ϵ + k_{ba}ϵ^{2}
b_{0}b_{1} = 1/2 + (j_{bb} + 1/2)ϵ + k_{bb}ϵ^{2}

for some integers j_{xy}, k_{xy} satisfying

0 ≤ j_{aa}, j_{ab}, j_{ba}, j_{bb} < \frac{n}{4}

|k_{aa}| , |k_{bb}| < n

|k_{ab}| , |k_{ba}| < \frac{n}{2}

and a_{0} ≠ b_{0}, a_{1} ≠ b_{1}.

Proof. From equation (5), we note that ϵ ≤ nϵ < 11/207 < 2^{−4}; we will use this trivial bound later without explicit comment.

From the proof of Theorem 1, we know that Case R4 must hold, i.e., there is no error introduced in the computation of the difference between a_{0} ⊗ a_{1} and b_{0} ⊗ b_{1}, and ulp(b_{0}b_{1}) = ulp(a_{0}a_{1}). From inequalities (2) and (4) above, this implies that

1/2 ≤ b_{0}b_{1} ≤ a_{0}a_{1} < 1
|ℜ(z_{2} − z_{0}z_{1})| ≤ |a_{0} ⊗ a_{1} − a_{0}a_{1}| + |b_{0} ⊗ b_{1} − b_{0}b_{1}| ≤ ϵ.

We can now obtain lower bounds on |z_{0}z_{1}| and |z_{2} − z_{0}z_{1}|, using the fact that (a_{0}a_{1})(b_{0}b_{1}) = (a_{0}b_{1})(b_{0}a_{1}):

|z_{0}z_{1}|^{2} = (a_{0}^{2} + b_{0}^{2}) (a_{1}^{2} + b_{1}^{2})

= (a_{0}a_{1})^{2} + (a_{0}b_{1})^{2} + (b_{0}a_{1})^{2} + (b_{0}b_{1})^{2}

≥ (1/2)^{2} + (a_{0}b_{1})^{2} + \frac{(1/2)^{4}}{(a_{0}b_{1})^{2}} + (1/2)^{2} ≥ 1

|z_{2} − z_{0}z_{1}|^{2} > |z_{0}z_{1}|^{2} ϵ^{2}(5 − nϵ) ≥ ϵ^{2}(5 − nϵ)

as well as an upper bound on |z_{0}z_{1}|:

|z_{0}z_{1}|^{2} · \frac{1024ϵ^{2}}{207} < |z_{2} − z_{0}z_{1}|^{2}
= |ℜ(z_{2} − z_{0}z_{1})|^{2} + |ℑ(z_{2} − z_{0}z_{1})|^{2}
< ϵ^{2} + (ϵ · (2a_{0}b_{1} + 2b_{0}a_{1}))^{2}
≤ ϵ^{2} + 4ϵ^{2} |z_{0}z_{1}|^{2}
|z_{0}z_{1}|^{2} < \frac{207}{196}

We now note that

(a_{0}b_{1})^{2} ≤ |z_{0}z_{1}|^{2} − (a_{0}a_{1})^{2} − (b_{0}b_{1})^{2}

≤ \frac{207}{196} − \frac{1}{4} − \frac{1}{4} = \frac{109}{196}

so b_{0}a_{1} ≤ a_{0}b_{1} ≤ \sqrt{109/196} < 1 and a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1} ≤ \sqrt{109/49} · (1 + ϵ) < 2; this implies that ulp(b_{0}a_{1}) ≤ ulp(a_{0}b_{1}) ≤ ulp(1/2) and ulp(a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1}) ≤ ulp(1), and therefore

|a_{0} ⊗ b_{1} − a_{0}b_{1}| ≤ ϵ/2
|b_{0} ⊗ a_{1} − b_{0}a_{1}| ≤ ϵ/2
|((a_{0} ⊗ b_{1}) ⊕ (b_{0} ⊗ a_{1})) − (a_{0} ⊗ b_{1} + b_{0} ⊗ a_{1})| ≤ ϵ
|ℑ(z_{2} − z_{0}z_{1})| ≤ ϵ/2 + ϵ/2 + ϵ = 2ϵ

which allows us to place upper bounds on |z_{2} − z_{0}z_{1}| and |z_{0}z_{1}|:

|z_{2} − z_{0}z_{1}|^{2} = |ℜ(z_{2} − z_{0}z_{1})|^{2} + |ℑ(z_{2} − z_{0}z_{1})|^{2} ≤ (ϵ)^{2} + (2ϵ)^{2} = 5ϵ^{2}

|z_{0}z_{1}|^{2} < \frac{|z_{2} − z_{0}z_{1}|^{2}}{ϵ^{2}(5 − nϵ)} ≤ \frac{5}{5 − nϵ}.

Combining the known lower bound ϵ^{2}(5 − nϵ) for |z_{2} − z_{0}z_{1}|^{2} with the upper bounds on the error contributed by each individual rounding step, we find that

ϵ/2 − (1 − \sqrt{1 − nϵ})ϵ < |a_{0} ⊗ a_{1} − a_{0}a_{1}| ≤ ϵ/2

ϵ/2 − (1 − \sqrt{1 − nϵ})ϵ < |b_{0} ⊗ b_{1} − b_{0}b_{1}| ≤ ϵ/2
ϵ/2 − (2 − \sqrt{4 − nϵ})ϵ < |a_{0} ⊗ b_{1} − a_{0}b_{1}| ≤ ϵ/2
ϵ/2 − (2 − \sqrt{4 − nϵ})ϵ < |b_{0} ⊗ a_{1} − b_{0}a_{1}| ≤ ϵ/2

and similarly, by combining the upper bound on |z_{0}z_{1}|^{2} with lower bound of 1/2 for each pairwise product, we obtain

1/2 ≤ b_{0}b_{1} ≤ a_{0}a_{1} ≤ \sqrt{\frac{5}{5 − nϵ} − \frac{3}{4}} = \sqrt{\frac{5 + 3nϵ}{20 − 4nϵ}}

1/2 ≤ b_{0}a_{1} ≤ a_{0}b_{1} ≤ \sqrt{\frac{5}{5 − nϵ} − \frac{3}{4}} = \sqrt{\frac{5 + 3nϵ}{20 − 4nϵ}}.

Now consider the possible values for a_{0}a_{1} which satisfy these restrictions. Since it is the product of two values which are expressible using t digits of significand, a_{0}a_{1} can be exactly represented using 2t digits of significand; but since 1/2 ≤ a_{0}a_{1} < 1, this implies that a_{0}a_{1} is an integer multiple of ϵ^{2}. There is therefore at least one pair of integers j_{aa}, k_{aa} with 0 ≤ j_{aa} < ϵ^{−1}/2, |k_{aa}| ≤ ϵ^{−1}/2 for which

a_{0}a_{1} = 1/2 + (j_{aa} + 1/2)ϵ + k_{aa}ϵ^{2}.

Since a_{0} ⊗ a_{1} is the closest multiple of ϵ to a_{0}a_{1}, this implies that

ϵ/2 − (1 − \sqrt{1 − nϵ})ϵ < |a_{0} ⊗ a_{1} − a_{0}a_{1}| = ϵ/2 − |k_{aa}| ϵ^{2}
|k_{aa}| ϵ < 1 − \sqrt{1 − nϵ} < 1 − (1 − nϵ) = nϵ

i.e., |k_{aa}| < n, and similarly

1/2 + j_{aa}ϵ ≤ a_{0}a_{1} ≤ \sqrt{\frac{5 + 3nϵ}{20 − 4nϵ}} < \sqrt{1/4 + nϵ/4} < 1/2 + \frac{nϵ}{4}

i.e., 0 ≤ j_{aa} < n/4.

Applying the same argument to a_{0}b_{1}, b_{0}a_{1}, and b_{0}b_{1} allows us to infer that they possess the same structure, as required. To complete the proof, we note that the rounding errors from the products a_{0}a_{1} and b_{0}b_{1} must be in opposite directions (in order that they accumulate when subtracted), while the rounding errors from the products a_{0}b_{1} and b_{0}a_{1} must be in the same direction (in order that they accumulate when added); consequently, we must have a_{0} ≠ b_{0} and a_{1} ≠ b_{1}. □

Corollary 1. Assume that the preconditions of Theorem 2 are satisfied, and assume further that

(6)  \frac{1}{2} ≤ a_{0} < 1

and n ≤ 2ϵ^{−1/2}. Then

\frac{1}{2} < a_{0}, b_{0}, a_{1}, b_{1} < 1.

Proof. Assume that a_{1} ≥ 1. Then we can write

a_{0} = 1/2 + Aϵ
a_{1} = 1 + 2Bϵ

for some 0 ≤ A, B < (2ϵ)^{−1}. From Theorem 2, we have

1/2 + (A + B)ϵ + 2ABϵ^{2} = a_{0}a_{1} = 1/2 + (j_{aa} + 1/2)ϵ + k_{aa}ϵ^{2}

for some 0 ≤ j_{aa} < n/4, |k_{aa}| < n.

As a result, we must have A + B ≤ n/4 ≤ 1/2 · ϵ^{−1/2}, and since 0 ≤ A, B this implies 0 ≤ 2ABϵ^{2} ≤ ϵ/8. However, by reducing the equation above modulo ϵ, we find that 2ABϵ^{2} ≡ ϵ/2 + k_{aa}ϵ^{2}, which contradicts our bounds on 2ABϵ^{2}. Consequently, we can conclude that a_{1} < 1. Now we note that a_{0}a_{1} > 1/2 and a_{0} < 1, so a_{1} > 1/2, and we have both of the bounds required for a_{1}.

Applying the same argument to the other products provides the same bounds for a_{0}, b_{0}, and b_{1}. □

Corollary 2. Assume that the preconditions of Corollary 1 are satisfied, and assume further that n ≤ ϵ^{−1/2} and ϵ ≤ 2^{−6}. Then

j_{aa} − j_{ab} − j_{ba} + j_{bb} = 0,
|a_{0} − b_{0}| · |a_{1} − b_{1}| < 3nϵ^{2}.

Proof. From Theorem 2, we obtain that

a_{0}(a_{1} − b_{1}) = a_{0}a_{1} − a_{0}b_{1} = (j_{aa} − j_{ab})ϵ + (k_{aa} − k_{ab})ϵ^{2}

where |j_{aa} − j_{ab}| < \frac{n}{4}, |k_{aa} − k_{ab}| < \frac{3n}{2}, and since a_{0} > \frac{1}{2} (from Corollary 1), we can conclude that |a_{1} − b_{1}| < \frac{n}{2}ϵ + 3nϵ^{2}. Since a_{1} and b_{1} are integer multiples of ϵ and 3nϵ^{2} < ϵ/2, we conclude that |a_{1} − b_{1}| ≤ \frac{n}{2}ϵ. Applying the same argument to the product a_{1}(a_{0} − b_{0}) provides the same bound for |a_{0} − b_{0}|.

We now note that

|(j_{aa} − j_{ab} − j_{ba} + j_{bb})ϵ + (k_{aa} − k_{ab} − k_{ba} + k_{bb})ϵ^{2}| = |a_{0} − b_{0}| · |a_{1} − b_{1}|

≤ (\frac{n}{2}ϵ)^{2}

< \frac{ϵ}{4}

from our assumed upper bound on n, and consequently we can conclude that j_{aa} − j_{ab} − j_{ba} + j_{bb} = 0. Finally, this allows us to write

|a_{0} − b_{0}| · |a_{1} − b_{1}| = |k_{aa} − k_{ab} − k_{ba} + k_{bb}| ϵ^{2}
< 3nϵ^{2}

as required. □

Corollary 3. Assume that the preconditions of Corollary 1 are satisfied, and assume further that n ≤ \frac{1}{4}ϵ^{−1/2}. Then

(a_{0} − b_{0})(a_{1} − b_{1}) = 2(j_{aa} − j_{ab})(j_{aa} − j_{ba})ϵ^{2}
(a_{0} − b_{0})(a_{1} − b_{1})k_{aa} = (k_{aa} − k_{ab})(k_{aa} − k_{ba})ϵ^{2}

Proof. For brevity and clarity, we will write (a_{0} − b_{0})(a_{1} − b_{1}) = xϵ^{2} and note that x is an integer between −3n and 3n, from Corollary 2. Then

xa_{0}a_{1} = \frac{x}{2} + x (j_{aa} + \frac{1}{2}) ϵ + xk_{aa}ϵ^{2},
xa_{0}a_{1} = \frac{a_{0}(a_{1} − b_{1})}{ϵ} · \frac{a_{1}(a_{0} − b_{0})}{ϵ}
= ((j_{aa} − j_{ab}) + (k_{aa} − k_{ab})ϵ) ((j_{aa} − j_{ba}) + (k_{aa} − k_{ba})ϵ)

= (j_{aa} − j_{ab})(j_{aa} − j_{ba}) + ((j_{aa} − j_{ab})(k_{aa} − k_{ba}) + (j_{aa} − j_{ba})(k_{aa} − k_{ab})) ϵ

+ (k_{aa} − k_{ab})(k_{aa} − k_{ba})ϵ^{2}.

Consequently,

x − 2(j_{aa} − j_{ab})(j_{aa} − j_{ba})

= (2(j_{aa} − j_{ab})(k_{aa} − k_{ba}) + 2(j_{aa} − j_{ba})(k_{aa} − k_{ab}) − (2j_{aa} + 1)x) ϵ

+ (2(k_{aa} − k_{ab})(k_{aa} − k_{ba}) − 2k_{aa}x) ϵ^{2},

|x − 2(j_{aa} − j_{ab})(j_{aa} − j_{ba})| ≤ (2\frac{n}{4} \frac{3n}{2} + 2\frac{n}{4} \frac{3n}{2} + 3 (\frac{n}{2} + 1) n) ϵ

+ (2\frac{3n}{2} \frac{3n}{2} + 6n^{2}) ϵ^{2}

= (3n^{2} + 3n) ϵ + \frac{21}{2}n^{2}ϵ^{2}

≤ \frac{3}{16} + \frac{3}{4}\sqrt{ϵ} + \frac{21}{32}ϵ < 1,

and since the only integer with absolute value less than one is zero, we can conclude that x = 2(j_{aa} − j_{ab})(j_{aa} − j_{ba}) as required.

We now consider xa_{0}a_{1}ϵ^{−2} modulo \frac{1}{2}ϵ^{−1}, and note that

xk_{aa} ≡ xa_{0}a_{1}ϵ^{−2}
≡ (k_{aa} − k_{ab})(k_{aa} − k_{ba})

and further that

|xk_{aa} − (k_{aa} − k_{ab})(k_{aa} − k_{ba})| ≤ 3n · n + \frac{3n}{2} · \frac{3n}{2}

= \frac{21n^{2}}{4} ≤ \frac{21ϵ^{−1}}{64} < \frac{1}{2}ϵ^{−1}

and therefore xk_{aa} = (k_{aa} − k_{ab})(k_{aa} − k_{ba}). □

Theorem 3. Let β = 2 and assume that z_{0} = a_{0} + b_{0}i, z_{1} = a_{1} + b_{1}i, and z_{2} = ((a_{0} ⊗ a_{1}) ⊖ (b_{0} ⊗ b_{1})) + ((a_{0} ⊗ b_{1}) ⊕ (b_{0} ⊗ a_{1}))i are such that

(1)  0 ≤ a_{0}, b_{0}, a_{1}, b_{1}
(2)  b_{0}b_{1} ≤ a_{0}a_{1}
(3)  b_{0}a_{1} ≤ a_{0}b_{1}
(4)  1/2 ≤ a_{0}a_{1} < 1
(6)  1/2 ≤ a_{0} < 1

and no overflow, underflow, or denormal values occur during the computation of z_{2}. Assume further that the results of arithmetic operations are correctly rounded to the nearest representable value, and that

(5)  \frac{|z_{2} − z_{0}z_{1}|}{|z_{0}z_{1}|} > ϵ\sqrt{5 − nϵ} > ϵ · max (\sqrt{1024/207}, \sqrt{32/7} + 2ϵ)

for some n < \frac{1}{4}ϵ^{−1/2} and ϵ ≤ 2^{−6}. Then there exist integers c_{0}, d_{0}, α_{0}, β_{0}, c_{1}, d_{1}, α_{1}, β_{1} satisfying

a_{0} = \frac{c_{0}}{d_{0}} (1 + α_{0}ϵ)  b_{0} = \frac{c_{0}}{d_{0}} (1 + β_{0}ϵ)
a_{1} = \frac{c_{1}}{d_{1}} (1 + α_{1}ϵ)  b_{1} = \frac{c_{1}}{d_{1}} (1 + β_{1}ϵ)
gcd(c_{0}, d_{0}) = 1  \frac{d_{0}}{2} ≤ c_{0} ≤ d_{0}
gcd(c_{1}, d_{1}) = 1  \frac{d_{1}}{2} ≤ c_{1} ≤ d_{1}
2c_{0}c_{1} = d_{0}d_{1} < 3n  \frac{1}{2} < a_{0},b_{0}, a_{1}, b_{1} < 1
α_{0} ≡ β_{0} ≡ −ϵ^{−1} (mod d_{0})  α_{0} ≠ β_{0}
α_{1} ≡ β_{1} ≡ −ϵ^{−1} (mod d_{1})  α_{1} ≠ β_{1}
min (α_{0}, β_{0}) + min (α_{1}, β_{1}) ≥ 0  max (|α_{0}| , |β_{0}|) · max (|α_{1}| , |β_{1}|) < n

Proof. Let the values j_{aa}, j_{ab}, j_{ba}, j_{bb}, k_{aa}, k_{ab}, k_{ba}, and k_{bb} be as constructed in Theorem 2, and further let g_{0} = gcd(j_{aa} − j_{ab}, (a_{1} − b_{1})/ϵ). From Corollary 1 we know that 1/2 < a_{1}, b_{1} < 1, so a_{1} and b_{1} are multiples of ϵ; consequently g_{0} must be an integer. By the same argument, g_{1} = gcd(j_{aa} − j_{ba}, (a_{0} − b_{0})/ϵ) is an integer.

Now note that

g_{0}|(a_{1} − b_{1})ϵ^{−1}|(a_{1} − b_{1})a_{0}ϵ^{−2} = (j_{aa} − j_{ab})ϵ^{−1} + (k_{aa} − k_{ab})

and since g_{0}|(j_{aa}−j_{ab}), we can conclude that g_{0}|(k_{aa}−k_{ab}). By the same argument, g_{1}|(k_{aa} − k_{ba}).

We now write

c_{0} = \frac{j_{aa} − j_{ab}}{g_{0}}  d_{0} = \frac{a_{1} − b_{1}}{g_{0}ϵ}  e_{0} = \frac{k_{aa} − k_{ab}}{g_{0}}
c_{1} = \frac{j_{aa} − j_{ba}}{g_{1}}  d_{1} = \frac{a_{0} − b_{0}}{g_{1}ϵ}  e_{1} = \frac{k_{aa} − k_{ba}}{g_{1}}

and note that these values are all integers; further, from Corollary 3 we have d_{0}d_{1}k_{aa} = e_{0}e_{1} and d_{0}d_{1} = 2c_{0}c_{1}, and since gcd(c_{0}, d_{0}) = gcd(c_{1}, d_{1}) = 1 by construction, this implies gcd(c_{0}, c_{1}) = 1.

We now observe that

a_{0} = \frac{a_{0}(a_{1} − b_{1})}{a_{1} − b_{1}} = \frac{c_{0}g_{0}ϵ + e_{0}g_{0}ϵ^{2}}{d_{0}g_{0}ϵ} = \frac{c_{0} + e_{0}ϵ}{d_{0}}

a_{1} = \frac{a_{1}(a_{0} − b_{0})}{a_{0} − b_{0}} = \frac{c_{1}g_{1}ϵ + e_{1}g_{1}ϵ^{2}}{d_{1}g_{1}ϵ} = \frac{c_{1} + e_{1}ϵ}{d_{1}}

and therefore

\frac{1}{2} + (j_{aa} + \frac{1}{2}) ϵ + k_{aa}ϵ^{2} = a_{0}a_{1}

= \frac{c_{0}c_{1}}{d_{0}d_{1}} + \frac{c_{0}e_{1} + e_{0}c_{1}}{d_{0}d_{1}} ϵ + \frac{e_{0}e_{1}}{d_{0}d_{1}} ϵ^{2}

= \frac{1}{2} + \frac{c_{0}e_{1} + e_{0}c_{1}}{d_{0}d_{1}} ϵ + k_{aa}ϵ^{2}

and thus (using d_{0}d_{1} = 2c_{0}c_{1})

c_{0}c_{1} (2j_{aa} + 1) = c_{0}e_{1} + e_{0}c_{1}.

Consequently c_{0}|e_{0}c_{1} and c_{1}|c_{0}e_{1}, and since gcd(c_{0}, c_{1}) = 1 it follows that c_{0}|e_{0} and c_{1}|e_{1}. Writing e_{0} = c_{0}α_{0}, e_{1} = c_{1}α_{1} for integers α_{0}, α_{1}, we now have

a_{0} = \frac{c_{0}}{d_{0}} (1 + α_{0}ϵ)  a_{1} = \frac{c_{1}}{d_{1}} (1 + α_{1}ϵ)

and taking β_{0} = α_{0} + 2c_{1}g_{1}, β_{1} = α_{1} + 2c_{0}g_{0}, we have

b_{0} = \frac{c_{0}}{d_{0}} (1 + β_{0}ϵ)  b_{1} = \frac{c_{1}}{d_{1}} (1 + β_{1}ϵ)

as required.

The remaining conditions can be obtained by remembering that a_{0}, b_{0}, a_{1}, and b_{1} are integer multiples of ϵ, and by using the bounds on j_{xy} and k_{xy} given in Theorem 2. □

Corollary 4. In IEEE 754 single-precision arithmetic (β = 2, t = 24, ϵ = 2^{−24}), using “nearest even” rounding mode, the values^{1}

a_{0} = \frac{3}{4}  b_{0} = \frac{3}{4}(1 − 4ϵ)  a_{1} = \frac{2}{3}(1 + 11ϵ)  b_{1} = \frac{2}{3}(1 + 5ϵ)

result in a relative error δ ≈ ϵ\sqrt{5 − 168ϵ} ≈ ϵ\sqrt{4.9999899864} in z_{2}, and δ is the largest possible relative error for for IEEE 754 single-precision inputs provided that overflow, underflow, and denormals do not occur.

Proof. Straightforward computation for the values given establishes that

a_{0}a_{1} = \frac{1}{2}(1 + 11ϵ)  a_{0} ⊗ a_{1} = \frac{1}{2}(1 + 12ϵ)
b_{0}b_{1} = \frac{1}{2}(1 + ϵ − 20ϵ^{2})  b_{0} ⊗ b_{1} = \frac{1}{2}
ℜ(z_{0}z_{1}) = 5ϵ + 10ϵ^{2}  ℜ(z_{2}) = 6ϵ
a_{0}b_{1} = \frac{1}{2}(1 + 5ϵ)  a_{0} ⊗ b_{1} = \frac{1}{2}(1 + 4ϵ)
b_{0}a_{1} = \frac{1}{2}(1 + 7ϵ − 44ϵ^{2})  b_{0} ⊗ a_{1} = \frac{1}{2}(1 + 6ϵ)
ℑ(z_{0}z_{1}) = 1 + 6ϵ − 22ϵ^{2}  ℑ(z_{2}) = 1 + 4ϵ
|z_{2} − z_{0}z_{1}|^{2} = ϵ^{2}(5 − 108ϵ + O(ϵ^{2}))
|z_{0}z_{1}|^{2} = 1 + 12ϵ + O(ϵ^{2})

and the ratio of these provides the error as stated.

To prove that this is the largest possible relative error for IEEE single-precision inputs, we note that the mappings z_{0} → z_{0}i, z_{1} → z_{1}i, (z_{0}, z_{1}) → (z̄_{0}, z̄_{1}), (z_{0}, z_{1}) → (z_{1}, z_{0}), z_{0} → z_{0} · 2^{j}, and z_{1} → z_{1} · 2^{k} do not affect the relative error in z_{2}; consequently, this allows us to assume without loss of generality that conditions (1-4) and (6) are satisfied by the worst-case inputs. Using the results of Theorem 3, an exhaustive computer search (taking about five minutes in MAPLE on the second author’s 1.4 GHz laptop) completes the proof. □

1Note that while \frac{2}{3} is not an IEEE 754 single-precision value, \frac{2}{3} (1 + 5ϵ) and \frac{2}{3} (1 + 11ϵ) are, since ϵ^{−1} + 5 ≡ ϵ^{−1} + 11 ≡ 0 (mod 3).

Corollary 5. In IEEE 754 double-precision arithmetic (β = 2, t = 53, ϵ = 2^{−53}), using “nearest even” rounding mode, the values

a_{0} = \frac{3}{4}(1 + 4ϵ)  b_{0} = \frac{3}{4}  a_{1} = \frac{2}{3}(1 + 7ϵ)  b_{1} = \frac{2}{3}(1 + 1ϵ)

result in a relative error in z_{2} of approximately ϵ\sqrt{5 − 96ϵ} ≈ ϵ\sqrt{4.9999999999999893}, and this is the worst possible provided that overflow, underflow, and denormals do not occur.

Proof. Straightforward computation for the values given establishes that

|z_{2} − z_{0}z_{1}|^{2} = ϵ^{2}(5 − 36ϵ + O(ϵ^{2}))
|z_{0}z_{1}|^{2} = 1 + 12ϵ + O(ϵ^{2})

and the ratio of these provides the error as stated.

As in Corollary 4, an exhaustive search using the results of Theorem 3 (again, taking just a few minutes) completes the proof. □

For β = 2 and t > 6, the constructions given in Corollaries 4 and 5 for a_{0}, b_{0}, a_{1}, b_{1} provide for even and odd t respectively relative errors of ϵ\sqrt{5 − 168ϵ + O(ϵ^{2})} and ϵ\sqrt{5 − 96ϵ + O(ϵ^{2})}. We believe that these are the worst-case inputs for all sufficiently large t when β = 2.

4. A note on methods

The existence of this paper serves a strong demonstration of the power of experimental mathematics. The initial result — the upper bound of \sqrt{5}ϵ — was discovered experimentally seven years ago, on the basis of testing a few million random single-precision products.

Experimental methods became even more important when it came to the results concerning worst-case inputs. Here the approach taken was to perform an exhaustive search, taking several hours on the second author’s laptop, of IEEE single-precision inputs, using only a few arguments from Theorem 1 to prune the search. Once the worst few sets of inputs had been enumerated, it became clear that they possessed the structure described in Theorem 3, and it was natural to conjecture that this structure would be satisfied by the worst-case inputs in any precision. As is common with such problems, once the required result was known, constructing a proof was fairly straightforward.

References

1. N.J. Higham, Accuracy and Stability of Numerical Algorithms, Second Edition, SIAM, 2002.

2. C. Percival, Rapid multiplication modulo the sum and difference of highly composite numbers, Math. Comp. 72 (2002), 387–395.

Mathematical Sciences Institute, Australian National University, Canberra, ACT 0200, Australia

   E-mail address: complex@rpbrent.com

IRMACS Centre, Simon Fraser University, Burnaby, BC, Canada

   E-mail address: cperciva@irmacs.sfu.ca

INRIA Lorraine/LORIA, 615 rue du Jardin Botanique, F-54602 Villers-lès-Nancy Cedex, France

   E-mail address: zimmerma@loria.fr