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-2014 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 if (ix > 1.0L && ix < 1.5L && rx < 0.5L)
  88     {
  89       if (rx < LDBL_EPSILON * LDBL_EPSILON)
  90         {
  91           long double ix2m1 = (ix + 1.0L) * (ix - 1.0L);
  92           long double s = __ieee754_sqrtl (ix2m1);
  93
  94           __real__ res = __log1pl (2.0L * (ix2m1 + ix * s)) / 2.0L;
  95           if (adj)
  96             __imag__ res = __ieee754_atan2l (rx, __copysignl (s, __imag__ x));
  97           else
  98             __imag__ res = __ieee754_atan2l (s, rx);
  99         }
 100       else
 101         {
 102           long double ix2m1 = (ix + 1.0L) * (ix - 1.0L);
 103           long double rx2 = rx * rx;
 104           long double f = rx2 * (2.0L + rx2 + 2.0L * ix * ix);
 105           long double d = __ieee754_sqrtl (ix2m1 * ix2m1 + f);
 106           long double dp = d + ix2m1;
 107           long double dm = f / dp;
 108           long double r1 = __ieee754_sqrtl ((dm + rx2) / 2.0L);
 109           long double r2 = rx * ix / r1;
 110
 111           __real__ res
 112             = __log1pl (rx2 + dp + 2.0L * (rx * r1 + ix * r2)) / 2.0L;
 113           if (adj)
 114             __imag__ res = __ieee754_atan2l (rx + r1, __copysignl (ix + r2,
 115                                                                    __imag__ x));
 116           else
 117             __imag__ res = __ieee754_atan2l (ix + r2, rx + r1);
 118         }
 119     }
 120   else if (ix == 1.0L && rx < 0.5L)
 121     {
 122       if (rx < LDBL_EPSILON / 8.0L)
 123         {
 124           __real__ res = __log1pl (2.0L * (rx + __ieee754_sqrtl (rx))) / 2.0L;
 125           if (adj)
 126             __imag__ res = __ieee754_atan2l (__ieee754_sqrtl (rx),
 127                                              __copysignl (1.0L, __imag__ x));
 128           else
 129             __imag__ res = __ieee754_atan2l (1.0L, __ieee754_sqrtl (rx));
 130         }
 131       else
 132         {
 133           long double d = rx * __ieee754_sqrtl (4.0L + rx * rx);
 134           long double s1 = __ieee754_sqrtl ((d + rx * rx) / 2.0L);
 135           long double s2 = __ieee754_sqrtl ((d - rx * rx) / 2.0L);
 136
 137           __real__ res = __log1pl (rx * rx + d + 2.0L * (rx * s1 + s2)) / 2.0L;
 138           if (adj)
 139             __imag__ res = __ieee754_atan2l (rx + s1,
 140                                              __copysignl (1.0L + s2,
 141                                                           __imag__ x));
 142           else
 143             __imag__ res = __ieee754_atan2l (1.0L + s2, rx + s1);
 144         }
 145     }
 146   else if (ix < 1.0L && rx < 0.5L)
 147     {
 148       if (ix >= LDBL_EPSILON)
 149         {
 150           if (rx < LDBL_EPSILON * LDBL_EPSILON)
 151             {
 152               long double onemix2 = (1.0L + ix) * (1.0L - ix);
 153               long double s = __ieee754_sqrtl (onemix2);
 154
 155               __real__ res = __log1pl (2.0L * rx / s) / 2.0L;
 156               if (adj)
 157                 __imag__ res = __ieee754_atan2l (s, __imag__ x);
 158               else
 159                 __imag__ res = __ieee754_atan2l (ix, s);
 160             }
 161           else
 162             {
 163               long double onemix2 = (1.0L + ix) * (1.0L - ix);
 164               long double rx2 = rx * rx;
 165               long double f = rx2 * (2.0L + rx2 + 2.0L * ix * ix);
 166               long double d = __ieee754_sqrtl (onemix2 * onemix2 + f);
 167               long double dp = d + onemix2;
 168               long double dm = f / dp;
 169               long double r1 = __ieee754_sqrtl ((dp + rx2) / 2.0L);
 170               long double r2 = rx * ix / r1;
 171
 172               __real__ res
 173                 = __log1pl (rx2 + dm + 2.0L * (rx * r1 + ix * r2)) / 2.0L;
 174               if (adj)
 175                 __imag__ res = __ieee754_atan2l (rx + r1,
 176                                                  __copysignl (ix + r2,
 177                                                               __imag__ x));
 178               else
 179                 __imag__ res = __ieee754_atan2l (ix + r2, rx + r1);
 180             }
 181         }
 182       else
 183         {
 184           long double s = __ieee754_hypotl (1.0L, rx);
 185
 186           __real__ res = __log1pl (2.0L * rx * (rx + s)) / 2.0L;
 187           if (adj)
 188             __imag__ res = __ieee754_atan2l (s, __imag__ x);
 189           else
 190             __imag__ res = __ieee754_atan2l (ix, s);
 191         }
 192       if (__real__ res < LDBL_MIN)
 193         {
 194           volatile long double force_underflow = __real__ res * __real__ res;
 195           (void) force_underflow;
 196         }
 197     }
 198   else
 199     {
 200       __real__ y = (rx - ix) * (rx + ix) + 1.0L;
 201       __imag__ y = 2.0L * rx * ix;
 202
 203       y = __csqrtl (y);
 204
 205       __real__ y += rx;
 206       __imag__ y += ix;
 207
 208       if (adj)
 209         {
 210           long double t = __real__ y;
 211           __real__ y = __copysignl (__imag__ y, __imag__ x);
 212           __imag__ y = t;
 213         }
 214
 215       res = __clogl (y);
 216     }
 217
 218   /* Give results the correct sign for the original argument.  */
 219   __real__ res = __copysignl (__real__ res, __real__ x);
 220   __imag__ res = __copysignl (__imag__ res, (adj ? 1.0L : __imag__ x));
 221
 222   return res;
 223 }