[Gllug] Basic numerical precision question

James Courtier-Dutton james.dutton at gmail.com
Sun Sep 26 17:13:47 UTC 2010


What your program is outputting is correct.
See the wiki page for an explanation.
http://en.wikipedia.org/wiki/Double_precision_floating-point_format

Summary: A Double gives you about 15 digits (base 10) to play with (53
bits in base 2), irrespective of where the decimal point is.

Also, if you change setprecision(40) you will see what is happening.
Each time you do an addition or a multiplication the amount of error
is increasing.
You only appear to need 10 digits, so the best way to reduce the error
values from accumulating is to round the values off to 10 digital
after each maths operation.
This will filter out the error values before they accumulate to cause problems.
In your particular calculation you can do a simple round after the
last calculation to filter out the 14th digit errors.

A very important consideration with floating point value calculations
is that comparing float values is very difficult, and particularly
nonsensical if you do if (a == b) because the answer will almost
always be false due to the introduced error values.
What "if (a == b)" really calculates is "if (a == (b + error))" where
"error" is some unknown, but likely to be small positive or negative
value.
So, the best you can do is do tests to see "if a is close enough to b".

Kind Regards

James



On 26 September 2010 00:18, Sanatan Rai <sanatan at gmail.com> 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
>
>
> --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
>
-- 
Gllug mailing list  -  Gllug at gllug.org.uk
http://lists.gllug.org.uk/mailman/listinfo/gllug




More information about the GLLUG mailing list