libquadmath/math/catanhq.c

   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
  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 }