accuracy of mathematical functions

Discussion in 'C Programming' started by Thomas Lumley, Jul 5, 2005.

  1. 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)
     
    Thomas Lumley, Jul 5, 2005
    #1
    1. Advertising

  2. Thomas Lumley wrote:

    > 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)).
     
    E. Robert Tisdale, Jul 5, 2005
    #2
    1. Advertising

  3. In article <> "Thomas Lumley" <> writes:
    > 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.)
    --
    dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
    home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/
     
    Dik T. Winter, Jul 6, 2005
    #3
  4. Thomas Lumley

    P.J. Plauger Guest

    "Thomas Lumley" <> wrote in message
    news:...

    > 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
     
    P.J. Plauger, Jul 6, 2005
    #4
  5. In article <daf11o$i4q$>,
    "E. Robert Tisdale" <> wrote:

    > 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, Jul 6, 2005
    #5
  6. In article <>, "Dik T. Winter" <>
    wrote:

    > 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.
     
    Christian Bau, Jul 6, 2005
    #6
  7. Thomas Lumley

    P.J. Plauger Guest

    "Christian Bau" <> wrote in message
    news:...

    > In article <daf11o$i4q$>,
    > "E. Robert Tisdale" <> wrote:
    >
    >> 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
     
    P.J. Plauger, Jul 6, 2005
    #7
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. =?Utf-8?B?Q2hhcmxlc0E=?=

    Accuracy and CSS

    =?Utf-8?B?Q2hhcmxlc0E=?=, Jan 16, 2006, in forum: ASP .Net
    Replies:
    5
    Views:
    672
  2. PercyP
    Replies:
    5
    Views:
    560
    PercyP
    Oct 20, 2004
  3. Replies:
    15
    Views:
    631
    David Harmon
    May 5, 2005
  4. Replies:
    21
    Views:
    664
  5. Takeshi NISHIMATSU

    lacking in mathematical functions

    Takeshi NISHIMATSU, May 6, 2009, in forum: Ruby
    Replies:
    0
    Views:
    117
    Takeshi NISHIMATSU
    May 6, 2009
Loading...

Share This Page