* [musl] Implementing csqrtl()
@ 2022-07-04 9:35 Nikolaos Chatzikonstantinou
2022-07-04 11:09 ` [musl] " Nikolaos Chatzikonstantinou
0 siblings, 1 reply; 6+ messages in thread
From: Nikolaos Chatzikonstantinou @ 2022-07-04 9:35 UTC (permalink / raw)
To: musl
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.
As a very basic test, I wrote a program that produces random
complex numbers in the square [0, N] x [0, N] for N=1,100 and
calculates csqrt{,f,l}() with my implementation, glibc and the
arbitrary precision mpc_sqrt() from MPC,
<https://www.multiprecision.org/mpc/>.
Glibc stays almost within 1 ulp in float and double, but my
implementation wasn't so good with float. The double
implementation seems to get the exact same results as glibc does.
I was not able to even test the long double version with this
method, because I did not write the code that produces random
long double complex numbers yet.
There's a few things that I don't quite understand here. One is, I'm
not sure why Kahan's implementation is accurate. For another, I don't
know how to do any sort of speed tests; I've read online that
microbenchmarking is not reproducible for math functions, and so
google's benchmark <https://github.com/google/benchmark> does not seem
to help. Of course I'm aware that the C implementation would not be
the one used in most systems, as the assembly implementations are
usually better. Finally, I don't know what the right way to test an
implementation for accuracy is: whether by using automation or writing
proofs. It seems the state of the art has evolved quite a bit from
1987, and yet I don't know where to look for information on this
topic, as it seems very specific to chips & the C std lib.
Feel free to provide any sort of criticism. I'm e-mailing this
implementation for the purpose of starting a discussion, but I'm
hoping to be able to contribute something in the near future.
Regards,
Nikolaos Chatzikonstantinou
^ permalink raw reply [flat|nested] 6+ messages in thread
* [musl] Re: Implementing csqrtl()
2022-07-04 9:35 [musl] Implementing csqrtl() Nikolaos Chatzikonstantinou
@ 2022-07-04 11:09 ` Nikolaos Chatzikonstantinou
2022-07-05 9:37 ` Szabolcs Nagy
0 siblings, 1 reply; 6+ messages in thread
From: Nikolaos Chatzikonstantinou @ 2022-07-04 11:09 UTC (permalink / raw)
To: musl
[-- 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
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [musl] Re: Implementing csqrtl()
2022-07-04 11:09 ` [musl] " Nikolaos Chatzikonstantinou
@ 2022-07-05 9:37 ` Szabolcs Nagy
2022-07-05 14:28 ` Nikolaos Chatzikonstantinou
0 siblings, 1 reply; 6+ messages in thread
From: Szabolcs Nagy @ 2022-07-05 9:37 UTC (permalink / raw)
To: Nikolaos Chatzikonstantinou; +Cc: musl
* 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
>
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [musl] Re: Implementing csqrtl()
2022-07-05 9:37 ` Szabolcs Nagy
@ 2022-07-05 14:28 ` Nikolaos Chatzikonstantinou
2022-07-05 15:35 ` Markus Wichmann
0 siblings, 1 reply; 6+ messages in thread
From: Nikolaos Chatzikonstantinou @ 2022-07-05 14:28 UTC (permalink / raw)
To: Nikolaos Chatzikonstantinou, musl
On Tue, Jul 5, 2022 at 9:37 AM Szabolcs Nagy <nsz@port70.net> wrote:
>
> * 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:
> > > 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.
>
> this will need a description
> and some testing.
Yes, I agree.How should it be tested?
> > + if ((isinf(x) || isinf(y)) && (isnan(r) || isinf(r))) {
>
> why do you need the && ?
> can r be other than inf or nan?
It's the case that x^2 + y^2 is inf for x,y finite.
> > + k = logbl(fmaxl(fabsl(x), fabsl(y)));
> > + x = scalbnl(x, -k);
> > + y = scalbnl(y, -k);
>
> k is unsigned
Good catch! I changed k to int and logbl to ilogbl.
Regards,
Nikolaos Chatzikonstantinou
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [musl] Re: Implementing csqrtl()
2022-07-05 14:28 ` Nikolaos Chatzikonstantinou
@ 2022-07-05 15:35 ` Markus Wichmann
2022-07-05 16:14 ` Nikolaos Chatzikonstantinou
0 siblings, 1 reply; 6+ messages in thread
From: Markus Wichmann @ 2022-07-05 15:35 UTC (permalink / raw)
To: musl; +Cc: Nikolaos Chatzikonstantinou
On Tue, Jul 05, 2022 at 02:28:32PM +0000, Nikolaos Chatzikonstantinou wrote:
> On Tue, Jul 5, 2022 at 9:37 AM Szabolcs Nagy <nsz@port70.net> wrote:
> >
> > * Nikolaos Chatzikonstantinou <nchatz314@gmail.com> [2022-07-04 11:09:44 +0000]:
> >
> > > + if ((isinf(x) || isinf(y)) && (isnan(r) || isinf(r))) {
> >
> > why do you need the && ?
> > can r be other than inf or nan?
>
> It's the case that x^2 + y^2 is inf for x,y finite.
>
Yeah, but if x or y is infinite then so is r. So the entire part in
front of the && is redundant, because it is contained in the second
part.
Unless I completely misunderstood how IEEE infinity works...
Ciao,
Markus
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [musl] Re: Implementing csqrtl()
2022-07-05 15:35 ` Markus Wichmann
@ 2022-07-05 16:14 ` Nikolaos Chatzikonstantinou
0 siblings, 0 replies; 6+ messages in thread
From: Nikolaos Chatzikonstantinou @ 2022-07-05 16:14 UTC (permalink / raw)
To: Markus Wichmann; +Cc: musl
On Tue, Jul 5, 2022 at 3:35 PM Markus Wichmann <nullplan@gmx.net> wrote:
>
> On Tue, Jul 05, 2022 at 02:28:32PM +0000, Nikolaos Chatzikonstantinou wrote:
> > On Tue, Jul 5, 2022 at 9:37 AM Szabolcs Nagy <nsz@port70.net> wrote:
> > >
> > > * Nikolaos Chatzikonstantinou <nchatz314@gmail.com> [2022-07-04 11:09:44 +0000]:
> > >
> > > > + if ((isinf(x) || isinf(y)) && (isnan(r) || isinf(r))) {
> > >
> > > why do you need the && ?
> > > can r be other than inf or nan?
> >
> > It's the case that x^2 + y^2 is inf for x,y finite.
> >
>
> Yeah, but if x or y is infinite then so is r. So the entire part in
> front of the && is redundant, because it is contained in the second
> part.
>
> Unless I completely misunderstood how IEEE infinity works...
Yes you are right, it should be if(isinf(x) || isinf(y)).
Regards,
Nikolaos Chatzikonstantinou
^ permalink raw reply [flat|nested] 6+ messages in thread
end of thread, other threads:[~2022-07-05 16:29 UTC | newest]
Thread overview: 6+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2022-07-04 9:35 [musl] Implementing csqrtl() Nikolaos Chatzikonstantinou
2022-07-04 11:09 ` [musl] " Nikolaos Chatzikonstantinou
2022-07-05 9:37 ` Szabolcs Nagy
2022-07-05 14:28 ` Nikolaos Chatzikonstantinou
2022-07-05 15:35 ` Markus Wichmann
2022-07-05 16:14 ` Nikolaos Chatzikonstantinou
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).