Maxima 浮点系数中的四舍五入为零行为
Round-off to zero behavior in Maxima float coefficients
所以我正在研究一个涉及一系列迭代的 Maxima 程序（具体来说是 Souriau-Frame Drazin 逆算法），每一步都会产生一个多项式。当多项式变为零（即所有系数变为零）时，我需要检查并停止我的迭代。
Maxima 似乎永远不会将小数字截断为零，直到达到 $10^{-323}$ 等荒谬的数字。
下面的代码片段说明了我需要什么：
(%i3) rat(1e-300);
rat: replaced 1.0E-300 by 1/999999999999999903803069407426113968898218766118103141789833949572356552411722264192305659040010509526872994217248819197070144216063125530186267630296136203765329090687113225440746189048800695790727969805197112921161540803823920273299782054992133678869364753954248541633605124057805104488924519071744 = 1.0E-300
(%o3)/R/ 1/9999999999999999038030694074261139688982187661181031417898339495723\
565524117222641923056590400105095268729942172488191970701442160631255301862676\
302961362037653290906871132254407461890488006957907279698051971129211615408038\
23920273299782054992
133678869364753954248541633605124057805104488924519071744
(%i4) rat(1e-323);
rat: replaced 9.0E-324 by^C
Maxima encountered a Lisp error:
SIMPLE-ERROR: Console interrupt.
Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.
(%i5) rat(1e-325);
rat: replaced 0.0 by 0/1 = 0.0
(%o5)/R/ 0
(%i6)
如您所见，不是 将 $10^{-300}$ 截断为零，它 挂起 （我不得不sigkill it) $10^{-323}$ 并且小于 $10^{-325}$ 的所有内容都设置为零。
我不知道这个324是从哪里来的。我想知道是否可以为我的代码减少这个。
编辑 1：这是我使用 rationalize 而不是 rat 时的输出：
(%i3) rationalize(1e-300);
(%o3) 6032057205060441/6032057205060440848842124543157735677050252251748505781\
796615064961622344493727293370973578138265743708225425014400837164813540499979\
063179105919597766951022193355091707896034850684039059079180396788349106095584\
290087446076413771468940477241550670753145517602931224392424029547429993824129\
889235158145614364972941312
(%i4) rationalize(1e-323);
(%o4) 1/1012011266536553091762476733594586535247783248820710591784506790137151\
697839976734459801918507185622475935389321584059556949043686928967384335066999\
703692549607587121382831806822334538710466081706198838392363725342810037417123\
463493090516778245797781704050282561793847761667073076152512660931637543230031\
31653853870546747392
(%i5) rationalize(1e-324);
(%o5) 0
编辑 2：这是“build_info();”的输出：
(%i6) build_info();
(%o6)
Maxima version: "5.43.2"
Maxima build date: "2020-02-21 05:22:38"
Host type: "x86_64-pc-linux-gnu"
Lisp implementation type: "GNU Common Lisp (GCL)"
Lisp implementation version: "GCL 2.6.12"
User dir: "/home/nidish/.maxima"
Temp dir: "/tmp"
Object dir: "/home/nidish/.maxima/binary/5_43_2/gcl/GCL_2_6_12"
Frontend: false
首先让我们退后一步。您希望找到的行为是什么？如果您需要将非常小的浮点数准确地转换为有理数，请尝试 rationalize 而不是 rat。这对 1e-323 是否正常工作？
如果您希望将小于公差的浮点数转换为零，我们需要采用不同的方法。我暂时不说了。
关于您观察到的具体行为，它似乎与实现有关；我在 Maxima + SBCL 中得到了不同的（仍然有问题）行为，它报告浮点溢出。 build_info(); 报告了什么？
我不知道这是否重要，但 1e-323 是所谓的非规范化浮点数 -- 它小于最小的规范化（全精度）浮点数，大约为 1e-308。
我了解到目标是用零替换小的（绝对值）浮点数。似乎没有内置功能。这是通过模式匹配机制实现的尝试。
首先定义一个规则来替换小浮点数，然后定义一个将规则应用于表达式的函数。
(%i4) matchdeclare(xx,floatnump) $
(%i5) defrule(squashing_rule,xx, if abs(xx) <= squashing_tolerance then 0 else xx);
(%o5) squashing_rule : xx -> (if abs(xx) <= squashing_tolerance then 0 else xx)
(%i6) squashing_tolerance:0.01 $
(%i7) squash_floats(expr):=applyb1(expr,squashing_rule) $
现在创建一个随机多项式。
(%i8) e:makelist(float((((2*random(2)-1)*(1+random(8)))/8) *10^-random(4)) *x^k,k,1,6);
2 3 4 5 6
(%o8) [- 3.75e-4 x, - 0.00625 x , - 0.05 x , 0.00625 x , 0.005 x , 0.5 x ]
(%i9) e1:apply("+",e);
6 5 4 3 2
(%o9) 0.5 x + 0.005 x + 0.00625 x - 0.05 x - 0.00625 x - 3.75e-4 x
对生成的多项式应用squash_floats。
(%i10) squash_floats(e1);
6 3
(%o10) 0.5 x - 0.05 x
更改压缩容差。
(%i11) squashing_tolerance:0.001;
(%i12) squash_floats(e1);
6 5 4 3 2
(%o12) 0.5 x + 0.005 x + 0.00625 x - 0.05 x - 0.00625 x
验证替换发生在嵌套表达式中。
(%i13) squash_floats(sin(1+1/e1));
1
(%o13) sin(----------------------------------------------------- + 1)
6 5 4 3 2
0.5 x + 0.005 x + 0.00625 x - 0.05 x - 0.00625 x
首先你说“我想知道多项式何时正好为零。”然后你说“如果多项式中的系数下降到阈值以下，我希望这些项完全从多项式中剔除”。所以你不希望多项式正好为零，你希望它在某个阈值（相对？绝对？）内为零。
恐怕我对Souriau-Frame Drazin算法不熟悉，但是看了Greville的论文，似乎所有的计算都是有理的（没有平方根等），所以我想知道如果使用完全精确的有理数而不是使用浮点数来执行计算是可行的。那么想必exact就是exact，完全不用担心阈值。
所以我正在研究一个涉及一系列迭代的 Maxima 程序（具体来说是 Souriau-Frame Drazin 逆算法），每一步都会产生一个多项式。当多项式变为零（即所有系数变为零）时，我需要检查并停止我的迭代。
Maxima 似乎永远不会将小数字截断为零，直到达到 $10^{-323}$ 等荒谬的数字。
下面的代码片段说明了我需要什么：
(%i3) rat(1e-300);
rat: replaced 1.0E-300 by 1/999999999999999903803069407426113968898218766118103141789833949572356552411722264192305659040010509526872994217248819197070144216063125530186267630296136203765329090687113225440746189048800695790727969805197112921161540803823920273299782054992133678869364753954248541633605124057805104488924519071744 = 1.0E-300
(%o3)/R/ 1/9999999999999999038030694074261139688982187661181031417898339495723\
565524117222641923056590400105095268729942172488191970701442160631255301862676\
302961362037653290906871132254407461890488006957907279698051971129211615408038\
23920273299782054992
133678869364753954248541633605124057805104488924519071744
(%i4) rat(1e-323);
rat: replaced 9.0E-324 by^C
Maxima encountered a Lisp error:
SIMPLE-ERROR: Console interrupt.
Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.
(%i5) rat(1e-325);
rat: replaced 0.0 by 0/1 = 0.0
(%o5)/R/ 0
(%i6)
如您所见，不是 将 $10^{-300}$ 截断为零，它 挂起 （我不得不sigkill it) $10^{-323}$ 并且小于 $10^{-325}$ 的所有内容都设置为零。
我不知道这个324是从哪里来的。我想知道是否可以为我的代码减少这个。
编辑 1：这是我使用 rationalize 而不是 rat 时的输出：
(%i3) rationalize(1e-300);
(%o3) 6032057205060441/6032057205060440848842124543157735677050252251748505781\
796615064961622344493727293370973578138265743708225425014400837164813540499979\
063179105919597766951022193355091707896034850684039059079180396788349106095584\
290087446076413771468940477241550670753145517602931224392424029547429993824129\
889235158145614364972941312
(%i4) rationalize(1e-323);
(%o4) 1/1012011266536553091762476733594586535247783248820710591784506790137151\
697839976734459801918507185622475935389321584059556949043686928967384335066999\
703692549607587121382831806822334538710466081706198838392363725342810037417123\
463493090516778245797781704050282561793847761667073076152512660931637543230031\
31653853870546747392
(%i5) rationalize(1e-324);
(%o5) 0
编辑 2：这是“build_info();”的输出：
(%i6) build_info();
(%o6)
Maxima version: "5.43.2"
Maxima build date: "2020-02-21 05:22:38"
Host type: "x86_64-pc-linux-gnu"
Lisp implementation type: "GNU Common Lisp (GCL)"
Lisp implementation version: "GCL 2.6.12"
User dir: "/home/nidish/.maxima"
Temp dir: "/tmp"
Object dir: "/home/nidish/.maxima/binary/5_43_2/gcl/GCL_2_6_12"
Frontend: false
首先让我们退后一步。您希望找到的行为是什么？如果您需要将非常小的浮点数准确地转换为有理数，请尝试 rationalize 而不是 rat。这对 1e-323 是否正常工作？
如果您希望将小于公差的浮点数转换为零，我们需要采用不同的方法。我暂时不说了。
关于您观察到的具体行为，它似乎与实现有关；我在 Maxima + SBCL 中得到了不同的（仍然有问题）行为，它报告浮点溢出。 build_info(); 报告了什么？
我不知道这是否重要，但 1e-323 是所谓的非规范化浮点数 -- 它小于最小的规范化（全精度）浮点数，大约为 1e-308。
我了解到目标是用零替换小的（绝对值）浮点数。似乎没有内置功能。这是通过模式匹配机制实现的尝试。
首先定义一个规则来替换小浮点数，然后定义一个将规则应用于表达式的函数。
(%i4) matchdeclare(xx,floatnump) $
(%i5) defrule(squashing_rule,xx, if abs(xx) <= squashing_tolerance then 0 else xx);
(%o5) squashing_rule : xx -> (if abs(xx) <= squashing_tolerance then 0 else xx)
(%i6) squashing_tolerance:0.01 $
(%i7) squash_floats(expr):=applyb1(expr,squashing_rule) $
现在创建一个随机多项式。
(%i8) e:makelist(float((((2*random(2)-1)*(1+random(8)))/8) *10^-random(4)) *x^k,k,1,6);
2 3 4 5 6
(%o8) [- 3.75e-4 x, - 0.00625 x , - 0.05 x , 0.00625 x , 0.005 x , 0.5 x ]
(%i9) e1:apply("+",e);
6 5 4 3 2
(%o9) 0.5 x + 0.005 x + 0.00625 x - 0.05 x - 0.00625 x - 3.75e-4 x
对生成的多项式应用squash_floats。
(%i10) squash_floats(e1);
6 3
(%o10) 0.5 x - 0.05 x
更改压缩容差。
(%i11) squashing_tolerance:0.001;
(%i12) squash_floats(e1);
6 5 4 3 2
(%o12) 0.5 x + 0.005 x + 0.00625 x - 0.05 x - 0.00625 x
验证替换发生在嵌套表达式中。
(%i13) squash_floats(sin(1+1/e1));
1
(%o13) sin(----------------------------------------------------- + 1)
6 5 4 3 2
0.5 x + 0.005 x + 0.00625 x - 0.05 x - 0.00625 x
首先你说“我想知道多项式何时正好为零。”然后你说“如果多项式中的系数下降到阈值以下，我希望这些项完全从多项式中剔除”。所以你不希望多项式正好为零，你希望它在某个阈值（相对？绝对？）内为零。
恐怕我对Souriau-Frame Drazin算法不熟悉，但是看了Greville的论文，似乎所有的计算都是有理的（没有平方根等），所以我想知道如果使用完全精确的有理数而不是使用浮点数来执行计算是可行的。那么想必exact就是exact，完全不用担心阈值。