[Gllug] Basic numerical precision question

[Sanatan Rai]

Hi,
    I am doing some numerical computations at home and see numerical
precision issues where there shouldn't be any.

   Briefly, I am subtracting two doubles that are numerically the
same and taking the square root. The difference is negative of
with magnitude 1e-14, which is greater than numerical_limits<double>::epsilon().
This is very strange, as I'd expect the difference to be less than
1e-31 or something.

  Or am I missing something very basic about numerical precision here?


  I am using g++ 4.3.2 on Debian Lenny. Sample code is appended here.

// BEGIN C++ Code
#include <iostream>
#include <iomanip>
#include <cmath>
#include <limits>
#include <vector>
using namespace std;

int main()
{
  double val[] = {29.1, 29.1, 29.1};
  double sum(0.0), sum2(0.0);
  for (int i = 0; i < 3; i++)
  {
    sum += val[i];
    sum2 += val[i]*val[i];
  }
  double avg = sum / 3.0;
  cout << setprecision(10)
       << "sum2 = " << sum2 << endl
       << "avgSum2 = " << sum2/3.0 << endl
       << "avg*avg = " << avg * avg << endl
       << "avgSum2 - avg*avg = " << (sum2/3.0 - avg * avg) << endl;
  return 0;
}
// END Code

Complie with g++ <filename>

My output:

sanat at fractal:~/src/test$ uname -as
Linux fractal 2.6.32-bpo.3-686 #1 SMP Wed Mar 17 14:31:18 UTC 2010
i686 GNU/Linux
sanat at fractal:~/src/test$ g++ -o prectest prectest.cpp
sanat at fractal:~/src/test$ ./prectest
sum2 = 2540.43
avgSum2 = 846.81
avg*avg = 846.81
avgSum2 - avg*avg = -1.924571613e-13 # This seems too big!!


On my system double is 8 bytes and numerical_limits<double>::epsilon()
yields 2.22045e-16. So am mystified as to why the difference
above is as large as -2e-13


--Sanatan
-- 
Sanatan Rai
3, Admirals Court,
30, Horselydown Lane,
London, SE1 2LJ.
+44-20-7403-2479.
-- 
Gllug mailing list  -  Gllug at gllug.org.uk
http://lists.gllug.org.uk/mailman/listinfo/gllug