One of the major concerns I have with this code, and codes derived from this, has to do with specific numerical underflow that can be seen to be happening before the series is summed.

1st:

  let pi_major = 3.1415927f32;
  let pi_minor = -0.00000008742278f32;

[...]

  (x - pi_major - pi_minor) *
  (x + pi_major + pi_minor) * p1 * x

Focus on the x, pi_minor and pi_major.

This is fine, pi_major is basically pi. pi_minor is a small constant.

One knows from basic numerical analytics as applied to computing systems floating point [0] to never do additions/subtractions along the lines of x +/- c*epsilon, and expect to have anything close to an accurate result.

What you are really doing is adding/subtracting numbers of such widely different scales. This is a bad thing in general, and the only real way to handle this is to increase precision. Otherwise you completely swamp the small number with a larger one, and you lose all precision from your smaller number.

Simple example in Julia (0.5.1 or higher)

julia> deltax = Float16(0.1)

Float16(0.1)

julia> sum = Float16(0)

Float16(0.0)

julia> for i=1:1000000 sum+=deltax end

julia> sum

Float16(256.0)

Yes, that is right, according to this, 1/10th as a Float16 summed 1 million times reports as 256. Which is wildly ... ridiculously ... wrong.

Then redo this as Float32

julia> deltax2 = Float32(0.1)

0.1f0

julia> sum2 = Float32(0)

0.0f0

julia> for i=1:1000000

       sum2+=deltax2

       end

julia> sum2

100958.34f0

Much closer, but you can still see the impact from the catastrophic loss of precision [1]. This is a very well known phenomenon, that people seem to forget about every now and then, until it comes back to bite you.

I illustrated this for an HPC/parallel programming class I taught a decade ago [2] using some simple C code. The math isn't the problem. Its the order of accumulation of error.

This should be fairly sobering to those whom are happy running f16 codes ... you need to be very careful about the order of your operations (far more so than with f32, or f64), lest you suffer a catastrophic loss of precision, and wind up with a code that "computes faster", or as we old HPC types like to say about codes that make these sorts of mistakes "generates garbage at a tremendous rate".

During my class a decade ago, I had an astronomer whom coded in Fortran complain to me bitterly that he's always done sums this way, and there was no way that this concept was correct. The C code was written to demonstrate the effect to him.

It is sobering how much garbage numerical code there is out there. It gets even harder to distinguish good numerical code from bad numerical code when people get lost in almost irrelevant things, like coding style, or language of implementation.

If you really want to freak someone out, and you have a calculator app, see if you can take 1/3 + 1/3 + 1/3 - 1 (or 1/7 + ... + 1/7 - -1). If they use less than f32, the 1/3 trick will get you a few ULP from zero. If they use something below f64, the 1/7 and a number of others will show this.

Basically do a forward and reverse function application of a non-zero non-pathological argument, and subtract the original value from this. A compiler may optimize this away. An interpreter wont. And you'll see the issue of precision right there in the non-zero result.

Note as well, some apps will try to be "smart" and handle that case (sort of like optimizing benchmark codes with specialized handling). This is more troubling than doing the right thing, but I've always been able to find ways to expose the precision limits with simple tests like this.

The limited precision isn't a problem, its how you make use of it, and use, compare, add/subtract numbers with it that matters.

[0] https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.h...

[1] https://en.wikipedia.org/wiki/Loss_of_significance

[2] https://scalability.org/2007/11/why-am-i-surprised-that-peop...

[edited for code clarity]