Mercurial > hg > CbC > CbC_gcc
comparison libquadmath/math/csqrtq.c @ 145:1830386684a0
gcc-9.2.0
author | anatofuz |
---|---|
date | Thu, 13 Feb 2020 11:34:05 +0900 |
parents | 04ced10e8804 |
children |
comparison
equal
deleted
inserted
replaced
131:84e7813d76e9 | 145:1830386684a0 |
---|---|
1 /* Complex square root of __float128 value. | 1 /* Complex square root of a float type. |
2 Copyright (C) 1997-2012 Free Software Foundation, Inc. | 2 Copyright (C) 1997-2018 Free Software Foundation, Inc. |
3 This file is part of the GNU C Library. | 3 This file is part of the GNU C Library. |
4 Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>. | 4 Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>. |
5 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. | 5 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. |
6 | 6 |
7 The GNU C Library is free software; you can redistribute it and/or | 7 The GNU C Library is free software; you can redistribute it and/or |
18 License along with the GNU C Library; if not, see | 18 License along with the GNU C Library; if not, see |
19 <http://www.gnu.org/licenses/>. */ | 19 <http://www.gnu.org/licenses/>. */ |
20 | 20 |
21 #include "quadmath-imp.h" | 21 #include "quadmath-imp.h" |
22 | 22 |
23 #ifdef HAVE_FENV_H | |
24 # include <fenv.h> | |
25 #endif | |
26 | |
27 | |
28 __complex128 | 23 __complex128 |
29 csqrtq (__complex128 x) | 24 csqrtq (__complex128 x) |
30 { | 25 { |
31 __complex128 res; | 26 __complex128 res; |
32 int rcls = fpclassifyq (__real__ x); | 27 int rcls = fpclassifyq (__real__ x); |
33 int icls = fpclassifyq (__imag__ x); | 28 int icls = fpclassifyq (__imag__ x); |
34 | 29 |
35 if (__builtin_expect (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE, 0)) | 30 if (__glibc_unlikely (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE)) |
36 { | 31 { |
37 if (icls == QUADFP_INFINITE) | 32 if (icls == QUADFP_INFINITE) |
38 { | 33 { |
39 __real__ res = HUGE_VALQ; | 34 __real__ res = HUGE_VALQ; |
40 __imag__ res = __imag__ x; | 35 __imag__ res = __imag__ x; |
41 } | 36 } |
42 else if (rcls == QUADFP_INFINITE) | 37 else if (rcls == QUADFP_INFINITE) |
43 { | 38 { |
44 if (__real__ x < 0.0Q) | 39 if (__real__ x < 0) |
45 { | 40 { |
46 __real__ res = icls == QUADFP_NAN ? nanq ("") : 0; | 41 __real__ res = icls == QUADFP_NAN ? nanq ("") : 0; |
47 __imag__ res = copysignq (HUGE_VALQ, __imag__ x); | 42 __imag__ res = copysignq (HUGE_VALQ, __imag__ x); |
48 } | 43 } |
49 else | 44 else |
50 { | 45 { |
51 __real__ res = __real__ x; | 46 __real__ res = __real__ x; |
52 __imag__ res = (icls == QUADFP_NAN | 47 __imag__ res = (icls == QUADFP_NAN |
53 ? nanq ("") : copysignq (0.0Q, __imag__ x)); | 48 ? nanq ("") : copysignq (0, __imag__ x)); |
54 } | 49 } |
55 } | 50 } |
56 else | 51 else |
57 { | 52 { |
58 __real__ res = nanq (""); | 53 __real__ res = nanq (""); |
59 __imag__ res = nanq (""); | 54 __imag__ res = nanq (""); |
60 } | 55 } |
61 } | 56 } |
62 else | 57 else |
63 { | 58 { |
64 if (__builtin_expect (icls == QUADFP_ZERO, 0)) | 59 if (__glibc_unlikely (icls == QUADFP_ZERO)) |
65 { | 60 { |
66 if (__real__ x < 0.0Q) | 61 if (__real__ x < 0) |
67 { | 62 { |
68 __real__ res = 0.0Q; | 63 __real__ res = 0; |
69 __imag__ res = copysignq (sqrtq (-__real__ x), | 64 __imag__ res = copysignq (sqrtq (-__real__ x), __imag__ x); |
70 __imag__ x); | |
71 } | 65 } |
72 else | 66 else |
73 { | 67 { |
74 __real__ res = fabsq (sqrtq (__real__ x)); | 68 __real__ res = fabsq (sqrtq (__real__ x)); |
75 __imag__ res = copysignq (0.0Q, __imag__ x); | 69 __imag__ res = copysignq (0, __imag__ x); |
76 } | 70 } |
77 } | 71 } |
78 else if (__builtin_expect (rcls == QUADFP_ZERO, 0)) | 72 else if (__glibc_unlikely (rcls == QUADFP_ZERO)) |
79 { | 73 { |
80 __float128 r; | 74 __float128 r; |
81 if (fabsq (__imag__ x) >= 2.0Q * FLT128_MIN) | 75 if (fabsq (__imag__ x) >= 2 * FLT128_MIN) |
82 r = sqrtq (0.5Q * fabsq (__imag__ x)); | 76 r = sqrtq (0.5Q * fabsq (__imag__ x)); |
83 else | 77 else |
84 r = 0.5Q * sqrtq (2.0Q * fabsq (__imag__ x)); | 78 r = 0.5Q * sqrtq (2 * fabsq (__imag__ x)); |
85 | 79 |
86 __real__ res = r; | 80 __real__ res = r; |
87 __imag__ res = copysignq (r, __imag__ x); | 81 __imag__ res = copysignq (r, __imag__ x); |
88 } | 82 } |
89 else | 83 else |
90 { | 84 { |
91 __float128 d, r, s; | 85 __float128 d, r, s; |
92 int scale = 0; | 86 int scale = 0; |
93 | 87 |
94 if (fabsq (__real__ x) > FLT128_MAX / 4.0Q) | 88 if (fabsq (__real__ x) > FLT128_MAX / 4) |
95 { | 89 { |
96 scale = 1; | 90 scale = 1; |
97 __real__ x = scalbnq (__real__ x, -2 * scale); | 91 __real__ x = scalbnq (__real__ x, -2 * scale); |
98 __imag__ x = scalbnq (__imag__ x, -2 * scale); | 92 __imag__ x = scalbnq (__imag__ x, -2 * scale); |
99 } | 93 } |
100 else if (fabsq (__imag__ x) > FLT128_MAX / 4.0Q) | 94 else if (fabsq (__imag__ x) > FLT128_MAX / 4) |
101 { | 95 { |
102 scale = 1; | 96 scale = 1; |
103 if (fabsq (__real__ x) >= 4.0Q * FLT128_MIN) | 97 if (fabsq (__real__ x) >= 4 * FLT128_MIN) |
104 __real__ x = scalbnq (__real__ x, -2 * scale); | 98 __real__ x = scalbnq (__real__ x, -2 * scale); |
105 else | 99 else |
106 __real__ x = 0.0Q; | 100 __real__ x = 0; |
107 __imag__ x = scalbnq (__imag__ x, -2 * scale); | 101 __imag__ x = scalbnq (__imag__ x, -2 * scale); |
108 } | 102 } |
109 else if (fabsq (__real__ x) < FLT128_MIN | 103 else if (fabsq (__real__ x) < 2 * FLT128_MIN |
110 && fabsq (__imag__ x) < FLT128_MIN) | 104 && fabsq (__imag__ x) < 2 * FLT128_MIN) |
111 { | 105 { |
112 scale = -(FLT128_MANT_DIG / 2); | 106 scale = -((FLT128_MANT_DIG + 1) / 2); |
113 __real__ x = scalbnq (__real__ x, -2 * scale); | 107 __real__ x = scalbnq (__real__ x, -2 * scale); |
114 __imag__ x = scalbnq (__imag__ x, -2 * scale); | 108 __imag__ x = scalbnq (__imag__ x, -2 * scale); |
115 } | 109 } |
116 | 110 |
117 d = hypotq (__real__ x, __imag__ x); | 111 d = hypotq (__real__ x, __imag__ x); |
118 /* Use the identity 2 Re res Im res = Im x | 112 /* Use the identity 2 Re res Im res = Im x |
119 to avoid cancellation error in d +/- Re x. */ | 113 to avoid cancellation error in d +/- Re x. */ |
120 if (__real__ x > 0) | 114 if (__real__ x > 0) |
121 { | 115 { |
122 r = sqrtq (0.5Q * (d + __real__ x)); | 116 r = sqrtq (0.5Q * (d + __real__ x)); |
123 s = 0.5Q * (__imag__ x / r); | 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); | |
124 } | 126 } |
125 else | 127 else |
126 { | 128 { |
127 s = sqrtq (0.5Q * (d - __real__ x)); | 129 s = sqrtq (0.5Q * (d - __real__ x)); |
128 r = fabsq (0.5Q * (__imag__ x / s)); | 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)); | |
129 } | 139 } |
130 | 140 |
131 if (scale) | 141 if (scale) |
132 { | 142 { |
133 r = scalbnq (r, scale); | 143 r = scalbnq (r, scale); |
134 s = scalbnq (s, scale); | 144 s = scalbnq (s, scale); |
135 } | 145 } |
136 | 146 |
147 math_check_force_underflow (r); | |
148 math_check_force_underflow (s); | |
149 | |
137 __real__ res = r; | 150 __real__ res = r; |
138 __imag__ res = copysignq (s, __imag__ x); | 151 __imag__ res = copysignq (s, __imag__ x); |
139 } | 152 } |
140 } | 153 } |
141 | 154 |