src/tests/wrongr.cpp

   1 #include "main.h"
   2 #include <cmath>
   3
   4 using namespace IPDF;
   5 /**
   6  * Example from Kahan's "Why is Floating-Point Computation so Hard to Debug when it Goes Wrong?"
   7  * Published on Kahan's website, March 2007  http://www.cs.berkeley.edu/~wkahan/WrongR.pdf
   8  * In this example we compute h(x) = |x| ... except it is not!
   9  * Note that the errors are NOT due to catastrophic cancellation (no subtraction) or division (no division) or accumulated rounding (only 128 operations)
  10  * http://www.cs.berkeley.edu/~wkahan/WrongR.pdf
  11  */
  12
  13 Real h(Real x)
  14 {
  15         Real y = abs(x);
  16         for (int k = 0; k < 128; ++k)
  17         {
  18                 y = sqrt(y);
  19         }
  20         Real w = y;
  21         for (int k = 0; k < 128; ++k)
  22         {
  23                 w = w*w;
  24         }
  25         return w;
  26 }
  27
  28 int main(int argc, char ** argv)
  29 {
  30         for (Real x = -2; x < 2; x+=1e-3)
  31         {
  32                 printf("%lf\t%lf\t%lf\n", x, Float(h(x)), abs(Float(x)));
  33         }
  34         return 0;
  35 }