* [musl] Implementing csqrtl() @ 2022-07-04 9:35 Nikolaos Chatzikonstantinou 2022-07-04 11:09 ` [musl] " Nikolaos Chatzikonstantinou 2024-05-28 20:48 ` Nikolaos Chatzikonstantinou 0 siblings, 2 replies; 13+ messages in thread From: Nikolaos Chatzikonstantinou @ 2022-07-04 9:35 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 13+ messages in thread
* [musl] Re: Implementing csqrtl() 2022-07-04 9:35 [musl] Implementing csqrtl() Nikolaos Chatzikonstantinou @ 2022-07-04 11:09 ` Nikolaos Chatzikonstantinou 2022-07-05 9:37 ` Szabolcs Nagy 2024-05-28 20:48 ` Nikolaos Chatzikonstantinou 1 sibling, 1 reply; 13+ messages in thread From: Nikolaos Chatzikonstantinou @ 2022-07-04 11:09 UTC (permalink / raw) To: musl [-- 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 ^ permalink raw reply [flat|nested] 13+ messages in thread
* Re: [musl] Re: Implementing csqrtl() 2022-07-04 11:09 ` [musl] " Nikolaos Chatzikonstantinou @ 2022-07-05 9:37 ` Szabolcs Nagy 2022-07-05 14:28 ` Nikolaos Chatzikonstantinou 0 siblings, 1 reply; 13+ messages in thread From: Szabolcs Nagy @ 2022-07-05 9:37 UTC (permalink / raw) To: Nikolaos Chatzikonstantinou; +Cc: musl * 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 > ^ permalink raw reply [flat|nested] 13+ messages in thread
* Re: [musl] Re: Implementing csqrtl() 2022-07-05 9:37 ` Szabolcs Nagy @ 2022-07-05 14:28 ` Nikolaos Chatzikonstantinou 2022-07-05 15:35 ` Markus Wichmann 0 siblings, 1 reply; 13+ messages in thread From: Nikolaos Chatzikonstantinou @ 2022-07-05 14:28 UTC (permalink / raw) To: Nikolaos Chatzikonstantinou, musl 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 ^ permalink raw reply [flat|nested] 13+ messages in thread
* Re: [musl] Re: Implementing csqrtl() 2022-07-05 14:28 ` Nikolaos Chatzikonstantinou @ 2022-07-05 15:35 ` Markus Wichmann 2022-07-05 16:14 ` Nikolaos Chatzikonstantinou 0 siblings, 1 reply; 13+ messages in thread From: Markus Wichmann @ 2022-07-05 15:35 UTC (permalink / raw) To: musl; +Cc: Nikolaos Chatzikonstantinou 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 ^ permalink raw reply [flat|nested] 13+ messages in thread
* Re: [musl] Re: Implementing csqrtl() 2022-07-05 15:35 ` Markus Wichmann @ 2022-07-05 16:14 ` Nikolaos Chatzikonstantinou 2024-05-29 1:38 ` Rich Felker 0 siblings, 1 reply; 13+ messages in thread From: Nikolaos Chatzikonstantinou @ 2022-07-05 16:14 UTC (permalink / raw) To: Markus Wichmann; +Cc: musl 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 ^ permalink raw reply [flat|nested] 13+ messages in thread
* Re: [musl] Re: Implementing csqrtl() 2022-07-05 16:14 ` Nikolaos Chatzikonstantinou @ 2024-05-29 1:38 ` Rich Felker 0 siblings, 0 replies; 13+ messages in thread From: Rich Felker @ 2024-05-29 1:38 UTC (permalink / raw) To: Nikolaos Chatzikonstantinou; +Cc: Markus Wichmann, musl 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 ^ permalink raw reply [flat|nested] 13+ messages in thread
* [musl] Re: Implementing csqrtl() 2022-07-04 9:35 [musl] Implementing csqrtl() Nikolaos Chatzikonstantinou 2022-07-04 11:09 ` [musl] " Nikolaos Chatzikonstantinou @ 2024-05-28 20:48 ` Nikolaos Chatzikonstantinou 2024-05-29 0:38 ` Gabriel Ravier 1 sibling, 1 reply; 13+ messages in thread From: Nikolaos Chatzikonstantinou @ 2024-05-28 20:48 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 13+ messages in thread
* Re: [musl] Re: Implementing csqrtl() 2024-05-28 20:48 ` Nikolaos Chatzikonstantinou @ 2024-05-29 0:38 ` Gabriel Ravier 2024-05-29 1:34 ` Rich Felker 0 siblings, 1 reply; 13+ messages in thread From: Gabriel Ravier @ 2024-05-29 0:38 UTC (permalink / raw) To: musl, Nikolaos Chatzikonstantinou 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 ^ permalink raw reply [flat|nested] 13+ messages in thread
* Re: [musl] Re: Implementing csqrtl() 2024-05-29 0:38 ` Gabriel Ravier @ 2024-05-29 1:34 ` Rich Felker [not found] ` <CAAQmekdsDO=KuJ4CuVjTnTU8Baz3++LdUr=TeKJ4bcLTfJ+1gA@mail.gmail.com> 0 siblings, 1 reply; 13+ messages in thread From: Rich Felker @ 2024-05-29 1:34 UTC (permalink / raw) To: Gabriel Ravier; +Cc: musl, Nikolaos Chatzikonstantinou 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 ^ permalink raw reply [flat|nested] 13+ messages in thread
[parent not found: <CAAQmekdsDO=KuJ4CuVjTnTU8Baz3++LdUr=TeKJ4bcLTfJ+1gA@mail.gmail.com>]
* Re: [musl] Re: Implementing csqrtl() [not found] ` <CAAQmekdsDO=KuJ4CuVjTnTU8Baz3++LdUr=TeKJ4bcLTfJ+1gA@mail.gmail.com> @ 2024-05-29 12:59 ` Rich Felker 2024-05-29 15:13 ` James Y Knight 0 siblings, 1 reply; 13+ messages in thread From: Rich Felker @ 2024-05-29 12:59 UTC (permalink / raw) To: Nikolaos Chatzikonstantinou; +Cc: musl 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 ^ permalink raw reply [flat|nested] 13+ messages in thread
* Re: [musl] Re: Implementing csqrtl() 2024-05-29 12:59 ` Rich Felker @ 2024-05-29 15:13 ` James Y Knight 2024-05-31 15:37 ` Szabolcs Nagy 0 siblings, 1 reply; 13+ messages in thread From: James Y Knight @ 2024-05-29 15:13 UTC (permalink / raw) To: musl; +Cc: Nikolaos Chatzikonstantinou [-- 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 --] ^ permalink raw reply [flat|nested] 13+ messages in thread
* Re: [musl] Re: Implementing csqrtl() 2024-05-29 15:13 ` James Y Knight @ 2024-05-31 15:37 ` Szabolcs Nagy 0 siblings, 0 replies; 13+ messages in thread From: Szabolcs Nagy @ 2024-05-31 15:37 UTC (permalink / raw) To: James Y Knight; +Cc: musl, Nikolaos Chatzikonstantinou * 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. ^ permalink raw reply [flat|nested] 13+ messages in thread
end of thread, other threads:[~2024-05-31 15:37 UTC | newest] Thread overview: 13+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2022-07-04 9:35 [musl] Implementing csqrtl() Nikolaos Chatzikonstantinou 2022-07-04 11:09 ` [musl] " Nikolaos Chatzikonstantinou 2022-07-05 9:37 ` Szabolcs Nagy 2022-07-05 14:28 ` Nikolaos Chatzikonstantinou 2022-07-05 15:35 ` Markus Wichmann 2022-07-05 16:14 ` Nikolaos Chatzikonstantinou 2024-05-29 1:38 ` Rich Felker 2024-05-28 20:48 ` Nikolaos Chatzikonstantinou 2024-05-29 0:38 ` Gabriel Ravier 2024-05-29 1:34 ` Rich Felker [not found] ` <CAAQmekdsDO=KuJ4CuVjTnTU8Baz3++LdUr=TeKJ4bcLTfJ+1gA@mail.gmail.com> 2024-05-29 12:59 ` Rich Felker 2024-05-29 15:13 ` James Y Knight 2024-05-31 15:37 ` Szabolcs Nagy
Code repositories for project(s) associated with this public inbox https://git.vuxu.org/mirror/musl/ This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for NNTP newsgroup(s).