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

Question

当使用具有不同优化选项的g++-11.1.1进行编译时，以下代码片段的行为是不同的。该代码段启用了浮点异常处理。我不认为算法细节是重要的(它是一个3D几何)。它所做的事情自然容易受到执行无效浮点操作的影响，但编写代码的目的是为了确保这些无效操作永远不会实际发生(例如，如果使用语句来确保非零分母)。

#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文件：

0 0 0
1 0 0
1 0 0
0 0 0

当像这样编译时：

g++-11 -O3 example.cpp

运行时会触发一个浮点异常，该异常会被捕获，然后程序通过signalHandler函数退出。

当像这样编译时：

g++-11 -O2 example.cpp

或者这样：

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

程序不会触发任何异常，正常退出。

这是g++-11.1.1中优化器中的一个bug吗？早期版本不显示此行为；所有优化选项都会导致可执行文件正常退出。或者，像O3支持的那些向量化优化是否不适合用于浮点异常处理？

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

Answer 1

这是Gcc-11.1中的一个错误。在11.2的发布候选版本中已经修复了这个问题。

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