Re: [musl] asinf pulls in sqrt (double version)

[Norbert Lange <nolange79@gmail.com>]

Am Do., 10. März 2022 um 18:05 Uhr schrieb Rich Felker <dalias@libc.org>:
>
> On Thu, Mar 10, 2022 at 04:44:37PM +0100, Norbert Lange wrote:
> > Hello,
> >
> > pulled some of the math routines for a float-only project,
> > what I found is that asinf pulls in sqrt, and I cant tell whether this
> > is necessary because of precision.
> >
> > The line is in src/math/asinf.c:
> > s = sqrt(z);
>
> It looks like it's intentional -- s has type double and is used in an
> expression where the additional precision would make a difference:
>
>         x = pio2 - 2*(s+s*R(z));
>
> I suspect it would be more work (and slower) to get suitable precision
> without the double intermediate here, but I have not done the analysis
> myself.

the extended precision is necessary *if* "2*(s+s*R(z))" is close to pio2,
*or* s is close to -1 * s*R(z), so that it would cancel out several of
the highest bits.

Given that this block results in 1 > |x| >= 0.5, which is atleast
0.5 away from pio2 and 0, my gut tells me that this is not the case.

So my argument would be that sqrtf might be viable without affecting
the result at all or more than one ulp.
I know that this doesn't mean much without a proof.

Any hints where the algorithm is coming from?

Norbert