[musl] Implementing csqrtl()

[musl] Implementing csqrtl()

[Nikolaos Chatzikonstantinou]

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.

Regards,
Nikolaos Chatzikonstantinou

[musl] Re: Implementing csqrtl()

[Nikolaos Chatzikonstantinou]

[-- Attachment #1: Type: text/plain, Size: 718 bytes --]

On Mon, Jul 4, 2022 at 9: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.

I forgot to attach the patch, but here it is.

Regards,
Nikolaos Chatzikonstantinou

[-- Attachment #2: 0001-add-csqrtl-implementation.patch --]
[-- Type: text/x-patch, Size: 2132 bytes --]

From 069a4165a743e217a73064d74e72f292ca5e2fe2 Mon Sep 17 00:00:00 2001
From: Nikolaos Chatzikonstantinou <nchatz314@gmail.com>
Date: Mon, 4 Jul 2022 18:07:57 +0900
Subject: [PATCH] add csqrtl() implementation
To: musl@lists.openwall.com

---
 src/complex/csqrtl.c | 67 +++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 63 insertions(+), 4 deletions(-)

diff --git a/src/complex/csqrtl.c b/src/complex/csqrtl.c
index 22539379..d28ec8e5 100644
--- a/src/complex/csqrtl.c
+++ b/src/complex/csqrtl.c
@@ -1,7 +1,66 @@
 #include "complex_impl.h"
+#include <fenv.h>
 
-//FIXME
-long double complex csqrtl(long double complex z)
-{
-	return csqrt(z);
+/* cssqsl() and csqrtl() taken from
+ * Kahan, W. (1987). Branch cuts for complex elementary functions.
+ */
+static inline long double complex _cssqsl(long double complex z) {
+#pragma STDC FENV_ACCESS ON
+  fenv_t env;
+  unsigned k = 0;
+  long double x, y, r;
+  int set_excepts;
+
+  feholdexcept(&env);
+  x = creal(z);
+  y = cimag(z);
+  r = x * x + y * y;
+  if ((isinf(x) || isinf(y)) && (isnan(r) || isinf(r))) {
+    r = INFINITY;
+  } else {
+    set_excepts = fetestexcept(FE_OVERFLOW | FE_UNDERFLOW);
+    if ((set_excepts & FE_OVERFLOW) ||
+        ((set_excepts & FE_UNDERFLOW) && isless(r, LDBL_MIN / LDBL_EPSILON))) {
+      k = logbl(fmaxl(fabsl(x), fabsl(y)));
+      x = scalbnl(x, -k);
+      y = scalbnl(y, -k);
+      r = x * x + y * y;
+    }
+  }
+  feupdateenv(&env);
+  return CMPLXL(r, k);
+}
+
+long double complex csqrtl(long double complex z) {
+  long double x, y, r, xi, eta;
+  unsigned k;
+
+  x = creal(z);
+  y = cimag(z);
+  z = _cssqsl(z);
+  r = creal(z);
+  k = cimag(z);
+  if (!isnan(x)) {
+    r = scalbnl(fabsl(x), -k) + sqrtl(r);
+  }
+  if (k & 1) {
+    k = (k - 1) / 2;
+  } else {
+    k = k / 2 - 1;
+    r *= 2;
+  }
+  r = scalbnl(sqrtl(r), k);
+  xi = r;
+  eta = y;
+  if (r != 0) {
+    if (!isinf(eta)) {
+      // TODO if eta underflowed, signal it
+      eta = (eta / r) / 2;
+    }
+    if (isless(x, 0)) {
+      xi = fabsl(eta);
+      eta = copysignl(r, y);
+    }
+  }
+  return CMPLX(xi, eta);
 }
-- 
2.36.1

Re: [musl] Re: Implementing csqrtl()

[Szabolcs Nagy]

* Nikolaos Chatzikonstantinou <nchatz314@gmail.com> [2022-07-04 11:09:44 +0000]:

> On Mon, Jul 4, 2022 at 9: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.
> 
> I forgot to attach the patch, but here it is.
> 
> Regards,
> Nikolaos Chatzikonstantinou

> From 069a4165a743e217a73064d74e72f292ca5e2fe2 Mon Sep 17 00:00:00 2001
> From: Nikolaos Chatzikonstantinou <nchatz314@gmail.com>
> Date: Mon, 4 Jul 2022 18:07:57 +0900
> Subject: [PATCH] add csqrtl() implementation
> To: musl@lists.openwall.com
> 

this will need a description
and some testing.

> ---
>  src/complex/csqrtl.c | 67 +++++++++++++++++++++++++++++++++++++++++---
>  1 file changed, 63 insertions(+), 4 deletions(-)
> 
> diff --git a/src/complex/csqrtl.c b/src/complex/csqrtl.c
> index 22539379..d28ec8e5 100644
> --- a/src/complex/csqrtl.c
> +++ b/src/complex/csqrtl.c
> @@ -1,7 +1,66 @@
>  #include "complex_impl.h"
> +#include <fenv.h>
>  
> -//FIXME
> -long double complex csqrtl(long double complex z)
> -{
> -	return csqrt(z);
> +/* cssqsl() and csqrtl() taken from
> + * Kahan, W. (1987). Branch cuts for complex elementary functions.
> + */
> +static inline long double complex _cssqsl(long double complex z) {
> +#pragma STDC FENV_ACCESS ON
> +  fenv_t env;
> +  unsigned k = 0;
> +  long double x, y, r;
> +  int set_excepts;
> +
> +  feholdexcept(&env);
> +  x = creal(z);
> +  y = cimag(z);
> +  r = x * x + y * y;
> +  if ((isinf(x) || isinf(y)) && (isnan(r) || isinf(r))) {

why do you need the && ?
can r be other than inf or nan?

> +    r = INFINITY;
> +  } else {
> +    set_excepts = fetestexcept(FE_OVERFLOW | FE_UNDERFLOW);
> +    if ((set_excepts & FE_OVERFLOW) ||
> +        ((set_excepts & FE_UNDERFLOW) && isless(r, LDBL_MIN / LDBL_EPSILON))) {
> +      k = logbl(fmaxl(fabsl(x), fabsl(y)));
> +      x = scalbnl(x, -k);
> +      y = scalbnl(y, -k);

k is unsigned

> +      r = x * x + y * y;
> +    }
> +  }
> +  feupdateenv(&env);
> +  return CMPLXL(r, k);
> +}
> +
> +long double complex csqrtl(long double complex z) {
> +  long double x, y, r, xi, eta;
> +  unsigned k;
> +
> +  x = creal(z);
> +  y = cimag(z);
> +  z = _cssqsl(z);
> +  r = creal(z);
> +  k = cimag(z);
> +  if (!isnan(x)) {
> +    r = scalbnl(fabsl(x), -k) + sqrtl(r);

k is unsigned

> +  }
> +  if (k & 1) {
> +    k = (k - 1) / 2;
> +  } else {
> +    k = k / 2 - 1;
> +    r *= 2;
> +  }
> +  r = scalbnl(sqrtl(r), k);
> +  xi = r;
> +  eta = y;
> +  if (r != 0) {
> +    if (!isinf(eta)) {
> +      // TODO if eta underflowed, signal it
> +      eta = (eta / r) / 2;
> +    }
> +    if (isless(x, 0)) {
> +      xi = fabsl(eta);
> +      eta = copysignl(r, y);
> +    }
> +  }
> +  return CMPLX(xi, eta);
>  }
> -- 
> 2.36.1
> 

Re: [musl] Re: Implementing csqrtl()

[Nikolaos Chatzikonstantinou]

On Tue, Jul 5, 2022 at 9:37 AM Szabolcs Nagy <nsz@port70.net> wrote:
>
> * Nikolaos Chatzikonstantinou <nchatz314@gmail.com> [2022-07-04 11:09:44 +0000]:
>
> > On Mon, Jul 4, 2022 at 9:35 AM Nikolaos Chatzikonstantinou
> > <nchatz314@gmail.com> wrote:
> > > 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.
> >
> > I forgot to attach the patch, but here it is.
>
> this will need a description
> and some testing.

Yes, I agree.How should it be tested?

> > +  if ((isinf(x) || isinf(y)) && (isnan(r) || isinf(r))) {
>
> why do you need the && ?
> can r be other than inf or nan?

It's the case that x^2 + y^2 is inf for x,y finite.

> > +      k = logbl(fmaxl(fabsl(x), fabsl(y)));
> > +      x = scalbnl(x, -k);
> > +      y = scalbnl(y, -k);
>
> k is unsigned

Good catch! I changed k to int and logbl to ilogbl.

Regards,
Nikolaos Chatzikonstantinou

Re: [musl] Re: Implementing csqrtl()

[Markus Wichmann]

On Tue, Jul 05, 2022 at 02:28:32PM +0000, Nikolaos Chatzikonstantinou wrote:
> On Tue, Jul 5, 2022 at 9:37 AM Szabolcs Nagy <nsz@port70.net> wrote:
> >
> > * Nikolaos Chatzikonstantinou <nchatz314@gmail.com> [2022-07-04 11:09:44 +0000]:
> >
> > > +  if ((isinf(x) || isinf(y)) && (isnan(r) || isinf(r))) {
> >
> > why do you need the && ?
> > can r be other than inf or nan?
>
> It's the case that x^2 + y^2 is inf for x,y finite.
>

Yeah, but if x or y is infinite then so is r. So the entire part in
front of the && is redundant, because it is contained in the second
part.

Unless I completely misunderstood how IEEE infinity works...

Ciao,
Markus

Re: [musl] Re: Implementing csqrtl()

[Nikolaos Chatzikonstantinou]

On Tue, Jul 5, 2022 at 3:35 PM Markus Wichmann <nullplan@gmx.net> wrote:
>
> On Tue, Jul 05, 2022 at 02:28:32PM +0000, Nikolaos Chatzikonstantinou wrote:
> > On Tue, Jul 5, 2022 at 9:37 AM Szabolcs Nagy <nsz@port70.net> wrote:
> > >
> > > * Nikolaos Chatzikonstantinou <nchatz314@gmail.com> [2022-07-04 11:09:44 +0000]:
> > >
> > > > +  if ((isinf(x) || isinf(y)) && (isnan(r) || isinf(r))) {
> > >
> > > why do you need the && ?
> > > can r be other than inf or nan?
> >
> > It's the case that x^2 + y^2 is inf for x,y finite.
> >
>
> Yeah, but if x or y is infinite then so is r. So the entire part in
> front of the && is redundant, because it is contained in the second
> part.
>
> Unless I completely misunderstood how IEEE infinity works...

Yes you are right, it should be if(isinf(x) || isinf(y)).

Regards,
Nikolaos Chatzikonstantinou

[musl] Re: Implementing csqrtl()

[Nikolaos Chatzikonstantinou]

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?

Regards,
Nikolaos Chatzikonstantinou

Re: [musl] Re: Implementing csqrtl()

[Gabriel Ravier]

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.


>
> Regards,
> Nikolaos Chatzikonstantinou

Re: [musl] Re: Implementing csqrtl()

[Rich Felker]

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.

Rich

Re: [musl] Re: Implementing csqrtl()

[Rich Felker]

On Tue, Jul 05, 2022 at 04:14:13PM +0000, Nikolaos Chatzikonstantinou wrote:
> On Tue, Jul 5, 2022 at 3:35 PM Markus Wichmann <nullplan@gmx.net> wrote:
> >
> > On Tue, Jul 05, 2022 at 02:28:32PM +0000, Nikolaos Chatzikonstantinou wrote:
> > > On Tue, Jul 5, 2022 at 9:37 AM Szabolcs Nagy <nsz@port70.net> wrote:
> > > >
> > > > * Nikolaos Chatzikonstantinou <nchatz314@gmail.com> [2022-07-04 11:09:44 +0000]:
> > > >
> > > > > +  if ((isinf(x) || isinf(y)) && (isnan(r) || isinf(r))) {
> > > >
> > > > why do you need the && ?
> > > > can r be other than inf or nan?
> > >
> > > It's the case that x^2 + y^2 is inf for x,y finite.
> > >
> >
> > Yeah, but if x or y is infinite then so is r. So the entire part in
> > front of the && is redundant, because it is contained in the second
> > part.
> >
> > Unless I completely misunderstood how IEEE infinity works...
> 
> Yes you are right, it should be if(isinf(x) || isinf(y)).

I don't think this is correct. It's possible for r (which actually
holds a value I would call r², perhaps r2 as var name?) to be infinite
even when x and y are finite. This sounds problematic. Don't you need
to work with actual r, guaranteed to be finite when x and y are, which
you could obtain from hypotl?

Rich

Re: [musl] Re: Implementing csqrtl()

[Rich Felker]

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

Re: [musl] Re: Implementing csqrtl()

[James Y Knight]

[-- Attachment #1: Type: text/plain, Size: 1240 bytes --]

On Wed, May 29, 2024 at 8:58 AM Rich Felker <dalias@libc.org> wrote:

> 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).


FWIW (not that it changes the recommendation in this instance), Clang has
put a bunch of work into addressing this and other floating-point
compliance issues over the past decade. It *is* considered an important
thing to implement correctly, and I believe it generally is working on the
common architectures, since a few years ago. Unfortunately it's not yet
implemented on every architecture.

There is no impact on benchmarks by implementing fenv correctly, because,
per the C standard, it is undefined behavior to access the floating-point
status flags or to execute any code with non-default floating-point control
modes in effect unless "fenv_access" is enabled -- either using "#pragma
STDC FENV_ACCESS ON" in the code, or via an implementation-defined
mechanism to set the default to on (in clang, the command-line option
"-ffp-model=strict"). Because it's opt-in, only the code which actually
requires it takes the performance hit.

[-- Attachment #2: Type: text/html, Size: 1888 bytes --]

Re: [musl] Re: Implementing csqrtl()

[Szabolcs Nagy]

* James Y Knight <jyknight@google.com> [2024-05-29 11:13:12 -0400]:
> There is no impact on benchmarks by implementing fenv correctly, because,
> per the C standard, it is undefined behavior to access the floating-point
> status flags or to execute any code with non-default floating-point control
> modes in effect unless "fenv_access" is enabled -- either using "#pragma
> STDC FENV_ACCESS ON" in the code, or via an implementation-defined
> mechanism to set the default to on (in clang, the command-line option
> "-ffp-model=strict"). Because it's opt-in, only the code which actually
> requires it takes the performance hit.

note: the libc *must* opt-in if it tries to conform to c99 annex f.
(at least for library calls that are expected to work in non-default
rounding mode and math code that is specified to set exception flags.)

so fenv correctness has an impact on libc and thus on code that calls
into libc.