我正在使用 @njaffa 的另一篇文章中描述的测试工具，试图制作一个更快、更准确的 erfc() 函数，这时我遇到了一个我仍然无法完全理解的特殊错误，并且需要一些帮助。我认为代码在突破某些界限的同时表现得非常好。
纯粹主义者可能会反对循环终止时计数器 argi 溢出，但这不是问题。当使用通常的 go fast 选项（特别是 /fp:fast）编译时，生成的代码是错误的。我添加了一些不影响问题的诊断并计算检测到的 nan 数量（但给出不一致的结果）。
NB 它仅在英特尔编译器上失败，并且必须设置为 /fp:fast 并发布（所有其他选项都可以正常工作，就像我尝试过的所有其他编译器一样）。 FP 严格和精确都可以在 Intel 上正常工作。
我尽可能地简化了他的代码，同时仍然保留了我在这里看到的错误行为，并且在下面以失败状态复制了它：
#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 下运行。
症状是，当 optimization 和 fp:fast 一起使用时，由于一些奇怪的编译器怪癖，absmaxerr 最终会出现 nan。我可以在某种程度上跟踪反汇编代码的执行，并用独特的负范数标记 max____ 变量，以使它们在 XMMn 寄存器集中脱颖而出。但我一直无法弄清楚它是如何发生的。我认为这可能是因为用 nans 执行库代码（测试所有值）。
代码的微小更改会导致生成的代码（和/或寄存器分配）发生根本性变化。特别是，必须至少执行一次循环分支才能显示问题。如果最后使用 while(false) ，编译器会生成完全不同（且有效）的代码，同样，如果防御测试或 printf 被放置在关键位置。
abserr is nan 有两项测试 - 一项在外循环级别处于活动状态，触发 0x7f800000
我的预感是，当 argi 是 nan 时，最大abs值更新逻辑和 maxabserr 的存储之间存在一些特殊的交互，但我无法确定它。 abserrloc 在第一个循环后不会改变并保持在 -4.0，但最初正确的 maxabserr 被践踏。感谢您的启发。
潜在问题：
>C 没有指定
if (abserr > maxabserr)为了更便携的功能而建议的更改：
if (isnan(abserr)) {
nancount++;
printf("NaN @ %8x ", argi);
}
// if (abserr > maxabserr) {
else if (abserr > maxabserr) {
abserrloc = arg;
maxabserr = abserr;
}