Computer Arithmetic and Self-Validating Numerical Methods

305 Pages·1990·10.99 MB·English

by  Christian Ullrich

Computer Arithmetic
and
Self-Validating
Numerical Methods
Edited by
Christian Ullrich
Institut für Informatik
Universität Basel
Basel, Switzerland
ACADEMIC PRESS, INC.
Harcourt Brace Jovanovich, Publishers
Boston San Diego New York
London Sydney Tokyo Toronto
This book is printed on acid-free paper. ®
Copyright © 1990 by Academic Press, Inc.
All rights reserved.
No part of this publication may be reproduced or
transmitted in any form or by any means, electronic
or mechanical, including photocopy, recording, or
any information storage and retrieval system, without
permission in writing from the publisher.
ACADEMIC PRESS, INC.
1250 Sixth Avenue, San Diego, CA 92101
United Kingdom Edition published by
ACADEMIC PRESS LIMITED
24-28 Oval Road, London NW1 7DX
Library of Congress Cataloging-in-Publication Data
Computer arithmetic and self-validating numerical methods /
edited by Christian Ullrich
p. cm—(Notes and reports in mathematics in science and
engineering; v. 7)
Proceedings of an international conference held Oct. 2-6,
1989, Basel, Switzerland.
Includes bibliographical references and index.
ISBN 0-12-708245-X (alk. paper)
1. Computer arithmetic—Congresses. 2. Numerical
analysis—Congresses. I. Ullrich, Christian. II. Series.
QA76.9.C62C663 1990
004'.01'51-dc20 90-816
CIP
Printed in the United States of America
90 91 92 93 9 8 7 6 5 4 3 2 1
Contributors
Numbers in parentheses indicate the pages on which the authors' contributions begin.
E. Adams (199), University of Karlsruhe, Institute for Applied Mathematics,
Kaiserstr. 12, D-7500 Karlsruhe 1, Federal Republic of Germany
G. Alefeld (55), University of Karlsruhe, Institute for Applied Mathematics, Kaiserstr.
12, D-7500 Karlsruhe 1, Federal Republic of Germany
H. Behnke (155), Technical University of Clausthal, Institute for Mathematics,
D-3392 Clausthal-Zellerfeld, Federal Republic of Germany
G. Bohlender (1), University of Karlsruhe, Institute for Applied Mathematics,
Kaiserstr. 12, D-7500 Karlsruhe 1, Federal Republic of Germany
L. Collatz (189), University of Hamburg, Institute for Applied Mathematics, Bun-
desstr. 55, D-2000 Hamburg 13, Federal Republic of Germany
G. F. Corliss (91), Marquette University, Department of Mathematics, Statistics, and
Computer Science, William Wehr Physics Building, Milwaukee, WI53233
F. Goerisch (137), Technical University of Clausthal, Institute for Mathematics,
D-3392 Clausthal-Zellerfeld, Federal Republic of Germany
K. Hafner (33), Siemens AG, Zentralbereich Forschung und Technik, Otto-Hahn-
Ring 6, D-8000 München 83, Federal Republic of Germany
Ζ. He (137), Technical University of Clausthal, Institute for Mathematics, D-3392
Clausthal-Zellerfeld, Federal Republic of Germany
E. Kaucher (269), University of Karlsruhe, Institute for Applied Mathematics,
Kaiserstr. 12, D-7500 Karlsruhe 1, Federal Republic of Germany
L. B. Rail (73), University of Wisconsin-Madison, Department of Mathematics, 610
Walnut Street, Madison, WI 53706
J. Schröder (247), University of Cologne, Institute for Mathematics, Weyertal 86-90,
D-5000 Cologne 41, Federal Republic of Germany
vii
viii Contributors
C. Schulz-Rinne (269), ΕΤΗ Zürich, Seminar for Applied Mathematics, ΕΤΗ Zen-
trum, CH-8092 Zurich, Switzerland
H. J. Stetter (171), Technical University of Vienna, Institute for Applied and Numer-
ical Mathematics, A-1040 Vienna, Austria
C. Ullrich (115), University of Basel, Institute for Informatics, Mittlere Str. 142,
CH-4056 Basel, Switzerland
Preface
Computer performance is expanding year by year. Increasingly expensive experi-
ments and development of prototypes can be replaced by mathematical models for
the simulation of technical processes on a computer. A good example is the large
field called numerical fluid dynamics. However, this process shifts the responsiblity
for reliable results to physical and mathematical models, and thus the automatic verifi-
cation of computed results becomes an essential tool. Automatic result verification
permits distinction between the effects of a mathematical model and computation-
al inaccuracies. Only when numerical errors are virtually eliminated can physical
and mathematical models be developed systematically.
Many workshops and tutorials have been devoted to this area during the last
decade. The most recent event, a first international conference on "Computer Arith-
metic and Self-Validating Numerical Methods," was held from October 2 to 6, 1989,
in Basel.
This volume is a collection of invited papers from this meeting. The contributed
papers will be published in a special issue of the IMACS Journal "Mathematics and
Computers in Simulation." Previous related volumes published by Academic Press
include Computer Arithmetic in Theory and Practice [Kulisch and Miranker, 1981],
Introduction to Interval Computations [Alefeld and Herzberger, 1983], Self-Validating
Numerics for Function Space Problems [Kaucher and Miranker, eds., 1984], A New
Approach to Scientific Computation [Kulisch and Miranker, eds., 1983], PASCAL-
SC: A Computer Language for Scientific Computation [Bohlender, Ullrich, Wolff
von Gundenberg, Rail, 1987] and Reliability in Computing [Moore, ed., 1988].
In the first part of this book, a number of papers give a tutorial introduction to
the following topics: computer arithmetic with operations of maximum accuracy,
differentiation arithmetic and enclosure methods, and programming languages for
self-validating numerical methods. In the second part, the authors discuss the deter-
mination of guaranteed bounds for eigenvalues by variational methods and the guaran-
ix
χ Preface
teed inclusion of solutions of differential equations. Finally, an appendix supplies
the IMACS-GAMM Resolution on Computer Arithmetic, which is intended to
influence and put pressure on manufacturers to implement computer arithmetic oper-
ations with necessary care.
Acknowledgments
Support for the symposium from the following sources is gratefully acknowledged:
the Institute for Informatics, the University of Basel, the Government of the Kanton
Basel-Stadt and the Nationalfonds of Switzerland; GAMM (Gesellschaft fur
Angewandte Mathematik und Mechanik); Gl (Gesellschaft fur Informatik); IMACS
(International Association for Mathematics and Computers in Simulation); F.
Hoffmann-La Roche & Co., Basel; IBM Schweiz; NCR Schweiz; Sandoz AG, Basel.
An excellent demonstration of software was made possible by efforts of Daniel
Hollenstein and others at the IBM offices at Basel and Zürich.
Many thanks are due to my collaborators Carlos Falco Korn, Stefan König, Roman
Reith, and others from the Institute for Informatics of the Basel University. Last but
not least, many thanks to our secretary, Agnes Mathys, who took care untiringly of
organizational work of the symposium.
xi
What Do We Need Beyond IEEE Arithmetic ?
G. Bohlender
Institut für Angewandte Mathematik
Universität Karlsruhe
Federal Republic of Germany
Abstract: While the four usual floating-point operations are the basis of real floating-
point arithmetic, the scalar product is the basis of the operations in higher numerical
spaces, such as matrices, vectors, etc. In addition, an exact scalar product is an in-
valuable tool for the verified solution of numerical problems by means of enclosure
methods. Therefore, computer arithmetic including such an exact scalar product is a
significant extension of IEEE arithmetic. In this paper, several algorithms for an im-
plementation are sketched, problems are discussed and solutions are suggested; finally,
typical designs and implementations are summarized and illustrated.
1 Introduction
Under unfavourable conditions, the small rounding error which is involved in each
floating-point operation can totally invalidate results. This may happen even after
only a few operations, as was demonstrated in [Rump 83a, Rump 83b, Ham 89, Ratz
89, Schu 89]. Let us study two simple examples:
50 50
1. The sum (1 -f 10 ) - 10 returns the wrong result 0 using ordinary floating-point
operations.
2. Consider a linear equation system A · χ = b with the matrix
and the vectors b = and χ = (x{). The true solution of this linear
system is x = 205117922,x = 83739041. But even in IEEE double precision
x 2
arithmetic, the following completely wrong results are computed
02 2 02 1 a =
X\ = Û&22/(ÛH " ~~ " a) 102558961
X2 = —a2i/(an · α 2 — αΐ2 · α ι) = 41869520.5
2 2
Note that the computation of x\ and x involves only four floating-point opera-
2
tions each!
Computer Arithmetic and Self-Validating 1 Copyright © 1990 by Academic Press, Inc.
Numerical Methods All rights of reproduction in any form reserved.
ISBN 0-12-708245-X
2 G. Bohlender
These examples show that ordinary floating-point arithmetic does not suffice to
compute reliable and precise results. While the first example can evidently be solved
by means of an exact scalar product or an exact sum, the second example requires
more subtle methods. For this purpose, enclosure methods have been developed which
compute a verified interval enclosure of a given problem. These enclosure methods are
again based on the exact evaluation of scalar products. For interval arithmetic see [Moo
66, Ale 74, Ale 83], for enclosure methods see [Kul 82, Kul 83a, Kau 84, Kau 87, Kul
88, Moo 88, Kul 89, Uli 90] and other papers in this volume. Enclosure methods are
typically programmed in computer languages like PASCAL-SC, FORTRAN-SC, etc.
which have special extensions for scientific and engineering computations; see [Nea 84,
Boh 86a, Boh 87a, Kul 87a, Kul 87b, Ble 87, Metz 88, Fal 89, Hus 88, Hus 89a, Hus
89b, Hahn 88, Hahn 89]. For an overview over systems which involve enclosure methods
see e.g. [Boh 89a].
In the following sections, let us concentrate on the implementation of an augmented
floating-point arithmetic and in particular on the optimal scalar product. Directed
roundings (which are needed for interval arithmetic) are required by the IEEE standard
for floating-point arithmetic [IEEE 85, IEEE 87] and are therefore available on many
modern floating-point processors. In section 2, the scalar product for floating-point
numbers is defined and its relation to the IEEE norm for floating-point arithmetic is
discussed. In section 3, several algorithms are roughly sketched for the computation
of scalar products with maximum accuracy. In section 4, some problems are described
which exist in the implementation of these algorithms in software or hardware and
several solutions for each problem are suggested. In section 5, some typical designs and
implementations of the scalar product in scientific research projects and in commercial
products are summarized including concepts for vector processors, parallel processors,
and VLSI implementations. In section 6, finally these ideas are summarized.
2 Scalar products and IEEE arithmetic
In the past twenty years, a uniform mathematical theory has been developed by Kulisch
and Miranker [Kul 76a, Kul 81]. It describes how an arithmetic operation * which is
defined on a superset R can be transferred on a subset F. The basic idea is to execute
the operation * in the superset R exactly and round the result into the subset F by
means of a so-called semimorphism • (a monotone projection with certain algebraic
properties). This method can be used to define operations for all customary spaces
of numerical computation, such as real and complex floating-point numbers, intervals,
matrices, vectors, etc.
Let us choose the set R = 1R of real numbers, and the subset F = F(6,p, e , e )
m t n max
of floating-point numbers which is characterized by four integer constants: base b > 2,
precision ρ > 1, minimum exponent e ; and maximum exponent e . The set F
m n m a x
contains floating-point numbers of the form
(-l)'.b<-(0.dd...d),
12 p
where s is the sign, e the exponent (e ; < e < e ), and d^ are base-ό digits of
m n m a x
the mantissa (0 < di < 6 — 1). (This representation differs slightly from the IEEE
What Do We Need Beyond IEEE Arithmetic? 3
representation where the decimal point is between the first and the second digit.) Let
* be any of the customary operations +,—,·,/ and let • be one of the roundings Ο
(round to nearest), V (round towards minus infinity), Δ (round towards plus infinity),
or Ο (truncate, round towards zero). Then the procedure described above agrees
τ
exactly with the definition of operations in the IEEE norms for floating-point arith-
metic ([IEEE 85] for binary arithmetic, [IEEE 87] for base-independent arithmetic).
Nevertheless, the two methods are clearly different. The aims of IEEE arithmetic are
standardization and reasonable definition of real floating-point arithmetic; a consider-
able effort was made to define computations with infinities, signed zeros, and "NaNs"
(for non-numeric data), and to provide exception handling, etc.
In contrast, the definition by semimorphism is not specialized on real floating-point
arithmetic but instead provides a uniform definition of arithmetic operations for all
customary higher numerical spaces as well. Therefore, it is less specific about data
formats, etc. While the term "vector processor" is usually defined in a technological
sense, here we can interpret the idea of a semimorphism as a mathematical definition
of a vector processor.
All these operations in the higher numerical spaces can be implemented efficiently,
if the basic operations for real floating-point numbers are defined by semimorphism
(or equivalently according to the IEEE norm), if the directed roundings toward plus or
minus infinity are available, and if a scalar product with maximum accuracy is available.
For all floating-point operands a;, 6, G F, such a scalar product shall compute the
correctly rounded result
• (a! · b + a · b + ... + a - b)
x 2 2 n n
as if it first produced an intermediate result correct to infinite precision and with un-
bounded range, and then rounded this intermediate result to the destination floating-
point system according to the selected rounding mode •. This contrasts with the
traditional computation of a scalar product which involves a rounding in each multi-
plication and addition. In the following sections, always such a scalar product with
maximum accuracy is meant, even if this is not explicitly mentioned.
Note that this definition conforms with the IMACS-GAMM Resolution on Com-
puter Arithmetic [IGR 89] about the elementary compound operations "accumulate"
for the components of a vector and "multiply and accumulate" for two vectors.
Let us conclude this section with a remark about scalar products and IEEE arith-
metic: in addition to the usual numeric values, a floating-point system according to the
IEEE norm contains special representations for signed infinity, signed zeros, and not-
a-number (which is intended to encode non-numeric information in the floating-point
format and is frequently used in connection with exception handling). These special
values can be handled in a scalar product as defined in the IEEE norm. Exceptions
and traps in a scalar product can again be handled as described in the IEEE norm; of
course the exceptions underflow, overflow and inexact can never occur during the com-
putation of the "infinitly precise" intermediate result, but only when this intermediate
result is rounded to the destination floating-point format.