[VA] "Introducing APRIL !" microchallenge

[Albert Chan]

(04-07-2022 09:09 AM)Werner Wrote:  

(04-07-2022 08:11 AM)J-F Garnier Wrote:  Yes, and HP also changed the formula to use 90 and [PI/2] instead of 45 and [PI/4]. Why?
...
Dividing by 90 may be less prone to error amplification than 45.

No, dividing by 90 and 45 has the same relative error ...

At first, I agree with Werner.
After a while, I switched side, and agree with J-F Garnier.

ULP error analysis is hard, because ULP is on a variable scale, depends on actual number.

It is even harder when number is close to boundary base^n, with different size ULP.
see The Shortest Decimal String That Round-Trips May Not Be The Nearest

We can assume truncation error, in terms of ULP, for both /90 or /45 about the same.
However, relative error depends on size of mantissa (range from 0.100 to 0.999...)

Say, we only consider mantissa ranged of 0.9 to 1
For x/90, it is 10% (mantissa(x) = 0.81 to 0.90)
For x/45, it is 5% (mantissa(x) = 0.405 to 0.45)

For the same ULP error, the bigger the mantissa, the smaller its relative error.
In terms of relative errors, we would expect x/90 twice as accurate, compared to x/45
(might not be twice as accurate, due to Benford's Law, leading digit is likely to be small)

Below count bad rad(x), with 10 million random n-digits numbers

Intermediate calculations chopped to 13 digits.
Final result rounded to 10 digits, then compared with true rad(x)
pi2 = 1.570796326795
pi4 = .7853981633975 (wrongly rounded pi/4 helps a little here)

Code:

        x*pi2/90    x*pi4/45    x/90*pi2    x/45*pi4
n=9     2500        2500        3240        6277
n=10    2175        2176        4053        5134

The winner is to do division last.