Keywords: IEEE 754, nombre flottant, multiplication complexe, erreur d'arrondi, analyse d'erreur

Bornes d'erreur pour la multiplication de nombres flottants complexes

Résumé : On considère une arithmétique flottante de t chiffres de précision en base β, où tous les calculs sont effectués avec arrondi au plus proche. Nous montrons que le produit de deux nombres complexes z 0 et z 1 peut être calculé avec une erreur absolue d'au plus

fournit des bornes de 2 -24 √ 5 et 2 -53 √ 5 respectivement pour l'erreur relative dans les formats simple et double précision du standard IEEE 754, en supposant qu'aucun dépassement de capacité ni nombre dénormalisé n'intervient. Nous donnons également les pires cas pour les formats simple et double précision du standard IEEE 754. Enfin, nous considérons des pires cas génériques, et nous évoquons brièvement la multiplication de Karatsuba.

In memory of Erin Brent (1947Brent ( -2005) )

Introduction

In an earlier paper [START_REF]Rapid multiplication modulo the sum and difference of highly composite numbers[END_REF], 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 √ 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 √ 8 which is commonly used [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] 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 = 1 2 ulp(1) = 1 2 β 1-t ; we also define ulp(0) = 0. We use the notations x ⊕ y, x y, and x ⊗ y to represent rounded floating-point addition, subtraction, and multiplication of the values x and y.

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 floating-point 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 √ 5 = 1 2 β 1-t √ 5.

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 0 ≤ a 0 , b 0 , a 1 , b 1 [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] b 0 b 1 ≤ a 0 a 1 [START_REF]Rapid multiplication modulo the sum and difference of highly composite numbers[END_REF] 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 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 )| ≤ 1 2 ulp(a 0 ⊗ b 1 + b 0 ⊗ a 1 ) ≤ 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.

INRIA 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 )| ≤ 1 2 ulp(a 0 a 1 ) + 1 2 ulp(b 0 b 1 ) + 1 2 ulp(a 0 ⊗ a 1 -b 0 ⊗ b 1 ) ≤ 1 2 ulp(a 0 ⊗ a 1 -b 0 ⊗ b 1 ) + 1 2 ulp(b 0 b 1 ) + 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 | = (z 2 -z 0 z 1 ) 2 + (z 2 -z 0 z 1 ) 2 < (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 ) ≤ 32 7 |z 0 z 1 | 2 - 4 7 (a 0 b 1 -b 0 a 1 ) 2 - 1 7 (2a 0 a 1 -5b 0 b 1 ) 2 + 2 2 |z 0 z 1 | ≤ 32/7 + 2 |z 0 z 1 | < √ 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) ≤ 1 2 ulp(y), we obtain | (z 2 -z 0 z 1 )| ≤ 1 2 ulp(a 0 a 1 ) + 1 2 ulp(b 0 b 1 ) + 1 2 ulp(a 0 ⊗ a 1 -b 0 ⊗ b 1 ) ≤ 7 8 ulp(a 0 a 1 ) ≤ • 7 4 a 0 a 1 RR n°6068
and therefore

|z 2 -z 0 z 1 | = (z 2 -z 0 z 1 ) 2 + (z 2 -z 0 z 1 ) 2 < 7 4 a 0 a 1 2 + (2a 0 b 1 + 2b 0 a 1 ) 2 = 1024 207 |z 0 z 1 | 2 - 196 207 (a 0 b 1 -b 0 a 1 ) 2 - 1 3312 (79a 0 a 1 -128b 0 b 1 ) 2 ≤ 1024/207 |z 0 z 1 | < √ 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 ) ≤ 1 β ulp(a 0 a 1 ) ≤ 1 2 ulp(a 0 a 1 )
so we have

| (z 2 -z 0 z 1 )| ≤ 1 2 ulp(a 0 a 1 ) + 1 2 ulp(b 0 b 1 ) ≤ 3 4 ulp(a 0 a 1 ) ≤ • 3 2 a 0 a 1 and consequently |z 2 -z 0 z 1 | = (z 2 -z 0 z 1 ) 2 + (z 2 -z 0 z 1 ) 2 < 3 2 a 0 a 1 2 + (2a 0 b 1 + 2b 0 a 1 ) 2 = 256 55 |z 0 z 1 | 2 - 36 55 (a 0 b 1 -b 0 a 1 ) 2 - 1 220 (23a 0 a 1 -32b 0 b 1 ) 2 ≤ 256/55 |z 0 z 1 | < √ 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 )

INRIA

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 | = (z 2 -z 0 z 1 ) 2 + (z 2 -z 0 z 1 ) 2 < (a 0 a 1 + b 0 b 1 ) 2 + (2a 0 b 1 + 2b 0 a 1 ) 2 = 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 ≤ √ 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 numbertheoretic 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 0 ≤ a 0 , b 0 , a 1 , b 1 (1) b 0 b 1 ≤ a 0 a 1 (2) b 0 a 1 ≤ a 0 b 1 (3) 1/2 ≤ a 0 a 1 < 1 (4)
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) |z 2 -z 0 z 1 | |z 0 z 1 | > √ 5 -n > • max 1024/207, 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 < n 4

|k aa | , |k bb | < n |k ab | , |k ba | < 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 2 0 + b 2 0 a 2 1 + b 2 1 = (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 + (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 |: 1), and therefore

|z 0 z 1 | 2 • 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 < 207 196 INRIA We now note that (a 0 b 1 ) 2 ≤ |z 0 z 1 | 2 -(a 0 a 1 ) 2 -(b 0 b 1 ) 2 ≤ 207 196 - 1 4 - 1 4 = 109 196 so b 0 a 1 ≤ a 0 b 1 ≤ 109/196 < 1 and a 0 ⊗ b 1 + b 0 ⊗ a 1 ≤ 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(
|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 < |z 2 -z 0 z 1 | 2 2 (5 -n ) ≤ 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 - √ 1 -n ) < |a 0 ⊗ a 1 -a 0 a 1 | ≤ /2 /2 -(1 - √ 1 -n ) < |b 0 ⊗ b 1 -b 0 b 1 | ≤ /2 /2 -(2 - √ 4 -n ) < |a 0 ⊗ b 1 -a 0 b 1 | ≤ /2 /2 -(2 - √ 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 ≤ 5 5 -n - 3 4 = 5 + 3n 20 -4n 1/2 ≤ b 0 a 1 ≤ a 0 b 1 ≤ 5 5 -n - 3 4 = 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 - √ 1 -n ) < |a 0 ⊗ a 1 -a 0 a 1 | = /2 -|k aa | 2 |k aa | < 1 - √ 1 -n < 1 -(1 -n ) = n RR n°6068
i.e., |k aa | < n, and similarly

1/2 + j aa ≤ a 0 a 1 ≤ 5 + 3n 20 -4n < 1/4 + n /4 < 1/2 + 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) 1 2 ≤ a 0 < 1 and n ≤ 2 -1/2 . Then 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 . INRIA 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 | < n 4 , |k aa -k ab | < 3n 2 , and since a 0 > 1 2 (from Corollary 1), we can conclude that |a 1 -b 1 | < 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 | ≤ 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 | ≤ n 2
2 < 4 from our assumed upper bound on n, and consequently we can conclude that j aa -j abj 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 ≤ 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 = x 2 + x j aa + 1 2 + xk aa 2 , xa 0 a 1 = a 0 (a 1 -b 1 ) • 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 , RR n°6068 |x -2(j aa -j ab )(j aa -j ba )| ≤ 2 n 4 3n 2 + 2 n 4 3n 2 + 3 n 2 + 1 n + 2 3n 2 3n 2 + 6n 2 2 = 3n 2 + 3n + 21 2 n 2 2 ≤ 3 16 + 3 4 √ + 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 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 + 3n 2 • 3n 2 = 21n 2 4 ≤ 21 -1 64 < 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 0 ≤ a 0 , b 0 , a 1 , b 1 (1) b 0 b 1 ≤ a 0 a 1 (2) b 0 a 1 ≤ a 0 b 1 (3) 1/2 ≤ a 0 a 1 < 1 (4) 1/2 ≤ a 0 < 1 (6)
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) |z 2 -z 0 z 1 | |z 0 z 1 | > √ 5 -n > • max 1024/207, 32/7 + 2 INRIA for some n < 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 = c 0 d 0 (1 + α 0 ) b 0 = c 0 d 0 (1 + β 0 ) a 1 = c 1 d 1 (1 + α 1 ) b 1 = c 1 d 1 (1 + β 1 ) gcd(c 0 , d 0 ) = 1 d 0 2 ≤ c 0 ≤ d 0 gcd(c 1 , d 1 ) = 1 d 1 2 ≤ c 1 ≤ d 1 2c 0 c 1 = d 0 d 1 < 3n 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 We now observe that

c 0 = j aa -j ab g 0 d 0 = a 1 -b 1 g 0 e 0 = k aa -k ab g 0 c 1 = j aa -j ba g 1 d 1 = a 0 -b 0 g 1 e 1 = k aa -k ba
a 0 = a 0 (a 1 -b 1 ) a 1 -b 1 = c 0 g 0 + e 0 g 0 2 d 0 g 0 = c 0 + e 0 d 0 a 1 = a 1 (a 0 -b 0 ) a 0 -b 0 = c 1 g 1 + e 1 g 1 2 d 1 g 1 = c 1 + e 1 d 1

RR n°6068

and therefore 1 2

+ j aa + 1 2 + k aa 2 = a 0 a 1 = c 0 c 1 d 0 d 1 + c 0 e 1 + e 0 c 1 d 0 d 1 + e 0 e 1 d 0 d 1 2 = 1 2 + 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 = c 0 d 0 (1 + α 0 ) a 1 = c 1 d 1 (1 + α 1 )
and taking

β 0 = α 0 + 2c 1 g 1 , β 1 = α 1 + 2c 0 g 0 , we have b 0 = c 0 d 0 (1 + β 0 ) b 1 = 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 values1

a 0 = 3 4 b 0 = 3 4 (1 -4 ) a 1 = 2 3 (1 + 11 ) b 1 = 2 3 (1 + 5 )
result in a relative error δ ≈ √ 5 -168 ≈ √ 4.9999899864 in z 2 , and δ is the worst possible provided that overflow, underflow, and denormals do not occur.

Proof. Straightforward computation for the values given establishes that

a 0 a 1 = 1 2 (1 + 11 ) a 0 ⊗ a 1 = 1 2 (1 + 12 ) b 0 b 1 = 1 2 (1 + -20 2 ) b 0 ⊗ b 1 = 1 2 (z 0 z 1 ) = 5 + 10 2 (z 2 ) = 6 a 0 b 1 = 1 2 (1 + 5 ) a 0 ⊗ b 1 = 1 2 (1 + 4 ) b 0 a 1 = 1 2 (1 + 7 -44 2 ) b 0 ⊗ a 1 = 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 worst possible relative error, we note that the mappings z 0 → z 0 i, z 1 → z 1 i, (z 0 , z 1 ) → ( z0 , z1 ), (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.

Corollary 5. In IEEE 754 double-precision arithmetic (β = 2, t = 53, = 2 -53 ), using "nearest even" rounding mode, the values

a 0 = 3 4 (1 + 4 ) b 0 = 3 4 a 1 = 2 3 (1 + 7 ) b 1 = 2 3 (1 + 1 )
result in a relative error in z 2 of approximately √ 5 -96 ≈ √ 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 5 -168 + O( 2) and 5 -96 + O( 2). We believe that these are the worst-case inputs for all sufficiently large t when β = 2.

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 √ 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. This may give a relative error of up to 8 . Note that there are different ways to compute p 1 +p 2 -p 3 , which may give different rounding errors: p 1 ⊕(p 2 p 3 ) as above, or (p 1 ⊕p 2 ) p 3 , or even (p 1 p 3 ) ⊕ p 2 .

The worst case we found with the above way of computing p 1 + p 2 -p 3 is z 1 = 260 + 278i, z 2 = 268 -278i, with a precision of 8 bits, which gives a ⊕ b = 536, c d = 544, p 1 = 290816, p 2 = -72192, p 3 = 74752, p 2 ⊕ p 3 = 2560, p 2 p 3 = -147456, and p 1 ⊕ (p 2 p 3 ) = 143360. The rounded product is thus z 2 = 143360 + 2560i, whereas the exact one is z 0 z 1 = 146964 + 2224i, with a relative error of about 6.30 . 4.3. Getting rid of the 2nd-order term in Case R1. The bound we get in Case R1 is the following:

|z 2 -z 0 z 1 | < 32/7 + 2 |z 0 z 1 | .
This is the only case among R1, . . . , R4 where we get a 2nd-order term 2 . We show how to get rid of that term in the case β = 2.

Lemma 1. Let x > 0 be rounded to a value y, with rounding to nearest. Then |y-x| ≤ 1+ x.

Proof. Remember = 1 2 ulp(1) = 1 2 β 1-t . Without loss of generality, one can assume 1 ≤ x < β. Then 1 ≤ y ≤ β, with an absolute error |y -x| ≤ .

If

x < 1 + , x is rounded to y = 1, thus the error is x -1 ≤ 1+ x. If 1 + ≤ x, the error is at most ≤ 1+ x.
We will now prove the following result, from which it follows |z 2 -z 0 z 1 | ≤ 32/7|z 0 z 1 |.

Lemma 2. Assume radix β = 2, precision t ≥ 2, and we are in Case R1, i.e., ulp(b 0 b 1 ) ≤ ulp(a 0 a 1 ) ≤ ulp(a 0 ⊗ a 1 -b 0 ⊗ b 1 ). Then:

| (z 2 -z 0 z 1 )| ≤ • (2a 0 a 1 -b 0 b 1 ).
Proof. Without loss of generality, we assume 1 ≤ a 0 a 1 < 2. Let k ≥ 0 be the exponent difference between a 0 a 1 and b 0 b 1 , i.e., ulp(a 0 a 1 ) = 2 k ulp(b 0 b 1 ).

First assume k = 0; then 1

≤ b 0 ⊗ b 1 ≤ 2. The only possibility that Case R1 holds is that a 0 ⊗ a 1 -b 0 ⊗ b 1 = 1, which implies a 0 ⊗ a 1 = 2 and b 0 ⊗ b 1 = 1. We thus have 2 -≤ a 0 a 1 ≤ 2 and 1 ≤ b 0 b 1 ≤ 1 + . The real error is bounded by 2 , and 2a 0 a 1 -b 0 b 1 ≥ 2(2 -) -(1 + ) ≥ 3(1 -) ≥ 2. Thus | (z 2 -z 0 z 1 )| ≤ (2a 0 a 1 -b 0 b 1 ). Now assume k ≥ 1, i.e., 2 -k ≤ b 0 b 1 < 2 1-k . Since b 0 ⊗ b 1 ≥ 2 -k , a 0 ⊗ a 1 -b 0 ⊗ b 1 is an integer multiple of 2 1-k ,
and is larger or equal to 1 by hypothesis. We distinguish three cases, depending on the position of a 0 ⊗ a 1 -b 0 ⊗ b 1 with respect to 1 + and 1 + 3 .

First assume 1 + 3 ≤ a 0 ⊗ a 1 -b 0 ⊗ b 1 .
The error when rounding a 0 a 1 is bounded by:

≤ 1 + 3 (a 0 ⊗ a 1 -b 0 ⊗ b 1 ).
The same bounds holds for the error when rounding a 0 ⊗ a 1 -b 0 ⊗ b 1 , and that on b 0 b 1 is bounded by 1+ b 0 b 1 from Lemma 1. The difference between (2a 0 a 1 -b 0 b 1 ) and the sum of those three bounds has the same sign as 4a 0 a 1 -7b

0 b 1 + (4a 0 a 1 -5b 0 b 1 ), which is non-negative since b 0 b 1 ≤ 1 2 a 0 a 1 . (Otherwise b 0 ⊗ b 1 ≥ 1 2 a 0 ⊗ a 1 , and a 0 ⊗ a 1 -b 0 ⊗ b 1 ≤ 1 2 a 0 ⊗ a 1 ≤ 1.) Now assume 1 + < a 0 ⊗ a 1 -b 0 ⊗ b 1 < 1 + 3 . Since a 0 ⊗ a 1 -b 0 ⊗ b 1 is an integer multiple of 2 1-k , we have 1 + + 2 1-k ≤ a 0 ⊗ a 1 -b 0 ⊗ b 1 ≤ 1 + 3 -2 1-k ,
and a 0 ⊗ a 1 -b 0 ⊗ b 1 is rounded to 1 + 2 , with an error at most (1 -2 1-k ). The real error is thus bounded by

+ (1 -2 1-k ) 1 + (a 0 ⊗ a 1 -b 0 ⊗ b 1 ) + 1 + b 0 b 1 ,
using again Lemma 1, and a 0 ⊗ a 1 -b 0 ⊗ b 1 > 1 + . In the case k = 1, the difference between (2a 0 a 1 -b 0 b 1 ) and the error bound has the same sign as (1 + )a 0 a 1 -(1 + 2 )b 0 b 1 , which is nonnegative since as above we have b 0 b 1 ≤ 1 2 a 0 a 1 . For k ≥ 2, using b 0 b 1 ≤ 2 1-k a 0 a 1 , the difference between (2a 0 a 1 -b 0 b 1 ) and the error bound has the same sign as 1 -

2 1-k + (2 1-k -2) ≥ 0. Now assume 1 ≤ a 0 ⊗ a 1 -b 0 ⊗ b 1 ≤ 1 + . In that case a 0 ⊗ a 1 -b 0 ⊗ b 1 is rounded down to 1. If a 0 a 1 is rounded down and b 0 b 1 is rounded up, then a 0 ⊗ a 1 -b 0 ⊗ b 1 ≤ a 0 a 1 -b 0 b 1 , thus the total error is bounded by 2 (a 0 a 1 -b 0 b 1 ) + b 0 b 1 ≤ (2a 0 a 1 -b 0 b 1 ).
If either a 0 a 1 is rounded up or b 0 b 1 is rounded down, the three errors have different signs, so the total error is bounded by the absolute sum of the two largest bounds, those for a 0 a 1 and a 0 ⊗a 1 -b 0 ⊗b 1 . Let a 0 ⊗a 1 -b 0 ⊗b 1 = 1+η, with η ≤ . The total error is thus bounded by + η. On the other side, we have

a 0 ⊗ a 1 = (a 0 ⊗ a 1 -b 0 ⊗ b 1 ) + b 0 ⊗ b 1 ≥ 1 + η + 2 -k , and thus a 0 a 1 ≥ 1 + η + 2 -k -. It follows 2a 0 a 1 -b 0 b 1 ≥ 2(1 + η -). Hence (2a 0 a 1 -b 0 b 1 ) -( + η) = (1 -2 )( -η) ≥ 0.
As a consequence, we can simplify Eq. ( 5) from Theorems 2 and 3 into:

|z 2 -z 0 z 1 | |z 0 z 1 | > √ 5 -n > • 1024/207.
4.4. Generic Worst Cases for Binary Precision t. We prove here the statement at the end of Section 3.

Theorem 4. For β = 2 and t ≥ 7, the worst cases are:

a 0 = 3 4 , b 0 = 3 4 (1 -4 ), a 1 = 2 3 (1 + 11 ), b 1 = 2 3
(1 + 5 ) for t even, and

a 0 = 3 4 (1 + 4 ), b 0 = 3 4 , a 1 = 2 3 (1 + 7 ), b 1 = 2 3 (1 + ) for t odd. INRIA • to 1/2[1 + (β 0 + α 1 + sign(β 0 α 1 )) ] otherwise. Thus a 0 ⊗ b 1 + b 0 ⊗ a 1 = 1/2[2 + (α 0 + β 1 + β 0 + α 1 + d ) ] where -2 ≤ d ≤ 2.
Since min(α 0 , β 0 ) + min(α 1 , β 1 ) ≥ 0, we have

α 0 + β 1 + β 0 + α 1 ≥ 2[min(α 0 , β 0 ) + min(α 1 , β 1 )] ≥ 0.
Moreover with the condition α 0 = β 0 and α 1 = β 1 , we have

α 0 + β 1 + β 0 + α 1 ≥ 2[min(α 0 , β 0 ) + min(α 1 , β 1 )] + 2. Thus α 0 + β 1 + β 0 + α 1 + d ≥ 0.
Since a 0 ⊗b 1 +b 0 ⊗a 1 is a multiple of 2 -p-1 and it rounds to a value ≥ 1, i.e., a multiple of 2 -p+1 , the rounding error is of the form k2 -p-1 with -2 ≤ k ≤ 2. Thus the total rounding error on the imaginary part is of the form 1/2de where -4 ≤ d ≤ 4.

Thus we have:

(z 2 ) = 1/2(α 0 + α 1 -β 0 -β 1 + c) (z 2 ) = 1/2[2 + (α 0 + β 1 + β 0 + α 1 + d) ]
where -2 ≤ c ≤ 2 and -4 ≤ d ≤ 4. The exact product z 0 z 1 has:

(z 0 z 1 ) = 1/2(1 + α 0 )(1 + α 1 ) -1/2(1 + β 0 )(1 + β 1 ) = 1/2[(α 0 + α 1 -β 0 -β 1 ) + (α 0 α 1 -β 0 β 1 ) 2 ] (z 0 z 1 ) = 1/2[2 + (α 0 + α 1 + β 0 + β 1 ) + (α 0 β 1 + β 0 α 1 ) 2 ]
The difference is then:

(z 2 -z 0 z 1 ) = 1/2[c -(α 0 α 1 -β 0 β 1 ) 2 ] (z 2 -z 0 z 1 ) = 1/2[d -(α 0 β 1 + β 0 α 1 ) 2 ]
with |c| ≤ 2 and |d| ≤ 4. Thus neglecting 2nd order terms we have |z

2 -z 0 z 1 | 2 ≈ 1/4[c 2 + d 2 ] 2 , and |z 0 z 1 | 2 ≈ 1, so the ratio is ≈ 1/4[c 2 + d 2 ] 2 .
To get an error near 5 2 , we need |c| = 2 and |d| = 4, which implies that:

• α 0 + α 1 is odd and > 0,

• β 0 + β 1 is odd and > 0,

• sign(α 0 α 1 ) = -sign(β 0 β 1 ) (if say α 0 is zero, then α 1 is odd, and the even-rule implies that we replace sign(α 0 α 1 ) by 1 if α 1 + 1 is a multiple of 4, and -1 otherwise), • α 0 + β 1 is odd and > 0, • β 0 + α 1 is odd and > 0,

• sign(α 0 β 1 ) = sign(β 0 α 1 ), with the same convention as above in case one is zero. If none is zero, the conditions sign(α 0 α 1 ) = -sign(β 0 β 1 ) and sign(α 0 β 1 ) = sign(β 0 α 1 ) are incompatible.

Thus we can assume without loss of generality that α 0 = 0, which gives:

• α 1 is odd and > 0, say 4k + l 1 with l 1 = 1 or 3,

• β 0 + β 1 is odd and > 0, INRIA • sign(l 1 -2) = -sign(β 0 β 1 ),
• β 1 is odd and > 0, say 4j + m 1 with m 1 = 1 or 3,

• β 0 + α 1 is odd and > 0,

• sign(m 1 -2) = sign(β 0 α 1 ).

This implies:

• m 1 = l 1 , thus α 1 -β 1 = 2 mod 4, • β 0 even, > 0 for l 1 = 1, < 0 for l 1 = 3, • α 1 > 0, odd, • β 1 > 0, odd.
We then get:

(z 2 -z 0 z 1 ) = 1/2[c + β 0 β 1 2 ] (z 2 -z 0 z 1 ) = 1/2[d -β 0 α 1 2 ]
with |c| = 2 and |d| = 4. Since c = sign(αα 1 ) -sign(β 0 β 1 ), and sign(d) = sign(β 0 α 1 ), we can write:

(z 2 -z 0 z 1 ) = 1/2c[ -1/2|β 0 β 1 | 2 ] (z 2 -z 0 z 1 ) = 1/2d[ -1/4|β 0 α 1 | 2 ] with |c| = 2 and |d| = 4. Thus |z 2 -z 0 z 1 | 2 = [ -1/2|β 0 β 1 | 2 ] 2 + [2 -1/2|β 0 α 1 | 2 ] 2 = 2 [5 -(|β 0 β 1 | + 2|β 0 α 1 |) + 1/4(|β 0 β 1 | 2 + |β 0 α 1 | 2 ) 2 ]. Now |z 0 | 2 = (c 0 /d 0 ) 2 [(1 + α 0 ) 2 + (1 + β 0 ) 2 ] = (c 0 /d 0 ) 2 [2 + 2(α 0 + β 0 ) + α 2 0 β 2 0 2 ]; |z 1 | 2 = (c 1 /d 1 ) 2 [(1 + α 1 ) 2 + (1 + β 1 ) 2 ] = (c 1 /d 1 ) 2 [2 + 2(α 1 + β 1 ) + α 2 1 β 2 1 2
]. Thus neglecting 2nd order terms and using α 0 = 0:

|z 0 z 1 | 2 ≈ 1 + (β 0 + α 1 + β 1 ) .
Finally |z 2 -z 0 z 1 | 2 /|z 0 z 1 | 2 is (neglecting 2nd order terms): Since α 1 , β 1 > 0, this simplifies to: thus we are looking for the minimum of:

E = |β 0 |(β 1 + 2α 1 ) + 5(β 0 + α 1 + β 1 )
under the contraints [remember α 0 = 0]:

• (a 1 ) either α 1 = 4k + 1 for k ≥ 0, β 1 = 4j + 3 for j ≥ 0, β 0 even > 0,

• (a 2 ) or α 1 = 4k + 3 for k ≥ 0, β 1 = 4j + 1 for j ≥ 0, β 0 even < 0,

• (b) and β 0 + β 1 > 0,

• (c) and β 0 + α 1 > 0.

Since α 0 = 0, and α 0 ≡ 2 t mod d 0 , this implies that d 0 is a power of two. For β 0 > 0, the expression E above simplifies to:

E = β 0 (β 1 + 2α 1 ) + 5(β 0 + α 1 + β 1 )
= β 0 (4j + 3 + 2(4k + 1)) + 0 + 4k + 1 + 4j + 3) = (4j + 8k + 10)β 0 + (20j + 20k + 20).

Since β 0 ≥ 4, this gives E ≥ 36j + 52k + 60. Since α 1 = 4k + 1 and β 1 = 4j + 3 must differ by at least 3, we cannot have both k = j = 0, thus the smallest value is obtained for j = 1 and k = 0, with E ≥ 96. For that case to apply, since α 1 ≡ -2 t mod d 1 , we need for d 1 = 3 that is t odd. (Otherwise if d 1 > 3, necessarily d 1 ≥ 5, thus α 1 and β 1 differ by at least 5. Then j = k = 0 does not work; same for j = 1, k = 0 which would give α 1 = 1, β 1 = 7 which would imply d 1 = 3; same for j = 0, k = 1 which would give α 1 = 5, β 1 = 3; same for j = k = 1 which would give α 1 = 5, β 1 = 7; thus the smallest solution would be j = 2, k = 1 which gives E ≥ 184.)

For β 0 < 0, say β 0 = -i, E simplifies to: E = i(β 1 + 2α 1 ) + 5(-i + α 1 + β 1 )

= i(4j + 1 + 2(4k + 3)) + 5(-i + 4k + 3 + 4j + 1) = (4j + 8k + 2)i + (20j + 20k + 20).

Since |β 0 | ≥ 4, i.e., i ≥ 4, this gives E ≥ 36j + 52k + 28. Since α 1 = 4k + 3 and β 1 = 4j + 1 must differ by at least 3, and we must have min(α 0 , β 0 ) + min(α 1 , β 1 ) ≥ 0, i.e. min(4k + 3, 4j + 1) ≥ 4, we need j, k ≥ 1. We cannot have j = k = 1 since this would give (α 1 , β 1 ) = (7, 5) which do not differ by ≥ 3, nor j = 2, k = 1 which would give (α 1 , β 1 ) = (7, 9), thus the smallest possible value is j = 1, k = 2 which gives E ≥ 168.

INRIA

g 1 and

 1 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.

2 5 -

 5 (|β 0 β 1 | + 2|β 0 α 1 |) 1 + (β 0 + α 1 + β 1 ) ≈ 2 [5 -(|β 0 β 1 | + 2|β 0 α 1 | + 5(β 0 + α 1 + β 1 )) ].

2 [ 5 -

 5 (|β 0 |(β 1 + 2α 1 ) + 5(β 0 + α 1 + β 1 )) ],

  If d 0 = 2, then c 0 = 1 (because d 0 /2 ≤ c 0 ≤ d 0 and gcd(c 0 , d 0 ) = 1), thus c 1 = d 1 (because 2c 0 c 1 = d 0 d 1 ), which contradicts gcd(c 1 , d 1 ) = 1. Thus d 0 is at least 4. Now the constraint β 0 ≡ 0 mod d 0 , together with β 0 = α 0 , implies that |β 0 | ≥ 4. d 1 cannot be 2, since otherwise c 0 and c 1 would be both powers of 2. (More generally d 1 cannot be divisible by 2, since 4 divides d 0 , then 2 divides c 1 , and gcd(c 1 , d 1 ) = 1.) Thus d 1 is at least 3. Now the constraint α 1 = β 1 mod d 1 implies that α 1 and β 1 differ by at least d 1 ≥ 3.

RR n°6068

Note that while

is not an IEEE 754 single-precision value, 2 3 (1 + 5 ) and 2 3 (1 + 11 ) are, since -1 + 5 ≡ -1 + 11 ≡ 0 (mod 3).INRIA

http://www.inria.fr ISSN 0249-6399

Appendix

The content of this report up to here will appear in Mathematics of Computation. This appendix gives some additional results that do not appear in the Mathematics of Computation article.

4.1.

A bit of history. In an article entitled Error Bounds for Polynomial Evaluation and Complex Arithmetic published in the IMA Journal of Numerical Analysis (1986), Frank W. J. Olver proves the bound of 16/3 for complex multiplication, and follows with this remark:

"Indeed, in unpublished work R. P. Brent has demonstrated that in base 2, for example, [the error term] can be reduced to √ 5 . . . "

Thus this "unpublished work" is now published, even if it took twenty years! 4.2. Karatsuba Multiplication. The √ 5 bound we obtained here holds for classical complex multiplication, with 4 real products. Karatsuba's algorithm enables one to perform a complex multiplication with only 3 real products, for example:

It is natural to ask what the bound becomes in that case. The worst case could come when all the errors accumulate in the real part and there is no error in the imaginary part. Take for example a ≈ b ≈ c ≈ -d ≈ 1 and assume that all rounding errors are maximal and contribute to the same direction:

We first compute the relative error for those cases, and then prove this is the largest possible error for a given value of t.

First assume t even. We have a 0 a 1 = 1 2 (1 + 11 ), which is rounded

(exact since we are in Case R4), and a 0 ⊗b 1 +b 0 ⊗a 1 = 1+5 is rounded to 1+4 . We thus have z 2 = 6 +i(1+4 ), whereas z 0 z 1 = (5 + 10 2 ) + i(1

and taking the ratio gives 2 (5 -168 + 2515 2 -. . .). Now assume t odd. We have a 0 a 1 = 1 2 (1 + 11 + 28 2 ), which is rounded

. Thus a 0 ⊗a 1 -b 0 ⊗b 1 = 6 (exact since we are in Case R4), and a 0 ⊗b 1 +b 0 ⊗a 1 = 1+7 is rounded to 1+8 . We thus have z 2 = 6 +i(1+8 ), whereas

and taking the ratio gives 2 (5 -96 + 1027 2 -. . .). Now assume t even, t ≥ 20. Since 1 4 -1/2 ≥ 256, we can apply Theorem 3 with n = 168. Assume we apply Theorem 3 with n = 97. We first assume that none of α 0 , β 0 , α 1 , β 1 is zero. We have

The 2nd order term (α 0 α 1 ) 2 < n 2 ≤ 97 2 < /2. Thus a 0 a 1 rounds:

and β 0 + β 1 are both even or ≤ 0, or three other cases, with a general form of 1/2(α 0 + α 1 -β 0 -β 1 + c) where -2 ≤ c ≤ 2; and a 0 a 1 -b 0 b 1 rounds to itself since we are in Case R4.

a 0 b 1 = 1/2(1 + α 0 )(1 + β 1 ) rounds:

• to 1/2[1 + (α 0 + β 1 ) ] if α 0 + β 1 is even or ≤ 0,

• to 1/2[1 + (α 0 + β 1 + sign(α 0 β 1 )) ] otherwise. b 0 a 1 = 1/2(1 + β 0 )(1 + α 1 ) rounds

• to 1/2[1 + (β 0 + α 1 ) ] if β 0 + α 1 is even or ≤ 0,