diff options
Diffstat (limited to 'src/complex/csqrt.c')
| -rw-r--r-- | src/complex/csqrt.c | 47 |
1 files changed, 46 insertions, 1 deletions
diff --git a/src/complex/csqrt.c b/src/complex/csqrt.c index f3a735ca..47b4c843 100644 --- a/src/complex/csqrt.c +++ b/src/complex/csqrt.c @@ -1,10 +1,55 @@ # define TGSOURCE "../complex/csqrt.c" #include "_tgmath.h" - #include <complex.h> +#include "math.h" +#include "fenv.h" TYPE complex TGFN(csqrt)(TYPE complex z) { + int classr = fpclassify(TGFN(creal)(z)); + int classi = fpclassify(TGFN(cimag)(z)); + int signr = signbit(TGFN(creal)(z)); + int signi = signbit(TGFN(cimag)(z)); + + if (classr == FP_ZERO && classi == FP_ZERO) { + return TGCMPLX(0.0, 0.0); + } + + if (classi == FP_INFINITE) { + /* for all reals, even NAN */ + return TGCMPLX(INFINITY, INFINITY); + } + + if (classr != FP_INFINITE && classi == FP_NAN) { + feraiseexcept(FE_INVALID); + return TGCMPLX(NAN, NAN); + } + + if (classr == FP_INFINITE && signr && classi != FP_INFINITE && !signi) { + return TGCMPLX(0.0, INFINITY); + } + + if (classr == FP_INFINITE && !signr && classi != FP_INFINITE && !signi) { + return TGCMPLX(INFINITY, 0.0); + } + + if (classr == FP_INFINITE && signr && classi == FP_NAN) { + return TGCMPLX(NAN, INFINITY); + } + + if (classr == FP_INFINITE && !signr && classi == FP_NAN) { + return TGCMPLX(INFINITY, NAN); + } + + if (classr == FP_NAN && classi != FP_INFINITE) { + feraiseexcept(FE_INVALID); + return TGCMPLX(NAN, NAN); + } + + if (classr == FP_NAN && classi == FP_NAN) { + return TGCMPLX(NAN, NAN); + } + return z; } |