Unexpected loss of precision when dividing doubles

问题:

I have a function getSlope which takes as parameters 4 doubles and returns another double calculated using this given parameters in the following way:

double QSweep::getSlope(double a, double b, double c, double d){
double slope;
slope=(d-b)/(c-a);
return slope;
}

The problem is that when calling this function with arguments for example:

getSlope(2.71156, -1.64161, 2.70413, -1.72219);

the returned result is:

10.8557

and this is not a good result for my computations. I have calculated the slope using Mathematica and the result for the slope for the same parameters is:

10.8452

or with more digits for precision:

10.845222072678331.

The result returned by my program is not good in my further computations. Moreover, I do not understant how does the program returns 10.8557 starting from 10.845222072678331 (supposing that this is the approximate result for the division)? How can I get the good result for my division?

thank you in advance, madalina

---

I print the result using the command line:

std::cout<<slope<<endl;

It may be that my parameters are maybe not good, as I read them from another program (which computes a graph; after I read this parameters fromt his graph I have just displayed them to see their value but maybe the displayed vectors have not the same internal precision for the calculated value..I do not know it is really strange. Some numerical errors appears..)

When the graph from which I am reading my parameters is computed, some numerical libraries written in C++ (with templates) are used. No OpenGL is used for this computation.

thank you, madalina

回答1:

I've tried with float instead of double and I get 10.845110 as a result. It still looks better than madalina result.

EDIT:

I think I know why you get this results. If you get a, b, c and d parameters from somewhere else and you print it, it gives you rounded values. Then if you put it to Mathemtacia (or calc ;) ) it will give you different result.

I tried changing a little bit one of your parameters. When I did:

double c = 2.7041304;

I get 10.845806. I only add 0.0000004 to c! So I think your "errors" aren't errors. Print a, b, c and d with better precision and then put them to Mathematica.

回答2:

The following code:

#include <iostream>
using namespace std;

double getSlope(double a, double b, double c, double d){
    double slope;
    slope=(d-b)/(c-a);
    return slope;
}

int main( ) {
    double s = getSlope(2.71156, -1.64161, 2.70413, -1.72219);
    cout << s << endl;
}

gives a result of 10.8452 with g++. How are you printing out the result in your code?

回答3:

Could it be that you use DirectX or OpenGL in your project? If so they can turn off double precision and you will get strange results.

You can check your precision settings with

std::sqrt(x) * std::sqrt(x)

The result has to be pretty close to x. I met this problem long time ago and spend a month checking all the formulas. But then I've found

D3DCREATE_FPU_PRESERVE

回答4:

The problem here is that (c-a) is small, so the rounding errors inherent in floating point operations is magnified in this example. A general solution is to rework your equation so that you're not dividing by a small number, I'm not sure how you would do it here though.

EDIT:

Neil is right in his comment to this question, I computed the answer in VB using Doubles and got the same answer as mathematica.

回答5:

The results you are getting are consistent with 32bit arithmetic. Without knowing more about your environment, it's not possible to advise what to do.

Assuming the code shown is what's running, ie you're not converting anything to strings or floats, then there isn't a fix within C++. It's outside of the code you've shown, and depends on the environment.

As Patrick McDonald and Treb brought both up the accuracy of your inputs and the error on a-c, I thought I'd take a look at that. One technique to look at rounding errors is interval arithmetic, which makes the upper and lower bounds which value represents explicit (they are implicit in floating point numbers, and are fixed to the precision of the representation). By treating each value as an upper and lower bound, and by extending the bounds by the error in the representation ( approx x * 2 ^ -53 for a double value x ), you get a result which gives the lower and upper bounds on the accuracy of a value, taking into account worst case precision errors.

For example, if you have a value in the range [1.0, 2.0] and subtract from it a value in the range [0.0, 1.0], then the result must lie in the range [below(0.0),above(2.0)] as the minimum result is 1.0-1.0 and the maximum is 2.0-0.0. below and above are equivalent to floor and ceiling, but for the next representable value rather than for integers.

Using intervals which represent worst-case double rounding:

getSlope(
 a = [2.7115599999999995262:2.7115600000000004144], 
 b = [-1.6416099999999997916:-1.6416100000000002357], 
 c = [2.7041299999999997006:2.7041300000000005888], 
 d = [-1.7221899999999998876:-1.7221900000000003317])
(d-b) = [-0.080580000000000526206:-0.080579999999999665783]
(c-a) = [-0.0074300000000007129439:-0.0074299999999989383218]

to double precision [10.845222072677243474:10.845222072679954195]

So although c-a is small compared to c or a, it is still large compared to double rounding, so if you were using the worst imaginable double precision rounding, then you could trust that value's to be precise to 12 figures - 10.8452220727. You've lost a few figures off double precision, but you're still working to more than your input's significance.

But if the inputs were only accurate to the number significant figures, then rather than being the double value 2.71156 +/- eps, then the input range would be [2.711555,2.711565], so you get the result:

getSlope(
 a = [2.711555:2.711565], 
 b = [-1.641615:-1.641605], 
 c = [2.704125:2.704135], 
 d = [-1.722195:-1.722185])
(d-b) = [-0.08059:-0.08057]
(c-a) = [-0.00744:-0.00742]

to specified accuracy [10.82930108:10.86118598]

which is a much wider range.

But you would have to go out of your way to track the accuracy in the calculations, and the rounding errors inherent in floating point are not significant in this example - it's precise to 12 figures with the worst case double precision rounding.

On the other hand, if your inputs are only known to 6 figures, it doesn't actually matter whether you get 10.8557 or 10.8452. Both are within [10.82930108:10.86118598].

回答6:

Better Print out the arguments, too. When you are, as I guess, transferring parameters in decimal notation, you will lose precision for each and every one of them. The problem being that 1/5 is an infinite series in binary, so e.g. 0.2 becomes .001001001.... Also, decimals are chopped when converting an binary float to a textual representation in decimal.

Next to that, sometimes the compiler chooses speed over precision. This should be a documented compiler switch.

回答7:

Patrick seems to be right about (c-a) being the main cause:

d-b = -1,72219 - (-1,64161) = -0,08058

c-a = 2,70413 - 2,71156 = -0,00743

S = (d-b)/(c-a)= -0,08058 / -0,00743 = 10,845222

You start out with six digits precision, through the subtraction you get a reduction to 3 and four digits. My best guess is that you loose additonal precision because the number -0,00743 can not be represented exaclty in a double. Try using intermediate variables with a bigger precision, like this:

double QSweep::getSlope(double a, double b, double c, double d)
{
    double slope;
    long double temp1, temp2;

    temp1 = (d-b);
    temp2 = (c-a);
    slope = temp1/temp2;

    return slope;
}

回答8:

While the academic discussion going on is great for learning about the limitations of programming languages, you may find the simplest solution to the problem is an data structure for arbitrary precision arithmetic.

This will have some overhead, but you should be able to find something with fairly guaranteeable accuracy.