Re: [PATCH] aarch64: Improve SVE sin polynomial

[Szabolcs Nagy <Szabolcs.Nagy@arm.com>]

The 08/04/2023 09:28, Paul Zimmermann wrote:
>        Hi Szabolcs,
> 
> > i'd say we want to minimize the maximum of
> > 
> >  (poly(x) - sin(x))/ulp(sin(x))
> > 
> > instead of
> > 
> >  (poly(x) - sin(x))/sin(x)
> 
> interesting question! I have asked the Sollya developers if this is possible.
> 
> Anyway, unless I'm wrong I find the largest ulp error with the polynomial
> from the patch is about 1.72 at x=0.252552902718105, and for the polynomial
> I obtained with Sollya it is about 1.43 at the same value of x.
> 
> See attached graph (blue = polynomial from the patch, red = polynomial
> I obtained with Sollya), which for 0 <= x <= pi/2 prints the absolute
> value of the ulp error as you define above.

neither the red nor the blue line has flat maximums, so they are not
optimal for |poly(x)-sin(x)|/ulp(sin(x)). but that only considers
approx error, not rounding.

as you can see the blue error goes down near pi/2, while it's high
around small values like 0.25. this is deliberate as at small values
there is tiny rounding error, while in [1,pi/2] there is cancellation
with huge errors (the order of x + c*x*x*x is smaller than x so the
terms are big an cancel each other after relatively big rounding
error).

so the actual worst case error after rounding is much bigger than
what sollya reports (approx error) and it is in [1,pi/2]. e.g.

_ZGVnN2v_sin(0x1.7a83373b951abp+0) got 0x1.fdd2e6d528af4p-1 want 0x1.fdd2e6d528af7p-1 - 0.234485 ulp

error: -2.7655 ulp

which shows that there is a point where rounding error takes over
and becomes dominant over approx error. the weighting function can
take this into account by roughly following ulp(sin(x)) for a while
and then in [1,pi/2] making it smaller (this is the weighting i
used, but it was put together manually, not computed precisely:
ideally we would compute rounding_error_bound(x) function based
on the poly eval and incorporate that into the weighting. for most
functions i don't expect this to help, but it matters for sin(x).)

i would expect simply minimizing the approx error (using the
ulp(sin(x)) weighting) would be better than your polynomial
(sin(x) weighting) after rounding but using our polynomial
(manually adjusted weighting) is even better (i verified it
experimentally that measured error near the approx error peaks
are roughly the same so the error is evenly spread out, but
i'm sure slightly better polynomial is possible to find with
a bit more research into the optimal design).