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