我正在玩@njaffa在另一篇帖子中描述的测试工具，试图创建更快、更准确的erfc()函数，这时我遇到了一个特殊的错误，我仍然无法完全理解它，希望得到一些帮助.我认为代码在突破某些界限的同时表现得非常好.
纯粹主义者可能会反对溢出计数器argi以终止循环，但这不是问题所在.当使用通常的GO FETHER选项编译时，特别是/fp:fast，生成的代码是错误的.我添加了一些不会影响问题的诊断，并计算了检测到的NAN的数量(但得出的结果不一致).
注意:只有在英特尔编译器上才会失败，必须将其设置为/fp:fast并释放(所有其他选项都能正常工作，我try 过的所有其他编译器也是如此).FP Strong和Precision在英特尔上都能正常工作.
我已经尽可能地简化了他的代码，同时仍然保留了我在这里看到的错误行为，它以失败的状态重现如下:
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
// njuffa's fast erfc(x) and test harness (stripped down as far as I can) see
// https://stackoverflow.com/questions/77741402/fast-implementation-of-complementary-error-function-with-5-digit-accuracy/
// It only fails on the latest Intel 2022 compiler when /fp:fast and /O2 /Ob2 /Ot
/* Fast computation of the complementary error function. For argument x > 0
erfc(x) = exp(-x*x) * P(x) / Q(x), where P(x) and Q(x) are polynomials.
If expf() is faithfully rounded, the following error bounds should hold:
Maximum relative error < 1.065e-5, maximum absolute error < 9.50e-6, and
maximum ulp error < 176.5
*/
float fast_erfcf (float x)
{
float a, c, e, p, q, r, s;
a = fabsf (x);
c = fminf (a, 10.5f);
s = -c * c;
e = expf (s);
q = 3.82346272e-1f; // 0x1.8785c8p-2f
p = -4.38243151e-5f; // -0x1.6fa000p-15
q = q * c + 1.30382371e+0f; // 0x1.4dc764p+0
p = p * c + 2.16852218e-1f; // 0x1.bc1d04p-3
q = q * c + 1.85278797e+0f; // 0x1.da5050p+0
p = p * c + 7.23953605e-1f; // 0x1.72aa0cp-1
q = q * c + 9.99991596e-1f; // 0x1.fffee6p-1
p = p * c + 1.00000000e+0f; // 0x1.000000p+0
r = e / q;
r = r * p;
if (x < 0.0f) r = 2.0f - r;
if (isnan(x)) r = x + x;
return r;
}
float uint32_as_float (uint32_t a)
{
float r;
memcpy (&r, &a, sizeof r);
return r;
}
uint32_t float_as_uint32 (float a)
{
uint32_t r;
memcpy (&r, &a, sizeof r);
return r;
}
uint64_t double_as_uint64 (double a)
{
uint64_t r;
memcpy (&r, &a, sizeof r);
return r;
}
double floatUlpErr (float res, double ref)
{
uint64_t refi, i, j, err;
int expoRef;
/* ulp error cannot be computed if either operand is NaN, infinity, zero */
if (isnan (res) || isnan (ref) || isinf(res) || isinf (ref) ||
(res == 0.0f) || (ref == 0.0)) {
return 0.0;
}
i = ((int64_t)float_as_uint32 (res)) << 32;
expoRef = (int)(((double_as_uint64 (ref) >> 52) & 0x7ff) - 1023);
refi = double_as_uint64 (ref);
if (expoRef >= 129) {
j = (refi & 0x8000000000000000ULL) | 0x7fffffffffffffffULL;
} else if (expoRef < -126) {
j = ((refi << 11) | 0x8000000000000000ULL) >> 8;
j = j >> (-(expoRef + 126));
j = j | (refi & 0x8000000000000000ULL);
} else {
j = ((refi << 11) & 0x7fffffffffffffffULL) >> 8;
j = j | ((uint64_t)(expoRef + 127) << 55);
j = j | (refi & 0x8000000000000000ULL);
}
err = (i < j) ? (j - i) : (i - j);
return err / 4294967296.0;
}
int main (void)
{
uint32_t argi, last_argi;
float arg, res, reff, abserrloc = NAN, ulperrloc = NAN;
double ref, abserr, ulperr;
double maxabserr = -0x1.2468Ap-127, maxulperr = -0x1.13579p-127; // tagged to spot in MMXreg
unsigned int nancount = 0, loopcount = 0;
// argi = 0; // original code to test entire range - uses XMM8h,XMM10h & fails
argi = float_as_uint32(-4.0); // this is a convenient hard failing value aka 0xc0800000 - uses XMM10h, XMM13h & fails
// argi = 0xff810000-0x007f0000; // this is the closest start value I found that fails with increment 0x007f0000
do {
if (argi == 0xff810000)
printf("This time it will fail\n");
arg = uint32_as_float (argi);
ref = erfc ((double)arg);
res = fast_erfcf (arg);
ulperr = floatUlpErr (res, ref);
if (ulperr > maxulperr) {
ulperrloc = arg;
maxulperr = ulperr;
}
abserr = fabs (res - ref);
if (isnan(abserr)) {
nancount++; // testing for nan here the bug persists
printf("NaN @ %8x ", argi); // this line will trigger breakpoint here
}
if (abserr > maxabserr) {
// if (isnan(abserr)) printf("abserr is Nan/n"); // testing for nan here makes results correct
abserrloc = arg;
maxabserr = abserr;
// printf("abserr %g @ %g ", abserr, arg); // this line also makes results correct
}
last_argi = argi;
// argi++; // original code with Nan (0xffffffff)
// argi += 0x007f0000; // the "nan" reported depends on this increment!
argi += 0x00400000;
// 0x007F0000 is the largest increment that is still a Nan (0x7fc100000)
// 0x00400000 is Nan (0xffc00000)
// 0x00200000 is Nan (0xffe00000)
// 0x00100000 is Nan (0xfff00000) etc.
loopcount++;
} while (argi>last_argi); (false); // changing this line to prevent looping generates correct code
// NB it is radically different to the looping code (that fails).
printf ("max abs err = %.6e %.6a %x @ % 15.8e %.6a\n"
"max ulp err = %.6e @ % 15.8e\n",
maxabserr, maxabserr, float_as_uint32(maxabserr), abserrloc, abserrloc,
maxulperr, ulperrloc);
printf("Nan count %i / %i\n", nancount, loopcount);
return EXIT_SUCCESS;
}
如果它按预期工作，则输出应为(舍入误差或取舍误差):
NaN @ ffc00000
max abs err = 1.541726e-08 0x1.08ddd1p-26 32846ee9 @ -4.00000000e+00 -0x1.000000p+2
max ulp err = 1.293293e-01 @ -4.00000000e+00
使用MS VC IDE的命令行选项包括:
/GS /W3 /Gy /Zi /O2 /Ob2 /fp:fast /D "NDEBUG" /D "_CONSOLE" /D "_UNICODE" /D "UNICODE"
/Qipo /Zc:forScope /arch:CORE-AVX2 /Oi /MD /Fa"x64\Release\" /EHsc /nologo
/Fo"x64\Release\" //fprofile-instr-use "x64\Release\" /Ot
/Fp"x64\Release\erfc_njaffa.pch"
因为它可能是相关的，CPU是i5-12600K AlderLake Family 6 Model 7 Step 2 Ext Model 97 Rev C0，我运行的是Win 11 Pro.
症状是，当优化和fp:fast一起使用时，由于一些奇怪的编译器怪癖，absmaxerr最终会在其中包含一个nan.我可以在一定程度上跟踪反汇编代码的执行，并使用独特的负反规范标记max_变量，使它们在XMMn寄存器集中脱颖而出.但我还没弄清楚这是怎么发生的.我认为这可能是因为使用nans(测试所有值)执行库代码.
对代码的微小更改会导致生成的代码(和/或寄存器分配)发生根本变化.特别是，必须至少使用一次循环分支才能显示问题.如果在最后使用While(FALSE)，编译器会生成完全不同的(和正常工作的)代码，同样，如果在关键位置放置防御性测试或printf，编译器也会生成完全不同的代码.
对于abserr is NaN有两个测试-一个在外循环级别活动，它触发0x7f800000<；argi<；=0xff800000<；argi<；=0xffffffff(当argi本身是NAN时，IOW).另一个被注释掉，将代码置于失败状态，并位于abserr；Maxabserr路径中.如果在那里启用，它永远不会触发，代码正常工作，而外部测试继续看到NaN的.
我的直觉是，当argi是NaN时，最大abs值更新逻辑和Maxabserr的存储之间存在一些特殊的交互作用，但我就是不能确定它.Abserrloc在第一次循环后不会改变，并保持在-4.0，但最初正确的Maxabserr会被践踏.谢谢你的启发.