Round-off to zero behavior in Maxima float coefficients

Question

So I am working on a Maxima program that involves a bunch of iterations (the Souriau-Frame Drazin Inverse Algorithm, to be specific), each step of which yields a polynomial. I need to check and stop my iterations when the polynomial goes to zero (i.e., all coefficients go to zero).

Maxima seems to never truncate small numbers to zero up until it reaches something absurd like $10^{-323}$ and so on.

The following code snippet gives an idea of what I need:

(%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) 

As one can see, it's not truncating $10^{-300}$ to zero, it hangs (and I had to sigkill it) for $10^{-323}$ and everything smaller than $10^{-325}$ is set to zero.

I don't know where this 324 is coming from. And I'd like to know if it's possible to reduce this for my code.

Edit 1: Here's the output if I used rationalize instead of 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

Edit 2: Here's the output to "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

Answer 1

I gather that the goal is to replace small (in absolute value) float with zero. There doesn't appear to be a built-in function for that. Here's an attempt at an implementation via the pattern matching machinery.

First define a rule to replace small floats, and define a function which applies the rule to an expression.

(%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) $

Now create a random polynomial.

(%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

Apply squash_floats to the generated polynomial.

(%i10) squash_floats(e1);
                             6         3
(%o10)                  0.5 x  - 0.05 x

Change the squashing tolerance.

(%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

Verify the replacement happens in nested expressions.

(%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

Answer 2

First let's step back a moment. What is the behavior you are hoping to find? If you need to convert very small floats to rational numbers accurately, try rationalize instead of rat. Does that work correctly for 1e-323?

If you want floats smaller than a tolerance to be converted to zero, we'll need to take a different approach. I'll hold off on that for the moment.

About the specific behavior you have observed, it appears to be implementation-dependent; I get a different (still buggy) behavior with Maxima + SBCL, which reports a floating point overflow. What does build_info(); report?

I don't know if it matters, but 1e-323 is a so-called denormalized float -- it is smaller than the smallest normalized (full precision) float, which is about 1e-308.

Answer 3

First you say "I want to know when a polynomial is exactly going to zero." And then you say "if a coefficient in a polynomial drops below a threshold, I want that terms to be completely thrown out of the polynomial". So you don't want the polynomial to be exactly zero, you want it to be zero within some threshold (relative? absolute?).

I'm afraid I'm not familiar with the Souriau-Frame Drazin algorithm, but looking at the Greville paper about it, it seems that all the calculations are rational (no square roots etc.), so I wonder if it's feasible to perform your calculations with completely exact rational numbers instead of using floating-point numbers. Then presumably exact means exact, and you don't need to worry about thresholds at all.