usage of signed zeros leads to wrong results with -ffast-math

found in files:

Eigen/src/Geometry/arch/Geometry_SIMD.h
Eigen/src/Core/arch/MSA/MathFunctions.h
Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
Eigen/src/Core/arch/ZVector/PacketMath.h
Eigen/src/Core/arch/AltiVec/PacketMath.h
Eigen/src/Core/arch/NEON/PacketMath.h
Eigen/src/LU/arch/InverseSize4.h

this was the buggy program that led to discovering this

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

using std::cout;
using std::endl;

int main() {
  Eigen::Matrix<double, 4, 4, Eigen::RowMajor> M;
  M << 1.275, 0, 0, 0, 0, 0.910, 0, -0, 0, -0.159, 0.923, 0, 0, -0.137, -0.137,
      0.934;

  Eigen::Matrix<double, 4, 4, Eigen::RowMajor> Minv = M.inverse();

  cout << "M:" << endl;
  cout << M << endl << endl;
  cout << "Minv:" << endl;
  cout << Minv << endl << endl;
  cout << "M*Minv:" << endl;
  cout << M * Minv << endl << endl;
}

the result ends up being wrong when compiled on g++-{9,10} -O3 -DNDEBUG -ffast-math
i believe this is due to the reliance on bit tricks for flipping the signs of floats in InverseSize4.h

    const double sign_mask1[2] = {0.0, -0.0};
    const double sign_mask2[2] = {-0.0, 0.0};
    const Packet2d sign_PN = pset<Packet2d>(sign_mask1);
    const Packet2d sign_NP = pset<Packet2d>(sign_mask2);
    d1 = pxor(rd, sign_PN);
    d2 = pxor(rd, sign_NP);
    // printing d1 and d2 at this point shows them to be equal
    // since the compiler can ignore the signs on zeros,
    // and ends up using `sign_PN` for both