```mailing list of musl libc
help / color / mirror / code / Atom feed```
```* [musl] Implementing csqrtl()
@ 2022-07-04  9:35 Nikolaos Chatzikonstantinou
2022-07-04 11:09 ` [musl] " Nikolaos Chatzikonstantinou
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
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()
2022-07-04  9:35 [musl] Implementing csqrtl() Nikolaos Chatzikonstantinou
@ 2022-07-04 11:09 ` Nikolaos Chatzikonstantinou
2022-07-05  9:37   ` Szabolcs Nagy
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

[-- 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
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()
2022-07-04 11:09 ` [musl] " Nikolaos Chatzikonstantinou
@ 2022-07-05  9:37   ` Szabolcs Nagy
2022-07-05 14:28     ` Nikolaos Chatzikonstantinou
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
>

```* 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
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

```* 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
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

```* Re: [musl] Re: Implementing csqrtl()
2022-07-05 15:35       ` Markus Wichmann
@ 2022-07-05 16:14         ` Nikolaos Chatzikonstantinou
0 siblings, 0 replies; 6+ 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

```end of thread, other threads:[~2022-07-05 16:29 UTC | newest]

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
```

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