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)