math/k_casinhl.c

   1 /* Return arc hyperbole sine for long double value, with the imaginary
   2    part of the result possibly adjusted for use in computing other
   3    functions.
   4    Copyright (C) 1997-2013 Free Software Foundation, Inc.
   5    This file is part of the GNU C Library.
   6
   7    The GNU C Library is free software; you can redistribute it and/or
   8    modify it under the terms of the GNU Lesser General Public
   9    License as published by the Free Software Foundation; either
  10    version 2.1 of the License, or (at your option) any later version.
  11
  12    The GNU C Library is distributed in the hope that it will be useful,
  13    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15    Lesser General Public License for more details.
  16
  17    You should have received a copy of the GNU Lesser General Public
  18    License along with the GNU C Library; if not, see
  19    <http://www.gnu.org/licenses/>.  */
  20
  21 #include <complex.h>
  22 #include <math.h>
  23 #include <math_private.h>
  24 #include <float.h>
  25
  26 /* To avoid spurious overflows, use this definition to treat IBM long
  27    double as approximating an IEEE-style format.  */
  28 #if LDBL_MANT_DIG == 106
  29 # undef LDBL_EPSILON
  30 # define LDBL_EPSILON 0x1p-106L
  31 #endif
  32
  33 /* Return the complex inverse hyperbolic sine of finite nonzero Z,
  34    with the imaginary part of the result subtracted from pi/2 if ADJ
  35    is nonzero.  */
  36
  37 __complex__ long double
  38 __kernel_casinhl (__complex__ long double x, int adj)
  39 {
  40   __complex__ long double res;
  41   long double rx, ix;
  42   __complex__ long double y;
  43
  44   /* Avoid cancellation by reducing to the first quadrant.  */
  45   rx = fabsl (__real__ x);
  46   ix = fabsl (__imag__ x);
  47
  48   if (rx >= 1.0L / LDBL_EPSILON || ix >= 1.0L / LDBL_EPSILON)
  49     {
  50       /* For large x in the first quadrant, x + csqrt (1 + x * x)
  51          is sufficiently close to 2 * x to make no significant
  52          difference to the result; avoid possible overflow from
  53          the squaring and addition.  */
  54       __real__ y = rx;
  55       __imag__ y = ix;
  56
  57       if (adj)
  58         {
  59           long double t = __real__ y;
  60           __real__ y = __copysignl (__imag__ y, __imag__ x);
  61           __imag__ y = t;
  62         }
  63
  64       res = __clogl (y);
  65       __real__ res += M_LN2l;
  66     }
  67   else if (rx >= 0.5L && ix < LDBL_EPSILON / 8.0L)
  68     {
  69       long double s = __ieee754_hypotl (1.0L, rx);
  70
  71       __real__ res = __ieee754_logl (rx + s);
  72       if (adj)
  73         __imag__ res = __ieee754_atan2l (s, __imag__ x);
  74       else
  75         __imag__ res = __ieee754_atan2l (ix, s);
  76     }
  77   else if (rx < LDBL_EPSILON / 8.0L && ix >= 1.5L)
  78     {
  79       long double s = __ieee754_sqrtl ((ix + 1.0L) * (ix - 1.0L));
  80
  81       __real__ res = __ieee754_logl (ix + s);
  82       if (adj)
  83         __imag__ res = __ieee754_atan2l (rx, __copysignl (s, __imag__ x));
  84       else
  85         __imag__ res = __ieee754_atan2l (s, rx);
  86     }
  87   else
  88     {
  89       __real__ y = (rx - ix) * (rx + ix) + 1.0;
  90       __imag__ y = 2.0 * rx * ix;
  91
  92       y = __csqrtl (y);
  93
  94       __real__ y += rx;
  95       __imag__ y += ix;
  96
  97       if (adj)
  98         {
  99           long double t = __real__ y;
 100           __real__ y = __copysignl (__imag__ y, __imag__ x);
 101           __imag__ y = t;
 102         }
 103
 104       res = __clogl (y);
 105     }
 106
 107   /* Give results the correct sign for the original argument.  */
 108   __real__ res = __copysignl (__real__ res, __real__ x);
 109   __imag__ res = __copysignl (__imag__ res, (adj ? 1.0L : __imag__ x));
 110
 111   return res;
 112 }