矢量化优化是否与浮点异常处理不兼容，或者 g++-11 中是否存在错误？

Question

当使用 g++-11.1.1 编译并使用不同的优化选项时，下面代码片段的行为是不同的。该代码段启用了浮点异常处理。我认为算法细节并不重要（它有点像 3D 几何）。它正在做的事情自然容易执行无效的浮点运算，但代码的编写使得这些无效运算永远不会实际发生（例如，if 语句用于确保非零分母）。

#include <array>
#include <cfenv>
#include <csignal>
#include <fstream>
#include <iostream>
#include <limits>
#include <cmath>
#include <vector>

// Vector structure -----------------------------------------------------------
struct Vector
{
    double xs[3];

    Vector(double x)
    {
        xs[0] = x;
        xs[1] = x;
        xs[2] = x;
    }

    Vector(double x, double y, double z)
    {
        xs[0] = x;
        xs[1] = y;
        xs[2] = z;
    }

    Vector(std::istream& is)
    {
        is >> xs[0] >> xs[1] >> xs[2];
    }
};

// Vector output stream operator ----------------------------------------------
inline std::ostream& operator<<(std::ostream& os, const Vector& v)
{   
    return os << '(' << v.xs[0] << ' ' << v.xs[1] << ' ' << v.xs[2] << ')';
}

// Vector geometry operators --------------------------------------------------
inline void operator+=(Vector& a, const Vector& b)
{
    a.xs[0] += b.xs[0];
    a.xs[1] += b.xs[1];
    a.xs[2] += b.xs[2];
}
inline void operator/=(Vector& a, const double b)
{
    a.xs[0] /= b;
    a.xs[1] /= b;
    a.xs[2] /= b;
}
inline Vector operator+(const Vector& a, const Vector& b)
{   
    return Vector(a.xs[0] + b.xs[0], a.xs[1] + b.xs[1], a.xs[2] + b.xs[2]);
}
inline Vector operator-(const Vector& a, const Vector& b)
{
    return Vector(a.xs[0] - b.xs[0], a.xs[1] - b.xs[1], a.xs[2] - b.xs[2]);
}
inline Vector operator*(const double& a, const Vector& b)
{
    return Vector(a*b.xs[0], a*b.xs[1], a*b.xs[2]);
}
inline Vector operator/(const Vector& a, const double& b)
{
    return Vector(a.xs[0]/b, a.xs[1]/b, a.xs[2]/b);
}
inline double operator&(const Vector& a, const Vector& b)
{
    return a.xs[0]*b.xs[0] + a.xs[1]*b.xs[1] + a.xs[2]*b.xs[2];
}
inline Vector operator^(const Vector& a, const Vector& b)
{
    return Vector
    (
        a.xs[1]*b.xs[2] - a.xs[2]*b.xs[1],
        a.xs[2]*b.xs[0] - a.xs[0]*b.xs[2],
        a.xs[0]*b.xs[1] - a.xs[1]*b.xs[0]
    );
}

// Polygon centre algorithm ---------------------------------------------------
template<class PointList>
typename PointList::value_type polygonCentre(const PointList& ps)
{
    typedef typename PointList::value_type Point;

    // Compute an estimate of the centre as the average of the points
    Point pAvg(0);
    for (typename PointList::size_type pi = 0; pi < ps.size(); ++ pi)
    {   
        pAvg += ps[pi];
    }
    pAvg /= ps.size();

    // Compute the polygon area normal and unit normal by summing up the
    // normals of the triangles formed by connecting each edge to the
    // point average.
    Point sumA(0);
    for (typename PointList::size_type pi = 0; pi < ps.size(); ++ pi)
    {
        const Point& p = ps[pi];
        const Point& pNext = ps[(pi + 1) % ps.size()];

        const Point a = (pNext - p)^(pAvg - p);
        
        sumA += a;
    }
    double sumAMagSqr = sumA & sumA; 
    const Point sumAHat = sumAMagSqr > 0 ? sumA/sqrt(sumAMagSqr) : Point(0);

    // Compute the area-weighted sum of the triangle centres
    double sumAn = 0;
    Point sumAnc(0);
    for (typename PointList::size_type pi = 0; pi < ps.size(); ++ pi)
    {
        const Point& p = ps[pi];
        const Point& pNext = ps[(pi + 1) % ps.size()];

        const Point a = (pNext - p)^(pAvg - p);
        const Point c = p + pNext + pAvg;

        const double an = a & sumAHat;

        sumAn += an;
        sumAnc += an*c;
    }

    // Complete calculating centres and areas. If the area is too small
    // for the sums to be reliably divided then just set the centre to
    // the initial estimate.
    if (sumAn > std::numeric_limits<double>::min())
    {   
        return (1.0/3.0)*sumAnc/sumAn;
    }
    else
    {
        return pAvg;
    }
}

// Signal handler -------------------------------------------------------------
void signalHandler(int signum)
{
    std::cout << "Signal " << signum << " caught." << std::endl;
    exit(1);
}

// Main routine ---------------------------------------------------------------
int main(int argc, char *argv[])
{
    feenableexcept(FE_INVALID);

    signal(SIGFPE, signalHandler);

    /*
    std::array<Vector, 4> ps
    ({
        Vector(0, 0, 0),
        Vector(1, 0, 0),
        Vector(1, 0, 0),
        Vector(0, 0, 0)
    });
    */

    std::ifstream is("example.dat");

    std::array<Vector, 4> ps
    ({
        Vector(is),
        Vector(is),
        Vector(is),
        Vector(is)
    });

    std::cout << "Centre = " << polygonCentre(ps) << std::endl;

    return 0;
}

具有以下数据example.dat文件：

像这样编译时：

g++-11 -O3 example.cpp

运行触发浮点数异常，被捕获，程序通过signalHandler函数退出

像这样编译时：

g++-11 -O2 example.cpp

或者这样：

g++-11 -O3 -fno-tree-slp-vectorize example.cpp

程序没有触发异常，正常退出

这是 g++-11.1.1 优化器中的错误吗？早期版本没有这种行为；所有优化选项都会导致可执行文件正常退出。或者，像 O3 启用的向量化优化是否不适合用于浮点异常处理？

编辑：从文件中读取点以确保优化不依赖于值

Answer 1

这是 Gcc-11.1 中的错误。它已在 11.2 的候选发布版本中修复。

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=101634

矢量化优化是否与浮点异常处理不兼容，或者 g++-11 中是否存在错误？

Are vectorize optimisations incompatible with floating point exception handling, or is there a bug in g++-11?

c++

g++

floating-point-exceptions