[libc-commits] [PATCH] D115828: [libc] Implement correctly rounded log2f based on RLIBM library.

[Tue Ly via Phabricator via libc-commits]

lntue added a comment.

In D115828#3199419 <https://reviews.llvm.org/D115828#3199419>, @zimmermann6 wrote:

> the new version applies cleanly to the main branch. I have tested it on x86_64 under Linux (haswell). I confirm it is CR for rounding to nearest, and I get 3 failures if I disable the 3 exceptional cases. For other rounding modes I get 8 failures for rounding towards zero (with the exceptional cases), 8 failures too for rounding towards -Inf, and 7 failures for rounding towards +Inf.
>
> I tried with a polynomial generated by Sollya and with this polynomial we need no exceptional cases, and the routine is CR for all rounding modes (please can someone confirm?):
>
>   double r = __llvm_libc::fputil::polyeval(
>       d, extra_factor, 0x1.71547652b7fefp+0, -0x1.715476500a42ep-1,
>       0x1.ec70917f77152p-2, -0x1.71482b204ea69p-2, 0x1.21da0eb07c659p-2);
>
> For the record this polynomial was obtained with the following input file (then run sollya log2f.sollya):
>
>   n = 5; /* polynomial degree */
>   P = 53; /* precision of the coefficients */
>   
>   pretty = proc(u) {
>     return ~(floor(u*1000)/1000);
>   };
>   
>   d = [0, 1/2^7];
>   f = log2(1+x);
>   w = 1;
>   p =  remez(f, n, d, w);
>   pf = fpminimax(log2(1+x), [|1,2,3,4,5|], [|P...|], d, absolute, floating, 0, p)\
>   ;
>   err_p = -log2(dirtyinfnorm(pf*w-f, d));
>   print (pf, pretty(err_p));
>
> Sollya is available from https://www.sollya.org/. Would you consider using the Sollya polynomial?



In D115828#3200306 <https://reviews.llvm.org/D115828#3200306>, @santoshn wrote:

> RLIBM's polynomial should also have 0 special case inputs and produce correctly rounded results for all inputs and for all rounding modes for log2f.
>
> I think the current polynomial  has special case inputs because the generated polynomial is generated with a different polynomial evaluation.
>
> Here is a revised polynomial that uses the current polynomial evaluation and produces correct results for all rounding modes and all inputs with zero violated inputs:
>
> Polynomial: y=-3.2945312494298684154217536821552091564491707891340621650044795387657359e-16 x^(0) + 1.4426950408890866217603843324468471109867095947265625000000000000000000e+00 x^(1) + -7.2134752022691861483849606884177774190902709960937500000000000000000000e-01 x^(2) + 4.8089833027421252653610395100258756428956985473632812500000000000000000e-01 x^(3) + -3.6069225263970772221711058591608889400959014892578125000000000000000000e-01 x^(4) + 2.8949201646411226729327381690382026135921478271484375000000000000000000e-01 x^(5)

Thanks Paul and Santosh for looking into this!  It looks like there are many polynomials that will make it correctly rounded for the entire range, which is a great news!

That also make me think that so actually the exception cases from logf (and maybe log10f in the near future) are coming from adding the extra factor: m*log(2) + log(f).  So if we can make that addition more accurate, we shouldn't have any exceptional cases.  Moreover, it's entirely possible that we can reduce the degree of our polynomials while maintaining the accuracy.

I've discussed with Christoph and see if he can get any update on those directions.


Repository:
  rG LLVM Github Monorepo

CHANGES SINCE LAST ACTION
  https://reviews.llvm.org/D115828/new/

https://reviews.llvm.org/D115828