Pascal-Sc. A Computer Language for Scientific Computation

298 Pages·1987·9.965 MB·

by  Gerd Bohlender, Christian Ullrich, Jürgen Wolff von Gudenberg, Werner Rheinboldt and Daniel Siewiorek (Auth.)

Pascal-SC
A Computer Language for Scientific Computation
Gerd Bohlender
Christian Ullrich
Jürgen Wolff von Gudenberg
Universität Karlsruhe
Institut für Angewandte Mathematik
Karlsruhe
Federal Republic of Germany
Louis B. Rail
Mathematics Research Center
University of Wisconsin—Madison
Madison, Wisconsin
ACADEMIC PRESS, INC.
Harcourt Brace Jovanovich, Publishers
Boston Orlando San Diego
New York Austin London Sydney
Tokyo Toronto
Copyright © 1987 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.
Orlando, Florida 32887
United Kingdom Edition published by
ACADEMIC PRESS INC. (LONDON) LTD
24-28 Oval Road, London NW1 7DX
Library of Congress Cataloging-in-Publication Data
Pascal-SC : a computer language for scientific computation.
(Perspectives in computing ; vol. 17)
Bibliography: p.
Includes index.
1. PASCAL-SC (Computer program language) 2. Computer arithmetic.
I. Bohlender, Gerd. II. Title: Scientific computation. III. Series.
QA76.73.P216P37 1987 005.13'3 87-14543
ISBN 0-12-111155-5
87 88 89 90 9 8 7 6 5 4 3 2 1
Printed in the United States of America
PREFACE
Pascal-SC is a powerful programming language developed as an extension
of standard Pascal. It offers numerous advantages and opens possibilities for
many new applications. The key new elements of this language are:
• Arithmetic operations with controlled rounding;
• An optimal scalar product;
• Functions with general result types;
• An operator concept for abstract data types;
• Overloading of procedures, functions, and operators;
• Strings;
• Dynamic arrays;
• Modules.
In addition, Pascal-SC provides a number of predefined standard modules
for calculation with complex numbers, real and complex intervals, and also
matrices and vectors with these component types.
Pascal-SC was developed over the last ten years at the Institute for Ap-
plied Mathematics at the University of Karlsruhe in West Germany. In the
course of this time, a number of publications and research reports appeared,
which in addition to the Authors of this book were produced by U. Allendoer-
fer, K. Braune, K. Gruener, E. Kaucher, R. Klatte, W. Kraemer, U. Kulisch,
and M. Neaga of the University of Karlsruhe, R. Kirchner and H.-W. Wip-
permann of the University of Kaiserslautern, H. Boehm and S. Rump of IBM
Germany, and W. L. Miranker of the IBM Thomas J. Watson Research Labo-
ratory, Yorktown Heights, U.S.A. The state of the development at present has
produced implementations of Pascal-SC for microcomputers with Z80, 8088,
and 68000 processors.
Following an informal introduction to standard Pascal, this book explains
the details of the new language elements of Pascal-SC, illustrated by examples
of numerical and non-numerical applications. The material is arranged so that
it is immediately suitable for use as a text. Even inexperienced programmers
will find that these new and important concepts will make it easier to write
v
VI PREFACE
clearer, better structured, and more powerful programs. In addition, this book
presents the basic fundamentals of the new computer arithmetic, and shows
how the tools provided by Pascal-SC can be used for controlled, maximally
accurate calculations for the solution of numerical problems.
A register and diagrams for the complete syntax of the language, descrip-
tions of its existing implementations, and a bibliography and indexes complete
the book.
We would like to thank all those who have contributed to the preparation
of this book by testing programs and the preparation of syntax diagrams or
typing the manuscript, particularly Mrs. K. Assisi and Mrs. U. Denninger,
Mr. D. Husung, Mr. W. Klein, Mr. M. Neaga, Mr. U. Nootbar in Germany,
and Mrs. Camilla Klyve in the U.S.A. The text of the American edition was
produced by the corresponding Author using TEX, but the syntax diagrams
still had to be done by hand and required the efforts of the capable staff of
Academic Press Boston.
Karlsruhe and Madison, Spring, 1987 The Authors
CHAPTER I
An Introduction to Pascal-SC
1.1. Why Pascal-SC?
The introduction of a new programming language requires careful justifica-
tion. The same applies to the modification or extension of an existing lan-
guage, particularly as successful and widely used as Pascal. However, the
present situation with regard to numerical computation is that existing pro-
gramming languages and the floating-point arithmetic they employ have not
bridged the gap between what engineers, scientists, and other users would like
to compute and what they can compute. This is particularly true with respect
to the amount of effort they have to devote to programming their problems.
For example, many problems lead to calculations with matrices or vec-
tors over the real or complex numbers. The architecture of the most powerful
computers available today is designed especially to perform these calcula-
tions. However, the most frequently used programming languages, such as
FORTRAN or Pascal, require that matrices and vectors be handled compo-
nentwise.
Since floating-point arithmetic only approximates the real operations, un-
avoidable rounding errors occur in the course of calculations. Except for Ada,
no widely available programming language makes any requirement regarding
the accuracy of floating-point arithmetic. For this reason, arithmetic has and
is being designed on the basis of how cheap and simple the hardware can be
made, or for fast execution time, particularly in the case of vector computers.
One hopes that the calculated results will lie in the vicinity of the solutions
sought, and the programmer responsible for the computation is burdened with
the task of a theoretical error analysis, which cannot be carried out in most
cases because the amount of data is too large.
In both numerical and non-numerical calculations, it is convenient to be
able to introduce abstract data types, that is, data types with the appropriate
operations.
Further, one would like to have the capability of separate translation of
program segments, and the use of libraries with the possibility of checking the
entries properly.
1
2 1. AN INTRODUCTION TO PASCAL-SC
In the following sections, we will describe in more detail how some of the
problems raised above are solved in Pascal-SC.
1.2. The propagation of roundoff error.
A basic requirement for floating-point arithmetic should be that at least the
four basic arithmetic operations +,—,*,/ are accurate in the sense that the
result of an operation is the closest floating-point number to the exact real
result. It will now be assumed that this simple requirement is satisfied. In
a single arithmetic operation, roundoff then effects only the last digit of the
result. The following examples show that in only a few successive operations,
roundoff error alters not only all digits of the result, but also its sign and
order of magnitude.
Example 1.2.1. Calculation of the value w = 9x4 — y4 + 2y2 three different
ways:
(a) w := 9*x*x*x*x - y*y*y*y + 2*y*y;
(b) w := 9*(x**4) - y**4 + 2*(y**2);
(c) w := (3*(x**2) - y**2)*(3*(x**2) + y**2) + 2*(y**2);
For x = 10864, y = 18817, these three expressions give, respectively,
(a) w= 1.58978 x 105,
(b) w = 8.41022 x 105,
(c) w = 1.000000000000.
In each case, 13-digit decimal arithmetic was used, and the result of each
operation was rounded to the closest floating-point number. For x — 18817
and y = 10864, all three evaluations give the same value w = 1.11442030725 x
1018.
Example 1.2.2. Calculate the value of the polynomial
(a) by multiplying out the powers,
(b) by the Horner scheme.
Using an accurate 53-digit binary arithmetic, which corresponds to about
16 decimal digits, the value p(0.707107) = 1.95541360881180 x 10"11 was
obtained for (a) and (b).
Example 1.2.3. By Gaussian elimination, the solutions of the linear system
of equations
41869520.5x - 102558961y = 0,
1.2. THE PROPAGATION OF ROUNDOFF ERROR 3
(a) using 12-digit decimal arithmetic are
x = 1.540378 ... x 105, y = 6.2885677... x 104;
(b) using 53-digit decimal arithmetic are
x = 8.656... x 107, y = 3.53... x 107.
Which of the above results are correct? How can the quality of a re-
sult be verified, without reevaluating the expression by hand? Is it possible
to give bounds within which the result can be guaranteed to lie? Can en-
tire expressions or algorithms be evaluated with the same accuracy as single
operations?
Answers to these questions should be expected from a programming lan-
guage for scientific computation. This language should also be easy to learn.
The assistance which Pascal-SC provides for the solutions of these prob-
lems is described in this book. The correct answers to the above examples
are:
Example 1.2.1: w = 1.0
Example 1.2.2: p[x) = -1.915273 ... x 10"11
Example 1.2.3: x = 2.05117922 x 108, y = 8.3739041 X 107
Without delving into details, the next few sections outline the theory
of computer arithmetic. For a full treatment, we recommend the book by
U. Kulisch and W. Mir anker, Computer Arithmetic in Theory and Practice.
1.3. Real floating-point arithmetic.
In a given floating-point system 5, two floating-point numbers u, v with u < v
are said to be adjacent if there is no floating-point number w such that u <
w < v. For x, y G S, the exact result x o y of an arithmetic operation o with
o G {+)—)*>/} is either a floating-point number or a real number w such
that u < w < V) where u, v are adjacent floating-point numbers. In the latter
case, w must be rounded to u or v in 5, since the result must be a floating-
point number. This property represents the basic requirement of reliability
of a floating-point arithmetic operation. Unfortunately, the floating-point
arithmetics implemented in software on many microcomputers or in hardware
on a number of mainframe computers do not satisfy this simplest requirement
of reliability. This situation has lead to the promulgation of standards for
floating-point arithmetic by the IEEE Computer Society.
An accurate rounding □ (x o y) of the exact result x o y gives either the
value u or v and minimizes the roundoff error
eps(a: o y) =□ (xo y) — x o y.
4 1. AN INTRODUCTION TO PASCAL-SC
If w is exactly halfway between u and v, then a definite, consistent rule is
applied to break the tie. In Pascal-SC, the value which is urthest from zero
is chosen as the answer. Thus, this rounding is antisymmetric, that is, □
(—x) = — □ (x), as required by the general theory.
The distance v — u between u and t; depends in particular on the precision,
or number of digits in the floating-point system used. This determines the
maximum roundoff error of a reliable operation. Of course, if the arithmetic
being used is not reliable, then roundoff error is not related in a simple way to
precision, and the attempt to "buy" greater accuracy by increasing precision
can be futile, as well as expensive.
There are other ways in which rounding to a reliable result can be carried
out, including:
(i) x is rounded upward to υ;
(ii) x is rounded downward to u.
The directed roundings (i) and (ii) are necessary to support interval
arithmetic, among other things. Rounding is said to be controllable if the
user can choose the way in which the result of a floating-point operation is
generated. Pascal-SC gives the user the choice of the rounding □ to the closest
floating-point number and the upward and downward directed roundings Δ,
V, for each of the four basic arithmetic operations +,—,*,/. Thus, a total of
twelve arithmetic operators are provided:
-f — * / (Rounding to the closest floating-point number);
+ > —> *> I > (Upward rounding);
+ < —< *< I < (Downward rounding).
Thus, the programmer can control the direction of rounding if desired, for
example to obtain guaranteed upper or lower bounds for the value of an
arithmetic expression. The reliability of Pascal-SC arithmetic also insures
that addition and multiplication are commutative, which is not true for a
number of existing floating-point arithmetics.
1.4. Complex floating-point arithmetic.
Numerical calculations frequently involve complex numbers. In contrast to
FORTRAN, standard Pascal does not provide the data type complex. Pascal-
SC provides a number of operators, functions, and procedures for computation
with complex numbers. The subroutines are organized into an arithmetic
module. In order to use these routines, the standard declaration
type complex = record re, im : real end;
is assumed. Thus, the standard representation of complex numbers in Pascal-
SC is in Cartesian coordinates. For input and output of complex numbers 2,
the standard format is (z. re, z. im).
1.4. COMPLEX FLOATING-POINT ARITHMETIC 5
In standard Pascal, procedures have to be used for complex arithmetic
and function evaluation. In Pascal-SC, operator overloading considerably sim-
plifies notation and hence programming. The following simple example shows
how to declare the operator +:
Example 1.4.1. Complex Addition
operator + (a,b: complex) res: complex;
var u: complex;
begin
u.re := a.re + b.re;
u.im := a.im + b.im;
res := u
end;
This declaration is essentially the same as for a function or procedure for
the same purpose. However, to add the two complex numbers v, w and assign
the result to z in the subsequent statement part of the program, one writes
only the simple statement
z := v + w;
All complex floating-point operations in Pascal-SC compute the real and
imaginary parts of their results with the highest possible accuracy. For addi-
tion and subtraction, operators similar to the example above are satisfactory;
however, multiplication and division require special algorithms to attain this
accuracy. Consider the following procedure, which does complex division ac-
cording to the ordinary formulas:
Example 1.4.2. Ordinary Complex Division
procedure cdiv(a,b: complex, var u: complex);
var denom: real;
begin
denom := b.re*b.re + b.im*b.im;
u.re := (a.re*b.re + a.im*b.im)/denom;
u.im := (a.im*b.re - a.re*b.im)/denom
end;
The Pascal-SC operator for complex division has a number of advantages
over the procedure cdiv. Among these are:
a. Accuracy. Because of the large number of rounding errors which
occur, ordinary Pascal functions and procedures for complex multiplication
or division, such as cdiv, cannot attain the accuracy of the complex Pascal-
SC operators, which give the best possible values for all results. In fact,
catastrophic cancellation can occur in ordinary calculations which renders
the computed results meaningless. For example, for
t; = (1.23456789,1.23456789),
w = (1.0000123e - 05,1.0000321e - 05),
6 1. AN INTRODUCTION TO PASCAL-SC
the result of the division is
v/w = (1.23454048308e + 05, -1.22216794612e + 00),
while the procedure cdiv(v,w,u) gives
u = (l.23454048308e + 05,-1.221673053),
where the incorrect digits of u are indicated in boldface. Even in this fairly
harmless looking example, the ordinary algorithm has lost five significant
digits in the imaginary part of the result of a single division, while all digits of
the Pascal-SC result are correct. To he sure, roundoff error will also increase
with successive use of Pascal-SC division, but at a slower and more predictable
rate.
b. Avoidance of overflow. In the Pascal-SC algorithms for complex multi-
plication and division overflow does not occur unless the real or imaginary part
of the result is greater in absolute value than the largest possible floating-point
number. By contrast, overflow can occur in the procedure cdiv during the cal-
culation of the intermediate results denom, u.re, or u.im, even though the final
result has real and imaginary parts which are representable as floating-point
numbers. For
t; = (3.0e + 99, -1.0e + 99), w = (l.Oe + 99, -1.0e + 99),
Pascal-SC complex division gives the result
v/w = (2.00000000000e + 00,1.00000000000e + 00),
while cdiv(v,w,u) overflows when trying to compute denom.
c. Ease of use. If u, υ, w are variables of type complex, then Pascal-SC
permits one to write
u := v/w;
in the source code of a program instead of
cdiv(v, w, u);
as would be required in standard Pascal. In case of complicated expressions
involving complex operands, the use of familiar mathematical notation with
operators instead of a sequence of procedure pays off. The program is more
likely to be written correctly in the first place, and will also be easier to
document and read.
1.5. Interval arithmetic.
Directed rounding makes it possible for the user to calculate bounds for the
values of arithmetic expressions. However, in the case of extensive expres-
sions, this is not very convenient and has to be done with great care. Interval