Subject: Comments on Elementary Functions from Abraham Ziv

From: Ronald F Boisvert <[email-masked]>
Date: Fri, 28 May 1999 11:22:29 -0400

Abraham Ziv of IBM Haifa has provided some comments on the minutes 
of the March Numerical Working Group subcommittee meeting regarding 
the computation of elementary functions.  Attached is the relevant
section of those minutes, along with his comments (delimited by
***BEGIN HAIFA*** and ***END HAIFA***).

============================================================================

4. Math library specification.

a The current specification for java.lang.Math says, roughly, to use the
  algorithms of fdlibm translated straightforwardly to Java. It is not clear
  who, if anyone, does this. In practice, Java programs which reference
  elementary functions are not producing the same results on all Java
  platforms.

b In JDK 1.2 Sun is distributing Java wrappers for fdlibm in C for
   java.lang.Math.

c Brophy announced that he has produced a complete translation of Sun's fdlibm
  in Java. This is posted at http://www.vni.com/corner/garage/grande/.
  Numerics Working Group members were invited to test this library. Moler
  agreed to test this implementation.

d Moreira announced that Abraham Ziv and colleagues at the IBM Haifa Labs have
  developed a suite of correctly rounded math functions for IEEE arithmetic
  in C.  These functions are claimed to cost only about 10% more than fdlibm
  in computing time. (This would be remarkable.) The Haifa researchers
  claim to have a proof of correctness for their algorithms, but it is
  proprietary. Several present expressed skepticism about proofs that are
  not public.

***BEGIN HAIFA***
Usually our routines are significantly faster than the
corresponding routines of fdlibm.
The proprietary parts are the details, not the main
methodology which was published in the open literature.
The proof of correctness (namely - yielding correctly
rounded results are already in the paper:

Abraham Ziv, "Fast Evaluation of Elementary
Mathematical Functions with Correctly Rounded Last Bit",
ACM Transactions on Mathematical Software,  Vol. 17,
(1991), pp. 410-423.

Only the error estimations are proprietary - in the same
way that libraries which are declared to yield 1 ULP do
not provide proof details.
***END HAIFA***

e Unfortunately, the footprint of the Haifa library seems large. The
  compiled binary library occupies about 250Kb (and compresses to 185 Kb). One
  of the reasons for its size is the inclusion of large tables. (The source
  form of the atan table is 500 kB, for example.)  It is probably possible
  to trade space for time in the algorithm, but it is not clear, say, how
  much slower a library of half the current size would run.

***BEGIN HAIFA***
The size of the atan tables is exaggerated here.
It is true that the file which includes the tables for the
atan is 500k bytes long. However, one should note that
each table is written in the file twice. Once in big
endian form and once in little endian form (because of
portability requirements. After compilation
only one version of the table is included in the code.
Moreover, the atan table is the largest in the library).
Also, note that there is no atan2 tables file. The same
file is used by both the atan and the atan2 routines.
So, it is not that terrible.
In view of the current price of computer memory
(which keeps dropping) memory space is not worth
trading for time.
***END HAIFA***

f Greg Tarsey's group at Sun (including KC Ng) are working on a similar
  library, which is not yet complete. It is based on slightly different
  algorithms than those used by the Haifa group. (The Haifa code works
  in fixed space, which their proofs show to be sufficient, while the Sun
  codes repeat computations in increasing precision as necessary until the
  result is satisfied.) The Sun codes run very roughly three times longer
  than fdlibm at present, but probably could be sped up considerably.

***BEGIN HAIFA***
It has been planned - from the very beginning -
to introduce into the routines a repeated loop which
will let the computations stop only when accuracy
conditions are fully met. This addition however was
planned to be introduced in a second version of the
library, and SUN's people were informed about this
intentions, - and this is still the case (we are waiting
for funding). Despite the above the probability to
obtain a non correctly rounded result with the present
version is lower than 2**(-600), i.e. much less
than the probability for a machine error.
***END HAIFA***

g Possible specification for strictfp mode: produce the correctly rounded
  result.
  Pros: This is the ideal definition.
  Cons: The Java code for this does not yet exist. Existing algorithms would be
        much too large for inclusion in the Java core. The resulting code
        might run too slowly (a factor of 2 is unacceptable). The code would
        be difficult to test.

***BEGIN HAIFA***
The code for correctly rounding library is not more
difficult to test than any other code. In both cases only
finite test files can be prepared (and we are willing to
supply them if required). However test files for a
correctly rounding library can be prepared independently,
by every one, applying known mathematical
principles. Such test files will be as correct and as good
as those coming from a 'standard' source. One should
be aware of the fact that the correctly rounding library
(like any other math library) must be compiled in an
appropriate way in order to achieve bit by bit
compatibility on various platforms. Thus, for instance,
no change of order of performance of arithmetic
operations, which is sometimes introduced by
optimising compilers, is permitted. So, a too high
optimisation level should be avoided.
***END HAIFA***

h Possible specification for strictfp mode: mandate use of VNI's native
  Java fdlibm.
  Pros: This code is available now. It is reasonably fast and has a small
        footprint.  It produces results that are usually within one unit in
        the last place (ulp) of the correctly rounded result.
  Cons: An operational definition like this would mandate incorrect results.
        Incorporating future imporvements in the algorithms yielding better
        results would change the behavior of existing Java class files.

***BEGIN HAIFA***
According to our tests, fdlibm routines do
not provide 1 ulp error bound (although it is close to this).
***END HAIFA***

i Possible specification for default mode: specify largest acceptable relative
  error in result for each function (e.g. 1 ulp), being sure that the bounds
  were satisfied by fdlibm.
  Pros: This allows flexibility of implementation. It allows (guarded) use of
        specialized hardware instructions for sqrt, sine, cosine, etc. It
        allows improvements in algorithms producing more accurate results
        to be trivially accommodated.
  Cons: Careful testing of fdlibm would have to be done to insure it satisfied
        the criteria; this might be difficult to do rigorously, but could be
        done satisfactorily by relying on comparisons with a multiprecision
        system like Maple or Mathematica. Exhaustive testing in pure Java
        would be difficult.

***BEGIN HAIFA***
Such a specification means that we give up the
principle of perfect compatibility among all platforms!
***END HAIFA***

j Tim Lindholm asked the Working Group to draft a proposal before April 9.

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