1

所以我正在研究一个涉及大量迭代的 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 它）$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
4

3 回答 3

2

我认为目标是用零替换小的（绝对值）浮点数。似乎没有内置功能。这是通过模式匹配机制实现的尝试。

首先定义一个规则来替换小浮点数，并定义一个将规则应用于表达式的函数。

(%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
于 2021-04-17T06:15:02.683 回答
0

首先让我们退后一步。您希望找到的行为是什么？如果您需要将非常小的浮点数准确地转换为有理数，请尝试rationalize使用rat. 这对 1e-323 是否正常工作？

如果您希望将小于容差的浮点数转换为零，我们需要采取不同的方法。我会暂时搁置。

关于您观察到的具体行为，它似乎取决于实现；我使用 Maxima + SBCL 得到不同的（仍然是错误的）行为，它报告浮点溢出。报告什么build_info();？

我不知道这是否重要，但 1e-323 是所谓的非规范化浮点数——它小于最小的规范化（全精度）浮点数，大约为 1e-308。

于 2021-04-16T03:58:36.627 回答
0

首先，您说“我想知道多项式何时准确归零。” 然后你说“如果多项式中的系数低于阈值，我希望将这些项完全排除在多项式之外”。因此，您不希望多项式恰好为零，而是希望它在某个阈值内为零（相对？绝对？）。

恐怕我对 Souriau-Frame Drazin 算法不熟悉，但是看一下关于它的 Greville 论文，似乎所有的计算都是有理的（没有平方根等），所以我想知道它是否可行您的计算使用完全精确的有理数，而不是使用浮点数。那么大概精确意味着精确，您根本不需要担心阈值。

于 2021-04-20T02:47:10.683 回答