Re: [PATCH] aarch64: Improve SVE sin polynomial

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

The 08/03/2023 15:21, Wilco Dijkstra wrote:
> Hi Paul,
> 
> > you can still improve the polynomial with the following one generated by
> > Sollya (https://www.sollya.org/) with relative error < 2^-52.454 on
> > [-pi/2, pi/2], instead of 2^-51.765 for the one you propose:
> 
> So the goal is to reduce the ULP error. Due to the way floating point works,
> errors are at their worst when the result is just above a power of 2. For the
> wide sin(x) polynomial it happens when sin(x) is close to but larger than 0.5.
> The poly is extremely accurate for small x. It's only problematic when x is large
> and we cross an exponent boundary.
> 
> Szabolcs can explain it better but what we ended up doing is to scale the
> sollya function so that errors at the end of the range are and reduced more.

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)

since ulp(x) is flat between 0.5 and 1 this can make a big
difference when designing a sin poly for -pi/2,pi/2.

e.g. you want to look at the abs error on [pi/6,pi/2], not
relative error when comparing the polynomials.

this is if you only consider approximation errors, but
there is an independent effect of large rounding errors
where operations in x + c*x*x*x cross powof2 boundaries
which is also near x = pi/2.

so for weighting func i used a polynomial that roughly
approximates ulp(sin(x)) but forces smaller approx error
where i expect big rounding error (this part was a bit
of black magic but i think the approx error part is
sound and covers most of the quality gain).

> 
> Note also this is not the final polynomial. The ULP is a bit high at 3.3. I had a
> more accurate one close to 2ULP but it could give results >1.0. So we need to
> teach sollya to estimate the end from below, and not allow it to overshoot.

we originally planned to use a poly with one more term,
but that made errors go in the wrong direction at x=pi/2,
and it could return 1.0 + 0x1p-52 which is not good for
sin(x). i could not fix this with weighting, but reducing
the degree the error goes into the other direction so
even though the error is bigger then it's probably better
in practice.

hope this helps.