The calculation of pi was accurate to 30 significant figures, and the calculation of k/pi was accurate to 30 significant figures, given our value of pi.
So far has performed as promised.
The computed value of k/pi is correct to 29 significant figures, 23 before the decimal place and 6 after.
So when we take the fractional part, we only have six significant figures, and that’s not enough.
That’s where things go wrong.
We get a value ⌊k/π⌋ that’s greater than 0.
5 in the seventh decimal place while the exact value is less than 0.
5 in the 25th decimal place.
We needed 25 – 6 = 19 more significant figures.
This is the core difficulty with floating point calculation: subtracting nearly equal numbers loses precision.
If two numbers agree to m decimal places, you can expect to lose m significant figures in the subtraction.
The error in the subtraction will be small relative to the inputs, but not small relative to the result.
Notice that we computed k/pi to 29 significant figures, and given that output we computed the fractional part exactly.
We accurately computed ⌊k/π⌋π, but lost precision when we computed we subtracted that value from k.
Our error in computing k – ⌊k/π⌋π was small relative to k, but not relative to the final result.
Our k is on the order of 1025, and the error in our subtraction is on the order of 10-7, but the result is on the order of 1.
There’s no bug in bc.
It carried out every calculation to the advertised accuracy, but it didn’t have enough working precision to produce the result we needed.
Related posts Quadratic formula and precision Computing the area of a polygon Math library functions that seem unnecessary.