Search code examples
c++mathboostquaternions

Wrong multiplication of quaternions


I made my quaternion realization and to test it I used boost. I test it in cycle and I always have errors in multiplication and division. Then I used Wolfram to find out who was wrong and I find out that boost was wrong.

I am not sure that the boost library was Wong, maybe some one can find out what is going on?

Console output example:

0
Operation : [ * ]
First vector =  [ -41.4168   -92.2373     -33.0126   -42.2364 ]
Second vector = [ -67.8087   -60.3523     58.8705   36.3265 ]
My multiplication =    [ 719.45   10041.3   5700.03   -6062.97 ]
Boost multiplication = [ 719.45   10041.3   10041.3   5700.03 ]

main.cpp

#include "quaternion.h"
#include <boost/math/quaternion.hpp>
#include <random>
#include <time.h>
#include <QDebug>

using boost::math::quaternion;

double fRand(double fMin, double fMax)
{
    double f = (double)rand() / RAND_MAX;
    return fMin + f * (fMax - fMin);
}

bool isEqual(double a, double b){
    return std::abs(a-b) < std::numeric_limits<double>::epsilon();
}

bool isEqual(Quaternion& quat1, quaternion<double> quat2){
    if (!isEqual(quat2.R_component_1(), quat1.real)) return false;
    if (!isEqual(quat2.R_component_2(), quat1.imagine.data.at(0))) return false;
    if (!isEqual(quat2.R_component_3(), quat1.imagine.data.at(1))) return false;
    if (!isEqual(quat2.R_component_4(), quat1.imagine.data.at(2))) return false;
    return true;
}

int main () {

    std::srand(std::time(nullptr));
    Quaternion compR;
    quaternion<double> boost;

    double a1, a2, a3, a4,
            b1, b2, b3, b4;
try {
    for (unsigned int i=0; i<10000000; ++i){
        a1 = fRand(-100, 100);
        a2 = fRand(-100, 100);
        a3 = fRand(-100, 100);
        a4 = fRand(-100, 100);

        b1 = fRand(-100, 100);
        b2 = fRand(-100, 100);
        b3 = fRand(-100, 100);
        b4 = fRand(-100, 100);

        Quaternion comp1{a1,a2,a3,a4};
        Quaternion comp2{b1,b2,b3,b4};
        quaternion<double> a(a1,a2,a3,a4);
        quaternion<double> b(b1,b2,b3,b4);


        //if (i%50000==0)
            qDebug() << i;

        compR = comp1+comp2;
        boost = a+b;
        if (!isEqual(compR, boost))
            throw std::runtime_error("+");

        compR = comp1-comp2;
        boost = a-b;
        if (!isEqual(compR, boost))
            throw std::runtime_error("-");

        compR = comp1*comp2;
        boost = a*b;
        if (!isEqual(compR, boost))
            throw std::runtime_error("*");

        compR = comp1/comp2;
        boost = a/b;
        if (!isEqual(compR, boost))
            throw std::runtime_error("/");

    }
    } catch (const std::runtime_error& error) {
        qDebug() << "Operation : [" << error.what() << "]";
        qDebug() << " First vector =  [" << a1 << " " << a2 << " " << " " << a3 << " " << a4 << "]\n" <<
                    "Second vector = [" << b1 << " " << b2 << " " << " " << b3 << " " << b4 << "]";
        qDebug() << " My multiplication =    [" << compR.real << " " << compR.imagine.data.at(0) << " " << compR.imagine.data.at(1)<< " " << compR.imagine.data.at(2) << "]\n" <<
                    "Boost multiplication = [" << boost.R_component_1() << " " << boost.R_component_2() << " " << boost.R_component_2() << " " << boost.R_component_3() << "]";
    }

    return 0;
}

Full project:https://github.com/Yegres5/quaternion


Solution

  • You have two issues:

    • you have a bug printing the boost quaternion (you print the index 2 twice)
    • use a larger epsilon. Using epsilon as it is can only be suitable, if your numbers around 1. Your result is around 10,000, so you need to use at least 10,000x larger epsilon. At this still may not be large enough, as quaternion multiplication uses several operations, each one can add some additional error. So, to be safe, multiply epsilon further by ~10 or so.