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

The behaviour of the code snippet below is different when compiled with g++-11.1.1 with different optimisation options. The snippet has floating point exception handling enabled. I don't think algorithmic details are important (it's a bit of 3D geometry). What it's doing is naturally susceptible to performing invalid floating point operations, but the code is written so that these invalid operations should never actually occur (e.g., if statements are used to ensure non-zero denominators).

#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;
}

With the following data example.dat file:

When compiled like this:

g++-11 -O3 example.cpp

Running triggers a floating point exception, which is caught, and the program exits through the signalHandler function.

When compiled like this:

g++-11 -O2 example.cpp

Or this:

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

The program does not trigger any exceptions and exits normally.

Is this a bug in the optimizer in g++-11.1.1? Earlier versions do not exhibit this behaviour; all optimisation options result in executables that exit normally. Alternatively, are vectorize optimisations like those enabled by O3 not suitable for use with floating point exception handling?

Edit: Read points from a file to ensure that the optimisation is not dependent on the values

Recommended topics

Hot tags