lib(ic) 的精确解
Exact solutions for lib(ic)
使用 ECLiPSe Prolog 的 lib(ic) I stumbled upon the following problem from David H. Bailey, "Resolving numerical anomalies in scientific computation." which I was referred to by the 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 页中的引文，我将其解释为数字被准确表示（现在 划过 ）。
Use a {3,5}, environment, the one that can simulate IEEE single precision. The input values are exactly representable. ...
{-1, 2}
...
That did the job, computing the exact answer with fewer than half the bits used by ...
否则第 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)
这里有几个问题共同造成了混乱：
示例中的三个常量与声明的不同
not 是否具有双浮点数的精确表示。
初始示例不涉及四舍五入是不正确的。
第一个例子看似正确的结果实际上是由于
幸运舍入误差。其他计算顺序给出不同的结果。
给出最接近的双精度浮点表示的准确结果
常数，确实不是零而是 2.2204460492503131e-16.
区间算法只有在输入时才能给出准确的结果
是准确的，这里不是这种情况。常数必须是
扩大到包含所需小数部分的区间。
关系算法就像 lib(ic) 提供的那样天生
不保证特定的评估顺序。出于这个原因，四舍五入
错误可能与功能评估期间遇到的错误不同。
然而，对于给定的常数，结果将是准确的。
下面更详细一些。正如我将展示的一些
使用 ECLiPSe 查询的点，事先简单介绍一下语法：
两个浮点数用双下划线隔开，如0.99__1.01
表示具有下限和上限的区间常数，在这种情况下
1附近的数字。
两个下划线分隔的整数，例如3_4
用分子和分母表示一个有理常数，在这个
case 四分之三.
为了演示第 (1) 点，将 float 表示形式转换为
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)
两个结果现在都包含零，很明显
结果的精度取决于评估顺序。
使用 ECLiPSe Prolog 的 lib(ic) I stumbled upon the following problem from David H. Bailey, "Resolving numerical anomalies in scientific computation." which I was referred to by the 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 页中的引文，我将其解释为数字被准确表示（现在 划过 ）。
Use a {3,5}, environment, the one that can simulate IEEE single precision. The input values are exactly representable. ...
{-1, 2}
...
That did the job, computing the exact answer with fewer than half the bits used by ...
否则第 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)
这里有几个问题共同造成了混乱：
示例中的三个常量与声明的不同 not 是否具有双浮点数的精确表示。
初始示例不涉及四舍五入是不正确的。
第一个例子看似正确的结果实际上是由于 幸运舍入误差。其他计算顺序给出不同的结果。
给出最接近的双精度浮点表示的准确结果 常数，确实不是零而是 2.2204460492503131e-16.
区间算法只有在输入时才能给出准确的结果 是准确的，这里不是这种情况。常数必须是 扩大到包含所需小数部分的区间。
关系算法就像 lib(ic) 提供的那样天生 不保证特定的评估顺序。出于这个原因，四舍五入 错误可能与功能评估期间遇到的错误不同。 然而，对于给定的常数，结果将是准确的。
下面更详细一些。正如我将展示的一些 使用 ECLiPSe 查询的点，事先简单介绍一下语法：
两个浮点数用双下划线隔开，如
0.99__1.01表示具有下限和上限的区间常数，在这种情况下 1附近的数字。两个下划线分隔的整数，例如
3_4用分子和分母表示一个有理常数，在这个 case 四分之三.
为了演示第 (1) 点，将 float 表示形式转换为 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)
两个结果现在都包含零，很明显 结果的精度取决于评估顺序。