blob: 77d61ec2d2322f87e9055123d273334b35d28993 (
plain)
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
|
#include <complex.h> // for complex, I, cabs, cimag, creal, csqrt
#include <math.h> // for fabs, sqrt
double complex csqrt(double complex z)
{
double complex w;
double x, y, r, t, scale;
x = creal(z);
y = cimag(z);
if (y == 0.0) {
if (x == 0.0) {
return 0.0 + y * (double complex)I;
}
r = fabs(x);
r = sqrt(r);
if (x < 0.0) {
return 0.0 + r * (double complex)I;
}
return r + y * (double complex)I;
}
if (x == 0.0) {
r = fabs(y);
r = sqrt(0.5 * r);
if (y > 0)
return r + r * (double complex)I;
return r - r * (double complex)I;
}
if ((fabs(x) > 4.0) || (fabs(y) > 4.0)) {
x *= 0.25;
y *= 0.25;
scale = 2.0;
} else {
x *= 1.8014398509481984e16; /* 2^54 */
y *= 1.8014398509481984e16;
scale = 7.450580596923828125e-9; /* 2^-27 */
}
w = x + y * (double complex)I;
r = cabs(w);
if (x > 0) {
t = sqrt(0.5 * r + 0.5 * x);
r = scale * fabs((0.5 * y) / t);
t *= scale;
} else {
r = sqrt(0.5 * r - 0.5 * x);
t = scale * fabs((0.5 * y) / r);
r *= scale;
}
if (y < 0)
return t - r * (double complex)I;
return t + r * (double complex)I;
}
|
