bfloat bug

[rtoy]

Imported from SourceForge on 2024-07-04 09:51:53 Created by nicolesharp100 on 2021-07-29 15:09:11 Original: https://sourceforge.net/p/maxima/bugs/3821

---

For the following calculation, I get different results depending on the value of fpprec for bfloat:

(%i112) bfloat(-(214826787042230210793264379160871222964571080869456289521*(8*%pi+3)^2*((394408937500000000000000000000*%i*(-8*%pi-3)*(4*%pi-3)*((8*%pi)/(4*%pi-3)+3/(4*%pi-3)))/(313860660122040131213690411963141883*(8*%pi+3)^2)+(394408937500000000000000000000*%i*(4*%pi-3)*(-(8*%pi)/(4*%pi-3)-3/(4*%pi-3)))/(313860660122040131213690411963141883*(-8*%pi-3))))/(30489448356056265625000000000000000000000000000000000000000000*(4*%pi-3)^2*((8*%pi)/(4*%pi-3)+3/(4*%pi-3)))),fpprec:32;

(%o112) 0.0b0

(%i109) bfloat(-(214826787042230210793264379160871222964571080869456289521*(8*%pi+3)^2*((394408937500000000000000000000*%i*(-8*%pi-3)*(4*%pi-3)*((8*%pi)/(4*%pi-3)+3/(4*%pi-3)))/(313860660122040131213690411963141883*(8*%pi+3)^2)+(394408937500000000000000000000*%i*(4*%pi-3)*(-(8*%pi)/(4*%pi-3)-3/(4*%pi-3)))/(313860660122040131213690411963141883*(-8*%pi-3))))/(30489448356056265625000000000000000000000000000000000000000000*(4*%pi-3)^2*((8*%pi)/(4*%pi-3)+3/(4*%pi-3)))),fpprec:64;

(%o109) 7.505588159636097600644900709241438622881609220057428301461265852b-76*%i

(%i110) bfloat(-(214826787042230210793264379160871222964571080869456289521*(8*%pi+3)^2*((394408937500000000000000000000*%i*(-8*%pi-3)*(4*%pi-3)*((8*%pi)/(4*%pi-3)+3/(4*%pi-3)))/(313860660122040131213690411963141883*(8*%pi+3)^2)+(394408937500000000000000000000*%i*(4*%pi-3)*(-(8*%pi)/(4*%pi-3)-3/(4*%pi-3)))/(313860660122040131213690411963141883*(-8*%pi-3))))/(30489448356056265625000000000000000000000000000000000000000000*(4*%pi-3)^2*((8*%pi)/(4*%pi-3)+3/(4*%pi-3)))),fpprec:128;

(%o110) 5.7015856673277740637677691905123095969144968122050015178699968321170963389954940908322433613327474741520172384511482893019560687b-140*%i

[rtoy]

Imported from SourceForge on 2024-07-04 09:51:54 Created by nicolesharp100 on 2021-07-30 18:13:45 Original: https://sourceforge.net/p/maxima/bugs/3821/#cf78

---

Still appears to be a bug, but I found a workaround by using

(%i*float(χμ) - %i*float(Z)^2*float(χε))/float(Z)^2;

instead of

float((%i*χμ - %i*Z^2*χε)/Z^2);

for the worksheet

α1: 137035999206*10^2*10^-11 $
i: %i $
lI: 1 $
mI: 1 $
np: 10^-12 $
qI: 1 $
tI: 1 $

α0: α1^-1 $
e: 1602176634*10^-19*10^-9*qI $
lm: lI $
vI: lI/tI $

c: 299792458*vI $
lpm: np*lm $
pI: mI*vI $

c2: c^2 $
LI: lI*pI $
r: 120*lpm $

h: 662607015*10^-34*10^-8*LI $
V: (4/3)*π*r^3 $

μ0: 2*α0*h/(c*e^2) $
Vα: V $

ε0: (c2*μ0)^-1 $
εr: -(8*π*Vα*V^-1 + 3)/(4*π*Vα*V^-1 - 3) $
ni: sqrt(-(8*π*Vα*V^-1/(4*π*Vα*V^-1 - 3)) - (3/(4*π*Vα*V^-1 - 3))) $

χε: εr*ε0 $
μr: ni^2/εr $

χμ: μr*μ0 $

Z: sqrt(χμ/χε) $

κω = (i*float(χμ) - i*float(Z)^2*float(χε))/float(Z)^2;

[rtoy]

Imported from SourceForge on 2024-07-04 09:51:57 Created by macrakis on 2021-08-04 23:32:20 Original: https://sourceforge.net/p/maxima/bugs/3821/#dd3c

---

Why do you think this is a bug? There are many floating-point calculations where different precisions will give different results. Can you identify a single arithmetic operation which is giving incorrectly rounded results?

[rtoy]

Imported from SourceForge on 2024-07-04 09:52:01 Created by nicolesharp100 on 2021-08-05 00:12:11 Original: https://sourceforge.net/p/maxima/bugs/3821/#dd3c/596d

---

I don't know if that is normal behavior or not, but bfloat gives wildly different results depending on fpprec, which is something I have not seen before: bfloat with fpprec 64 is about 8E-76*%i, bfloat with fpprec 128 is about 6E-140*%i, and bfloat with fpprec 256 is about 7E-268*%i, when it appears that the correct value might be approximately 4E-27*%i. I don't see how that can be a rounding error and it invalidates the ability to use bfloat to compute the result without using previous floats.

[rtoy]

Imported from SourceForge on 2024-07-04 09:52:05 Created by macrakis on 2021-08-05 00:30:01 Original: https://sourceforge.net/p/maxima/bugs/3821/#1a0a

---

No, it is not a rounding error. It is catastrophic cancellation.

Such errors do not "invalidate the ability to use bfloat to compute the result". They just mean that you need to do more careful numerical analysis. I don't know what you mean by "without using previous floats".

Calculating symbolically (apply ratsimp to your original expression), the exact result is actually zero, not 4e-27.

If you calculate the two subexpressions

(394408937500000000000000000000*%i*((-8*%pi)-3)*(4*%pi-3)*((8*%pi)/
  (4*%pi-3)+3/(4*%pi-3)))/(313860660122040131213690411963141883*(8*%pi+3)^2)

(394408937500000000000000000000*%i*(4*%pi-3)*((-(8*%pi)/(4*%pi-3))-3/
  (4*%pi-3)))/(313860660122040131213690411963141883*((-8*%pi)-3))

you will see that they evaluate to roughly plus and minus 1.25e-6. The calculatedsum of these two numbers at various values of fpprec is much smaller than 10^-fpprec.

[rtoy]

Imported from SourceForge on 2024-07-04 09:52:08 Created by nicolesharp100 on 2021-08-05 00:32:26 Original: https://sourceforge.net/p/maxima/bugs/3821/#bc62

---

Here is something even worse, I get a value opposite in sign when I use bfloat instead of float:

(%i330) (i*bfloat(χμ) - i*bfloat(Z)^2*bfloat(χε))/bfloat(Z)^2,fpprec:32;
(%o330) -6.0892565348022981111491567369641b-44*%i

(%i331) (i*bfloat(χμ) - i*bfloat(Z)^2*bfloat(χε))/bfloat(Z)^2,fpprec:64;
(%o331) 1.50111763192721952012898014184828772457632184401148566029225317b-75*%i

Not really sure what is going on there or if it is supposed to do that. If there is something wrong with the input, then a warning message to the user that the result is unreliable might be a good feature to add.

[rtoy]

Imported from SourceForge on 2024-07-04 09:52:12 Created by nicolesharp100 on 2021-08-05 00:35:58 Original: https://sourceforge.net/p/maxima/bugs/3821/#1a0a/ccce

---

Thank you for that explanation. If float or bfloat can detect these kinds of instances, I think that a warning message to the user would be very helpful.

[rtoy]

Imported from SourceForge on 2024-07-04 09:52:15 Created by macrakis on 2021-08-05 00:36:04 Original: https://sourceforge.net/p/maxima/bugs/3821/#bc62/09bd

---

Doing and interpreting floating-point calculations is surprisingly hard. My guess is that you haven't studied numerical analysis yet. To get a taste of the kinds of issues that come up, I recommend George Forsythe's very old, but still very valid, How do you solve a quadratic equation.

[rtoy]

Imported from SourceForge on 2024-07-04 09:52:19 Created by robert_dodier on 2021-08-05 00:36:19 Original: https://sourceforge.net/p/maxima/bugs/3821/#3fae

---

Hi Nicole, it's not clear to me what the bug is. In the original example, it appears that the argument for bfloat is actually identically equal to zero, and in fact ratsimp applied to that expression returns 0. However, bfloat (likewise float) evaluates expressions in just the form in which they are given to bfloat (or float); bfloat (and float) doesn't try to reorder expressions in any way. So the exact numerical result depends on the order of operations, and the result for floating point operations can therefore be nonzero due to roundoff error. The error is smaller (but still nonzero) for larger values of fpprec.

Since that's consistent with what's shown by the example, it looks to me like the behavior of Maxima here is to be expected, and it's not a bug. Is there something here that I've overlooked that is unexpected?

[rtoy]

Imported from SourceForge on 2024-07-04 09:52:22 Created by macrakis on 2021-08-05 00:39:08 Original: https://sourceforge.net/p/maxima/bugs/3821/#1a0a/ccce/9f4d

---

I'm not sure what you have in mind. How is the bfloat system supposed to know that two numbers that are slightly different "should be" exactly the same? Or that two numbers that are the same "should be" a little different?

There are techniques for automatically generating error bounds (interval arithmetic), but they have their own issues.