1/* Compute complex base 10 logarithm.
2 Copyright (C) 1997-2020 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 <https://www.gnu.org/licenses/>. */
19
20#include <complex.h>
21#include <math.h>
22#include <math_private.h>
23#include <math-underflow.h>
24#include <float.h>
25
26/* log_10 (2). */
27#define LOG10_2 M_LIT (0.3010299956639811952137388947244930267682)
28
29/* pi * log10 (e). */
30#define PI_LOG10E M_LIT (1.364376353841841347485783625431355770210)
31
32CFLOAT
33M_DECL_FUNC (__clog10) (CFLOAT x)
34{
35 CFLOAT result;
36 int rcls = fpclassify (__real__ x);
37 int icls = fpclassify (__imag__ x);
38
39 if (__glibc_unlikely (rcls == FP_ZERO && icls == FP_ZERO))
40 {
41 /* Real and imaginary part are 0.0. */
42 __imag__ result = signbit (__real__ x) ? PI_LOG10E : 0;
43 __imag__ result = M_COPYSIGN (__imag__ result, __imag__ x);
44 /* Yes, the following line raises an exception. */
45 __real__ result = -1 / M_FABS (__real__ x);
46 }
47 else if (__glibc_likely (rcls != FP_NAN && icls != FP_NAN))
48 {
49 /* Neither real nor imaginary part is NaN. */
50 FLOAT absx = M_FABS (__real__ x), absy = M_FABS (__imag__ x);
51 int scale = 0;
52
53 if (absx < absy)
54 {
55 FLOAT t = absx;
56 absx = absy;
57 absy = t;
58 }
59
60 if (absx > M_MAX / 2)
61 {
62 scale = -1;
63 absx = M_SCALBN (absx, scale);
64 absy = (absy >= M_MIN * 2 ? M_SCALBN (absy, scale) : 0);
65 }
66 else if (absx < M_MIN && absy < M_MIN)
67 {
68 scale = M_MANT_DIG;
69 absx = M_SCALBN (absx, scale);
70 absy = M_SCALBN (absy, scale);
71 }
72
73 if (absx == 1 && scale == 0)
74 {
75 __real__ result = (M_LOG1P (absy * absy)
76 * ((FLOAT) M_MLIT (M_LOG10E) / 2));
77 math_check_force_underflow_nonneg (__real__ result);
78 }
79 else if (absx > 1 && absx < 2 && absy < 1 && scale == 0)
80 {
81 FLOAT d2m1 = (absx - 1) * (absx + 1);
82 if (absy >= M_EPSILON)
83 d2m1 += absy * absy;
84 __real__ result = M_LOG1P (d2m1) * ((FLOAT) M_MLIT (M_LOG10E) / 2);
85 }
86 else if (absx < 1
87 && absx >= M_LIT (0.5)
88 && absy < M_EPSILON / 2
89 && scale == 0)
90 {
91 FLOAT d2m1 = (absx - 1) * (absx + 1);
92 __real__ result = M_LOG1P (d2m1) * ((FLOAT) M_MLIT (M_LOG10E) / 2);
93 }
94 else if (absx < 1
95 && absx >= M_LIT (0.5)
96 && scale == 0
97 && absx * absx + absy * absy >= M_LIT (0.5))
98 {
99 FLOAT d2m1 = M_SUF (__x2y2m1) (absx, absy);
100 __real__ result = M_LOG1P (d2m1) * ((FLOAT) M_MLIT (M_LOG10E) / 2);
101 }
102 else
103 {
104 FLOAT d = M_HYPOT (absx, absy);
105 __real__ result = M_SUF (__ieee754_log10) (d) - scale * LOG10_2;
106 }
107
108 __imag__ result = M_MLIT (M_LOG10E) * M_ATAN2 (__imag__ x, __real__ x);
109 }
110 else
111 {
112 __imag__ result = M_NAN;
113 if (rcls == FP_INFINITE || icls == FP_INFINITE)
114 /* Real or imaginary part is infinite. */
115 __real__ result = M_HUGE_VAL;
116 else
117 __real__ result = M_NAN;
118 }
119
120 return result;
121}
122
123declare_mgen_alias (__clog10, clog10)
124