Re: [musl] Re: Implementing csqrtl()

[Szabolcs Nagy <nsz@port70.net>]

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