当我运行fma优化的角方案多项式计算(用于余弦近似)时，尽管缺乏- fast-math (GCC)，但它在FX8150上的误差为0.161 ulps，而在godbolt.org服务器上的误差为0.154 ulps。
如果这是由硬件引起的，并且每个硬件的精度不同，那么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
这个实例(我不知道是xeon还是epyc)进一步提高到平均0.152秒。
关于c++语言，没有很强的要求，它主要是实现定义的，正如@Maxpm在评论中指出的前面的答案所述。
浮点精度的主要标准是IEEE-754。现在大多数供应商(至少几乎所有最近的主流x86-64 cpu和大多数主流gpu)都正确地实现了它。这不是c++标准所要求的，但您可以通过std::numeric_limits<T>::is_iec559检查。
IEEE-754标准要求正确计算操作(即。误差小于1 ULP)，使用正确的舍入方法。标准支持不同的舍入方法，但最常见的是舍入到最接近的舍入。该标准还要求一些操作，如FMA，以相同的要求来实施。因此，不能期望以优于1 ULP的精度计算每个操作的结果使用此标准(四舍五入可能有助于平均达到0.5 ULP，甚至对于使用的实际算法更好)。
在实践中，兼容ieee -754硬件供应商的计算单元在内部使用更高的精度，以便满足所提供输入的要求。但是，当结果存储在内存中时，它们需要按照IEEE-754的方式进行四舍五入。在x86-64处理器上，像SSE、AVX和AVX-512这样的SIMD寄存器有一个众所周知的固定大小。对于浮点操作，每个通道都是16位(半浮点)，32位(浮点)或64位(双精度)。每个指令都应采用符合IEEE-754标准的舍入。虽然处理器理论上可以实现聪明的优化，如融合两个FP指令在一个(只要精度是<1 ULP)， AFAIK还没有这样做(虽然融合一些指令，如条件分支)。
IEEE-754平台之间的差异可能是由于编译器或硬件供应商的FP单元配置所致。
对于编译器，优化可以在满足IEEE-754要求的同时提高精度。例如，在代码中使用FMA指令是一种优化，可以提高结果的精度，但在x86-64平台上编译器并不强制这样做(事实上，并非所有x86-64处理器都支持它)。编译器可能出于某些原因使用单独的乘法+加法指令(Clang有时会这样做)。编译器可以使用比目标处理器更好的精度预计算一些常量(例如，GCC以更高的精度操作FP数以生成编译时常量)。此外，可以使用不同的舍入方法来计算常数。
对于硬件供应商，由于默认舍入模式可以从一个平台更改到另一个平台。在你的情况下，非常小的差异可能是由于这个原因。四舍五入模式可以是"四舍五入，四舍五入"。在一个平台上，"轮到最近，从零开始";在另一个平台上，产生了非常小但明显的差异。您可以使用答案中提供的C代码设置舍入模式。还要注意，异常数字有时在某些平台上被禁用，因为它们的开销非常高(请参阅此以获取更多信息)，尽管它使结果不符合IEEE-754。你应该检查一下是不是这样。
简而言之，两个符合ieee -754标准的平台之间的差异<1 ULP是完全正常的，实际上在非常不同的平台之间是相当频繁的(例如;Clang on POWER vs GCC on x86-64)。