[<prev] [next>] [thread-next>] [day] [month] [year] [list]
Date: Mon, 4 Jul 2022 09:35:05 +0000
From: Nikolaos Chatzikonstantinou <nchatz314@...il.com>
To: musl@...ts.openwall.com
Subject: Implementing csqrtl()
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
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.
