Re: [Maxima-discuss] numerical issues with bfloat

[Michel T. <ta...@lp...>]

Le 01/10/2025 à 10:12, Igor Pesando a écrit :
>
> Now if you try to do the same computation with R using the same mpfr
> library (I don't know whether the same version) you get ....

Consider, following Stavros idea, the following:

Maxima 5.48.1 https://maxima.sourceforge.io
using Lisp SBCL 2.5.8

(%i1) 
qq:((42328762422919168*3^(19/2)-1020404800499791299075*sqrt(2))*%pi)/(5*2^(19/2))
  +((114866054188452371105939465*sqrt(3)-249900741459*2^(99/2))*%pi)/(5242880*3^(9/2))
  +((44912031264987141343855529*sqrt(3)-781679147445*2^(93/2))*%pi)/(13107200*3^(3/2))
  +((1439185410761439252972725*sqrt(3)-1603107148923*2^(81/2))*%pi)/(1638400*sqrt(3))
  -((101083151603075712162029*sqrt(3)-900771348897*2^(75/2))*%pi)/(81920*3^(3/2))-
  (397992477574879*%pi)/16511297126400$

(%i2) qr:radcan(qq);
                                      31/2
(%o2) - ((1391611950183735955030016 3
51/2
                 - 24455840869029408432698763712577 sqrt(2)) 
%pi)/(492075 2    )
(%i3) qf:factor(qq);
(%o3) ((24455840869029408432698763712577 sqrt(2) sqrt(3)
51/2  19/2
                       - 59904331359825180393845645377536) %pi)/(25 
2     3    )
(%i4) q0: -inpart(qf,4,2) /inpart(qf,4,1)$

(%i5) fpprec:16;
(%o5)                                 16
(%i6) bfloat(qq);
(%o6)                         4.188123926016877b3
(%i7) bfloat(qr);
(%o7)                        - 2.044971377673459b3
(%i8) bfloat(q0);
(%o8)                                1.0b0
(%i9) fpprec:32;
(%o9)                                 32
(%i10) bfloat(qq);
(%o10)               1.0940966878634238343622068201596b-5
(%i11) bfloat(qr);
(%o11)               1.0940967459989769199596357848314b-5
(%i12) bfloat(q0);
(%o12)                1.0000000000000000000000023512804b0
(%i13) fpprec:64;
(%o13)                                64
(%i14) bfloat(qq);
(%o14) 1.094096767106645616816035040876573195641179441484968568346303529b-5
(%i15) bfloat(qr);
(%o15) 1.094096767106645616816035040876573195640339186572047842998702713b-5
(%i16) bfloat(q0);
(%o16) 1.000000000000000000000002351280569685031554276286551780025082575b0
(%i17)

You see that qr and qf are much simpler than qq, computing with 16 
digits of precision gives totally wrong results for qq, and qr. With 32 
digits the result is basically correct , you see that qq and qr are 
close, computing with 64 bits confirms that. Note that the 32 digits 
computation provides only 7 correct digits, by comparison with 64.

On the other hand, the results computed by R are obviously wrong. 
Conclusion, qq is hard to compute and in such case one needs to increase 
fpprec, until the results stabilize. Probably the knob you have used 
with R is incorrect.



-- 
Michel Talon