Re: correct rounding or not?

[Carlos O'Donell <codonell at redhat dot com>]

On 1/29/20 11:36 AM, paul zimmermann wrote:
>        Hi,
> 
> in sysdeps/ieee754/dbl-64 one can find code coming from IBM's
> Accurate Mathematical Library (libultim), with comments claiming
> for correct rounding, for example in s_sin.c one can read around
> lines 194-197:
> 
> /*******************************************************************/
> /* An ultimate sin routine. Given an IEEE double machine number x   */
> /* it computes the correctly rounded (to nearest) value of sin(x)  */
> /*******************************************************************/
> ...
> __sin (double x)
> 
> However, it appears that the __sin function is not correctly rounded,
> for example for x=8.5546900000006210e-01 it gives 7.5487859961632653e-01
> whereas the correct rounding is 7.5487859961632664e-01 (according to mpfr).

Confirmed.

The correct infinite precision value is (computed by Wolfram Alpha as an 
independent check against mpfr):

0.7548785996163265 87068141643701359060493582844968198559261...

With the differences being:

5.7068141643701359060493582844968198559261 × 10^-17

-5.2931858356298640939506417155031801440739 × 10^-17

Thus 7.5487859961632664e-01 (0x1.827f72a39ad45p-1) is closer than...

The result of 7.5487859961632653e-01 (0x1.827f72a39ad44p-1).

Note the answers is off only by 1 ULP.

> Maybe the original code from libultim was modified, in particular the
> code that guarantees correct rounding (rounding test + slow path) was
> removed?
> 
> Anyway, I believe the comments should be modified to accurately
> describe what the current code does.

I agree with you and Szabolcs, the comment should be removed.

Do you have a copyright assignment in place with the FSF for glibc?

-- 
Cheers,
Carlos.