From mboxrd@z Thu Jan 1 00:00:00 1970 X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on inbox.vuxu.org X-Spam-Level: X-Spam-Status: No, score=-1.0 required=5.0 tests=MAILING_LIST_MULTI, T_SCC_BODY_TEXT_LINE autolearn=ham autolearn_force=no version=3.4.4 Received: (qmail 20170 invoked from network); 5 Jul 2022 09:37:20 -0000 Received: from second.openwall.net (193.110.157.125) by inbox.vuxu.org with ESMTPUTF8; 5 Jul 2022 09:37:20 -0000 Received: (qmail 17820 invoked by uid 550); 5 Jul 2022 09:37:17 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: List-ID: Reply-To: musl@lists.openwall.com Received: (qmail 17788 invoked from network); 5 Jul 2022 09:37:16 -0000 Date: Tue, 5 Jul 2022 11:37:04 +0200 From: Szabolcs Nagy To: Nikolaos Chatzikonstantinou Cc: musl@lists.openwall.com Message-ID: <20220705093704.GY1320090@port70.net> Mail-Followup-To: Nikolaos Chatzikonstantinou , musl@lists.openwall.com References: MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: Subject: Re: [musl] Re: Implementing csqrtl() * Nikolaos Chatzikonstantinou [2022-07-04 11:09:44 +0000]: > On Mon, Jul 4, 2022 at 9:35 AM Nikolaos Chatzikonstantinou > wrote: > > > > Hello list, > > > > I wanted to implement some function from > > > > 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 > 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 > > -//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 >