使用ECLiPSe Prolog的lib(ic),我偶然发现了来自David H. Bailey, "Resolving numerical anomalies in scientific computation."的以下问题,该问题由Unum book引用。实际上,它只是其中的一部分。首先,让我用(is)/2表示方程式。另外,请注意,所有这些十进制数字在基数2浮点数(包括IEEE)中都有确切的表示形式:

ECLiPSe Constraint Logic Programming System [kernel]
...
Version 6.2development #21 (x86_64_linux), Wed May 27 20:58 2015
[eclipse 1]: lib(ic).
...
Yes (0.36s cpu)
[eclipse 2]: X= -1, Y = 2, Null is 0.80143857*X+1.65707065*Y-2.51270273.

X = -1
Y = 2
Null = 0.0
Yes (0.00s cpu)


因此,这实际上是0.0(根本没有舍入)。但是现在用$=代替is一样:

[eclipse 3]: X= -1, Y = 2, Null $= 0.80143857*X+1.65707065*Y-2.51270273.

X = -1
Y = 2
Null = 2.2204460492503131e-16__2.2204460492503131e-16
Yes (0.00s cpu)


此间隔不包含0.0。我知道间隔算术通常有点过于近似:

[eclipse 4]: 1 $= sqrt(1).

Delayed goals:
    0 $= -1.1102230246251565e-16__2.2204460492503131e-16
Yes (0.00s cpu)


但是至少方程式成立!但是,在第一种情况下,不再包含零。显然我还不了解。我也尝试了eval/1,但无济于事。

[eclipse 5]: X= -1, Y = 2, Null $= eval(0.80143857*X+1.65707065*Y-2.51270273).

X = -1
Y = 2
Null = 2.2204460492503131e-16__2.2204460492503131e-16
Yes (0.00s cpu)


Null不包含0.0的原因是什么?



(在@jschimpf令人惊讶的答案之后编辑)

这是我在187页上的引文,我认为这是数字的准确表示(现在是逐句显示)。


使用{3,5}环境,该环境可以模拟IEEE单精度。输入值可精确表示。 ... {-1,2} ...就完成了这项工作,用少于...所用位的一半来计算确切的答案


否则,语句页面184包含:



...

0.80143857 x + 1.65707065 y = 2.51270273


这些方程式看起来足够纯真。假设精确的十进制输入,则通过x = -1和y = 2精确求解该系统。


这是使用SICStus的library(clpq)重新检查的:

| ?- {X= -1,Y=2,
     A = 80143857/100000000,
     B = 165707065/100000000,
     C = 251270273/100000000,
     Null = A*X+B*Y-C}.
X =  -1,
Y = 2,
A = 80143857/100000000,
B = 33141413/20000000,
C = 251270273/100000000,
Null = 0 ?
yes


因此-1,2是确切的解决方案。



精确的配方

这是一个在输入系数中不存在舍入问题的公式,解决方案仍然是-∞... +∞。因此很正确,但不可用。

[eclipse 2]: A = 25510582, B = 52746197, U = 79981812,
C = 80143857, D = 165707065, V = 251270273,
A*X+B*Y$=U,C*X+D*Y$=V.

A = 25510582
B = 52746197
U = 79981812
C = 80143857
D = 165707065
V = 251270273
X = X{-1.0Inf .. 1.0Inf}
Y = Y{-1.0Inf .. 1.0Inf}


Delayed goals:
    52746197 * Y{-1.0Inf .. 1.0Inf} + 25510582 * X{-1.0Inf .. 1.0Inf} $= 79981812
    80143857 * X{-1.0Inf .. 1.0Inf} + 165707065 * Y{-1.0Inf .. 1.0Inf} $= 251270273
Yes (0.00s cpu)

最佳答案

以下几个问题共同造成了混乱:


除了要求之外,示例中的三个常量
没有确切的表示形式为double float。
最初的示例不涉及舍入是不正确的。
第一个示例中看似正确的结果实际上是由于
幸运的舍入错误。其他计算顺序给出不同的结果。
给出最精确的结果,给出
常数确实不是零,而是2.2204460492503131e-16。
区间算术只能在输入时给出准确的结果
是准确的,这里不是这种情况。常数必须是
扩展为包含所需小数部分的间隔。
就像lib(ic)提供的关系算术一样,它本质上就是
不保证特定的评估顺序。因此,舍入
错误可能不同于功能评估期间遇到的错误。
但是,对于给定的常数,结果将是准确的。


以下将更详细地介绍。我将示范一些
使用ECLiPSe查询指向要点,以下是有关语法的简短提示:


由双下划线分隔的两个浮点,例如0.99__1.01
表示带有下限和上限的区间常数
1附近的数字
用单个下划线分隔的两个整数,例如3_4
用分子和分母表示一个有理常数
案例四分之三。


为了演示点(1),请转换以下内容的浮点表示形式:
0.80143857成理性。这给出了精确的分数
3609358445212343/4503599627370496,虽然接近但不完全相同,
到预期的十进制小数80143857/100000000。浮点数
因此,表示形式不准确:

?- F is rational(0.80143857), F =\= 80143857_100000000.
F = 3609358445212343_4503599627370496
Yes (0.00s cpu)


下面显示了结果如何取决于评估顺序
(以上第3点;请注意,我已通过以下方式简化了原始示例:
摆脱无关的乘法):

?- Null is -0.80143857 + 3.3141413 - 2.51270273.
Null = 0.0
Yes (0.00s cpu)

?- Null is -2.51270273 + 3.3141413 - 0.80143857.
Null = 2.2204460492503131e-16
Yes (0.00s cpu)


顺序相关性证明发生舍入误差(点2)。对于那些熟悉浮点运算的人来说,很容易看到
添加-0.80143857 + 3.3141413时,0.80143857的精度为两位
在调整操作数的指数时迷路。实际上是
这个幸运的舍入错误使OP看起来看似正确!

实际上,第二个结果相对于
常数的浮点表示。我们可以证明这一点
通过使用精确的有理算法重复计算:

?- Null is rational(-0.80143857) + rational(3.3141413) - rational(2.51270273).
Null = 1_4503599627370496
Yes (0.00s cpu)

?- Null is rational(-2.51270273) + rational(3.3141413) - rational(0.80143857).
Null = 1_4503599627370496
Yes (0.00s cpu)


由于加法是使用精确的理性完成的,因此现在的结果是
与订单无关,并且由于1_4503599627370496 =:= 2.2204460492503131e-16
这确认了上面获得的非零浮点结果(点4)。

区间算术如何在这里提供帮助?通过与
包含真实值的时间间隔,这样结果将始终
输入方面要准确。因此,拥有
包含的输入间隔(ECLiPSe术语中的有界实数)
所需的真实值。可以通过编写它们来获得
显式向下,例如0.80143856__0.80143858;
通过转换精确数字,例如合理使用
breal(80143857_100000000);或指示解析器自动
将所有浮点数扩展为有界实数区间,如下所示:

?- set_flag(syntax_option, read_floats_as_breals).
Yes (0.00s cpu)

?- Null is -0.80143857 + 3.3141413 - 2.51270273.
Null = -8.8817841970012523e-16__1.3322676295501878e-15
Yes (0.00s cpu)

?- Null is -2.51270273 + 3.3141413 - 0.80143857.
Null = -7.7715611723760958e-16__1.2212453270876722e-15
Yes (0.00s cpu)


这两个结果现在都包含零,并且很明显
结果的精度取决于评估顺序。

08-16 06:50