Hi everyone,

Below is a simple code snippet that I isolated from a buggy project. It computes a division of a float zero by a double (the weird value of the double is taken from a case where the bug shows up; the value is around 2.86346e-45). The division is performed both with primitive data type float (variable "a"), and with an Eigen 1x1 mtrix templated by type float (variable "c"). I compiled it with GCC 4.2.1 on a Mac, and its results are a = 0 (correct) but c = nan (wrong). Can you explain it? Maybe it has something to do with mixing float and double in a division and somehow Eigen isn't handling the type promotion correctly?

Thank you very much!

Code: Select all

#include <iostream>
#include <Eigen/Core>

int main()
{
    unsigned char buf[8] = { 0x4d, 0x27, 0xd6, 0x46, 0xf4, 0x58, 0xb0, 0x36 };  // reverse the order if your machine is big endian
    double b;
    b = *((double *)buf);

    float a;
    a = 0;
    a /= b;

    Eigen::Matrix<float, 1, 1> c;
    c[0] = 0;
    c /= b;

    std::cout << "b = " << b << " a = " << a << " c = " << c << std::endl;

    return 0;
}


Here is a simplified version of the same problem, without any doubles or type conversion. The error only appears if you compile without optimization:

Code: Select all

#include <iostream>
#include <Eigen/Core>

int main()
{
    float a = 0, b = 2.8e-45;
    Eigen::Matrix<float, 1, 1> c;
    c[0] = 0;

    a /= b;
    c /= b;

    std::cout << "b = " << b << " a = " << a << " c = " << c << std::endl;
    std::cout << "1/b = " << 1/b << std::endl;

    return 0;
}


I think that the issue here is that Eigen evaluates "c /= b" as "c *= 1/b". Normally, these are the same. But in this case, b is one of the smallest numbers that fits in a float (it is twice the smallest subnormal positive number), and 1/b overflows, as the last statement shows. So Eigen computes 0 times infinity, which is NaN.

For Eigen developers: This is in SelfCwiseBinaryOp.h:188. Why do we do this? Performance? It's been there for a long time; I traced it back to revision 87a668aa3764 in November 2009 but it may well have been introduced before that.

The issue disappears if we change it to

Code: Select all

template<typename Derived>
inline Derived& DenseBase<Derived>::operator/=(const Scalar& other)
{
  typedef typename Derived::PlainObject PlainObject;
  SelfCwiseBinaryOp<internal::scalar_quotient_op<Scalar>, Derived, typename PlainObject::ConstantReturnType> tmp(derived());
  tmp = PlainObject::Constant(rows(),cols(),other);
  return derived();
}


Yes, it's been like this since the beginning. There is the same optimization in Functors.h:493. I guess this is an optimization that we could leave to the user: mat * (1/a)

Thank you very much for the insight.

jitseniesen wrote:Here is a simplified version of the same problem, without any doubles or type conversion. The error only appears if you compile without optimization:

Code: Select all

#include <iostream>
#include <Eigen/Core>

int main()
{
    float a = 0, b = 2.8e-45;
    Eigen::Matrix<float, 1, 1> c;
    c[0] = 0;

    a /= b;
    c /= b;

    std::cout << "b = " << b << " a = " << a << " c = " << c << std::endl;
    std::cout << "1/b = " << 1/b << std::endl;

    return 0;
}


I think that the issue here is that Eigen evaluates "c /= b" as "c *= 1/b". Normally, these are the same. But in this case, b is one of the smallest numbers that fits in a float (it is twice the smallest subnormal positive number), and 1/b overflows, as the last statement shows. So Eigen computes 0 times infinity, which is NaN.

For Eigen developers: This is in SelfCwiseBinaryOp.h:188. Why do we do this? Performance? It's been there for a long time; I traced it back to revision 87a668aa3764 in November 2009 but it may well have been introduced before that.

The issue disappears if we change it to

Code: Select all

template<typename Derived>
inline Derived& DenseBase<Derived>::operator/=(const Scalar& other)
{
  typedef typename Derived::PlainObject PlainObject;
  SelfCwiseBinaryOp<internal::scalar_quotient_op<Scalar>, Derived, typename PlainObject::ConstantReturnType> tmp(derived());
  tmp = PlainObject::Constant(rows(),cols(),other);
  return derived();
}


Why is this an optimization? Does a*(1/b) compute faster than a/b?

ggael wrote:Yes, it's been like this since the beginning. There is the same optimization in Functors.h:493. I guess this is an optimization that we could leave to the user: mat * (1/a)

not for a single division, but when 'a' is a vector or matrix doing one 1/b plus n '*' with n>1 is quite faster than n '/'.

I see.