1 | /* Compute complex base 10 logarithm for complex __float128. |
---|---|

2 | Copyright (C) 1997-2012 Free Software Foundation, Inc. |

3 | This file is part of the GNU C Library. |

4 | Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. |

5 | |

6 | The GNU C Library is free software; you can redistribute it and/or |

7 | modify it under the terms of the GNU Lesser General Public |

8 | License as published by the Free Software Foundation; either |

9 | version 2.1 of the License, or (at your option) any later version. |

10 | |

11 | The GNU C Library is distributed in the hope that it will be useful, |

12 | but WITHOUT ANY WARRANTY; without even the implied warranty of |

13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

14 | Lesser General Public License for more details. |

15 | |

16 | You should have received a copy of the GNU Lesser General Public |

17 | License along with the GNU C Library; if not, see |

18 | <http://www.gnu.org/licenses/>. */ |

19 | |

20 | #include "quadmath-imp.h" |

21 | |

22 | |

23 | /* log_10 (2). */ |

24 | #define M_LOG10_2q 0.3010299956639811952137388947244930267682Q |

25 | |

26 | |

27 | __complex128 |

28 | clog10q (__complex128 x) |

29 | { |

30 | __complex128 result; |

31 | int rcls = fpclassifyq (__real__ x); |

32 | int icls = fpclassifyq (__imag__ x); |

33 | |

34 | if (__builtin_expect (rcls == QUADFP_ZERO && icls == QUADFP_ZERO, 0)) |

35 | { |

36 | /* Real and imaginary part are 0.0. */ |

37 | __imag__ result = signbitq (__real__ x) ? M_PIq : 0.0Q; |

38 | __imag__ result = copysignq (__imag__ result, __imag__ x); |

39 | /* Yes, the following line raises an exception. */ |

40 | __real__ result = -1.0Q / fabsq (__real__ x); |

41 | } |

42 | else if (__builtin_expect (rcls != QUADFP_NAN && icls != QUADFP_NAN, 1)) |

43 | { |

44 | /* Neither real nor imaginary part is NaN. */ |

45 | __float128 absx = fabsq (__real__ x), absy = fabsq (__imag__ x); |

46 | int scale = 0; |

47 | |

48 | if (absx < absy) |

49 | { |

50 | __float128 t = absx; |

51 | absx = absy; |

52 | absy = t; |

53 | } |

54 | |

55 | if (absx > FLT128_MAX / 2.0Q) |

56 | { |

57 | scale = -1; |

58 | absx = scalbnq (absx, scale); |

59 | absy = (absy >= FLT128_MIN * 2.0Q ? scalbnq (absy, scale) : 0.0Q); |

60 | } |

61 | else if (absx < FLT128_MIN && absy < FLT128_MIN) |

62 | { |

63 | scale = FLT128_MANT_DIG; |

64 | absx = scalbnq (absx, scale); |

65 | absy = scalbnq (absy, scale); |

66 | } |

67 | |

68 | if (absx == 1.0Q && scale == 0) |

69 | { |

70 | __float128 absy2 = absy * absy; |

71 | if (absy2 <= FLT128_MIN * 2.0Q * M_LN10q) |

72 | __real__ result |

73 | = (absy2 / 2.0Q - absy2 * absy2 / 4.0Q) * M_LOG10Eq; |

74 | else |

75 | __real__ result = log1pq (absy2) * (M_LOG10Eq / 2.0Q); |

76 | } |

77 | else if (absx > 1.0Q && absx < 2.0Q && absy < 1.0Q && scale == 0) |

78 | { |

79 | __float128 d2m1 = (absx - 1.0Q) * (absx + 1.0Q); |

80 | if (absy >= FLT128_EPSILON) |

81 | d2m1 += absy * absy; |

82 | __real__ result = log1pq (d2m1) * (M_LOG10Eq / 2.0Q); |

83 | } |

84 | else if (absx < 1.0Q |

85 | && absx >= 0.75Q |

86 | && absy < FLT128_EPSILON / 2.0Q |

87 | && scale == 0) |

88 | { |

89 | __float128 d2m1 = (absx - 1.0Q) * (absx + 1.0Q); |

90 | __real__ result = log1pq (d2m1) * (M_LOG10Eq / 2.0Q); |

91 | } |

92 | else if (absx < 1.0Q && (absx >= 0.75Q || absy >= 0.5Q) && scale == 0) |

93 | { |

94 | __float128 d2m1 = __quadmath_x2y2m1q (absx, absy); |

95 | __real__ result = log1pq (d2m1) * (M_LOG10Eq / 2.0Q); |

96 | } |

97 | else |

98 | { |

99 | __float128 d = hypotq (absx, absy); |

100 | __real__ result = log10q (d) - scale * M_LOG10_2q; |

101 | } |

102 | |

103 | __imag__ result = M_LOG10Eq * atan2q (__imag__ x, __real__ x); |

104 | } |

105 | else |

106 | { |

107 | __imag__ result = nanq (""); |

108 | if (rcls == QUADFP_INFINITE || icls == QUADFP_INFINITE) |

109 | /* Real or imaginary part is infinite. */ |

110 | __real__ result = HUGE_VALQ; |

111 | else |

112 | __real__ result = nanq (""); |

113 | } |

114 | |

115 | return result; |

116 | } |

117 |