[Gllug] Basic numerical precision question

['lesleyb']

On Sun, Sep 26, 2010 at 12:18:10AM +0100, Sanatan Rai wrote:
> 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


On the basis that machine epsilon is the smallest positive number 
distinguishable from 0, and that the absolute value of your result 
is larger than epsilon, your problem is where exactly?

But I am still suspicious - 
 (a + b) * (a + b) = a^2 + b^2 +2ab > a^2+b^2 when a,b > 0.

sum2 is equivalent to a^2 + b^2 
sum  is equivalent to (a+b)*(a+b)

Thus sum > sum2 always.

(You should disregard the sign of the mantissa in your result.)

Regards

L.

-- 
Gllug mailing list  -  Gllug at gllug.org.uk
http://lists.gllug.org.uk/mailman/listinfo/gllug