我在 CI 期间注意到此代码
#include <stdio.h>
#include <math.h>
int main()
{
const double d = 3.81219767052986080458;
fprintf(stderr, "%.20f\n", sin(d));
return 0;
}
编译后,gcc test.c -lm -o test
对于我的两个测试平台,返回的结果略有不同。
对于 gcc 版本 13.2.0 libc 2.39-0ubuntu8.3 Kubuntu 24.04
,结果为 -0.621460098738919 27544 ;
对于 gcc 版本 12.2.0 libc 2.36-9+deb12u7+ci202405171200+astra5+b1 Astra Linux 1.8
,结果为-0.621460098738919 38646
我想知道是什么影响了它?libm?gcc?可能是 CPU?
这种细微的差异可能是在紧密循环的集成过程中累积起来的。
有没有一种简单的方法可以摆脱它而不会影响性能?例如,gcc 编译器参数?
11
最佳答案
2
我想知道什么会影响它?libm?
对于相同输入,的差异sin
通常是由标准 C 库中的 的不同实现引起的sin
。 的实现sin
可能位于名为 的“数学库”中libm
,也可能位于其他文件中,具体取决于平台。
(另一种可能性是,输入数字3.81219767052986080458
在不同的编译器版本中被解析得不同,但这种可能性较小。)
理想情况下,sin
实现会返回等于其输入操作数的数学正弦值的数字,该数字四舍五入为浮点格式可表示的最接近值。1这将产生您显示的第一个结果,或者更准确地说,-.62146009873891927544065083566238172352313995361328125。但是,C 标准并不要求这样做;它允许实现产生未正确四舍五入的结果。
对于某些数学库函数,人类尚不清楚如何在已知界限的执行时间内计算出具有正确舍入的整个函数。我们可以使用扩展精度算法将函数近似到任何所需的精度,并且可以编写迭代算法来计算越来越多的精度,直到结果非常准确,以至于我们知道应该如何舍入。但是,某些函数的某些结果可能非常接近舍入变化的点(可表示值的中间值),因此需要很高的精度。我们尚不知道所有函数都需要多少精度。因此,我们可以编写一个在有足够时间的情况下有效的迭代算法,但是,对于某些函数,我们不知道需要多少时间来覆盖所有情况。
许多用户无法容忍数学库例程在某些情况下运行时间过长。他们希望例程能在较短的固定或有限的时间内运行。同样,目前人类不可能实现所有数学库函数,因为我们不知道所有数学库函数在最坏情况下需要多少精度。
float
对于 IEEE-754 binary32 和 binary64 格式(通常用于和)的正弦函数,情况并非如此double
。它的所有最坏情况都已完全映射。因此,可以实现正确舍入的正弦函数,并且可以以合理的速度完成,尽管有些人更喜欢没有正确舍入的更快实现。
正在致力于进一步开发正确舍入的实现。
即使有人选择不实现正确舍入的sin
,3.81219767052986080458 中的差异仍然会让人略感意外,因为它太小了。对于周期函数,通常的设计是将基区间之外的参数减少到基区间,然后求一个近似于基区间内函数的多项式。对于正弦,实现通常至少使用一点额外的精度来减少基区间并进行多项式求值,对于这种情况,一点额外的精度就足够了。我确实注意到,双精度最接近的 3.81219767052986080458 的正弦距离两个最接近的可表示值之间的中点仅约 0.2%。因此,八九个额外的位应该足够了。但也许您的某个正弦实现没有这样做。
有没有简单的方法可以摆脱它而不影响性能?
如果问题出在sin
实现上,则必须使用不同的数学库。其性能可能会有所不同。
脚注
1对于一般函数,舍入会打破平局,取其有效数字中最低位数为偶数(二进制浮点数的位)的最接近可表示值。这与正弦无关,因为其唯一可用有理正弦表示的操作数是零,因此所有操作数均不涉及平局。
|
此类调查以十六进制表示,%a
可以阐明:
int main() {
double x;
x = 3.81219767052986080458;
printf("x: % -25a % .21g\n", x, x);
double y = sin(x);
printf("sin(x): % -25a % .21g\n", y, y);
long double yl = sinl(x);
printf("sinl(x): % -25La % .21Lg\n", yl, yl);
y = -0.62146009873891927544;
printf("13.2.0 : % -25a % .21g\n", y, y);
y = -0.62146009873891938646;
printf("12.2.0: % -25a % .21g\n", y, y);
return 0;
}
输出:
x: 0x1.e7f617e06814dp+1 3.81219767052986080458
sin(x): -0x1.3e30049fb486ap-1 -0.621460098738919386463
sinl(x): -0x1.3e30049fb48697f8p-1 -0.621460098738919330735
13.2.0 : -0x1.3e30049fb4869p-1 -0.621460098738919275441
12.2.0: -0x1.3e30049fb486ap-1 -0.621460098738919386463
由此我们可以看出
- 代码x =得出的结果非常接近。超过 15 位有效小数的代码值值得怀疑。这个结果非常接近 – 准确到 21 位。
3.81219767052986080458;
double
x
包括sin()
我自己的这 2 个系统,答案的差异都在最后一个单位 (ULP) 以内。
使用long double
数学运算,我们得到sinl()
-0x1.3e30049fb4869_7f8p-1。请注意,7f8 几乎位于sin()
OP 报告的两个系统较低精度结果的中间。-0x1.3e30049fb4869p-1 大约接近 100 分之 4,因此是更好的答案。
我想知道是什么影响了它?libm?gcc?可能是 CPU?
许多因素可能会影响这一点,包括运行时数学(FLT_EVAL_METHOFD 的报告值)、优化设置、使用的库甚至 CPU。
有没有一种简单的方法可以摆脱它而不会影响性能?例如 gcc 编译器参数
不。
—
[更多的]
虽然只有几行,但可能有 3 处不准确之处。
double d = 3.81219767052986080458;
21 位十进制常数被转换为double
。尽管从十进制 FP 常数到 的转换double
并不精确,但在这里已经足够接近,并且对最终结果没有明显的影响。
x before: 0x1.e7f617e06814_cp+1 3.81219767052986_036049...
x: 0x1.e7f617e06814_dp+1 3.81219767052986_0804577992894337512552738189697265625
code: 3.81219767052986_080458
x after: 0x1.e7f617e06814_ep+1 3.81219767052986_124867...
sin(d)
是超越函数。即使x
精确,结果也几乎总是不是有理数(因此不能精确地表示为double
)。好的将返回与正确答案相差约 1double sin(double)
的结果。非常好的将返回与正确答案结果。对于常见的,我们可以看到 OP 找到的 2 个答案都与更精确的答案相差近 0.5 ULP:一个在上方,一个在下方。sin (x) 13.2.0比更好的答案低约 0.46 ULP,而sin(x) 12.2.0比更好的答案高约 0.54 ULP。double sin(double)
double
// x from above.
sin(x) 13.2.0 : -0x1.3e30049fb4869000p-1 -0.621460098738919_275441...
sinl(x): -0x1.3e30049fb48697f8p-1 -0.621460098738919_330735...
sin(x) 12.2.0: -0x1.3e30049fb486a000p-1 -0.621460098738919_386463...
我怀疑这是因为 13.2.0 中的库函数更好。然而,即使这个例子有点糟糕,但新库的平均水平可能更好。
fprintf(stderr, "%.20f\n", sin(d));
-0x1.3e30049fb4869000p-1
和的转换-0x1.3e30049fb486a000p-1
都打印出四舍五入的结果。然而,由于要求的精度足够,转换也做得很好,四舍五入的输出结果并没有比预期的更差。但我仍然觉得使用输出"%a"
更有参考价值。
这种细微的差别可能会在紧密循环的集成过程中累积起来。
然而,OP 发现的细微差别是实现中的 1 位差异,而不是sin(x)
与数学函数sine(x)之间的数学差异。
实施差异sin(x)
可以凸显一致性问题 – OP 发现了这一点。然而,OP 应该关注的是数学差异– 而不是解决。这需要更大规模的分析和大量代码才能评估。
如果sin()
实现差异产生了实际结果问题(除了精确的一致性),对此我表示怀疑(需要查看更大的应用程序),那么 OP 的代码可能需要更高精度的类型。
我强烈怀疑任何一个结果都OP 的需求。
2
-
它似乎取决于范围缩减的具体方式。我可以通过直接计算 sin(x) 或使用 pi 的标准双精度表示计算 sin(pi()-x) 在 Excel 中获取任一值。大多数系统数学库的范围缩减效果都比这更好。正如您所说,真实值相当接近两者的中间值。
– -
@MartinBrown 下次使用 π 而不是 pi() 😉。这可能是一个减少的问题,但我还怀疑减少的计算存在差异
x
。IAC,它只是一个接近一半情况(+/- 0.04 ULP)的数据点。
–
|
–
double
有 53 位尾数,相当于大约 15 位有效十进制数字。超过这个数字的一切都毫无意义。如果您使用的数值算法非常的差异就会导致其发散,那么您可能需要找到更好的算法。–
double
可以显示超过 15 位十进制数字,但精度最高只能达到 15 位;超出此范围的任何数字都不可信赖。由于您的输入d
超过 15 位,因此不能保证sin
在两个平台上都收到完全相同的输入 + 结果的前 15 位数字相同,这基本上就是您可以进行有意义比较的全部内容。–
d
以确保不同的编译器生成相同的值d
?–
printf("%d\n", FLT_EVAL_METHOD);
针对每个平台进行报告以更深入地了解该问题。–
|