Re: correct rounding or not?

Re: correct rounding or not?

[Wilco Dijkstra]

Hi Paul,

/*******************************************************************/
/* An ultimate sin routine. Given an IEEE double machine number x   */
/* it computes the correctly rounded (to nearest) value of sin(x)  */
/*******************************************************************/

> Anyway, I believe the comments should be modified to accurately
> describe what the current code does.

Yes that comment is out  of date - I have posted a patch to fix that [1].
This math code is ancient, extremely slow and buggy, so the plan is to
rewrite like we've done with exp(2)(f)/pow(f)/log(f)/sinf/cosf/sincosf.

There is ongoing work to remove the slow paths first [2], since that causes
the huge 100x slowdowns that most people are complaining about.
Rewriting the actual math code then gives another 2-5x gain.

Cheers,
Wilco

[1] https://www.sourceware.org/ml/libc-alpha/2020-01/msg00647.html
[2] https://www.sourceware.org/ml/libc-alpha/2020-01/msg00574.html

Re: correct rounding or not?

[paul zimmermann]

       Dear Wilco,

for [1] I don't know what "it computes *the* rounded value of sin(x)" means.
The header of the file says that the error is ~0.55ulp, that would be clearer.
Maybe the original libultim source contained an error bound for the fast path,
but unfortunately I do not have it any more.

However, in revision e4d8276, I see the largest 'corr' factor is 1.07,
which corresponds to a relative error of 0.07/2^55, thus to a maximum
error of 0.5175 ulps. But maybe the code was changed since then.

Note also that we can still find "correctly rounded" in e_atan2.c, s_atan.c,
and s_tan.c.

For [2] I fully support this patch. I checked independently with mpcheck
(with glibc 2.29, and with the current development version with your patch)
and found no difference (I only give the output for the development version
since both are identical):

zimmerma@tomate:/localdisk/zimmerma/glibc/build$ ./testrun.sh /tmp/mpcheck-double --num=14400000 --seed=1
/tmp/mpcheck-double --num=14400000 --seed=1
GCC: 9.2.1 20200123
GNU libc version: 2.30.9000
GNU libc release: development
MPFR 4.0.2
[precision=53, seed=1, emin=-1021, emax=1024]
Testing function asin for precision 53, exponent 0 [seed=1]
Max. errors for asin [exp. 0]: 0 (rndn), 5.59e-3 (rndz), 1.00 (rndu), 1.00 (rndd)
Testing function asin for precision 53, exponent -10 [seed=1]
Max. errors for asin [exp. -10]: 0 (rndn), 0 (rndz), 1.84e-5 (rndu), 2.04e-5 (rndd)
Testing function acos for precision 53, exponent 0 [seed=1]
Max. errors for acos [exp. 0]: 0 (rndn), 2.91e-2 (rndz), 2.76e-1 (rndu), 1.38e-2 (rndd)
Testing function acos for precision 53, exponent -10 [seed=1]
Max. errors for acos [exp. -10]: 0 (rndn), 0 (rndz), 2.76e-1 (rndu), 0 (rndd)
Max. errors: 0 (rndn), 2.91e-2 (rndz), 1.00 (rndu), 1.00 (rndd) [seed=1]
Incorrect roundings: 14487298 (basic 0)
Wrong side of directed rounding: 14487298

Also, "make bench" on my machine gives the following results:

    function & fast path & slow path \\ \hline
    acos     & 168 & 160894 \\ % 235 after patch [2]
    asin     & 160 & 163083 \\ % 216 after patch [2]
    atan     & 127 & 29389 \\ % 144 bits
    sin      & 92 & 198 \\ % 768 bits
    cos      & 115 & 211 \\ % 768 bits
    tan      & 189 & 172531 \\ % 768 bits
    exp      & 25  & 26 \\ % 144 bits
    log      & 20  &    \\
    pow      & 35 & 67 \\ % 240 bits

Note that atan and tan still have slow "slow paths".

Also I wonder, now that the "slow path" for asin takes only 216 cycles,
whereas the fast path takes 160 cycles, what would you obtain if you use
the slow path only?

Best regards,
Paul

> From: Wilco Dijkstra <Wilco.Dijkstra@arm.com>
> Date: Thu, 30 Jan 2020 13:05:26 +0000
> 
> Hi Paul,
> 
> /*******************************************************************/
> /* An ultimate sin routine. Given an IEEE double machine number x   */
> /* it computes the correctly rounded (to nearest) value of sin(x)  */
> /*******************************************************************/
> 
> > Anyway, I believe the comments should be modified to accurately
> > describe what the current code does.
> 
> Yes that comment is out  of date - I have posted a patch to fix that [1].
> This math code is ancient, extremely slow and buggy, so the plan is to
> rewrite like we've done with exp(2)(f)/pow(f)/log(f)/sinf/cosf/sincosf.
> 
> There is ongoing work to remove the slow paths first [2], since that causes
> the huge 100x slowdowns that most people are complaining about.
> Rewriting the actual math code then gives another 2-5x gain.
> 
> Cheers,
> Wilco
> 
> [1] https://www.sourceware.org/ml/libc-alpha/2020-01/msg00647.html
> [2] https://www.sourceware.org/ml/libc-alpha/2020-01/msg00574.html=

Re: correct rounding or not?

[Wilco Dijkstra]

Hi Paul,

> for [1] I don't know what "it computes *the* rounded value of sin(x)" means.
> The header of the file says that the error is ~0.55ulp, that would be clearer.
> Maybe the original libultim source contained an error bound for the fast path,
> but unfortunately I do not have it any more.

Well we don't have good error bounds for sin/cos, the ULP estimates are based
on the original code with the slow paths.

> However, in revision e4d8276, I see the largest 'corr' factor is 1.07,
> which corresponds to a relative error of 0.07/2^55, thus to a maximum
> error of 0.5175 ulps. But maybe the code was changed since then.

__sin used to have this:

  else if (k < 0x3feb6000)
    {
      res = do_sin (x, 0, &cor);
      retval = (res == res + 1.096 * cor) ? res : slow1 (x);

From this I estimated the worst-case ULP as < 0.55. It could well be a bit
better in reality.

> Note also that we can still find "correctly rounded" in e_atan2.c, s_atan.c,
> and s_tan.c.

Yes, those still have slow paths which give correctly rounded results.

> /tmp/mpcheck-double --num=14400000 --seed=1

That looks like a useful tool - we did something similar for random testing
using long double math to test double functions.

> Also, "make bench" on my machine gives the following results:

    function & fast path & slow path \\ \hline

    acos     & 168 & 160894 \\ % 235 after patch [2]
    asin     & 160 & 163083 \\ % 216 after patch [2]
    atan     & 127 & 29389 \\ % 144 bits
    sin      & 92 & 198 \\ % 768 bits
    cos      & 115 & 211 \\ % 768 bits
    tan      & 189 & 172531 \\ % 768 bits
    exp      & 25  & 26 \\ % 144 bits
    log      & 20  &    \\
    pow      & 35 & 67 \\ % 240 bits

> Note that atan and tan still have slow "slow paths".

Yes the 1000x slowdown to try to round correctly is completely ridiculous!

It's also interesting to see that the new exp/log worst case is 5x faster than
sin/cos fast case...

> Also I wonder, now that the "slow path" for asin takes only 216 cycles,
> whereas the fast path takes 160 cycles, what would you obtain if you use
> the slow path only?

The fast path is likely already more than accurate enough. It's unlikely to be
worth it though - rewriting these functions from scratch is the way forward.

Cheers,
Wilco

Re: correct rounding or not?

[paul zimmermann]

       Dear Wilco,

> __sin used to have this:
> 
>   else if (k < 0x3feb6000)
>     {
>       res = do_sin (x, 0, &cor);
>       retval = (res == res + 1.096 * cor) ? res : slow1 (x);

this would correspond to a maximum error of 0.5+0.096/4 = 0.524 ulp.
This corresponds to what I see after about 2.3e9 random tests (5.23e-1).

> >From this I estimated the worst-case ULP as < 0.55. It could well be a bit
> better in reality.

Best regards,
Paul

Re: correct rounding or not?

[Carlos O'Donell]

On 1/30/20 4:13 AM, paul zimmermann wrote:
>        Dear Carlos,
> 
>> I agree with you and Szabolcs, the comment should be removed.
>>
>> Do you have a copyright assignment in place with the FSF for glibc?
> 
> not yet. I requested the forms on January 16, but did not receive them yet.
> I just asked assign@gnu.org again.
>
Paul,

Thanks! If this stalls too much please ping me and I can talk to the FSF
directly.

-- 
Cheers,
Carlos.

Re: correct rounding or not?

[paul zimmermann]

       Dear Carlos,

> I agree with you and Szabolcs, the comment should be removed.
> 
> Do you have a copyright assignment in place with the FSF for glibc?

not yet. I requested the forms on January 16, but did not receive them yet.
I just asked assign@gnu.org again.

Paul

Re: correct rounding or not?

[Carlos O'Donell]

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.

Re: correct rounding or not?

[Szabolcs Nagy]

On 29/01/2020 16:36, 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).
> 
> 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.

yes the comments were likely left untouched
when the slow paths were removed, i agree
that this should be fixed.

correct rounding or not?

[paul zimmermann]

       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).

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.

Best regards,
Paul Zimmermann