1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
|
#include <complex.h> // for complex, I, cimagl, creall, csqrtl
#include <float.h> // for LDBL_MAX
#include <math.h> // for fabsl, copysignl, hypotl, sqrtl, isinf, INFINITY
#include <stdbool.h> // for bool, false, true
#define THRESH (LDBL_MAX / 2.414213562373095048801688724209698L)
#define cpackl(r, i) ((r) + (i) * (long double complex)I)
long double complex csqrtl(long double complex z)
{
long double complex result;
long double a, b;
long double t;
bool scale;
a = creall(z);
b = cimagl(z);
/* Handle special cases. */
if (z == 0.0L)
return (cpackl((long double)0.0L, b));
if (isinf(b))
return (cpackl((long double)INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return (cpackl(a, t)); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
* csqrt(inf + NaN i) = inf + NaN i
* csqrt(inf + y i) = inf + 0 i
* csqrt(-inf + NaN i) = NaN +- inf i
* csqrt(-inf + y i) = 0 + inf i
*/
if (signbit(a))
return (cpackl(fabsl(b - b), copysignl(a, b)));
return (cpackl(a, copysignl(b - b, b)));
}
/*
* The remaining special case (b is NaN) is handled just fine by
* the normal code path below.
*/
/* Scale to avoid overflow. */
if (fabsl(a) >= THRESH || fabsl(b) >= THRESH) {
a *= 0.25L;
b *= 0.25L;
scale = true;
} else {
scale = false;
}
/* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0L) {
t = sqrtl((a + hypotl(a, b)) * 0.5L);
result = cpackl(t, b / (2.0L * t));
} else {
t = sqrtl((-a + hypotl(a, b)) * 0.5L);
result = cpackl(fabsl(b) / (2.0L * t), copysignl(t, b));
}
/* Rescale. */
if (scale)
return (result * 2.0L);
return (result);
}
|
