[ut.na] NA Digest Volume 88 : Issue 15

krj@csri.toronto.edu (Ken Jackson) (04/13/88)

NA Digest   Tuesday, April 12, 1988   Volume 88 : Issue 15

Today's Editor: Cleve Moler

Today's Topics:
 
     Fortran Code Generation by Computer Algebra Systems
     Temporary Address Change for Hans Mittlemann
     Models of the Human Head or Ear
     Transcendental Functions on the 80387
     Inverse Square Root of a Matrix
     International Conference on Applications of Supercomputers
     Floating-Point Indoctrination by Kahan

-------------------------------------------------------

From: Heikki Apiola <APIOLA%FINFUN.BITNET@forsythe.stanford.edu>
Date: Sun, 3 Apr 88 20:37 O
Subject: Fortran Code Generation by Computer Algebra Systems

I would like to know about existing computer algebra based (especially
MACSYMA) Fortran code generation packages and their availability.
There have been articles about GENTRAN by Paul S. Wang  in J. of Symbolic
Computing.

Please reply to:

Heikki Apiola
e-mail: apiola@finfun.bitnet

Scientific Computing Centre
State Computer Centre
PO-BOX 40
SF-02101 Espoo
Finland


------------------------------

From: Hans D. Mittelmann  <aihdm@asuacad.bitnet, na.mittelmann>
Date: 8 Apr 88   19:35 MST
Subject: Temporary Address Change for Hans Mittlemann

My address from 1 May 1988 until 31 July 1988 will be:

    Institut fuer Angewandte Mathematik
    Universitaet Erlangen-Nuernberg
    Martensstrasse 3
    8520 Erlangen
    West Germany

    HJW@DERRZE0.BITNET   or   NA.MITTELMANN@NA-NET.STANFORD.EDU

   --Hans Mittelmann

------------------------------

From: Jim Buchanan <ACAD8005%RYERSON.BITNET@forsythe.stanford.edu>
Date: Tue, 12 Apr 88 16:40:24 EST
Subject: Models of the Human Head or Ear

I am looking for information on mathematical models of the human head
in general, or of the ear in particular.  Does anyone have any
information in this area?   Thanks in advance.
  --Jim Buchanan


------------------------------

From: Warren E Ferguson <smu!ferguson@uunet.UU.NET>
Date: Tue, 12 Apr 88 15:31:27 CDT
Subject: Transcendental Functions on the 80387

I am thinking of having students in one of my classes create a partial
emulator of the INTEL 80387 chip.  After looking at the spec sheets for the
80387, it appears that INTEL is using CORDIC's to compute their special
functions.  Does anyone know of literature describing which CORDIC's INTEL is
using for these functions?  Did members of the IEEE Binary Standard ever
suggest how the Standard could be used to compute the basic transcendental
functions?

[Editor's Note: See the announcement of Kahan's lectures at Sun
 later in this Digest.]

Warren Ferguson, Southern Methodist University 
na.ferguson
ferguson%smu.uucp@uunet


------------------------------

From: Cleve Moler <na.moler@na-net.stanford.edu>
Date: Tue, 12 Apr 88 22:14:26 PDT
Subject: Inverse Square Root of a Matrix

I was asked recently if I had any software to compute the inverse
square root of a symmetric, positive definite matrix.  All I had to
suggest was EISPACK and the obvious diagonalization algorithm.
Has anybody got anything faster?  The setting involves lots of
matrices of small order, say less than 10-by-10.
  --Cleve


------------------------------

From: [Three levels of forwarding]
Date: Fri, 08 Apr 88 14:05:16 -0800
Subject: International Conference on Applications of Supercomputers


Announcement of ASE 89

The International Conference on Applications of Supercomputers in
Engineering will be held 5-7 September 1989 (yes 1989) at Southampton
University, United Kingdom.  This International Conference is convened
to increase the awareness of the potential of the new supercomputers 
amongst scientists and engineers.

Papers are invited on the following topics, and the abstracts should reach
the Organizers by 30 December 1988.  Topics include:
 	* computer architecture and current supercomputers
	* use and development of transputers
 	* case studies in
 		hydrodynamics
		groundwater 
		environmental sciences, etc.
	* new parallel algorithms for finite differences
		and finite and boundary element methods
	* optimization problems 
	* structural and dynamic analysis including crash and impact 
		simulation applications in automotive engineering
	* use of supercomputers in CAD and CIM
	* benchmarks and comparison between different processors
 	* fluid dynamics and atmospheric science applications
	* expert system applications 
	* future trends.

For further information, please contact:

	Liz Newman
	Conference Secretary
	Computational Mechanics Institute
	Ashurst Lodge
	Ashurt
	Southampton, S04 2AA
	England

	Telephone: 	(042129) 3223
	Telex:		47388
	Attn:		COMPMECH
	Fax:		(042129) 2853



------------------------------

From: David Hough <dgh@Sun.COM>
Date: Fri, 8 Apr 88 17:47:25 PDT
Subject: Floating-Point Indoctrination by Kahan


     Sun Microsystems presents a public series of lectures:

 Computer System Support for Scientific and Engineering Computation

                           Prof. W. Kahan
                      University of California

 The environment for computing
 Comparisons of floating-point arithmetics in computers
 Programmers' models of floating-point arithmetic
 Comparison of computers' floating-point register architectures
 Real elementary transcendental functions
 Complex arithmetic
 Exception handling
 Numerical lapses in standard programming languages

     A detailed syllabus follows this announcement.

FORMAT

     The lecture series will start at 8:00 AM on Tuesday May
3, 1988, in the Training Room at Sun Microsystems building
7, 2700 Coast Avenue, Mountain View, CA. From US 101, take
the Amphitheater Parkway exit in Mountain View and turn left
on Garcia.  Turn right on Marine Way, then right on Coast
Avenue to its end.

     The first lecture will provide an overview of the sub-
ject matter to be discussed in the course, and the last lec-
ture will summarize the material presented previously and
survey advanced topics and directions for further research.
Approximately 24 lectures will be presented Tuesdays and
Thursdays through approximately July 29.  Each lecture will
start at 8:00 AM and last for two hours.

     Lectures are open to the public subject to space avail-
able, and persons with experience in some of the lecture
topics are encouraged to attend and contribute to the dis-
cussion. The amount of space available will be determined by
the number of pre-registrants, so persons desiring to insure
a place should register in advance with the Sun Training and
Development Department using the form below.

     Persons with training or experience comparable to that
represented by a B. S. degree in Computer Science should be
able to understand most of the conclusions. The mathematical
background required for a deep understanding of the lecture
topics varies considerably, especially for some of the
proofs. Chapter 4 of Knuth's Seminumerical Algorithms exem-
plifies the presentation level - accessible to undergradu-
ates, but containing interesting mathematics at all levels.
Note that the scope of the course is much broader than
Knuth's book.

     Written lecture notes will be distributed to regis-
trants about a week after each lecture.  These notes may
subsequently be commercially published to supersede the 1973
Implementation of Algorithms lecture notes distributed by
NTIS.

     The first and last lectures will be videotaped; regis-
trants will receive a copy of the first for their personal
use.

     This lecture series does not provide academic credit.
Registrants interested in obtaining formal recognition of
participation may present an acceptable individual or group
project.  Outstanding project reports will be included in
the published course lecture notes. An example of such a
project is to reimplement some of the 4.3 BSD libm functions
using fast table-driven techniques.  Prof. Kahan will be
available after lectures for consultation on projects.

     Questions pertaining to the technical scope and depth
of the course may be addressed to David Hough,
dhough@sun.com, (415) 691-2268.  Questions pertaining to
registration procedures should be directed to Janet Rath,
jrath@sun.com, Sun Training and Development, (415) 494-6058
or 494-8008 x 6058.

VENDOR LITERATURE

     Vendors of semiconductor devices, calculators, comput-
ers, and software products oriented toward scientific and
engineering computation are encouraged to contribute techni-
cal literature describing their products for distribution to
attendees.  Please contact David Hough after 22 April to
determine quantities required.

REGISTRATION - general

     The registration fee guarantees a place and provides
lecture notes and a videotape of the initial lecture of 3
May.  Registration by 22 April 1988 is recommended, but late
registrants will be accommodated subject to available space
through at least the second class meeting, Thursday 5 May.
The registration fee is $1000, subject to reduction if more
participants register than projected.

     Send a check for the registration fee, payable to Sun
Microsystems, with the Registration Form below.

     Full-time students at accredited institutions may
register without fee; they will receive no videotape.  Such
students should submit current evidence of their status in
lieu of a check.

REGISTRATION FORM

     Please send the following information to:

        Sun Microsystems Training and Development
                Attention: Janet Rath - Floating-Point Lectures
        MS A6-05
        2550 Garcia
        Mountain View, CA 94043

Your name:


USENET or ARPANET address for electronic mail, if any:


Mailing address:




FULL-TIME STUDENTS ONLY:  Your institution, department, next degree and projected date:



Telephone number during business hours:


Highest academic degree and major field:


Circle any of the following languages of which you have a working knowledge:

Ada  APL  BASIC  C  C++  Fortran  Modula-2  Pascal  PL/I

Name computers, if any, of which you have working knowledge of
assembly-language programming:


Name computers, if any, of which you have a working knowledge of the floating-point
architecture, such as number representations, rounding methods, etc.


What do you want to get out of this course?









SYLLABUS

     The following details indicate the scope of the lec-
tures and point to references.  While the focus will be on
general principles, particular areas of interest expressed
by attendees will also be discussed.

                               Lecture Plan

Lecture number                             Topic

1                Prospective overview from an elementary level (videotaped)
2-21             See details below
22-23            Student project reports
24               Retrospective overview from an advanced level (videotaped)


1.  Environment for computing
Gresham's Law: "the bad drives out the good"; similarly the fast
        drives out the slow even if the fast is wrong.
Computer Users: sophisticated but numerically naive.
        Many layers of software from diverse sources.
Portable software in the public domain: LINPACK, EISPACK, ...
Commercial Libraries: IMSL, NAG, ... [Rice 1983]
Standards: for programming languages but little else.
Proprietary rights, copyright, patent, trade secrets.
        L. J. Comrie's trick
        Execute-only code, copy protection.
Liability and exculpation: law of torts
        Reliability means conformity to reasonable expectations
        Marketing claims vs. absolute and negligent liabilities
        Obligation to correct errors; change vs. compatibility

2.  Comparisons of floating-point arithmetics in computers
Conventional floating-point arithmetics: [Sterbenz 1975]
        Radix, Word size, Range/Precision trade-off
                Why is binary best?  When is decimal best?
Commercially significant machines:
        IBM 7090/7094, 360, 370, PC's
        Univac 11xx
        CDC 6600, Cyber 1xx, 205
        Burroughs 6500
        DEC PDP-11, VAX
        Crays
        Hewlett-Packard 3000
        Calculators by Hewlett-Packard and Texas Instruments
IEEE Standard 754:
        Chips: Intel 80x87, Motorola 6888X, ATT 32106
                National 32081, Weitek 1164/5, AMD 29027,
                Fairchild Clipper, Analog Devices 3212/22
                BIT, Inmos T800, ...
        Micros: IBM PC's, Sun-3 and Sun-4, Apple Macintosh II, ...
        Mainframes: ELXSI 6400, HP Spectrum
Nonconventional representations of real numbers:
        Logarithmic (Radix 1)
        Variable exponent width - why are they all bad ideas? Demmel 1987
                Matsui-Iri 1983, Hamada 1987
                Symmetric level-index arithmetic, Clenshaw-Olver 1987
        Floating-slash rational arithmetic: Kornerup-Matula 1987
Interval arithmetic: Moore, Herzberger and Alefeld
        what it is and what good it is, briefly
Multi-precision schemes:
        Conventional double and extended precisions
        Numerical Turing language - Hull 1987
        MP package for Fortran - Brent
        MACSYMA Bigfloat - Fateman
        Kulisch-Miranker paradigm, IBM ACRITH package 1983
        Kahan's indefinitely extended IEEE 754/854
Implementation of algebraic operations +-*/ rem sqrt:
        Fast carry propagation: fast prefix calculation, VLSI
        Floating-point add/subtract, shifters, normalizers
        Multiplication: iterative, Booth re-encoding,
                Wallace tree of carry-save adders
        Division: higher radix [Taylor]; preconditioned
                iterative using fast multiply, can be grubby [Cray]
                can be cleaned up and tested [Kahan 1987]
                Why keep (division time)/(multiplication time) < 4
                Interruptible remainder
        Square root by hardware, by software, or without division
        Rounding in theory: guard-round-sticky bits
                rerounding without double rounding (C. Lee 1988)
                carry-free roundings (chop, jam, ROM)
        Rounding as practiced on commercially significant machines
        DPROD, or multiply-add primitives
Conversions: of precisions, float-integer, BCD-binary
Exceptions: Traps/Faults, Flags, Defaults, Presubstitution,
        Invalid operation, Division by zero, Over/Underflow,
        as practiced on commercially significant machines
                partial overflow on Cray
                partial underflow on CDC Cyber 1xx
                gradual underflow in IEEE 754/854
                NaNs, infinities, and reserved operands
        Range extension via over/underflow tallies [Kahan 1966]
Precision doubling, and further extensions [Pichat 1972]
Ultra-wide multiplication (Karatsuba, Cook, Schoenhage)

3.  Models of floating-point arithmetic for programmers
Purpose of models: promote correctness proofs and portability
Axioms: non-categorical axiomatization is too loose
        32 by Wijngaarden 1966
        20 by Brown 1981
Axioms: over-axiomatization is too slow
        Kulisch's ringoids, semi-morphism, optimality
Descriptions:
        crude inequalities, Wilkinson 1959
        parameterization, Sterbenz 1973
        environmental inquiries in Fortran, ADA, ...
Prescriptions:
        Computer manufacturer's manuals
        IEEE standards 754/854
Tests for correct floating-point arithmetic:
        Test batteries: the IEEE 754 test-vectors tape
        Test patterns: Schryer 1980, NAG 1987
        Inquiries and consistency tests: MACHAR Cody 1980, PARANOIA Kahan 1986
        P-adic tests for correctly rounded sqrt: Kahan 1966
        P-adic tests for correctly rounded division: Kahan 1987
Identities that persist despite roundoff
        Weak monotonicity
        Cancellation:
                p-q is exact if 1/2 <= p/q <= 2
                When does integer m = (m/n)*n rounded?
                When does x = (x+y) - y rounded ?
                When does sqrt(x*x) = |x| ?
        Avoiding drift:
                during conversion, Matula 1968
                during iteration; Knuth-Reiser proof
                Harry Diamond's theorem
Error analysis:
        Conventional "forward" analysis (Hotelling)
        "Backward" analysis (Turing, Givens, Bauer, Wilkinson)
        Combined approaches: the web of relationships
        Misconceptions corrected:
                Rounding errors are NOT infinitesimal
                Rounding errors are NOT random variables
                Unstable computations are often free from cancellation
                Cancellation is often a GOOD thing in equation solving,
                        extrapolation, divided differences
                Simulation is easy, prediction difficult - removable singularities
                Backward error analysis is an explanation, not a license:
                        accuracy not transitive under composition
                Precision does not limit accuracy (Dekker, Pichat)
Why double+ precision is a good rule of thumb:
        double zeros, normal equations, divided differences
What "unstable", "ill-conditioned", "ill-posed" mean; how to
        help an ill-posed problem get well
Interval arithmetic vs. Running Error-Analysis
Examples:
        When arbitrarily high precision might not help
        Triangle with given side-lengths
        Quadratic and Cubic equations: betrayed by books
        Polynomial equations in general
        Linear systems of equations and eigenproblems
        Numerical quadrature and infinite series
        Trajectory problems and boundary value problems

4.  Comparison of computers' floating-point register archi-
tectures
One accumulator: IBM 7094, DEC PDP 10, GE/Honeywell 635
Stack-top: B5500, HP 3000
Small stack: HP calculators, Inmos T800, Intel 80x87
        indefinitely long stack intended for 8087
Effect on procedures' argument passing and register saving
Extended width:
        internally in HP calculators
        Kulisch-Miranker's hidden Super-Accumulator
        in registers: IBM 7094, DEC PDP-10, GE 635
        IEEE 754: 8087, 68881, WE 32106
Orthogonal registers: IBM 370, CDC 6600, NS 32081, Cray,
        DEC VAX, ELXSI 6400, Weitek, AMD, Analog Devices, etc.
Vector registers: Cray; simulated vectors in CDC Cyber 205
        Fast vector arithmetic without vector registers.
Influences on programming languages' expression evaluation:
        Precisions available
        Mixed-precision expressions:
                Fundamentalist's Fortran for IBM 370, DEC VAX
                Widest available for C on PDP-11; 8087; 68881
                "scan for widest" seems best [Corbett]
Conflict between compiler's pseudo-machine and actual registers

5.  Real elementary transcendental functions
Argument reduction for log, exponential, and trig functions
        How to simulate 2/pi to infinite precision, or not.
CORDIC algorithms = argument reduction carried to its limit
Polynomial, rational, and other "best" approximation
        The Remez algorithm, with compensation for roundoff
Interpolation in BIG tables: IBM 370, DEC VAX VMS
        J. C. P. Miller's trick to suppress roundoff: Gal
"Kernel function" approach; monotonicity achieved thereby
Special properties for special functions:
        Weak monotonicity: EXP(X+Y**2) >= EXP(X)
        Cardinal values: cos(0) = 1, |cos(x)| <= 1 .
        Huge arguments: cos(3083975227) = -.00000000007486714572...
        Removable singularities: does 0**0 = 1 ?
        Unremovable singularities: ln(0)=-infinity; ln(-1) is NaN.
Accuracy claims and proofs: the Table-Maker's Dilemma
Testing: Cody 1980, Liu 1987

6.  Complex Arithmetic
Mixed Real and Complex expressions
Multiply/Divide despite exceptions in concurrent environments
Different approaches in infinity
The sign of zero: functions like sqrt on slitted domains
Accurate elementary transcendental functions
Peter Tang's complex Remez algorithms

7.  Exception handling
An Exception is not an Error unless handled badly
Why we must neither always ABORT nor always merely CONTINUE
Retrospective Diagnostics resolve the dilemma.
Let's not get enmeshed in the operating system's kernel.
Concurrent/parallel/pipelines machines restrict the options:
        how can we pause and resume during debugging?
Infinity and NaN/Indefinite/Reserved operand allow CONTINUE option.
Presubstitution, and some of its applications.
Flags permit recovery from ill-presubstituted exceptions.
Circumventing over/underflow in extended products/quotients:
        Kahan's 1966 implementations on 7094 and 360.
        Barnett's 1987 implementations on VAX and Sun-3.
The problem of Scope:
        solved without language support by Apple SANE 1986
        with minimal language support:
                APL's localization of system variables
        proper language support achieves minimal run-time cost:
                diminish incidence of forgotten SAVE/RESTOREs
                avoid unnecessary SAVE/RESTORE in leaf subprograms.
Comparisons with:
        "ON condition" in PL/I and Basic
        "drop-through" in ADA.
        "Error-Handlers" in DEC VMS Fortran.
        Signals thrown and caught by Unix.
        Augmented data-structures for Data-Flow.

8.  Numerical lapses in standard programming languages:
Misguided attempts:
        to define the semantics of arithmetic
                Brown's model in ADA, Fortran 8x, ANSI C
                Ambiguous environmental inquiries
        to define occurrences like 1.0/0.0 and 0.0**0.0 as Errors.
        to ignore integer overflow
Running roughshod over parentheses (C)
Moving code out of loops despite side-effects.
Other would-be optimizations like
        misuse of register variables causing schizophrenia
        misuse of "0-x" instead of "-x", or
                "-(x-y)" instead of "y-x" , and
                " x-y > 0" instead of "x > y".
Pascal's fixed array-dimensions, lack of double precision,
        inadequate provision for precompiled libraries
C's awkward multi-subscripted arrays inside procedures
Fortran's unprotected argument lists
Inability to overload arithmetic operators to introduce:
        complex arithmetic
        matrix arithmetic
        interval arithmetic
Awkward treatment of
        "Eureka" problem (Break)
        "Lost search-party" problem (compare Exception Handling)
        Mixed-precision expressions (Corbett)
        Re-running a program with higher precision
        Mixed real-complex expressions
        Inexact constants, inaccurate conversions
        Output too wide for its field
What else should compiler writers know about floating point:
        Farnum's 1988 paper


------------------------------

End of NA Digest
**************************
-------

Reposted by

-- 
Kenneth R. Jackson,                   krj@csri.toronto.edu (csnet)
Department of Computer Science,       uunet!csri.toronto.edu!krj (uucp)
University of Toronto,                krj@csri.toronto.cdn (ean x.400)
Toronto, Canada  M5S 1A4              krj%csri.toronto.edu@relay.cs.net (arpa)
(416) 978-7075                        krj@csri.utoronto (bitnet)