lib(ic)的精确解决方案
使用ECLiPSe Prolog的lib(ic)
我偶然发现了David H. Bailey的以下问题,“解决科学计算中的数字异常”。 这是我在Unum书中提到的。 实际上,它只是其中的一部分。 首先,让我用(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页的书中引用的意思,我解释为这些数字是完全表示的(现在通过抚摸)。
使用可以模拟IEEE单精度的{3,5}环境。 输入值可以精确表示。 ...
{-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)
有几个问题在这里引起混淆:
除了声明之外,示例中的三个常量没有双精度浮点数的精确表示。
最初的例子不涉及四舍五入。
第一个例子看似正确的结果实际上是由于一个幸运的四舍五入错误。 其他计算顺序给出不同的结果。
给定最接近常数的双浮点表示的确切结果实际上不是零,而是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)
现在两个结果都包含零,结果的精确度取决于评估顺序。
链接地址: http://www.djcxy.com/p/85799.html