Search code examples
precisionmaximawxmaxima

Round-off to zero behavior in Maxima float coefficients


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

Solution

  • 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