1 /* Compute complex natural logarithm. 2 Copyright (C) 1997-2018 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 __complex128 23 clogq (__complex128 x) 24 { 25 __complex128 result; 26 int rcls = fpclassifyq (__real__ x); 27 int icls = fpclassifyq (__imag__ x); 28 29 if (__glibc_unlikely (rcls == QUADFP_ZERO && icls == QUADFP_ZERO)) 30 { 31 /* Real and imaginary part are 0.0. */ 32 __imag__ result = signbitq (__real__ x) ? (__float128) M_PIq : 0; 33 __imag__ result = copysignq (__imag__ result, __imag__ x); 34 /* Yes, the following line raises an exception. */ 35 __real__ result = -1 / fabsq (__real__ x); 36 } 37 else if (__glibc_likely (rcls != QUADFP_NAN && icls != QUADFP_NAN)) 38 { 39 /* Neither real nor imaginary part is NaN. */ 40 __float128 absx = fabsq (__real__ x), absy = fabsq (__imag__ x); 41 int scale = 0; 42 43 if (absx < absy) 44 { 45 __float128 t = absx; 46 absx = absy; 47 absy = t; 48 } 49 50 if (absx > FLT128_MAX / 2) 51 { 52 scale = -1; 53 absx = scalbnq (absx, scale); 54 absy = (absy >= FLT128_MIN * 2 ? scalbnq (absy, scale) : 0); 55 } 56 else if (absx < FLT128_MIN && absy < FLT128_MIN) 57 { 58 scale = FLT128_MANT_DIG; 59 absx = scalbnq (absx, scale); 60 absy = scalbnq (absy, scale); 61 } 62 63 if (absx == 1 && scale == 0) 64 { 65 __real__ result = log1pq (absy * absy) / 2; 66 math_check_force_underflow_nonneg (__real__ result); 67 } 68 else if (absx > 1 && absx < 2 && absy < 1 && scale == 0) 69 { 70 __float128 d2m1 = (absx - 1) * (absx + 1); 71 if (absy >= FLT128_EPSILON) 72 d2m1 += absy * absy; 73 __real__ result = log1pq (d2m1) / 2; 74 } 75 else if (absx < 1 76 && absx >= 0.5Q 77 && absy < FLT128_EPSILON / 2 78 && scale == 0) 79 { 80 __float128 d2m1 = (absx - 1) * (absx + 1); 81 __real__ result = log1pq (d2m1) / 2; 82 } 83 else if (absx < 1 84 && absx >= 0.5Q 85 && scale == 0 86 && absx * absx + absy * absy >= 0.5Q) 87 { 88 __float128 d2m1 = __quadmath_x2y2m1q (absx, absy); 89 __real__ result = log1pq (d2m1) / 2; 90 } 91 else 92 { 93 __float128 d = hypotq (absx, absy); 94 __real__ result = logq (d) - scale * (__float128) M_LN2q; 95 } 96 97 __imag__ result = atan2q (__imag__ x, __real__ x); 98 } 99 else 100 { 101 __imag__ result = nanq (""); 102 if (rcls == QUADFP_INFINITE || icls == QUADFP_INFINITE) 103 /* Real or imaginary part is infinite. */ 104 __real__ result = HUGE_VALQ; 105 else 106 __real__ result = nanq (""); 107 } 108 109 return result; 110 } 111