From: Nikolaos Chatzikonstantinou <nchatz314@gmail.com>
To: musl@lists.openwall.com
Subject: [musl] Re: Implementing csqrtl()
Date: Mon, 4 Jul 2022 11:09:44 +0000 [thread overview]
Message-ID: <CAAQmekc76iSjRots+cQuWMMsBPnPEV5emiCH8wAsmVxRH+OpEA@mail.gmail.com> (raw)
In-Reply-To: <CAAQmekd_LMvH+3sTu8E8LNB4zop_aB9BgKRe2M+HKtd6nERTzg@mail.gmail.com>
[-- 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
next prev parent reply other threads:[~2022-07-04 11:31 UTC|newest]
Thread overview: 6+ messages / expand[flat|nested] mbox.gz Atom feed top
2022-07-04 9:35 [musl] " Nikolaos Chatzikonstantinou
2022-07-04 11:09 ` Nikolaos Chatzikonstantinou [this message]
2022-07-05 9:37 ` [musl] " Szabolcs Nagy
2022-07-05 14:28 ` Nikolaos Chatzikonstantinou
2022-07-05 15:35 ` Markus Wichmann
2022-07-05 16:14 ` Nikolaos Chatzikonstantinou
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=CAAQmekc76iSjRots+cQuWMMsBPnPEV5emiCH8wAsmVxRH+OpEA@mail.gmail.com \
--to=nchatz314@gmail.com \
--cc=musl@lists.openwall.com \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
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).