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