Unexpected loss of precision when dividing doubles

Question

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

Answer 1

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

Answer 2

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.

Answer 3

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?

Answer 4

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

Answer 5

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.

Answer 6

yes this is what I did. I listed the parameters with precision up to 15 digits and so my former parameters where just some approximations.

well..if my parameters are chopped than I am guessing my computations will be approximate for sure. thank you, madalina

Answer 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;
}

Answer 8

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.

Answer 9

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].

Answer 10

Okay so this is an old question, but I have the exact same scenario: finding the slope in regions which have arbitrarily close points.

What brought to stackoverflow was the idea of finding a method to rejig the equation. I was hoping to verify that:

float get_slope(float dXa, float dXb, float dYa, float dYb) {
    return (dXa - dXb)/(dYa - dYb);
}

could be improved with:

float get_slope(float dXa, float dXb, float dYa, float dYb) {
    return  dXa/(dYa - dYb) - dXb(dYa - dYb);
}

My questions then are:

1. Is there a better way to rework the equation?
2. Should there be multiple equations (paths) in the get_slope() function?

Answer 11

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.