|
145
|
1 /* Complex square root of a float type.
|
|
|
2 Copyright (C) 1997-2018 Free Software Foundation, Inc.
|
|
111
|
3 This file is part of the GNU C Library.
|
|
|
4 Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>.
|
|
|
5 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
|
|
|
6
|
|
|
7 The GNU C Library is free software; you can redistribute it and/or
|
|
|
8 modify it under the terms of the GNU Lesser General Public
|
|
|
9 License as published by the Free Software Foundation; either
|
|
|
10 version 2.1 of the License, or (at your option) any later version.
|
|
|
11
|
|
|
12 The GNU C Library is distributed in the hope that it will be useful,
|
|
|
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
|
15 Lesser General Public License for more details.
|
|
|
16
|
|
|
17 You should have received a copy of the GNU Lesser General Public
|
|
|
18 License along with the GNU C Library; if not, see
|
|
|
19 <http://www.gnu.org/licenses/>. */
|
|
|
20
|
|
|
21 #include "quadmath-imp.h"
|
|
|
22
|
|
|
23 __complex128
|
|
|
24 csqrtq (__complex128 x)
|
|
|
25 {
|
|
|
26 __complex128 res;
|
|
|
27 int rcls = fpclassifyq (__real__ x);
|
|
|
28 int icls = fpclassifyq (__imag__ x);
|
|
|
29
|
|
145
|
30 if (__glibc_unlikely (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE))
|
|
111
|
31 {
|
|
|
32 if (icls == QUADFP_INFINITE)
|
|
|
33 {
|
|
|
34 __real__ res = HUGE_VALQ;
|
|
|
35 __imag__ res = __imag__ x;
|
|
|
36 }
|
|
|
37 else if (rcls == QUADFP_INFINITE)
|
|
|
38 {
|
|
145
|
39 if (__real__ x < 0)
|
|
111
|
40 {
|
|
|
41 __real__ res = icls == QUADFP_NAN ? nanq ("") : 0;
|
|
|
42 __imag__ res = copysignq (HUGE_VALQ, __imag__ x);
|
|
|
43 }
|
|
|
44 else
|
|
|
45 {
|
|
|
46 __real__ res = __real__ x;
|
|
|
47 __imag__ res = (icls == QUADFP_NAN
|
|
145
|
48 ? nanq ("") : copysignq (0, __imag__ x));
|
|
111
|
49 }
|
|
|
50 }
|
|
|
51 else
|
|
|
52 {
|
|
|
53 __real__ res = nanq ("");
|
|
|
54 __imag__ res = nanq ("");
|
|
|
55 }
|
|
|
56 }
|
|
|
57 else
|
|
|
58 {
|
|
145
|
59 if (__glibc_unlikely (icls == QUADFP_ZERO))
|
|
111
|
60 {
|
|
145
|
61 if (__real__ x < 0)
|
|
111
|
62 {
|
|
145
|
63 __real__ res = 0;
|
|
|
64 __imag__ res = copysignq (sqrtq (-__real__ x), __imag__ x);
|
|
111
|
65 }
|
|
|
66 else
|
|
|
67 {
|
|
|
68 __real__ res = fabsq (sqrtq (__real__ x));
|
|
145
|
69 __imag__ res = copysignq (0, __imag__ x);
|
|
111
|
70 }
|
|
|
71 }
|
|
145
|
72 else if (__glibc_unlikely (rcls == QUADFP_ZERO))
|
|
111
|
73 {
|
|
|
74 __float128 r;
|
|
145
|
75 if (fabsq (__imag__ x) >= 2 * FLT128_MIN)
|
|
111
|
76 r = sqrtq (0.5Q * fabsq (__imag__ x));
|
|
|
77 else
|
|
145
|
78 r = 0.5Q * sqrtq (2 * fabsq (__imag__ x));
|
|
111
|
79
|
|
|
80 __real__ res = r;
|
|
|
81 __imag__ res = copysignq (r, __imag__ x);
|
|
|
82 }
|
|
|
83 else
|
|
|
84 {
|
|
|
85 __float128 d, r, s;
|
|
|
86 int scale = 0;
|
|
|
87
|
|
145
|
88 if (fabsq (__real__ x) > FLT128_MAX / 4)
|
|
111
|
89 {
|
|
|
90 scale = 1;
|
|
|
91 __real__ x = scalbnq (__real__ x, -2 * scale);
|
|
|
92 __imag__ x = scalbnq (__imag__ x, -2 * scale);
|
|
|
93 }
|
|
145
|
94 else if (fabsq (__imag__ x) > FLT128_MAX / 4)
|
|
111
|
95 {
|
|
|
96 scale = 1;
|
|
145
|
97 if (fabsq (__real__ x) >= 4 * FLT128_MIN)
|
|
111
|
98 __real__ x = scalbnq (__real__ x, -2 * scale);
|
|
|
99 else
|
|
145
|
100 __real__ x = 0;
|
|
111
|
101 __imag__ x = scalbnq (__imag__ x, -2 * scale);
|
|
|
102 }
|
|
145
|
103 else if (fabsq (__real__ x) < 2 * FLT128_MIN
|
|
|
104 && fabsq (__imag__ x) < 2 * FLT128_MIN)
|
|
111
|
105 {
|
|
145
|
106 scale = -((FLT128_MANT_DIG + 1) / 2);
|
|
111
|
107 __real__ x = scalbnq (__real__ x, -2 * scale);
|
|
|
108 __imag__ x = scalbnq (__imag__ x, -2 * scale);
|
|
|
109 }
|
|
|
110
|
|
|
111 d = hypotq (__real__ x, __imag__ x);
|
|
|
112 /* Use the identity 2 Re res Im res = Im x
|
|
|
113 to avoid cancellation error in d +/- Re x. */
|
|
|
114 if (__real__ x > 0)
|
|
|
115 {
|
|
|
116 r = sqrtq (0.5Q * (d + __real__ x));
|
|
145
|
117 if (scale == 1 && fabsq (__imag__ x) < 1)
|
|
|
118 {
|
|
|
119 /* Avoid possible intermediate underflow. */
|
|
|
120 s = __imag__ x / r;
|
|
|
121 r = scalbnq (r, scale);
|
|
|
122 scale = 0;
|
|
|
123 }
|
|
|
124 else
|
|
|
125 s = 0.5Q * (__imag__ x / r);
|
|
111
|
126 }
|
|
|
127 else
|
|
|
128 {
|
|
|
129 s = sqrtq (0.5Q * (d - __real__ x));
|
|
145
|
130 if (scale == 1 && fabsq (__imag__ x) < 1)
|
|
|
131 {
|
|
|
132 /* Avoid possible intermediate underflow. */
|
|
|
133 r = fabsq (__imag__ x / s);
|
|
|
134 s = scalbnq (s, scale);
|
|
|
135 scale = 0;
|
|
|
136 }
|
|
|
137 else
|
|
|
138 r = fabsq (0.5Q * (__imag__ x / s));
|
|
111
|
139 }
|
|
|
140
|
|
|
141 if (scale)
|
|
|
142 {
|
|
|
143 r = scalbnq (r, scale);
|
|
|
144 s = scalbnq (s, scale);
|
|
|
145 }
|
|
|
146
|
|
145
|
147 math_check_force_underflow (r);
|
|
|
148 math_check_force_underflow (s);
|
|
|
149
|
|
111
|
150 __real__ res = r;
|
|
|
151 __imag__ res = copysignq (s, __imag__ x);
|
|
|
152 }
|
|
|
153 }
|
|
|
154
|
|
|
155 return res;
|
|
|
156 }
|
