Re: correct rounding or not?

[paul zimmermann <Paul dot Zimmermann at inria dot fr>]

       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=