#4106 sign bug involving large integers and %pi

[Robert Dodier]

From Stackoverflow: https://stackoverflow.com/questions/75536022/cas-maxima-problem-with-relational-operator-in-function-definition

Looks like the observed behavior is due to the following sign bug, as seen by the incorrect results in %o8, %o10, %o15, and %o16.

(%i2) fpprec:100$

(%i3) fpprintprec: 8 $

(%i4) xx: 11468322278445318 - 36028797018963968/%pi;
                                  36028797018963968
(%o4)         11468322278445318 - -----------------
                                         %pi
(%i5) yy: expand (%pi*xx);
(%o5)       11468322278445318 %pi - 36028797018963968
(%i6) zz: expand (2*yy);
(%o6)       22936644556890636 %pi - 72057594037927936
(%i7) bfloat (xx);
(%o7)                     7.0899905b-1
(%i8) is (xx > 1/2);
(%o8)                         false
(%i9) bfloat (yy);
(%o9)                      2.2273862b0
(%i10) is (yy > %pi/2); 
(%o10)                        false
(%i11) bfloat (zz);
(%o11)                     4.4547724b0
(%i12) is (zz > %pi);
(%o12)                        true
(%i13) [xx0, yy0, zz0]: expand ([xx - 1/2, yy - %pi/2, zz - %pi]);
        22936644556890635   36028797018963968
(%o13) [----------------- - -----------------, 
                2                  %pi
22936644556890635 %pi
--------------------- - 36028797018963968, 
          2
22936644556890635 %pi - 72057594037927936]
(%i14) bfloat ([xx0, yy0, zz0]);
(%o14)      [2.0899905b-1, 6.5658987b-1, 1.3131797b0]
(%i15) is (xx0 > 0);
(%o15)                        false
(%i16) is (yy0 > 0);
(%o16)                        false
(%i17) is (zz0 > 0);
(%o17)                        true
(%i18) float ([xx, yy, zz]);
(%o18)                   [0.0, 0.0, 0.0]
(%i19) float ([xx0, yy0, zz0]);
(%o19)                 [- 2.0, - 4.0, 0.0]

With fewer digits involved, the bug goes away, from what I can tell. But the numbers as shown are exact -- integers and %pi. There are heuristics in the sign code involving float or bigfloat evaluation; maybe that's where it is going wrong.

[Robert Dodier]

Some investigation shows SIGN1 calls NUMER, which calls $FLOAT on 22936644556890635/2 - 36028797018963968/%pi. This can't work right in general; I suspect it would be easy to make up further examples which fail.

The obvious thing to do is to call $BFLOAT with a large fpprec, but I wonder what other effects that's going to have.

[Stavros Macrakis]

Some parts of Maxima work hard to get correct results for cases like this.
For example, floor and ceiling (written by Barton) evaluate at three different bfloat precisions. His pretty-good-floor-or-ceiling could easily be extended for the sign case.

Fateman correctly points out that you don't want loop-termination conditions to be expensive to evaluate. But they won't be if the calculation is already a number. It's only when it's an expression with constants like %e and %pi or functions of exact numbers that you run into this problem. And if you write something perverse involving these things in your loop condition, your code deserves to run slowly!

[Robert Dodier]

Yes, probably something like the floor/ceiling code is needed -- maybe some parts can be repurposed for SIGN.

The expense of loop termination isn't relevant -- first of all, Maxima needs to get it right. Anyway the stuff about if x > 1/2 then ... isn't a loop termination test.

[Richard Fateman]

Just to be clear,
float (11468322278445318 - 36028797018963968/%pi); --> returns 0.0, and 0.0>1/2 is false.
How much work do you want Maxima to do, for example in testing for
inequalities for termination of a loop?
Here's some more info
(%i12) for i: 11468322278445315 thru
11468322278445320 do print([i , float (i - 36028797018963968/%pi)]);
[11468322278445315,-2.0]
[11468322278445316,-2.0]
[11468322278445317,-2.0]
[11468322278445318,0.0]
[11468322278445319,2.0]
[11468322278445320,2.0]

(%o12) done
...
with fpprintprec:4, using bfloat instead of float, the answers are
[11468322278445315,-2.291b0]
[11468322278445316,-1.291b0]
[11468322278445317,-2.91b-1]
[11468322278445318,7.09b-1]
[11468322278445319,1.709b0]
[11468322278445320,2.709b0]

[Robert Dodier]

The expressions in question are exact, so the intrusion of fixed precision floats is a surprise.

Maxima can't figure out everything, but this is a serious limitation. One can create large integers and do a lot of stuff with them, but Maxima falls down on on a random subproblem? Really? I don't think so.