accuracy of mathematical functions

[Thomas Lumley]

What does the standard guarantee about the accuracy of eg the
trigonometric functions? There is obviously an implementation-dependent
upper bound on the accuracy, since the answer is stored in a double,
but is this bound actually achieved?

First, a simple example

Suppose I want the arc cosine of -1, and use acos(-1.0) to compute it.
The correct answer to 25 hexadecimal digits is
-3.243F6A8885A308D313198A2E. This falls between two exactly
representable doubles, call them A and B. Does the standard guarantee
that I get one of A and B as the answer (or even better, than I get the
one that is closer to the right answer)?

A more difficult example.

Suppose I want sin(exp(100)). The value of exp(100) is not exactly
representable in a double, and in fact, all numbers in the range
exp(100)-pi to exp(100)+pi have the same closest representation as a
double. Does this authorize the implementation to give any answer it
likes (in the interval [-1,1]), or is it required to give a value close
to the sine of the number it uses to represent exp(100)? Requiring the
right answer here would be a lot of work.


Thomas Lumley (thomas at drizzle dot net)

 

[E. Robert Tisdale]

What does the standard guarantee
about the accuracy of e.g. the trigonometric functions?

Almost nothing.
But other standards such as IEEE 754 may govern accuracy.

There is obviously an implementation-dependent upper bound on the accuracy
since the answer is stored in a double,
but is this bound actually achieved?

First, a simple example

Suppose [that] I want the arc cosine of -1, and use acos(-1.0) to compute it.
The correct answer to 25 hexadecimal digits is

-3.243F6A8885A308D313198A2E.
>
This falls between two exactly representable doubles, call them A and B.
Does the standard guarantee that I get one of A and B as the answer?
(or even better, than I get the one that is closer to the right answer)?

A more difficult example.

Suppose I want sin(exp(100)).
The value of exp(100) is not exactly representable in a double and,
in fact, all numbers in the range exp(100)-pi to exp(100)+pi
have the same closest representation as a double.
Does this authorize the implementation to give any answer it likes
(in the interval [-1,1]), or is it required to give a value close
to the sine of the number it uses to represent exp(100)?
Requiring the right answer here would be a lot of work.

The typical implementation [of IEEE floating-point arithmetic]
gets floating-point addition, subtraction, multiplication, division
and square root to an accuracy of 1/2 unit last place (ulp)!

In general, the trigonometric and transcendental functions
are [practically] impossible to get to better than 1 ulp.
The accuracy of trigonometric and transcendental function
is usually specified only for a small domain (i.e. -pi to +pi)
but, in fact, most implementations give astounding accuracy
over a *much* larger domain.
I don't think that any standard specifies accuracy for sin(exp(100)).

 

[Dik T. Winter]

> What does the standard guarantee about the accuracy of eg the
> trigonometric functions?

Nothing. It is a quality of implementation issue.

> Suppose I want the arc cosine of -1, and use acos(-1.0) to compute it.
> The correct answer to 25 hexadecimal digits is
> -3.243F6A8885A308D313198A2E. This falls between two exactly
> representable doubles, call them A and B. Does the standard guarantee
> that I get one of A and B as the answer (or even better, than I get the
> one that is closer to the right answer)?

No such guarantee.

> Suppose I want sin(exp(100)). The value of exp(100) is not exactly
> representable in a double, and in fact, all numbers in the range
> exp(100)-pi to exp(100)+pi have the same closest representation as a
> double. Does this authorize the implementation to give any answer it
> likes (in the interval [-1,1]), or is it required to give a value close
> to the sine of the number it uses to represent exp(100)? Requiring the
> right answer here would be a lot of work.

It may give any answer it likes. I may note that I know of implementations
of the sine function that can on occasion give a result that is slightly
outside the range [-1,1]. This is also allowed. (A naive implementation
of the Cody-Waite algorithm might do that.)

 

[P.J. Plauger]

What does the standard guarantee about the accuracy of eg the
trigonometric functions? There is obviously an implementation-dependent
upper bound on the accuracy, since the answer is stored in a double,
but is this bound actually achieved?

I know of microprocessors that get 0 ulp (units in the least significant
place) essentially all the time. They do so by computing polynomial
approximations to an extra 12 bits in firmware. But those of us who
write portable software-only math functions consider 1 ulp the gold
standard. Your typical C library has math functions that are *much*
worse.

First, a simple example

Suppose I want the arc cosine of -1, and use acos(-1.0) to compute it.
The correct answer to 25 hexadecimal digits is
-3.243F6A8885A308D313198A2E. This falls between two exactly
representable doubles, call them A and B. Does the standard guarantee
that I get one of A and B as the answer (or even better, than I get the
one that is closer to the right answer)?

The C Standard doesn't. Supplemental standards, such as IEEE 754
(aka IEC 60559 are more ambitious. In the best of all worlds,
every function should round to the nearest representable value
(or to the nearest even value if the true value is exactly in
the middle).

A more difficult example.

Suppose I want sin(exp(100)). The value of exp(100) is not exactly
representable in a double, and in fact, all numbers in the range
exp(100)-pi to exp(100)+pi have the same closest representation as a
double. Does this authorize the implementation to give any answer it
likes (in the interval [-1,1]), or is it required to give a value close
to the sine of the number it uses to represent exp(100)? Requiring the
right answer here would be a lot of work.

The C Standard doesn't say you can punt for sin/cos/tan of a large
argument, so in principle you have to return the nearest internal
representation corresponding to the function evaluated as if the
argument were exact. We do this with our latest library (in the
process of being packaged), as does Sun's latest C library. But
you're right, it's a lot of work; and many people don't think it's
worth the bother. We do the work anyway, for that handful of
people who know what they're doing and care about truly accurate
math functions.

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com

 

[Christian Bau]

In general, the trigonometric and transcendental functions
are [practically] impossible to get to better than 1 ulp.

Actually, getting to within 0.6 ulp for the trigonometric and
transcendental functions is not very difficult. Check any Java
implementation.

 

[Christian Bau]

It may give any answer it likes. I may note that I know of implementations
of the sine function that can on occasion give a result that is slightly
outside the range [-1,1]. This is also allowed. (A naive implementation
of the Cody-Waite algorithm might do that.)

An interesting mathematical problem is this: Implement sin (x) and cos
(x) in such a way that for all values of x,

double s = sin (x)
double c = cos (x)
double t = s*s + c*c;

produces a value with t <= 1.0, while maintaining highest possible
precision for sin and cos.

 

[P.J. Plauger]

In general, the trigonometric and transcendental functions
are [practically] impossible to get to better than 1 ulp.

Actually, getting to within 0.6 ulp for the trigonometric and
transcendental functions is not very difficult. Check any Java
implementation.

Which is almost certainly written in C, BTW.

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com