Precision required in accuracy assessment

[blapie]

Hello,

There has been a number of inaccuracies reported in issues. Some of them reporting errors slightly above 1ULP on *_u10 routines. In #570 we seem to be in a case where MPFR calculates an error of 1 + eps where eps < unit roundoff, naturally eps gets rounded out once we display the final error in double precision.

I can understand that we wish to represent errors in the working precision, in this case double precision, as we have to set a limit to our final representation of errors. It seems like an acceptable compromise, although some user might expect ULP errors to be strictly inferior to threshold, in this case error < 1ULP.

What I'm wondering though is whether we actually need MPFR computation to be that precise.
Other libraries like Arm Optimized Routines only use an MPFR extended precision of 96 bits, which does not seem to limit accuracy assessment. I can dig to see what other libraries use.

So I have these 2/3 questions:

• What is the theoretical motivation for setting the extended precision to max exponent + number of bits + small offset, e.g. ~1023+53 in double precision ?
• Can we reduce that at least a little bit - say 256bits. This would greatly accelerate testing.
• Are there corner cases that definitely require the current formula? Can this be done on a limited number of routines/accuracies?

[shibatch]

In current testers, the accuracy of MPFR is already set close to the minimum required.

The reason why testing functions such as sinpi requires accuracy of more than 1000 bits is because that many digits of PI are needed to process the reduction of argument. However, since the latest version of MPFR includes sinpi and other functions, it is not necessary to specify a precision of more than 1000 bits to test those functions. Since only a few test items are calculated with high precision, modifying the code for those functions will not reduce the test time that much.

[blapie]

Thanks for the insight. Sounds like I missed a few of these aspects. That answers the last 2 questions.
As for my first one, where does the formula (max exponent + number of bits + small offset) comes from?

[shibatch]

That's basically the number of bits required to represent DBL_MAX * PI.

[shibatch]

To explain it a little better, for the calculation of sin function, the argument must first be reduced to a range between -PI/2 and PI/2.
To do so, we need to calculate remainder(x, PI/2) for argument x.
Since x reaches up to DBL_MAX, we need PI to a number of digits that can represent DBL_MAX to the decimal point + 53 bits for this computation.