			   Updates for the
			     Cody & Waite
	  Fortran Elementary Function Test Package (ELEFUNT)

			    [24-May-2002]

				  by

			  Nelson H. F. Beebe
		   Center for Scientific Computing
			  University of Utah
		  Department of Mathematics, 110 LCB
			 155 S 1400 E RM 233
		    Salt Lake City, UT 84112-0090
				 USA
	      Email: beebe@math.utah.edu, beebe@acm.org,
	    beebe@computer.org, beebe@ieee.org (Internet)
	       WWW URL: http://www.math.utah.edu/~beebe
		      Telephone: +1 801 581 5254
		FAX: +1 801 585 1640, +1 801 581 4148


=======
Preface
=======

To get started immediately, with minimal reading, skip down to the
section ``How to run the ELEFUNT tests''.


============
Introduction
============

This directory contains the Cody & Waite Fortran Elementary Function
Test Package (ELEFUNT) (see below for references).

There are three main programs:

	talldp (test all double precision)
	tallqp (test all quadruple precision)
	tallsp (test all single precision)


=======================================================
Updates for quadruple precision and IEEE 754 arithmetic
=======================================================

The tallqp program is new (23-May-2002): it was developed by applying
the author's Fortran precision translator package (program dtoq from
the distribution in

	ftp://ftp.math.utah.edu/pub/misc/dtosstod.tar.gz
	ftp://ftp.math.utah.edu/pub/misc/qtoddtoq.tar.gz

), and then fixing minor problems, extending FORMAT item widths,
changing comments and FORMAT statement text, and in a couple of
places, introducing higher-precision magic constants.  The random
number generator, QRAN(), has been updated to generate up to 116 bits
of randomness, sufficient for the 113-bit precision of 128-bit
IEEE-754-like implementations.

The original subroutine MACHAR in earlier tall*.f files has been
replaced by the 1988 version from reference Cody:1988:AMS below, and
available in

	ftp://netlib.bell-labs.com/netlib/toms/665.gz

That newer version updated MACHAR to handle IEEE 754 arithmetic, for
which previous versions failed.  However, that 1988 code is old enough
that it does not correctly deal with modern compilers that keep
expressions in registers: on Intel IA-32, HP/Intel IA-64, and Motorola
68000, these registers are wider than memory words, and offer higher
precision and wider range, foiling Cody's careful coding in MACHAR.
The solution has been to force storage into memory by wrapping 35
expressions in calls to a function, STORE().


====================================
Results for various computer systems
====================================

The logs/tall*.lst* files record results for several current
systems.  The names reflect the system, compiler, and where relevant,
optimization levels.  For example,

	tallqp.lst-intel-ia32-ifc-g

was produced on an Intel IA-32 architecture (Pentium III 600MHz dual
CPUs) with the Intel Fortran Compiler (ifc), with optimization level
-g.  Compilers on IA-32 may behave quite differently with -g and -O
optimizations: with -g, library routines are used for the elementary
functions, while with -O, inline instructions are generated to access
the Intel 80x87 elementary function floating-point instructions.


==========================================
Quadruple-precision background and caveats
==========================================

The quadruple-precision tests need some background explanation:

(1) IEEE 754 defines four floating-point formats:
	* 32-bit (single-precision),
	* single-extended precision (size unspecified, but > 32, and
	  in practice, not implemented in desktop processors),
	* 64-bit (double-precision), and
	* extended precision (size unspecified, but > 64: only three
          architecture families implement it: Intel IA-32 (80-bit),
	  HP/Intel IA-64 (82-bit), and Motorola 68000 (80-bit)).

(2) IEEE 754 does not define a quadruple-precision format, but several
    vendors IBM S/390 (hardware), Hal/Fujitsu SPARC64-GP (hardware),
    Sun Solaris SPARC (software), Lahey/Fujitsu lf95 on Intel IA-32
    (software), and Intel ifc on Intel IA-32 and HP/Intel IA-64
    (software)) have implemented identical representations:

	1 sign bit, 15 exponent bits biased by 16383, 1 hidden integer
	bit, 112 stored fraction bits

(3) IBM RS/6000 AIX and SGI MIPS IRIX compilers implement quadruple
    precision from paired doubles: both implementations have very
    serious problems that are described in a paper by Nelson
    H. F. Beebe and James S. Ball under review at ACM Transactions on
    Mathematical Software.

(4) HP/Compaq/DEC Alpha OSF/1 does not offer any quadruple-precision
    support at all.  The Alpha architecture defines 128-bit IEEE 754
    style arithmetic, called `X-floating', but no compilers for OSF/1
    or GNU Linux on Alpha implement it; I don't know whether any do on
    OpenVMS.

(5) None of the GNU compilers for Fortran, C, C++, or other languages
    provide any quadruple-precision support.

(6) Java defines a crippled subset of 32-bit and 64-bit IEEE 754, but
    offers no other floating-point formats. Java has BigInteger and
    BigDecimal packages for arbitrary-precision integer and
    fixed-decimal arithmetic: see

        http://www.math.utah.edu/~beebe/software/java/fibonacci/

    for details, and comments on the floating-point support in more
    than 30 programming languages.

(7) IBM S/360 hexadecimal floating-point arithmetic has been available
    on mainframes since 1967, implemented in hardware in the larger
    models, and in software in smaller ones.

(8) DEC VAX arithmetic defines a 128-bit quadruple-precision form,
    similar to IEEE 754, but lacking NaN, Infinity, and subnormals.
    Although there are 128-bit floating-point instructions defined in
    the VAX instruction set, I am unaware of any hardware
    implementations: all compilers that I used under DEC VAX VMS in
    the 1980s emulated those instructions by software trap handlers.

(9) None of the ANSI/ISO Fortran standards up to Fortran 95 define
    quadruple-precision support, so Fortran compiler vendors usually
    follow venerable IBM S/360 practice, with some variations, for the
    naming of elementary functions, the representation of numbers using
    the Q exponent letter, and the Q FORMAT items.

(10) C99 (ISO/IEC 9899:1999(E) Programming languages -- C) defines
    support for long double: number representation (L suffix),
    *printf()/*scanf() I/O format items, and elementary functions.
    However, as of the time of writing this [24-May-2002], none of the
    more than 50 C compilers implements C99, although several provide
    reasonably common support for long double.  However, C99 only
    recognizes float, double, and long double, with the expectation
    that these map to 32-bit, 64-bit, and either 80-bit or 128-bit
    representations.  A vendor such as HP, which has Fortran and C
    compiler support for four precisions (32-, 64-, 80-/82-, and
    128-bit) then has to invent new language keywords to distinguish
    between the two longer precisions.

Thus, tallqp.f assumes S/360 quadruple-precision conventions.  The
only systems currently available to the author on which it can be
tested are Sun Solaris SPARC (f77, f90, f95), and Intel IA-32 ifc.
Our Lahey/Fujitsu lf95 compiler, which also has 128-bit support in
software, stopped working after an operating system upgrade, and has
yet to be updated to a working version.  My guest account on IBM S/390
GNU/Linux expired a few months ago, so that mainframe access
disappeared.


==========
What's new
==========

Details of the changes can be found in the companion ChangeLog file,
but can be briefly summarized as follows:

	* Support for quadruple-precision (tallqp.f) as well as the
	  normal single- (tallsp.f) and double (talldp.f) tests.

	* Changes to the environmental inquiry routines, machar(), to
          support IEEE 754 arithmetic, including those systems with
          registers that are wider than memory words.

	* Source code can be compiled with Fortran 77, Fortran 90,
	  Fortran 95, and High-Performance Fortran compilers.

	* Output listings for several current systems are available in
          ./logs.

	* New top-level Makefile to control building and testing.

	* New helper script (./elefunt.sh) to run the tests with
	  several different Fortran compilers, giving the output
	  listing distinctive names that identify the system,
	  compiler, and optimization level


============================
How to run the ELEFUNT tests
============================

On most systems, you can build and test the programs like these
examples. The default target, all, runs only single- and
double-precision tests: use the all3 target to also include
quadruple-precision tests.

	# Build and run tallsp and talldp with default Fortran and -g
	# optimization level:
	make

	# Build and run just tallqp with defaults:
	make tallqp

	# Build with a different compiler (GNU g77):
	make FC=g77

	# Need special option for Sun Fortran 90 (and 95) to avoid
	# premature process termination from intentional arithmetic
	# exceptions:
	make all3 FC=f90 XFFLAGS='-ftrap=%none'

	# Build with Intel Fortran Compiler, suppressing the mass of
	# complaints about Fortran 66 and Fortran 77 old-style syntax
	# (-w means no warnings):
	make all3 FC=ifc XFFLAGS=-w


====================
Result summarization
====================

Because the output files, tall*.lst, are voluminous, a small data
reduction program is provided, to be used like this:

	awk -f elefunt.awk tallsp.lst

Its header comments show how to modify its filtering with additional
command-line options.  By default, it reduces a 1450-line output file
to just a few lines of results from tests with larger-than-expected
errors.


======================
The origins of ELEFUNT
======================

The original ELEFUNT code was developed for the book

@String{pub-PH                  = "Pren{\-}tice-Hall"}
@String{pub-PH:adr              = "Upper Saddle River, NJ 07458, USA"}

@Book{Cody:SME80,
  author =       "William J. {Cody, Jr.} and William Waite",
  title =        "Software Manual for the Elementary Functions",
  publisher =    pub-PH,
  address =      pub-PH:adr,
  pages =        "x + 269",
  year =         "1980",
  ISBN =         "0-13-822064-6",
  LCCN =         "QA331 .C635 1980",
  bibdate =      "Tue Dec 14 23:28:38 1993",
}

and made available via netlib.  There are several technical reports
describing ELEFUNT results for various platforms listed in
W. J. Cody's complete publication bibliography at

	http://www.math.utah.edu/pub/bibnet/authors/c/cody-william-j.bib

One of ELEFUNT's routines was later described, with updates for
support of IEEE 754 arithmetic, in:

@String{j-TOMS                  = "ACM Transactions on Mathematical Software"}

@Article{Cody:1988:AMS,
  author =       "W. J. Cody",
  title =        "{Algorithm 665}: {MACHAR}: {A} Subroutine to
                 Dynamically Determine Machine Parameters",
  journal =      j-TOMS,
  volume =       "14",
  number =       "4",
  pages =        "303--311",
  month =        dec,
  year =         "1988",
  CODEN =        "ACMSCU",
  ISSN =         "0098-3500",
  bibdate =      "Sun Sep 04 22:33:58 1994",
  URL =          "http://www.acm.org/pubs/citations/journals/toms/1988-14-4/p303-cody/",
  acknowledgement = ack-nhfb,
  keywords =     "algorithms",
  subject =      "{\bf G.4}: Mathematics of Computing, MATHEMATICAL
                 SOFTWARE, Portability. {\bf G.1.0}: Mathematics of
                 Computing, NUMERICAL ANALYSIS, General, Computer
                 arithmetic.",
}

That was followed five years later by a description of a package for
testing the complex elementary functions:

@Article{Cody:1993:ACP,
  author =       "W. J. Cody",
  title =        "{Algorithm 714}: {CELEFUNT}: {A} Portable Test Package
                 for Complex Elementary Functions",
  journal =      j-TOMS,
  volume =       "19",
  number =       "1",
  pages =        "1--21",
  month =        mar,
  year =         "1993",
  CODEN =        "ACMSCU",
  ISSN =         "0098-3500",
  bibdate =      "Mon Sep 05 09:15:25 1994",
  URL =          "http://www.acm.org/pubs/citations/journals/toms/1993-19-1/p1-cody/",
  abstract =     "This paper discusses CELEFUNT, a package of Fortran
                 programs for testing complex elementary functions.",
  acknowledgement = ack-nhfb,
  keywords =     "algorithms; measurement; performance",
  subject =      "{\bf G.4}: Mathematics of Computing, MATHEMATICAL
                 SOFTWARE, Certification and testing. {\bf G.1.0}:
                 Mathematics of Computing, NUMERICAL ANALYSIS, General,
                 Numerical algorithms.",
}
