How to deal with large sums of floating points (C++)

Question

I'm trying to understand why computing determinant of a 3x3 matrix in this way gives broken results:

#include <iostream>
#include <array>
#include <vector>

int main(int argc, char const *argv[])
{
    const int dim = 3;
    std::vector<std::vector<float>> data;
    data.resize(dim, std::vector<float>(dim, 0));
    std::array<float, dim*dim> myArray = {
        -1000.0, -1000.0, -1000.0, 
        -1000.0, -1000.0, -1000.0, 
        -1000.0, -1000.0, -999.0
    };

    for(int i = 0; i < dim; i++){
        for(int j = 0; j < dim; j++){
            data[i][j] = myArray[dim*i + j];
        }
    }
    float det = 
    data[0][0] * data[1][1] * data[2][2] + 
    data[0][1] * data[1][2] * data[2][0] + 
    data[0][2] * data[2][1] * data[1][0] - 
    data[2][0] * data[1][1] * data[0][2] - 
    data[0][0] * data[1][2] * data[2][1] - 
    data[1][0] * data[0][1] * data[2][2];

    float anotherDet = 0;
    for(int i = 0; i < 2; i++){
        anotherDet = anotherDet + 1000 * 1000 * 1000;
    }
    for(int i = 0; i < 2; i++){
        anotherDet = anotherDet - 1000 * 1000 * 1000;
    }
    anotherDet = anotherDet - 1000 * 1000 * 999;
    anotherDet = anotherDet + 1000 * 1000 * 999;

    std::cout << det << " " << anotherDet << std::endl;
    return 0;
}

The output of this code is -64 0, whereas the real determinant of such matrix should be zero. My suspicion is that the error is related to the floating points precision, but I still don't get the logic for why such discrepancies happen.

I tried debugging, and indeed I found that there is a rounding error that gets carried on. -2999000060 is the value of det in the middle of that large sum.

But if I try the following:

float det2 = 
    1000 * 1000 * 1000 + 
    1000 * 1000 * 1000 + 
    1000 * 1000 * 999 - 
    1000 * 1000 * 1000 - 
    1000 * 1000 * 1000 - 
    1000 * 1000 * 999;

It gives the correct value of 0, but only because its an arithmetic operation, and thus there are not rounding errors.

Edit: I understand now how operations with floats that big could carry some rounding error.

And when you used your debugger to run your program, what did you see? This is precisely what a debugger is for. If you don't know how to use a debugger this is a good opportunity to learn how to use it to run your program one line at a time, monitor all variables and their values as they change, and analyse your program's logical execution flow. Knowing how to use a debugger is a required skill for every C++ developer, no exceptions. With your debugger's help you should be able to quickly find all bugs in this and all future programs you write, without having to ask anyone for help. — Sam Varshavchik, Jun 18 '20 at 16:50
Does this answer your question? [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) — EOF, Jun 18 '20 at 16:56
@SamVarshavchik I tried using a debugger, but because ```float det``` is a one-liner, I can't really see whats going on in that sum. — hao123, Jun 18 '20 at 17:00
@EOF If I try to explicitly write the determinant expression with numeric values, the result is 0 instead of -64. I still don't get why — hao123, Jun 18 '20 at 17:03
Then split that big sum in multiple individual steps, using some temporaries; and check the value at every step. — peppe, Jun 18 '20 at 17:04
@hao123 I have no idea what you mean. Maybe you should add your observations to the question, including a [mre]? — EOF, Jun 18 '20 at 17:08
Unrelated: Here's [a way to define a matrix](https://stackoverflow.com/a/2076668/4581301) that is significantly faster than the `vector` of `vector`s approach for small matrices. You'll recognize the trick it uses because it's exactly how `myArray` works. — user4581301, Jun 18 '20 at 17:12
Right! Understood now! Because it is a sum of integers, instead of floats — hao123, Jun 18 '20 at 17:21
@hao123 Your example of `float det2 = 1000 * 1000 ...` is evaluated **as an integer type**. You won't see rounding errors there. To see this, change the first literal *to `float`*: `float det2 = 1000.f * 1000 ...` — EOF, Jun 18 '20 at 17:22
What is your question here? What do you not understand about floating point error? — eesiraed, Jun 18 '20 at 17:39

How to deal with large sums of floating points (C++)

0 Answers0