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 | } |

144 |