[<prev] [next>] [<thread-prev] [thread-next>] [day] [month] [year] [list]
Date: Tue, 5 Jul 2022 11:37:04 +0200
From: Szabolcs Nagy <nsz@...t70.net>
To: Nikolaos Chatzikonstantinou <nchatz314@...il.com>
Cc: musl@...ts.openwall.com
Subject: Re: Re: Implementing csqrtl()
* Nikolaos Chatzikonstantinou <nchatz314@...il.com> [2022-07-04 11:09:44 +0000]:
> On Mon, Jul 4, 2022 at 9:35 AM Nikolaos Chatzikonstantinou
> <nchatz314@...il.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@...il.com>
> Date: Mon, 4 Jul 2022 18:07:57 +0900
> Subject: [PATCH] add csqrtl() implementation
> To: musl@...ts.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
>
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
