Vector addition in C++ is not coming out exactly right

Question

I am trying to add one vector to another vector in my C++ program and addition is coming out inexact.

This is the vector class (using cmath library):

class Vec{
    float dir, mag;
    public:
        Vec(float dir, float mag){
            this->dir = dir;
            this->mag = mag;
        }
        float getX(){
            return cos(dir)*mag;
        }
        Vec operator + (Vec v2){
            float triangleX = cos(dir)*mag+cos(v2.dir)*v2.mag;
            float triangleY = sin(dir)*mag+sin(v2.dir)*v2.mag;
            return Vec(atan2(triangleY, triangleX), sqrt(pow(triangleX,2)+pow(triangleY,2)));
    }
};

And this is the main function:

int main(){
    Vec v1(0, 3); // 0º
    Vec v2(3.14159265/2, 3); // 90º
    Vec v3(3.14159265, 3); // 180º
    std::cout.precision(15);
    std::cout<<"v1: "<<v1.getX()<<std::endl;
    std::cout<<"v1+v2: "<<(v1+v2).getX()<<std::endl;
    std::cout<<"v1+v3: "<<(v1+v3).getX()<<std::endl;
    return 0;
}

And this is the output:

v1: 3
v1+v2: 2.99999976158142
v1+v3: 1.98007097372034e-014

As you can see, the first output v1 is fine.

The second output is the addition of 0 degrees and 90 degrees (an angle that was not supposed to affect the x component), his x component is close to 3, but not exactly 3.

The third output is the addition of 2 opposite vectors with the same magnitude, it was supposed to be 0, but it is not what shows here.

What is causing these strange additions and how do I make them exact?

using double is recommended, not only because it's more accurate but most modern architectures are optimized for double, not float — phuclv, May 04 '14 at 13:14
@LưuVĩnhPhúc Do you know of a platform where using `double` is faster than `float`? Afaik all platforms support both, most are still twice as fast with `float`. And the memory footprint of a `float` is also only half that of `double` so that memory bound computations will perform faster with `float` on all platforms. Still, I'm entirely with you that `double` should be used as a default. — cmaster - reinstate monica, May 04 '14 at 13:26
@cmaster I heard many people say that they see double is faster than float. Actually I haven't done a benchmark before, but most answers on this sites state that doubles on modern platforms are at least just as fast as floats. http://stackoverflow.com/questions/1074474/should-i-use-double-or-float http://stackoverflow.com/questions/4584637/double-or-float-which-is-faster — phuclv, May 07 '14 at 04:15
@LưuVĩnhPhúc Reading the answers on the link, I can't find anybody claiming that doubles **are** faster than float. Rightly so, because that would be wrong. On all platforms I know, the opposite is true. But I don't know all platforms, so I can't rule out the possibility that there is this wild CPU that crunches doubles faster than floats. So, what is probably true is that doubles **might** be faster on some (esoteric) platforms. But that was precisely my question: Do you **know** of a specific platform where doubles are faster? — cmaster - reinstate monica, May 08 '14 at 05:15

score 4 · Answer 1 · edited May 23 '17 at 12:33

4

v1 + v3 is almost 0.0, the code is working correctly.

Why is it not exactly 0.0?

because some numbers can not be represented exactly as doubles. see: C/C++: 1.00000 <= 1.0f = False for some explanations.

Also: Pi is an irrational number, and can not be represented exactly in any base that is a natural number.

Thus any calculation involving Pi is never completely precise.

And: sin, cos, sqrt are usually all implemented as algorithms that don't return completely exact results, for example as approximative numeric algorithms.

edited May 23 '17 at 12:33

Community

1
1

answered May 04 '14 at 13:12

alain

11,939
2
31
51

I've intended to comment about that just right before you edit the post. Also, the use of literals like that instead of PI constant is not recommended – phuclv May 04 '14 at 13:17
Correct me if I'm wrong, but I thought that `sin` & co. are only approximative when you turn on the fast-math compiler flag? – cmaster - reinstate monica May 04 '14 at 13:28
@cmaster I don't know in which libraries/with which flags they are approximative, but one of the faster implementations is a table with linear interpolation. (That's actually not approximative, but not completely precise too) – alain May 04 '14 at 13:36
@cmaster I edited my answer to be a bit more accurate. Thanks for the input. – alain May 04 '14 at 13:45
@alain: To be pedantic, one could exactly represent pi as 1 base pi. – Edward May 04 '14 at 13:53
@Edward I think you can never be too pedantic with mathematics ;-) I'll change that, thanks. – alain May 04 '14 at 13:58

score 1 · Accepted Answer · edited May 23 '17 at 12:00

The basic problem you're having is one of limited precision of float and the value of pi that you're using. Moving to a double will help, and since you're already including <cmath> you should use the value there which is M_PI and is accurate to at least the precision of double. With that said, you're still not going to get exact answers with this approach. (The answer from alain does a good job of explaining why.)

There are some improvements that could be made. One is a neat trick using a C++11 feature which is "user-defined string literals." If you add this definition to your code:

constexpr long double operator"" _deg(long double deg) {
    return deg*M_PI/180;
}

You can now append _deg to any long double literal and it will automatically be converted to radians at compile time. Your main function would then look like this:

int main(){
    Vec v1(0.0_deg, 3); 
    Vec v2(90.0_deg, 3); 
    Vec v3(180.0_deg, 3); 
    // ...
}

The next thing you could do would be to store the x and y coordinates and only do the trigonometric manipulations when needed. That version of Vec might look like this:

class Vec{
    double x,y;
    public:
        Vec(double dir, double mag, bool cartesian=false) : x(dir), y(mag) {
            if (!cartesian) {
                x = mag*cos(dir);
                y = mag*sin(dir);
            }
        }
        double getX() const {
            return x;
        }
        Vec operator + (const Vec &v2){
            return Vec(x+v2.x, y+v2.y, true);
        }
}

Note that I've created a bool value for the constructor which tells whether the input is to be a magnitude and direction or an x and y value. Also note that the getX() is declared const because it doesn't alter the Vec and that the argument to operator+ is also a const reference for the same reason. When I make those changes on my machine (a 64-bit machine), I get the following output:

v1: 3
v1+v2: 3
v1+v3: 0

Vector addition in C++ is not coming out exactly right

2 Answers2