1 /* Return arc hyperbolic tangent for a complex float type.
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
catanhq(__complex128 x)23 catanhq (__complex128 x)
24 {
25   __complex128 res;
26   int rcls = fpclassifyq (__real__ x);
27   int icls = fpclassifyq (__imag__ x);
28 
29   if (__glibc_unlikely (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE))
30     {
31       if (icls == QUADFP_INFINITE)
32 	{
33 	  __real__ res = copysignq (0, __real__ x);
34 	  __imag__ res = copysignq (M_PI_2q, __imag__ x);
35 	}
36       else if (rcls == QUADFP_INFINITE || rcls == QUADFP_ZERO)
37 	{
38 	  __real__ res = copysignq (0, __real__ x);
39 	  if (icls >= QUADFP_ZERO)
40 	    __imag__ res = copysignq (M_PI_2q, __imag__ x);
41 	  else
42 	    __imag__ res = nanq ("");
43 	}
44       else
45 	{
46 	  __real__ res = nanq ("");
47 	  __imag__ res = nanq ("");
48 	}
49     }
50   else if (__glibc_unlikely (rcls == QUADFP_ZERO && icls == QUADFP_ZERO))
51     {
52       res = x;
53     }
54   else
55     {
56       if (fabsq (__real__ x) >= 16 / FLT128_EPSILON
57 	  || fabsq (__imag__ x) >= 16 / FLT128_EPSILON)
58 	{
59 	  __imag__ res = copysignq (M_PI_2q, __imag__ x);
60 	  if (fabsq (__imag__ x) <= 1)
61 	    __real__ res = 1 / __real__ x;
62 	  else if (fabsq (__real__ x) <= 1)
63 	    __real__ res = __real__ x / __imag__ x / __imag__ x;
64 	  else
65 	    {
66 	      __float128 h = hypotq (__real__ x / 2, __imag__ x / 2);
67 	      __real__ res = __real__ x / h / h / 4;
68 	    }
69 	}
70       else
71 	{
72 	  if (fabsq (__real__ x) == 1
73 	      && fabsq (__imag__ x) < FLT128_EPSILON * FLT128_EPSILON)
74 	    __real__ res = (copysignq (0.5Q, __real__ x)
75 			    * ((__float128) M_LN2q
76 			       - logq (fabsq (__imag__ x))));
77 	  else
78 	    {
79 	      __float128 i2 = 0;
80 	      if (fabsq (__imag__ x) >= FLT128_EPSILON * FLT128_EPSILON)
81 		i2 = __imag__ x * __imag__ x;
82 
83 	      __float128 num = 1 + __real__ x;
84 	      num = i2 + num * num;
85 
86 	      __float128 den = 1 - __real__ x;
87 	      den = i2 + den * den;
88 
89 	      __float128 f = num / den;
90 	      if (f < 0.5Q)
91 		__real__ res = 0.25Q * logq (f);
92 	      else
93 		{
94 		  num = 4 * __real__ x;
95 		  __real__ res = 0.25Q * log1pq (num / den);
96 		}
97 	    }
98 
99 	  __float128 absx, absy, den;
100 
101 	  absx = fabsq (__real__ x);
102 	  absy = fabsq (__imag__ x);
103 	  if (absx < absy)
104 	    {
105 	      __float128 t = absx;
106 	      absx = absy;
107 	      absy = t;
108 	    }
109 
110 	  if (absy < FLT128_EPSILON / 2)
111 	    {
112 	      den = (1 - absx) * (1 + absx);
113 	      if (den == 0)
114 		den = 0;
115 	    }
116 	  else if (absx >= 1)
117 	    den = (1 - absx) * (1 + absx) - absy * absy;
118 	  else if (absx >= 0.75Q || absy >= 0.5Q)
119 	    den = -__quadmath_x2y2m1q (absx, absy);
120 	  else
121 	    den = (1 - absx) * (1 + absx) - absy * absy;
122 
123 	  __imag__ res = 0.5Q * atan2q (2 * __imag__ x, den);
124 	}
125 
126       math_check_force_underflow_complex (res);
127     }
128 
129   return res;
130 }
131