不同CPU的FMA指令是否有不同的中间精度？如果是，那么编译器如何平衡浮点行为？

Question

当我运行 fma 优化的 horner 方案多项式计算（用于余弦近似）时，它在 FX8150 上产生 0.161 ulps 错误，但在 godbolt.org 服务器上产生 0.154 ulps 尽管缺少 -ffast-数学（海湾合作委员会）。

如果这是由硬件引起的，并且如果每个硬件的精度不同，C++ 编译器如何保持不同机器之间的浮点精度？

编程语言规范是否只有最低精度要求，以便任何 cpu 供应商都可以根据需要提高精度？

最小可重现样本：

#include<iostream>
        // only optimized for [-1,1] input range
        template<typename Type, int Simd>
        inline
        void cosFast(Type * const __restrict__ data, Type * const __restrict__ result) noexcept
        {
            alignas(64)
            Type xSqr[Simd];
            
            for(int i=0;i<Simd;i++)
            {
                xSqr[i] =   data[i]*data[i];
            }   
            for(int i=0;i<Simd;i++)
            {
                result[i] =     Type(2.425144155360214881511638e-05);
            }
            for(int i=0;i<Simd;i++)
            {
                result[i] =     result[i]*xSqr[i] + Type(-0.001388599083010255696990498);
            }
            for(int i=0;i<Simd;i++)
            {
                result[i] =     result[i]*xSqr[i] + Type(0.04166657759826541962411284);
            }       
            for(int i=0;i<Simd;i++)
            {
                result[i] =     result[i]*xSqr[i] + Type(-0.4999999436679569697616898);
            }       
            for(int i=0;i<Simd;i++)
            {
                result[i] =     result[i]*xSqr[i] + Type(0.9999999821855363180134191);
            }


        }


#include<cstring>
template<typename T>
uint32_t GetUlpDifference(T a, T b)
{
    uint32_t aBitValue;
    uint32_t bBitValue;
    std::memcpy(&aBitValue,&a,sizeof(T));
    std::memcpy(&bBitValue,&b,sizeof(T));
    return (aBitValue > bBitValue) ?
           (aBitValue - bBitValue) :
           (bBitValue - aBitValue);
}
#include<vector>
template<typename Type>
float computeULP(std::vector<Type> real, std::vector<Type> approximation)
{
    int ctr = 0;
    Type diffSum = 0;
    for(auto r:real)
    {
        Type diff = GetUlpDifference(r,approximation[ctr++]);
        diffSum += diff;
    }
    return diffSum/ctr;
}

template<typename Type>
float computeMaxULP(std::vector<Type> real, std::vector<Type> approximation)
{
    int ctr = 0;
    Type mx = 0;
    int index = -1;
    Type rr = 0;
    Type aa = 0;
    for(auto r:real)
    {
        Type diff = GetUlpDifference(r,approximation[ctr++]);
        if(mx<diff)
        {
            mx = diff;
            rr=r;
            aa=approximation[ctr-1];
            index = ctr-1;
        }
    }
    std::cout<<"("<<index<<":"<<rr<<"<-->"<<aa<<")";
    return mx;
}
#include<cmath>
void test()
{
    constexpr int n = 8192*64;
    std::vector<float> a(n),b(n),c(n);
    for(int i=0;i<n;i++)
        a[i]=(i-(n/2))/(float)(n/2);

    // approximation
    for(int i=0;i<n;i+=16)
        cosFast<float,16>(a.data()+i,b.data()+i);

    // exact
    for(int i=0;i<n;i++)
        c[i] = std::cos(a[i]);
    
    std::cout<<"avg. ulps: "<<computeULP(b,c)<<std::endl;
    std::cout<<"max. ulps: "<<computeMaxULP(b,c)<<std::endl;
}

int main()
{
    test();
    return 0;
}

使用 FMA 的证明：

https://godbolt.org/z/Y4qYMoxcn

.L23:
    vmovups ymm3, YMMWORD PTR [r12+rax]
    vmovups ymm2, YMMWORD PTR [r12+32+rax]
    vmulps  ymm3, ymm3, ymm3
    vmulps  ymm2, ymm2, ymm2
    vmovaps ymm1, ymm3
    vmovaps ymm0, ymm2
    vfmadd132ps     ymm1, ymm7, ymm8
    vfmadd132ps     ymm0, ymm7, ymm8
    vfmadd132ps     ymm1, ymm6, ymm3
    vfmadd132ps     ymm0, ymm6, ymm2
    vfmadd132ps     ymm1, ymm5, ymm3
    vfmadd132ps     ymm0, ymm5, ymm2
    vfmadd132ps     ymm1, ymm4, ymm3
    vfmadd132ps     ymm0, ymm4, ymm2
    vmovups YMMWORD PTR [r13+0+rax], ymm1
    vmovups YMMWORD PTR [r13+32+rax], ymm0
    add     rax, 64
    cmp     rax, 2097152
    jne     .L23

这个实例（我不知道它是至强还是epyc）将它进一步提高到平均0.152 ulps。

Answer 1

对于C++语言，没有什么强烈的要求，主要是implementation-defined，正如@Maxpm在评论中指出的

floating-point 精度的主要标准是 IEEE-754。它通常被当今大多数供应商正确实施（至少几乎所有最近的主流 x86-64 CPU 和大多数主流 GPU）。它不是 C++ 标准所要求的，但您可以使用 std::numeric_limits<T>::is_iec559.

进行检查

IEEE-754 标准要求使用正确的舍入方法正确计算运算（即误差小于 1 ULP）。有不同的舍入方法 supported by the norm 但最常见的是舍入到最近的舍入。该标准还要求以相同的要求实施 FMA 等操作。因此，您不能期望使用此标准计算的结果的精度优于每次操作 1 ULP（四舍五入可能有助于平均达到 0.5 ULP 甚至更好的实际算法使用过）。

实际上，符合 IEEE-754 标准的硬件供应商的计算单元在内部使用更高的精度，以便无论提供什么输入都能满足要求。尽管如此，当结果存储在内存中时，它们需要按照 IEEE-754 的方式正确舍入。在 x86-64 处理器上，SIMD 寄存器（如 SSE、AVX 和 AVX-512 之一）具有 well-known 固定大小。对于 floating-point 操作，每个通道都是 16 位（half-float）、32 位（浮点）或 64 位（双精度）。符合 IEEE-754 的舍入应适用于每条指令。虽然处理器理论上可以实现巧妙的优化，例如将两条 FP 指令融合为一条（只要精度 <1 ULP），但 AFAIK none 做到了这一点（尽管对某些指令（如条件分支）进行了融合）。

IEEE-754 平台之间的差异可能是由于编译器或硬件供应商的 FP 单元配置所致。

关于编译器，优化可以提高精度，同时符合 IEEE-754。例如，在您的代码中使用 FMA 指令是一种提高结果精度的优化，但编译器在 x86-64 平台上并不强制执行此操作（事实上，并非所有 x86-64 处理器都支持它） .出于某些原因，编译器可能会使用单独的乘法+加法指令（Clang 有时会这样做）。编译器可以使用比目标处理器更高的精度 pre-compute 一些常量（例如，GCC 以更高的精度对 FP 数字进行操作以生成 compile-time 常量）。此外，可以使用不同的舍入方法来计算常量。

关于硬件供应商，默认舍入模式可以从一个平台更改为另一个平台。在您的情况下，非常小的差异可能是由于此。舍入模式在一个平台上可能是“舍入到最近，与偶数相等”，而在另一个平台上可能是“舍入到最近，从零开始舍入”，导致非常小但可见的差异。您可以使用 this answer. Note also that denormal numbers are sometimes disabled on some platforms because of their very-high overhead (see this 中提供的 C 代码设置舍入模式以获取更多信息），尽管这会使结果不符合 IEEE-754 标准。你应该检查一下是否是这种情况。

简而言之，两个 IEEE-754 兼容平台之间的差异 <1 ULP 是完全正常的，并且实际上在非常不同的平台之间非常频繁（例如，POWER 上的 Clang 与 x86-64 上的 GCC）。

不同CPU的FMA指令是否有不同的中间精度？如果是，那么编译器如何平衡浮点行为？

Do FMA instructions of different CPUs have different intermediate accuracy? If yes, then how does a compiler equalize the floating-point behavior?

c++

gcc

floating-point

precision

floating-accuracy