13

Do the C and C++ standards require the math operations in math.h on floating points (i.e. sqrt, exp, log, sin, ...) to return numerically best solution?

For a given (exact and valid) input there can obviously in general not be an exact floating point output from these functions. But is the output required to be the representable value nearest to the mathematically exact one?

If not, are there any requirements on precision whatsoever (possibly platform-specific / in other standards ?), so that I am able to make worst-case estimates of calculation errors in my code? What are typical limits on numerical errors of modern implementations?

  • 6
    Basically no. IEEE-754 doesn't even specify most of these. There are some limits on the number of representable digits but there was a time when wonky proprietary floating-point implementations proliferated, and as a result the language leaves most everything about the floating-point results unspecified. That is the C spirit after all – doynax Jan 06 '14 at 08:35
  • 1
    @doynax: +1 An implementation needn't even follow IEEE 754 specification for floats/doubles. – legends2k Jan 06 '14 at 08:38
  • @doynax: What do you mean IEEE 754 does not specify most of these? IEEE 754-2008 table 9.1 recommends `exp`, `log`, `sin`, and more. – Eric Postpischil Jan 06 '14 at 20:17
  • @EricPostpischil: It would appear I am out-of-date, having only read IEEE 754-1985. From a cursory reading of the 2008 revision it seems that the accuracy of these functions is left unspecified aside from a few identities and the curious requirement to signal inexact results if-and-only-if they are inexact. This seems needlessly expensive to handle. Surely no-one would expect exact results from the transcendental functions anyway, with the possible exception of the exponential function? – doynax Jan 06 '14 at 21:12
  • @doynax: Section 9.2, where the table appears, says the functions should be correctly rounded. That means the error must be the minimum possible given the rounding mode and the format; in round-to-nearest mode, the nearest representable value must be returned (breaking ties with the usual rule). I agree, these results should not generally be expected except for those that have been demonstrated to be feasible (as by the [CRlibm project](http://lipforge.ens-lyon.fr/www/crlibm/). See my comment with [MSalter’s answer](http://stackoverflow.com/a/20946263/298225); this was changed post-committee. – Eric Postpischil Jan 06 '14 at 21:29

1 Answers1

4

No, and for good reason. In general, you'd need an infinite precision (and infinite time) to determine the exact mathematical result. Now most of the times you need only a few extra iterations to determine sufficient bits for rounding, but this number of extra bits depend on the exact result (simply put: when the result is close to .5 ULP). Even determining the extra number of iterations required is highly non-trivial. As a result, requiring exact results is far, far slower than a pragmatic approach.

MSalters
  • 173,980
  • 10
  • 155
  • 350
  • +1 Nice answer. Note: Agree with functions like `exp(), log(), sin()`, but I think `sqrt()` can always be answered to the nearest to the mathematically exact one. Wonder what IEEE-754 says about `sqrt()`? – chux - Reinstate Monica Jan 06 '14 at 17:25
  • 3
    @chux: Correctly rounded implementations of `exp` and `log` with reasonable performance are available in [CRlbim](http://lipforge.ens-lyon.fr/www/crlibm/). `sin` is available correctly rounded over [-π, +π]. IEEE 754 recommends that all these and some other functions (in table 9.1) be implemented with correct rounding. In committee, I opposed that because no viable correclty-rounded implementations existed (and still do not) for some of them, and the committee passed the standard without saying they should be correctly rounded, but it was changed in post-committee editing or balloting. – Eric Postpischil Jan 06 '14 at 20:13
  • There is a typo in my previous comment; “CRlbim” should be “CRlibm”. It is a very nice project to implement mathematical functions with correct rounding and good performance, and with proofs. – Eric Postpischil Jan 06 '14 at 21:31
  • Thanks - an impressive reference and table. My reference version of table 9.1 implies a recommended correctly rounded domain for `sin()/cos()` of (−∞, +∞). So not withstanding K.C. Ng's "Good to the Last Bit" work, it is still (2014) only [-π, +π]? – chux - Reinstate Monica Jan 06 '14 at 21:52
  • To clarify what Eric says about IEEE-754 slightly, the intent is not to recommend that the `math.h` functions be correctly rounded, but rather that correctly rounded versions be provided, possibly as well as versions that are not correctly rounded. – Stephen Canon Jan 10 '14 at 17:18