|
111
|
1 /* Complex square root of __float128 value.
|
|
|
2 Copyright (C) 1997-2012 Free Software Foundation, Inc.
|
|
|
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 #ifdef HAVE_FENV_H
|
|
|
24 # include <fenv.h>
|
|
|
25 #endif
|
|
|
26
|
|
|
27
|
|
|
28 __complex128
|
|
|
29 csqrtq (__complex128 x)
|
|
|
30 {
|
|
|
31 __complex128 res;
|
|
|
32 int rcls = fpclassifyq (__real__ x);
|
|
|
33 int icls = fpclassifyq (__imag__ x);
|
|
|
34
|
|
|
35 if (__builtin_expect (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE, 0))
|
|
|
36 {
|
|
|
37 if (icls == QUADFP_INFINITE)
|
|
|
38 {
|
|
|
39 __real__ res = HUGE_VALQ;
|
|
|
40 __imag__ res = __imag__ x;
|
|
|
41 }
|
|
|
42 else if (rcls == QUADFP_INFINITE)
|
|
|
43 {
|
|
|
44 if (__real__ x < 0.0Q)
|
|
|
45 {
|
|
|
46 __real__ res = icls == QUADFP_NAN ? nanq ("") : 0;
|
|
|
47 __imag__ res = copysignq (HUGE_VALQ, __imag__ x);
|
|
|
48 }
|
|
|
49 else
|
|
|
50 {
|
|
|
51 __real__ res = __real__ x;
|
|
|
52 __imag__ res = (icls == QUADFP_NAN
|
|
|
53 ? nanq ("") : copysignq (0.0Q, __imag__ x));
|
|
|
54 }
|
|
|
55 }
|
|
|
56 else
|
|
|
57 {
|
|
|
58 __real__ res = nanq ("");
|
|
|
59 __imag__ res = nanq ("");
|
|
|
60 }
|
|
|
61 }
|
|
|
62 else
|
|
|
63 {
|
|
|
64 if (__builtin_expect (icls == QUADFP_ZERO, 0))
|
|
|
65 {
|
|
|
66 if (__real__ x < 0.0Q)
|
|
|
67 {
|
|
|
68 __real__ res = 0.0Q;
|
|
|
69 __imag__ res = copysignq (sqrtq (-__real__ x),
|
|
|
70 __imag__ x);
|
|
|
71 }
|
|
|
72 else
|
|
|
73 {
|
|
|
74 __real__ res = fabsq (sqrtq (__real__ x));
|
|
|
75 __imag__ res = copysignq (0.0Q, __imag__ x);
|
|
|
76 }
|
|
|
77 }
|
|
|
78 else if (__builtin_expect (rcls == QUADFP_ZERO, 0))
|
|
|
79 {
|
|
|
80 __float128 r;
|
|
|
81 if (fabsq (__imag__ x) >= 2.0Q * FLT128_MIN)
|
|
|
82 r = sqrtq (0.5Q * fabsq (__imag__ x));
|
|
|
83 else
|
|
|
84 r = 0.5Q * sqrtq (2.0Q * fabsq (__imag__ x));
|
|
|
85
|
|
|
86 __real__ res = r;
|
|
|
87 __imag__ res = copysignq (r, __imag__ x);
|
|
|
88 }
|
|
|
89 else
|
|
|
90 {
|
|
|
91 __float128 d, r, s;
|
|
|
92 int scale = 0;
|
|
|
93
|
|
|
94 if (fabsq (__real__ x) > FLT128_MAX / 4.0Q)
|
|
|
95 {
|
|
|
96 scale = 1;
|
|
|
97 __real__ x = scalbnq (__real__ x, -2 * scale);
|
|
|
98 __imag__ x = scalbnq (__imag__ x, -2 * scale);
|
|
|
99 }
|
|
|
100 else if (fabsq (__imag__ x) > FLT128_MAX / 4.0Q)
|
|
|
101 {
|
|
|
102 scale = 1;
|
|
|
103 if (fabsq (__real__ x) >= 4.0Q * FLT128_MIN)
|
|
|
104 __real__ x = scalbnq (__real__ x, -2 * scale);
|
|
|
105 else
|
|
|
106 __real__ x = 0.0Q;
|
|
|
107 __imag__ x = scalbnq (__imag__ x, -2 * scale);
|
|
|
108 }
|
|
|
109 else if (fabsq (__real__ x) < FLT128_MIN
|
|
|
110 && fabsq (__imag__ x) < FLT128_MIN)
|
|
|
111 {
|
|
|
112 scale = -(FLT128_MANT_DIG / 2);
|
|
|
113 __real__ x = scalbnq (__real__ x, -2 * scale);
|
|
|
114 __imag__ x = scalbnq (__imag__ x, -2 * scale);
|
|
|
115 }
|
|
|
116
|
|
|
117 d = hypotq (__real__ x, __imag__ x);
|
|
|
118 /* Use the identity 2 Re res Im res = Im x
|
|
|
119 to avoid cancellation error in d +/- Re x. */
|
|
|
120 if (__real__ x > 0)
|
|
|
121 {
|
|
|
122 r = sqrtq (0.5Q * (d + __real__ x));
|
|
|
123 s = 0.5Q * (__imag__ x / r);
|
|
|
124 }
|
|
|
125 else
|
|
|
126 {
|
|
|
127 s = sqrtq (0.5Q * (d - __real__ x));
|
|
|
128 r = fabsq (0.5Q * (__imag__ x / s));
|
|
|
129 }
|
|
|
130
|
|
|
131 if (scale)
|
|
|
132 {
|
|
|
133 r = scalbnq (r, scale);
|
|
|
134 s = scalbnq (s, scale);
|
|
|
135 }
|
|
|
136
|
|
|
137 __real__ res = r;
|
|
|
138 __imag__ res = copysignq (s, __imag__ x);
|
|
|
139 }
|
|
|
140 }
|
|
|
141
|
|
|
142 return res;
|
|
|
143 }
|
