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
  68     {
  69       __real__ y = (rx - ix) * (rx + ix) + 1.0;
  70       __imag__ y = 2.0 * rx * ix;
  71
  72       y = __csqrtl (y);
  73
  74       __real__ y += rx;
  75       __imag__ y += ix;
  76
  77       if (adj)
  78         {
  79           long double t = __real__ y;
  80           __real__ y = __copysignl (__imag__ y, __imag__ x);
  81           __imag__ y = t;
  82         }
  83
  84       res = __clogl (y);
  85     }
  86
  87   /* Give results the correct sign for the original argument.  */
  88   __real__ res = __copysignl (__real__ res, __real__ x);
  89   __imag__ res = __copysignl (__imag__ res, (adj ? 1.0L : __imag__ x));
  90
  91   return res;
  92 }