N1884=05-0144
2005-08-26

Further Restrictions on Special Math Functions

P.J. Plauger & Chris Walker
Dinkumware, Ltd.


INTRODUCTION

We've encountered several additional problems in implementing the special
math functions for TR1. You can argue that that's our problem, as
implementers, but we believe that the issue is not so simple. What we've
learned is that the problem areas we've encountered are also problem
areas for potential users of the functions, for reasons we explain below.
Thus, we recommend that the problematic domains, or whole functions in
two cases, be removed from the special math requirements.

The unifying principle at work here is what we've come to call the
*sensitivity* of functions. We follow the conventional rule that a
library function should treat each input value x as if it were exact.
(You don't know where it's been, in any case.) It is the job of the
implementers to return the nearest internal representation corresponding
to the actual value of the function for a given x. Nevertheless, it's
worth asking what happens to the output value if you change x by the
smallest possible amount -- conventionally 1 ULP (unit in the least
significant place). We take the magnitude of the larger of the two
changes, also expressed in ULP, as the *sensitivity* of the function
for the input value x.

It is well understood that large sensitivity causes problems for the
programmer. A small error in the input value can lead to surprisingly
large errors in the function value. It doesn't help if a library does
a perfect job if the program cannot guarantee perfect function inputs --
and most programs can't. Large sensitivity is not necessarily an excuse
for the implementers to do a sloppy job, however; if there's any chance
that the input is correct, then the most meaningful output should be
generated. But it should come as no surprise that the biggest
challenges to implementers typically occur where sensitivity is high.
So the community is not well served by demanding great accuracy for
domains that are not likely to be usable in a meaningful program.

The commonest cause of high sensitivity is a very large ratio of
first derivative to function value, |f'(x)/f(x)|. tgamma(x), for
example, is notoriously sensitive for large function values. But the
peculiarities of floating-point arithmetic create other cases as
well:

-- If a function has a zero value for an irrational argument value,
such as expint (one) or lgamma (lots), the region near the zero has
high sensitivity. The implementer can patch a small number of such
zeros, but an infitite number requires heroic measures that are
probably not worthwhile.

-- A cyclic function, such as sin(x) or a Bessel, gets progressively
more sensitive as the argument grows larger. Put simply, each ULP
change in the argument makes ever larger leaps around the unit circle,
so the sequence of function values becomes erratic. Techniques exist
for reducing arguments for many of the cyclic functions; whether
they're worth employing beyond a certain sensitivity is debatable.

-- A really fast moving function, such as assoc_laguerref for large
order, can actually toggle between plus and minus infinity with 1
ULP changes in the argument value. Determining the proper sign for
each function value can be both difficult and of only academic
interest.

The most useful specification for a math library would probably be
one that identifies regions of "high" sensitivity -- and provides
a defensible definition of "high" in the bargain. We at Dinkumware
may eventually produce such a document, but not right now. The
best we can offer at the moment is a few more specification changes
for the special math functions, to eliminate known problem areas.


SUGGESTED CHANGES

assoc_laguerre -- The m argument should also be implementation
defined for m >= 128.

assoc_lgauerref -- If n >= 75 AND m >= 75, IEEE 32-bit is pretty
meaningless (toggles between plus and minus infinity). For example,
the value of assoc_laguerref(127, 127, 41.85287) is -0.64511883e44,
with a sensitivity of 24920640. We can probably get the signs right
by storing a table of roots, but it's hard to imagine a program that
would make good use of the function values.

assoc_legendre -- The m argument should also be implementation
defined for m >= 128.

assoc_legendref -- If n >= 35 and m >= 35, IEEE 32-bit is pretty
meaningless (toggles between plus and minus infinity). For example,
the value of assoc_legendre(50, 35, 0.17249578) is -0.33007636e52,
with a sensitivity of 27800896.

conf_hyperg -- This function has so many sensitive subdomains that
it is essentially useless computationally (however pretty it is as
a unifying concept mathematically). It should be removed from TR1.

hermitef -- If n >= 47 and |x| < 10, IEEE 32-bit is pretty
meaningless (toggles between plus and minus infinity). For example,
the value of hermitef(127, 0.19673838) is -0.13458131e119,
with a sensitivity of 31215745. We can probably get the signs right
by storing a table of roots, but it's hard to imagine a program that
would make good use of the function values.

hyperg -- This function has so many sensitive subdomains that
it is essentially useless computationally (however pretty it is as
a unifying concept mathematically). It should be removed from TR1.

sph_legendre -- The m argument should also be implementation defined
for m >= 128.