From b95fd30aafb3f47d2bfc8aa8f75600db4fdb995c Mon Sep 17 00:00:00 2001 From: Sam Moore Date: Fri, 2 May 2014 00:21:11 +0800 Subject: [PATCH] Test from one of Kahan's articles "Why is Floating-Point Computation so Hard to Debug when it Goes Wrong?" This is a MUCH better example than massive rounding errors as in calculatepi.test --- src/tests/wrongr.cpp | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 src/tests/wrongr.cpp diff --git a/src/tests/wrongr.cpp b/src/tests/wrongr.cpp new file mode 100644 index 0000000..a7f5ec0 --- /dev/null +++ b/src/tests/wrongr.cpp @@ -0,0 +1,35 @@ +#include "main.h" +#include + +using namespace IPDF; +/** + * Example from Kahan's "Why is Floating-Point Computation so Hard to Debug when it Goes Wrong?" + * Published on Kahan's website, March 2007 http://www.cs.berkeley.edu/~wkahan/WrongR.pdf + * In this example we compute h(x) = |x| ... except it is not! + * Note that the errors are NOT due to catastrophic cancellation (no subtraction) or division (no division) or accumulated rounding (only 128 operations) + * http://www.cs.berkeley.edu/~wkahan/WrongR.pdf + */ + +Real h(Real x) +{ + Real y = abs(x); + for (int k = 0; k < 128; ++k) + { + y = sqrt(y); + } + Real w = y; + for (int k = 0; k < 128; ++k) + { + w = w*w; + } + return w; +} + +int main(int argc, char ** argv) +{ + for (Real x = -2; x < 2; x+=1e-3) + { + printf("%lf\t%lf\t%lf\n", x, Float(h(x)), abs(Float(x))); + } + return 0; +} -- 2.20.1