10 ms·

Intel Underestimates Error Bounds by 1.3 Quintillion

424 points by ghusbands 12y ago

CurtHagenlocher 12y ago
Of possible interest, I researched a related issue once and here's what I found: The particular approximation used by the FPU instructions fsin and fcos is valid in the domain [-π /4, π /4]. Inputs outside this domain are reduced by taking the remainder when dividing by π /4 and then using trigonometric identities to fix up values in the other three quadrants (ie sin(x) = sin(π /2 – x) for x in [π /4, π /2]. The Intel FPU is “off” because it uses an approximation for π /4 when performing the reduction. In the mid 90s, some engineers at AMD figured out how to do a more precise reduction and implemented this reduction in the K5 FPU. AMD went back to being bit-for-bit compatible with the x87 behavior in the very next generation of their CPUs, presumably because it either broke or appeared to break too many applications.
- CurtHagenlocher 12y ago
  I just noticed that Bruce has moved from Valve to Google... . Interesting!
- sillysaurus3 12y ago
  It's not possible to implement a more precise fsin() without breaking apps? One scenario I can imagine an app breaking is if someone executed fsin(x), pasted the result into their code, and then tested some other float against it using the == operator. I suppose a more common breakage scenario would be if someone did something like floor(x + fsin(x)), since that'd very slightly change where the floor "triggers" with respect to the input. But how could that cause an app to break? Every scenario I can think of comes back to someone using == to test floats, which everyone knows is a no-no for dynamic results, such as the output of a function like fsin(). I guess a programmer wouldn't really expect the result of "x + fsin(x)" to ever change for a certain value of x, so maybe they executed it with that certain value, stored the result, and did an equality comparison later, which would always work unless the implementation of fsin() changed, and the thought process went something like, "This should be reliable because realistically the implementation of fsin() won't change." Can anyone think of some obvious way that making fsin() more accurate could cause serious breakage in an app that isn't testing floats with the equality operator?
  - Guvante 12y ago
    > One scenario I can imagine an app breaking is if someone executed fsin(x), pasted the result into their code, and then tested some other float against it using the == operator. While this sounds grandious, thanks to aggressive inlining this could be the case even if you wrote in your logic: sin(x) == sin(constant) Or even the following depending on your choice of delta. abs(sin(x) - sin(constant) < delta
    - sillysaurus3 12y ago
      Is it even valid for a compiler to transform sin(constant) into a hardcoded constant embedded into the executable? I guess since sin(x) is computed by an FPU in hardware, then it's inherently "dynamically linked." But in that case, why is it valid for a compiler to transform sin(constant) into a hardcoded result? Is there some standard somewhere which says that the implementation of sin/cos/etc can't ever change, and therefore it's valid for the compiler to do that? Or do compilers just make an assumption that the implementation of sin(x) won't ever change?
      - darkmighty 12y ago
        It breaks binary compatibility. The compiler substitutes some value (e.g. the broken value), so if the sin(x) evaluation changes (to the correct one), then binary execution will fail.
        DannyBee 12y ago
        Why would the compiler substitute the broken value unless the compiler was broken? Compilers generally don't use specialized local architecture instructions to compute their results, because they can't (It would break cross-compiling or targeting different platforms in a lot of cases) In fact, it's more likely the other way around: Compiler substitutes proper value. actual calls on sin(x) on the local platform give wrong answer.
      - DannyBee 12y ago
        "but if it's dynamic, is it even valid for a compiler to transform sin(constant) into a hardcoded constant embedded into the executable?" Yes. I honestly don't remember all the rules here, but it does in fact do so, legally.
luu 12y ago
This is one of those things that's well known by people who spend a lot of time doing low-level nonsense that really ought to get wider press. You really don't want to use x86 hardware transcendental instructions. I've seen people try to measure the effectiveness of their approximation by comparing against built-in hardware, but their approximation is probably better than the hardware! If you're not writing your own assembly, you'll be ok on most modern platforms if you just link with "-lm" and use whatever math library happens to be lying around as the default, but it's still possible to get caught by this. On obscure platforms, it's basically a coin flip. I used to work for another CPU vendor, and when we implemented more accurate sin/cos functions, some benchmarks would fail us for getting the wrong result. Turns out those benchmarks hardcode a check for a couple results, and that check is based on what Intel has been doing for ages. My recollection is that we had a switch to enable the extra accuracy, but that it wasn't enabled by default because it was too much of a risk to break compatibility with Intel. If that sounds too risk averse, there's a lot of code out there that depends on your processor precisely matching Intel's behavior on undefined flags and other things that code has no business depending on. It's basically the same thing Raymond Chen is always talking about, but down one level of abstraction. I've seen what happens if you just implement something that matches the "spec" (the Intel manual) and do whatever you want for "undefined" cases. You get a chip that's basically useless because existing software will crash on it.
- JoeAltmaier 12y ago
  The problem is not so much "Intel stuff is old and broken" as it is "Intel wrote the standard and there's no process to move it forward"?
  - sillysaurus3 12y ago
    OpenGL is a great example of how painful it is to write apps that are capable of not breaking when the standard it relies on changes. You basically end up with capability checks or version checks which cause your code to become complicated. I suppose there could've been a function that returned the implementation version of fsin/fcos/etc, which was incremented whenever Intel changed the implementation. That way benchmarks could test against that particular version. But it'd be hard to come up with a precise, consistent way to implement that. For example, do you increment the version whenever any transcendental function changes? If you do that, then you get benchmark code riddled with special cases for very specific versions, which someone has to constantly look up in documentation somewhere. You'd basically have to memorize the effect of every version increment. On the other hand, you could implement a more elaborate scheme, like having a version number for each of the transcendental functions, but then you'd either need to hardcode the number of version numbers returned by the API, or expose a way to add to the total number of version numbers in case someone wanted to add another function, which is extra API complexity. I'm not necessarily arguing against having some kind of process by which apps could have the implementation of sin/cos/etc changed, just explaining how it gets complicated pretty quickly.
conistonwater 12y ago
There is a mathematical reason for why this is not completely horrible, even though this would catch you off guard. When x is approximately pi, the function x -> sin(x) is ill-conditioned. [0] So the relative error between the computed value ybar=fl(sin(x)) and the true value y=sin(x) is (ybar-y)/y and is quite large. However: you can find a number xbar that is very close to x, in the sense that (xbar-x)/x is very small, such that the computed "inaccurate" value ybar=fl(sin(x)) matches sin(xbar) exactly. Thus this procedure is backward-stable, but not forward-stable. For "typical" numerically-stable algorithms, this might not be such an issue. Very interesting writeup, anyway. I also found this post by James Gosling [1] describing the same issue, and the slower software workaround implemented in JVM; the post is from 2005. The paper that explain K5's implementation is "The K5 transcendental functions" by Lynch, Ahmed, Schulte, Callaway, Tisdale, but it is behind an academic paywall. [0] http://en.wikipedia.org/wiki/Condition_number http://en.wikipedia.org/wiki/Condition_number [1] https://blogs.oracle.com/jag/entry/transcendental_meditation https://blogs.oracle.com/jag/entry/transcendental_meditation
acqq 12y ago
https://software.intel.com/en-us/forums/topic/289702 https://software.intel.com/en-us/forums/topic/289702 "Shane Story (Intel) Sun, 06/27/2010 - 11:44 (...) rigorous error bounds for most of the x87 transcendental functions (f2xm1, fsincos, fyl2x...) were provided when the Pentium processor released. I recall an appendix in the orginal Pentium manual which showed a collection of ulp (unit in the last place) charts - don't recall the upper bound for all of them but it was something like 1.8 ulps. This would be for double-extended 80 bit fp inputs. Results for the fsin/fcos were for arguments < pi/4 because for the reduction for larger inputs, a 66-bit value of pi (or pi/2) was used." It's really a documentation bug. Sometimes some important details get lost after the documentation rewrites.
msie 12y ago
I'm not worthy, I'm not worthy... Edit: The intel engineer was probably smart, but then there's this guy. And then there's the commenters here. Lesson: There's always someone out there smarter than you.
- bsder 12y ago
  The lesson is actually "This was designed 35+ years ago in different constraints." Memory and ROM were small and expensive and processors were in-order and slow. Putting the transcendental in microcode made sense, back then. Now, memory is huge; processors are superscalar and speculative execution; vector units exist; and everybody ignores the built-in transcendentals, but Intel can't dump them for backwards compatibility.
  - msie 12y ago
    Yes, but from the text it seemed like fsin was broken to begin with: It is surprising that fsin is so inaccurate. I could perhaps forgive it for being inaccurate for extremely large inputs (which it is) but it is hard to forgive it for being so inaccurate on pi which is, ultimately, a very ‘normal’ input to the sin() function.
  - thaumasiotes 12y ago
    It's hard for the constraints of old to justify lying about your accuracy in the documentation you produce, though...?
Aardwolf 12y ago
> However this technique relies on printing the value of the double-precision approximation of pi to 33 digits, and up to VC++ 2013 this doesn’t work How can an IDE prevent printing a number? If all else fails, you can reinterpret_cast it to integer and handle the bits yourself?
- azernik 12y ago
  First, a nitpick: VC++ is a compiler toolchain (including libraries), not just an IDE. Also, the next sentence: So, you either need *custom printing code* or you need to wait for Dev 14 which has improved printing code.
colanderman 12y ago
Note that sin(double(pi)) should not return zero. sin(pi) is, mathematically, zero, but since float/double/long-double can’t represent pi it would actually be incorrect for them to return zero when passed an approximation to pi. Huh? sin(pi ± epsilon) better damn well return 0 if epislon is closer to 0 than the minimum representable floating-point value. Otherwise the processor is rounding incorrectly. (I realize epsilon in this case is likely to be much larger than the minimum representable float, but it's not guaranteed and thus the article's reasoning is flawed.)
- conistonwater 12y ago
  No, it is not pi that must get rounded, but sin(pi). So in sin(double(pi)), first double(pi) is computed, which is not exactly pi, then an approximation to sin(double(pi)) is computed, which is not exactly zero. It would be absolutely wrong to round a non-zero number to zero.
  - colanderman 12y ago
    Yes, we're both talking about sin(pi) (hence why I claimed "sin(pi ± epsilon) better damn well return 0"). I take issue with this assertion: It would be absolutely wrong to round a non-zero number to zero. Why? I know IEEE floats have subnormals and subnormals can be VERY small, but if the error is less than half the minimum representable float, why should it be rounded away from zero? Yes, zero is "special" because you can't divide by it, but for stability, my understanding is that you shouldn't be dividing by anything that can be so close to zero as to potentially round to it anyway.
    - conistonwater 12y ago
      sin(pi+eps) is a quantity on the order of eps. Not only is it non-zero, it is also so far from zero that it can be represented quite accurately with a floating-point number (which is a separate concern from the accuracy with which it is computed). The relative error introduced in rounding 1e-16 down to 0 is 100%, or approximately 1e16 ulp, which is why that would be unacceptable. Floating point arithmetic cares about relative errors much much more than absolute. Also, while dividing by zero is a problem, dividing by 1e-16 should not be a problem, at all. 1e-16 is a perfectly ordinary representable number.
pedrow 12y ago
Possible solution: revive the Indiana Pi Bill[0] (assuming 3.2 has an exact representation in binary?) [0]:http://en.wikipedia.org/wiki/Indiana_Pi_Bill http://en.wikipedia.org/wiki/Indiana_Pi_Bill
- leni536 12y ago
  assuming 3.2 has an exact representation in binary? It does not have an exact representation in binary.
  - thaumasiotes 12y ago
    Easy solution: amend the value to 3.25, or 3.125
leoh 12y ago
I wonder what a great standard for this is. Incidentally, do TI or HP calculators give good results?
- StefanKarpinski 12y ago
  Wolfram Alpha (and Mathematica) are generally completely accurate at the cost of being less efficient. I use those to check most things I want to verify.
- fnordfnordfnord 12y ago
  Nowhere near double precision.
preillyme 12y ago
Wow, when doing argument reduction from the 80-bit long-double format there will be almost complete cancellation? Damn.
nroets 12y ago
This is a storm in a tea cup. "The absolute error ... was fairly constant at about 4e-21, which is quite small. For many purposes this will not matter" Good programmers know that most floating point calculations introduce small errors and find ways to make their programs cope.
- lvh 12y ago
  That's pretty bad for a double that you're about to do a bunch of operations with. It's pretty easy to throw numerical precision away.
jhallenworld 12y ago
Related: here is a very good article and comments on estimating transcendental functions: http://www.coranac.com/2009/07/sines/ http://www.coranac.com/2009/07/sines/ I'm annoyed that one of the links in the comments is broken: there was a page where you entered your function and it provided a chebyshev polynomial to generate it to a particular precision. EDIT: well here is another one: http://goddard.net.nz/files/js/chebyshev/cheby.html http://goddard.net.nz/files/js/chebyshev/cheby.html
- edwintorok 12y ago
  There is some interesting info here as well on how to compute most of libm (including transcedental functions) in a correctly-rounded/precise way: http://lipforge.ens-lyon.fr/frs/download.php/153/crlibm-1.0beta3.pdf http://lipforge.ens-lyon.fr/frs/download.php/153/crlibm-1.0b...
  - jhallenworld 12y ago
    Also, the old standard reference for this was "Computer Approximations" by Hart and Cheney

424 points by ghusbands 12y ago

Of possible interest, I researched a related issue once and here's what I found: The particular approximation used by the FPU instructions fsin and fcos is valid in the domain [-π /4, π /4]. Inputs outside this domain are reduced by taking the remainder when dividing by π /4 and then using trigonometric identities to fix up values in the other three quadrants (ie sin(x) = sin(π /2 – x) for x in [π /4, π /2]. The Intel FPU is “off” because it uses an approximation for π /4 when performing the reduction. In the mid 90s, some engineers at AMD figured out how to do a more precise reduction and implemented this reduction in the K5 FPU. AMD went back to being bit-for-bit compatible with the x87 behavior in the very next generation of their CPUs, presumably because it either broke or appeared to break too many applications.

CurtHagenlocher 12y ago
I just noticed that Bruce has moved from Valve to Google... . Interesting!
sillysaurus3 12y ago
It's not possible to implement a more precise fsin() without breaking apps? One scenario I can imagine an app breaking is if someone executed fsin(x), pasted the result into their code, and then tested some other float against it using the == operator. I suppose a more common breakage scenario would be if someone did something like floor(x + fsin(x)), since that'd very slightly change where the floor "triggers" with respect to the input. But how could that cause an app to break? Every scenario I can think of comes back to someone using == to test floats, which everyone knows is a no-no for dynamic results, such as the output of a function like fsin(). I guess a programmer wouldn't really expect the result of "x + fsin(x)" to ever change for a certain value of x, so maybe they executed it with that certain value, stored the result, and did an equality comparison later, which would always work unless the implementation of fsin() changed, and the thought process went something like, "This should be reliable because realistically the implementation of fsin() won't change." Can anyone think of some obvious way that making fsin() more accurate could cause serious breakage in an app that isn't testing floats with the equality operator?
- Guvante 12y ago
  > One scenario I can imagine an app breaking is if someone executed fsin(x), pasted the result into their code, and then tested some other float against it using the == operator. While this sounds grandious, thanks to aggressive inlining this could be the case even if you wrote in your logic: sin(x) == sin(constant) Or even the following depending on your choice of delta. abs(sin(x) - sin(constant) < delta
  - sillysaurus3 12y ago
    Is it even valid for a compiler to transform sin(constant) into a hardcoded constant embedded into the executable? I guess since sin(x) is computed by an FPU in hardware, then it's inherently "dynamically linked." But in that case, why is it valid for a compiler to transform sin(constant) into a hardcoded result? Is there some standard somewhere which says that the implementation of sin/cos/etc can't ever change, and therefore it's valid for the compiler to do that? Or do compilers just make an assumption that the implementation of sin(x) won't ever change?
    - darkmighty 12y ago
      It breaks binary compatibility. The compiler substitutes some value (e.g. the broken value), so if the sin(x) evaluation changes (to the correct one), then binary execution will fail.
      - DannyBee 12y ago
        Why would the compiler substitute the broken value unless the compiler was broken? Compilers generally don't use specialized local architecture instructions to compute their results, because they can't (It would break cross-compiling or targeting different platforms in a lot of cases) In fact, it's more likely the other way around: Compiler substitutes proper value. actual calls on sin(x) on the local platform give wrong answer.
    - DannyBee 12y ago
      "but if it's dynamic, is it even valid for a compiler to transform sin(constant) into a hardcoded constant embedded into the executable?" Yes. I honestly don't remember all the rules here, but it does in fact do so, legally.

luu 12y ago

This is one of those things that's well known by people who spend a lot of time doing low-level nonsense that really ought to get wider press. You really don't want to use x86 hardware transcendental instructions. I've seen people try to measure the effectiveness of their approximation by comparing against built-in hardware, but their approximation is probably better than the hardware! If you're not writing your own assembly, you'll be ok on most modern platforms if you just link with "-lm" and use whatever math library happens to be lying around as the default, but it's still possible to get caught by this. On obscure platforms, it's basically a coin flip. I used to work for another CPU vendor, and when we implemented more accurate sin/cos functions, some benchmarks would fail us for getting the wrong result. Turns out those benchmarks hardcode a check for a couple results, and that check is based on what Intel has been doing for ages. My recollection is that we had a switch to enable the extra accuracy, but that it wasn't enabled by default because it was too much of a risk to break compatibility with Intel. If that sounds too risk averse, there's a lot of code out there that depends on your processor precisely matching Intel's behavior on undefined flags and other things that code has no business depending on. It's basically the same thing Raymond Chen is always talking about, but down one level of abstraction. I've seen what happens if you just implement something that matches the "spec" (the Intel manual) and do whatever you want for "undefined" cases. You get a chip that's basically useless because existing software will crash on it.

JoeAltmaier 12y ago
The problem is not so much "Intel stuff is old and broken" as it is "Intel wrote the standard and there's no process to move it forward"?
- sillysaurus3 12y ago
  OpenGL is a great example of how painful it is to write apps that are capable of not breaking when the standard it relies on changes. You basically end up with capability checks or version checks which cause your code to become complicated. I suppose there could've been a function that returned the implementation version of fsin/fcos/etc, which was incremented whenever Intel changed the implementation. That way benchmarks could test against that particular version. But it'd be hard to come up with a precise, consistent way to implement that. For example, do you increment the version whenever any transcendental function changes? If you do that, then you get benchmark code riddled with special cases for very specific versions, which someone has to constantly look up in documentation somewhere. You'd basically have to memorize the effect of every version increment. On the other hand, you could implement a more elaborate scheme, like having a version number for each of the transcendental functions, but then you'd either need to hardcode the number of version numbers returned by the API, or expose a way to add to the total number of version numbers in case someone wanted to add another function, which is extra API complexity. I'm not necessarily arguing against having some kind of process by which apps could have the implementation of sin/cos/etc changed, just explaining how it gets complicated pretty quickly.

conistonwater 12y ago

There is a mathematical reason for why this is not completely horrible, even though this would catch you off guard. When x is approximately pi, the function x -> sin(x) is ill-conditioned. [0] So the relative error between the computed value ybar=fl(sin(x)) and the true value y=sin(x) is (ybar-y)/y and is quite large. However: you can find a number xbar that is very close to x, in the sense that (xbar-x)/x is very small, such that the computed "inaccurate" value ybar=fl(sin(x)) matches sin(xbar) exactly. Thus this procedure is backward-stable, but not forward-stable. For "typical" numerically-stable algorithms, this might not be such an issue. Very interesting writeup, anyway. I also found this post by James Gosling [1] describing the same issue, and the slower software workaround implemented in JVM; the post is from 2005. The paper that explain K5's implementation is "The K5 transcendental functions" by Lynch, Ahmed, Schulte, Callaway, Tisdale, but it is behind an academic paywall. [0] http://en.wikipedia.org/wiki/Condition_number http://en.wikipedia.org/wiki/Condition_number [1] https://blogs.oracle.com/jag/entry/transcendental_meditation https://blogs.oracle.com/jag/entry/transcendental_meditation

acqq 12y ago

https://software.intel.com/en-us/forums/topic/289702 https://software.intel.com/en-us/forums/topic/289702 "Shane Story (Intel) Sun, 06/27/2010 - 11:44 (...) rigorous error bounds for most of the x87 transcendental functions (f2xm1, fsincos, fyl2x...) were provided when the Pentium processor released. I recall an appendix in the orginal Pentium manual which showed a collection of ulp (unit in the last place) charts - don't recall the upper bound for all of them but it was something like 1.8 ulps. This would be for double-extended 80 bit fp inputs. Results for the fsin/fcos were for arguments < pi/4 because for the reduction for larger inputs, a 66-bit value of pi (or pi/2) was used." It's really a documentation bug. Sometimes some important details get lost after the documentation rewrites.

msie 12y ago

I'm not worthy, I'm not worthy... Edit: The intel engineer was probably smart, but then there's this guy. And then there's the commenters here. Lesson: There's always someone out there smarter than you.

bsder 12y ago
The lesson is actually "This was designed 35+ years ago in different constraints." Memory and ROM were small and expensive and processors were in-order and slow. Putting the transcendental in microcode made sense, back then. Now, memory is huge; processors are superscalar and speculative execution; vector units exist; and everybody ignores the built-in transcendentals, but Intel can't dump them for backwards compatibility.
- msie 12y ago
  Yes, but from the text it seemed like fsin was broken to begin with: It is surprising that fsin is so inaccurate. I could perhaps forgive it for being inaccurate for extremely large inputs (which it is) but it is hard to forgive it for being so inaccurate on pi which is, ultimately, a very ‘normal’ input to the sin() function.
- thaumasiotes 12y ago
  It's hard for the constraints of old to justify lying about your accuracy in the documentation you produce, though...?

Aardwolf 12y ago

> However this technique relies on printing the value of the double-precision approximation of pi to 33 digits, and up to VC++ 2013 this doesn’t work How can an IDE prevent printing a number? If all else fails, you can reinterpret_cast it to integer and handle the bits yourself?

azernik 12y ago
First, a nitpick: VC++ is a compiler toolchain (including libraries), not just an IDE. Also, the next sentence: So, you either need *custom printing code* or you need to wait for Dev 14 which has improved printing code.

colanderman 12y ago

Note that sin(double(pi)) should not return zero. sin(pi) is, mathematically, zero, but since float/double/long-double can’t represent pi it would actually be incorrect for them to return zero when passed an approximation to pi. Huh? sin(pi ± epsilon) better damn well return 0 if epislon is closer to 0 than the minimum representable floating-point value. Otherwise the processor is rounding incorrectly. (I realize epsilon in this case is likely to be much larger than the minimum representable float, but it's not guaranteed and thus the article's reasoning is flawed.)

conistonwater 12y ago
No, it is not pi that must get rounded, but sin(pi). So in sin(double(pi)), first double(pi) is computed, which is not exactly pi, then an approximation to sin(double(pi)) is computed, which is not exactly zero. It would be absolutely wrong to round a non-zero number to zero.
- colanderman 12y ago
  Yes, we're both talking about sin(pi) (hence why I claimed "sin(pi ± epsilon) better damn well return 0"). I take issue with this assertion: It would be absolutely wrong to round a non-zero number to zero. Why? I know IEEE floats have subnormals and subnormals can be VERY small, but if the error is less than half the minimum representable float, why should it be rounded away from zero? Yes, zero is "special" because you can't divide by it, but for stability, my understanding is that you shouldn't be dividing by anything that can be so close to zero as to potentially round to it anyway.
  - conistonwater 12y ago
    sin(pi+eps) is a quantity on the order of eps. Not only is it non-zero, it is also so far from zero that it can be represented quite accurately with a floating-point number (which is a separate concern from the accuracy with which it is computed). The relative error introduced in rounding 1e-16 down to 0 is 100%, or approximately 1e16 ulp, which is why that would be unacceptable. Floating point arithmetic cares about relative errors much much more than absolute. Also, while dividing by zero is a problem, dividing by 1e-16 should not be a problem, at all. 1e-16 is a perfectly ordinary representable number.

pedrow 12y ago

Possible solution: revive the Indiana Pi Bill[0] (assuming 3.2 has an exact representation in binary?) [0]:http://en.wikipedia.org/wiki/Indiana_Pi_Bill http://en.wikipedia.org/wiki/Indiana_Pi_Bill

leni536 12y ago
assuming 3.2 has an exact representation in binary? It does not have an exact representation in binary.
- thaumasiotes 12y ago
  Easy solution: amend the value to 3.25, or 3.125

leoh 12y ago

I wonder what a great standard for this is. Incidentally, do TI or HP calculators give good results?

StefanKarpinski 12y ago
Wolfram Alpha (and Mathematica) are generally completely accurate at the cost of being less efficient. I use those to check most things I want to verify.
fnordfnordfnord 12y ago
Nowhere near double precision.

preillyme 12y ago

Wow, when doing argument reduction from the 80-bit long-double format there will be almost complete cancellation? Damn.

nroets 12y ago

This is a storm in a tea cup. "The absolute error ... was fairly constant at about 4e-21, which is quite small. For many purposes this will not matter" Good programmers know that most floating point calculations introduce small errors and find ways to make their programs cope.

lvh 12y ago
That's pretty bad for a double that you're about to do a bunch of operations with. It's pretty easy to throw numerical precision away.

jhallenworld 12y ago

Related: here is a very good article and comments on estimating transcendental functions: http://www.coranac.com/2009/07/sines/ http://www.coranac.com/2009/07/sines/ I'm annoyed that one of the links in the comments is broken: there was a page where you entered your function and it provided a chebyshev polynomial to generate it to a particular precision. EDIT: well here is another one: http://goddard.net.nz/files/js/chebyshev/cheby.html http://goddard.net.nz/files/js/chebyshev/cheby.html

edwintorok 12y ago
There is some interesting info here as well on how to compute most of libm (including transcedental functions) in a correctly-rounded/precise way: http://lipforge.ens-lyon.fr/frs/download.php/153/crlibm-1.0beta3.pdf http://lipforge.ens-lyon.fr/frs/download.php/153/crlibm-1.0b...
- jhallenworld 12y ago
  Also, the old standard reference for this was "Computer Approximations" by Hart and Cheney