Re: [musl] Re: Implementing csqrtl()

[Rich Felker <dalias@libc.org>]

On Tue, May 28, 2024 at 11:42:38PM -0400, Nikolaos Chatzikonstantinou wrote:
> On Tue, May 28, 2024 at 9:33 PM Rich Felker <dalias@libc.org> wrote:
> >
> > On Wed, May 29, 2024 at 02:38:22AM +0200, Gabriel Ravier wrote:
> > > On 5/28/24 10:48 PM, Nikolaos Chatzikonstantinou wrote:
> > > >On Mon, Jul 4, 2022 at 5:35 AM Nikolaos Chatzikonstantinou
> > > ><nchatz314@gmail.com> wrote:
> > > >>Hello list,
> > > >>
> > > >>I wanted to implement some function from
> > > >><https://wiki.musl-libc.org/open-issues.html#Complex-math>
> > > >>which is an open issue in the wiki.
> > > >>
> > > >>One of the missing complex functions is csqrtl(), the long double
> > > >>version of complex square root. I was able to find a 1987 article from
> > > >>W. Kahan, titled "Branch cuts for complex elementary functions." that
> > > >>contained an implementation for complex square root for arbitrary
> > > >>floating-point numbers. In this e-mail you'll find an attached git
> > > >>patch with the implementation.
> > > >>
> > > >>As a very basic test, I wrote a program that produces random
> > > >>complex numbers in the square [0, N] x [0, N] for N=1,100 and
> > > >>calculates csqrt{,f,l}() with my implementation, glibc and the
> > > >>arbitrary precision mpc_sqrt() from MPC,
> > > >><https://www.multiprecision.org/mpc/>.
> > > >>
> > > >>Glibc stays almost within 1 ulp in float and double, but my
> > > >>implementation wasn't so good with float. The double
> > > >>implementation seems to get the exact same results as glibc does.
> > > >>
> > > >>I was not able to even test the long double version with this
> > > >>method, because I did not write the code that produces random
> > > >>long double complex numbers yet.
> > > >>
> > > >>There's a few things that I don't quite understand here. One is, I'm
> > > >>not sure why Kahan's implementation is accurate. For another, I don't
> > > >>know how to do any sort of speed tests; I've read online that
> > > >>microbenchmarking is not reproducible for math functions, and so
> > > >>google's benchmark <https://github.com/google/benchmark> does not seem
> > > >>to help. Of course I'm aware that the C implementation would not be
> > > >>the one used in most systems, as the assembly implementations are
> > > >>usually better. Finally, I don't know what the right way to test an
> > > >>implementation for accuracy is: whether by using automation or writing
> > > >>proofs. It seems the state of the art has evolved quite a bit from
> > > >>1987, and yet I don't know where to look for information on this
> > > >>topic, as it seems very specific to chips & the C std lib.
> > > >>
> > > >>Feel free to provide any sort of criticism. I'm e-mailing this
> > > >>implementation for the purpose of starting a discussion, but I'm
> > > >>hoping to be able to contribute something in the near future.
> > > >Hello, I want to revive the discussion on this message. Last time I
> > > >attached the patch I was ignored, apart from some nitpicks on the C
> > > >source in the patch, which was a draft anyway. Missing complex math
> > > >functions are still an open issue. I can't help but think that last
> > > >time I was ignored apart from some nitpick on some signed/unsigned
> > > >type, as if getting schooled on C is the purpose of the mailing list
> > > >-- bizarre! the patch was a draft to start a discussion on resolving
> > > >the complex math issue. What I concluded from the responses I received
> > > >is that the musl project does not have a maintainer who is capable of
> > > >understanding the underlying issues with the complex math algorithms
> > > >in libc. Is that still the case and is there any interest in resolving
> > > >it perhaps by finding a maintainer who can work on this issue?
> > >
> > >
> > > I cannot read people's minds to determine exactly why it is that
> > > they did not thoroughly review the implementation but my guess would
> > > be that, unable to understand the mathematical details of how the
> > > code worked, they only reviewed the parts of it they were more
> > > familiar with (e.g. the "nitpick on some signed/unsigned type") - I
> > > guess it makes sense it would seem odd as the only thing to get
> > > reviewed at all, but normally those parts would also get reviewed
> > > regardless even if people were better able to understand the
> > > mathematics of the patch.
> >
> > Indeed, I didn't intend to ignore it, just left review to folks whose
> > focus is that part of the code. I'm sorry it looks like discussion
> > trailed off.
> >
> > With that said, the unsigned thing was not a nitpick afaict; it was
> > one reason the code would not work as written. I'll follow up with
> > some comments/review of my own on the code.
> 
> No, don't follow up with a review. I'm not looking for a review of the
> code right now. It's Kahan's algorithm from 1987, the implementation
> is irrelevant at this point in time. I included it just as a proof of
> concept, not as an actual candidate. Perhaps I shouldn't had included
> it at all.

I don't see any reason it's not a good candidate, but as implemented
it had bugs that were never resolved.

> From my original e-mail, I'm asking
> 
> (1) Why is Kahan's implementation accurate?

I'm not sure.

> (2) What sort of speed benchmarking tests would musl be satisfied with?

Not really needed unless the proposed change is intended to improve
something that's slow. If the existing code is incorrect, there's no
burden to *measure* performance of the correct replacement, just to
avoid doing things that are conceptually slow.

> (3) What sort of accuracy benchmarking tests would musl be satisfied with?

I'm not aware of whether there's a requirement on the complex sqrt
functions to be correctly rounded like for the real ones. If not, just
being a lot better than the wrong one we have now (in particular, not
turning finite inputs into infinities!) would be very welcome.

Paul Zimmermann is regularly testing and comparing different libcs'
math implementations, but I don't recall whether complex functions are
included. libc-test also has quite a few math tests, but no complex.
This is a potential area for enhancement, but not an acceptance
condition.

> (4) Where can I find more information on the topics above relevant to
> the requirements that musl places on a submission for complex math
> functions?

In short, musl does not have any "fixes must come with tests" policy.
The fix just needs to be defensibly correct, which depends on being
read by (and comprehensible to) the people familiar with that part of
the code.

One thing to point out from the old patch is that you can't depend on
fenv to tell you things like you did there. Not all archs (in
particular, not soft-float archs) have fenv available, and existing
compilers all just wrongly reorder things with respect to fenv access
and break it and don't care to fix this (probably because it would
hurt their benchmarks to fix it the right way).

So like, if you need to check for overflow, you can't perform the
operations then check fenv. You need to do like you would for integer
overflow, and check the conditions that would lead to overflow before
doing them.

> Obviously the nature of mailing lists is that conversations evolve and
> branch off, but I'm feeling that I'm again being ignored. When I have
> answers to the questions I'm asking, the next step would be for me to
> put effort in trying to get an actual patch going for musl. If for
> nothing else, the project can at least update the page mentioning the
> missing complex functions with the additional information that arises
> from this conversation.

Sounds good.

Rich