msh*_*ldt 1 c floating-point unit-testing sse floating-accuracy
我正在尝试为一些简单的向量数学函数编写单元测试,这些函数在单精度浮点数的数组上运行.这些函数使用SSE内在函数,当在32位系统上运行测试时,我会得到误报(至少我认为)(测试通过64位).当操作通过数组运行时,我累积越来越多的舍入错误.这是一个单元测试代码和输出的片段(我的实际问题如下):
测试设置:
static const int N = 1024;
static const float MSCALAR = 42.42f;
static void setup(void) {
input = _mm_malloc(sizeof(*input) * N, 16);
ainput = _mm_malloc(sizeof(*ainput) * N, 16);
output = _mm_malloc(sizeof(*output) * N, 16);
expected = _mm_malloc(sizeof(*expected) * N, 16);
memset(output, 0, sizeof(*output) * N);
for (int i = 0; i < N; i++) {
input[i] = i * 0.4f;
ainput[i] = i * 2.1f;
expected[i] = (input[i] * MSCALAR) + ainput[i];
}
}
然后我的主要测试代码调用要测试的函数(用于生成expected数组的相同计算)并根据expected上面生成的数组检查其输出.检查是为了接近(在0.0001内)不平等.
样本输出:
0.000000 0.000000 delta: 0.000000
44.419998 44.419998 delta: 0.000000
...snip 100 or so lines...
2043.319946 2043.319946 delta: 0.000000
2087.739746 2087.739990 delta: 0.000244
...snip 100 or so lines...
4086.639893 4086.639893 delta: 0.000000
4131.059570 4131.060059 delta: 0.000488
4175.479492 4175.479980 delta: 0.000488
...etc, etc...
我知道我有两个问题:
所以我的问题是,为浮点数据的数学运算编写有意义且可移植的单元测试的正确方法是什么？
*便携式我的意思是应该传递32位和64位架构.
根据评论,我们看到正在测试的函数基本上是:
for (int i = 0; i < N; ++i)
D[i] = A[i] * b + C[i];
where A[i],b和C[i],D[i]都有类型float.当提到一个迭代的数据,我会用a,c和d为A[i],C[i],和D[i].
下面是对测试此函数时可用于容错的内容的分析.首先,我想指出我们可以设计测试,以便没有错误.我们可以选择的值A[i],b,C[i],并D[i]让所有的结果,无论最终和中间结果,是完全表示,没有舍入误差.显然,这不会测试浮点运算,但这不是目标.目标是测试函数的代码:它是否执行计算所需函数的指令？只需选择能够揭示使用正确数据的任何失败,添加,乘法或存储到正确位置的值,就足以揭示函数中的错误.我们相信硬件正确地执行浮点并且没有测试它; 我们只想测试函数是否正确编写.为了实现这一点,我们可以,例如,设置b为2的幂,A[i]各种小整数,以及C[i]各种小整数乘以b.如果需要,我可以更准确地详细说明这些值的限制.那么所有结果都是准确的,任何比较允许容差的需要都会消失.
除此之外,让我们继续进行错误分析.
目标是找到函数实现中的错误.为此,我们可以忽略浮点运算中的小错误,因为我们所寻求的错误种类几乎总是会导致大错误:使用了错误的操作,使用了错误的数据,或者结果没有存储在期望的位置,因此实际结果几乎总是与预期结果非常不同.
现在的问题是我们应该容忍多少错误？由于错误通常会导致较大的错误,因此我们可以将容差设置得相当高.但是,在浮点数中,"高"仍然是相对的; 与数万亿的值相比,100万的误差很小,但是当输入值在1中时,发现错误太高.所以我们至少应该做一些分析来决定水平.
正在测试的功能将使用SSE内在函数.这意味着,对于i上面的循环中的每一个,它将执行浮点乘法和浮点加法,或者将执行融合浮点乘法加法.后者的潜在错误是前者的一个子集,所以我将使用前者.浮点运算用于a*b+c进行一些舍入,以便计算大约为a•b + c的结果(解释为精确的数学表达式,而不是浮点数).如果所有值都在浮点格式的正常范围内,我们可以写出(a•b•(1+e0)+c)•(1+e1)针对某些误差e0和e1 计算的精确值,其大小为2到24.(2 -24是在IEEE-754 32位二进制浮点格式中以舍入到最接近模式的任何正确舍入的基本浮点运算中可能发生的最大相对误差.舍入到舍入到最接近的模式将数学值最多改为有效数中最低有效位的一半,即最高有效位以下23位.)
接下来,我们考虑测试程序为其预期值产生的值.它使用C代码d = a*b + c;.(我已将问题中的长名称转换为较短的名称.)理想情况下,这也会计算IEEE-754 32位二进制浮点的乘法和加法.如果确实如此,那么结果将与正在测试的功能相同,并且不需要允许任何容差进行比较.但是,C标准允许实现在执行浮点运算时具有一定的灵活性,并且存在比标准允许更多自由的不符合实现.
一种常见的行为是表达式的计算精度高于其标称类型.一些编译器可以a*b + c使用double或long double算术计算.C标准要求将结果转换为演员表或作业中的名义类型; 必须丢弃额外的精度.如果C实现使用额外的精度,那么计算将继续:a*b以额外的精度计算,精确地产生a•b,因为double并且long double具有足够的精度来精确地表示任何两个float值的乘积.然后,AC实现可以将此结果舍入到float.这不太可能,但无论如何我还是允许的.但是,我也忽略它,因为它将预期结果移动到更接近被测函数的结果,我们只需要知道可能发生的最大错误.所以我会继续,更糟糕(更遥远)的情况,到目前为止的结果是a•b.然后c添加,产生(a•b + c)•(1 + e2)的某些e2,其幅度最多为2-53(64位二进制格式中正常数字的最大相对误差).最后,将此值转换float为赋值给d(a•b + c)•(1 + e2)•(1 + e3),对于某些幅度最大为2 -24的 e3 .
现在我们有正确运算函数计算的精确结果的表达式,(a•b•(1 + e0)+ c)•(1 + e1),对于测试代码计算的精确结果,(a•b) + c)•(1 + e2)•(1 + e3),我们可以计算它们可以有多大差异的界限.简单代数告诉我们确切的区别是a•b•(e0 + e1 + e0•e1-e2-e3-e2•e3)+ c•(e1-e2-e3-e2•e3).这是e0,e1,e2和e3的简单函数,我们可以看到它的极值发生在e0,e1,e2和e3的潜在值的端点.由于值的符号可能性之间的相互作用,存在一些复杂性,但是对于最坏的情况,我们可以简单地允许一些额外的错误.差值的最大值的界限是| a•b |•(3•2 -24 +2 -53 +2 -48)+ | c |•(2•2 -24 +2 -53 +2 -77).
因为我们有足够的空间,所以我们可以简化这一点,只要我们在使价值更大的方向上这样做.例如,使用| a•b |•3.001•2 -24 + | c |•2.001•2 -24可能很方便.该表达式应该足以允许在浮点计算中进行舍入,同时检测几乎所有的实现错误.
请注意,表达式与最终值不成比例,a*b+c由正在测试的函数或测试程序计算得出.这意味着,通常,使用相对于被测试函数或测试程序计算的最终值的公差的测试是错误的.适当的测试形式应该是这样的:
double tolerance = fabs(input[i] * MSCALAR) * 0x3.001p-24 + fabs(ainput[i]) * 0x2.001p-24;
double difference = fabs(output[i] - expected[i]);
if (! (difference < tolerance))
// Report error here.
总之,这给了我们一个容差,它比浮点舍入更大的任何可能的差异,所以它永远不应该给我们一个误报(报告测试函数不会被破坏).但是,与我们想要检测的错误导致的错误相比,它非常小,所以它应该很少给我们一个假阴性(无法报告实际错误).
(注意,还有计算公差的舍入误差,但它们小于我在系数中使用.001所允许的斜率,所以我们可以忽略它们.)
(另请注意,这! (difference < tolerance)不等同difference >= tolerance.如果函数产生NaN,由于错误,任何比较都会产生错误:both difference < tolerance和difference >= toleranceyield都会生成false,但会! (difference < tolerance)产生true.)