Home | History | Annotate | Line # | Download | only in math 
      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 (at) 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 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