math/k_casinhf.c

   1 /* Return arc hyperbole sine for float value, with the imaginary part
   2    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 /* Return the complex inverse hyperbolic sine of finite nonzero Z,
  27    with the imaginary part of the result subtracted from pi/2 if ADJ
  28    is nonzero.  */
  29
  30 __complex__ float
  31 __kernel_casinhf (__complex__ float x, int adj)
  32 {
  33   __complex__ float res;
  34   float rx, ix;
  35   __complex__ float y;
  36
  37   /* Avoid cancellation by reducing to the first quadrant.  */
  38   rx = fabsf (__real__ x);
  39   ix = fabsf (__imag__ x);
  40
  41   if (rx >= 1.0f / FLT_EPSILON || ix >= 1.0f / FLT_EPSILON)
  42     {
  43       /* For large x in the first quadrant, x + csqrt (1 + x * x)
  44          is sufficiently close to 2 * x to make no significant
  45          difference to the result; avoid possible overflow from
  46          the squaring and addition.  */
  47       __real__ y = rx;
  48       __imag__ y = ix;
  49
  50       if (adj)
  51         {
  52           float t = __real__ y;
  53           __real__ y = __copysignf (__imag__ y, __imag__ x);
  54           __imag__ y = t;
  55         }
  56
  57       res = __clogf (y);
  58       __real__ res += (float) M_LN2;
  59     }
  60   else if (rx >= 0.5f && ix < FLT_EPSILON / 8.0f)
  61     {
  62       float s = __ieee754_hypotf (1.0f, rx);
  63
  64       __real__ res = __ieee754_logf (rx + s);
  65       if (adj)
  66         __imag__ res = __ieee754_atan2f (s, __imag__ x);
  67       else
  68         __imag__ res = __ieee754_atan2f (ix, s);
  69     }
  70   else if (rx < FLT_EPSILON / 8.0f && ix >= 1.5f)
  71     {
  72       float s = __ieee754_sqrtf ((ix + 1.0f) * (ix - 1.0f));
  73
  74       __real__ res = __ieee754_logf (ix + s);
  75       if (adj)
  76         __imag__ res = __ieee754_atan2f (rx, __copysignf (s, __imag__ x));
  77       else
  78         __imag__ res = __ieee754_atan2f (s, rx);
  79     }
  80   else if (ix > 1.0f && ix < 1.5f && rx < 0.5f)
  81     {
  82       if (rx < FLT_EPSILON * FLT_EPSILON)
  83         {
  84           float ix2m1 = (ix + 1.0f) * (ix - 1.0f);
  85           float s = __ieee754_sqrtf (ix2m1);
  86
  87           __real__ res = __log1pf (2.0f * (ix2m1 + ix * s)) / 2.0f;
  88           if (adj)
  89             __imag__ res = __ieee754_atan2f (rx, __copysignf (s, __imag__ x));
  90           else
  91             __imag__ res = __ieee754_atan2f (s, rx);
  92         }
  93       else
  94         {
  95           float ix2m1 = (ix + 1.0f) * (ix - 1.0f);
  96           float rx2 = rx * rx;
  97           float f = rx2 * (2.0f + rx2 + 2.0f * ix * ix);
  98           float d = __ieee754_sqrtf (ix2m1 * ix2m1 + f);
  99           float dp = d + ix2m1;
 100           float dm = f / dp;
 101           float r1 = __ieee754_sqrtf ((dm + rx2) / 2.0f);
 102           float r2 = rx * ix / r1;
 103
 104           __real__ res
 105             = __log1pf (rx2 + dp + 2.0f * (rx * r1 + ix * r2)) / 2.0f;
 106           if (adj)
 107             __imag__ res = __ieee754_atan2f (rx + r1, __copysignf (ix + r2,
 108                                                                    __imag__ x));
 109           else
 110             __imag__ res = __ieee754_atan2f (ix + r2, rx + r1);
 111         }
 112     }
 113   else if (ix == 1.0f && rx < 0.5f)
 114     {
 115       if (rx < FLT_EPSILON / 8.0f)
 116         {
 117           __real__ res = __log1pf (2.0f * (rx + __ieee754_sqrtf (rx))) / 2.0f;
 118           if (adj)
 119             __imag__ res = __ieee754_atan2f (__ieee754_sqrtf (rx),
 120                                              __copysignf (1.0f, __imag__ x));
 121           else
 122             __imag__ res = __ieee754_atan2f (1.0f, __ieee754_sqrtf (rx));
 123         }
 124       else
 125         {
 126           float d = rx * __ieee754_sqrtf (4.0f + rx * rx);
 127           float s1 = __ieee754_sqrtf ((d + rx * rx) / 2.0f);
 128           float s2 = __ieee754_sqrtf ((d - rx * rx) / 2.0f);
 129
 130           __real__ res = __log1pf (rx * rx + d + 2.0f * (rx * s1 + s2)) / 2.0f;
 131           if (adj)
 132             __imag__ res = __ieee754_atan2f (rx + s1,
 133                                              __copysignf (1.0f + s2,
 134                                                           __imag__ x));
 135           else
 136             __imag__ res = __ieee754_atan2f (1.0f + s2, rx + s1);
 137         }
 138     }
 139   else if (ix < 1.0f && rx < 0.5f)
 140     {
 141       if (ix >= FLT_EPSILON)
 142         {
 143           if (rx < FLT_EPSILON * FLT_EPSILON)
 144             {
 145               float onemix2 = (1.0f + ix) * (1.0f - ix);
 146               float s = __ieee754_sqrtf (onemix2);
 147
 148               __real__ res = __log1pf (2.0f * rx / s) / 2.0f;
 149               if (adj)
 150                 __imag__ res = __ieee754_atan2f (s, __imag__ x);
 151               else
 152                 __imag__ res = __ieee754_atan2f (ix, s);
 153             }
 154           else
 155             {
 156               float onemix2 = (1.0f + ix) * (1.0f - ix);
 157               float rx2 = rx * rx;
 158               float f = rx2 * (2.0f + rx2 + 2.0f * ix * ix);
 159               float d = __ieee754_sqrtf (onemix2 * onemix2 + f);
 160               float dp = d + onemix2;
 161               float dm = f / dp;
 162               float r1 = __ieee754_sqrtf ((dp + rx2) / 2.0f);
 163               float r2 = rx * ix / r1;
 164
 165               __real__ res
 166                 = __log1pf (rx2 + dm + 2.0f * (rx * r1 + ix * r2)) / 2.0f;
 167               if (adj)
 168                 __imag__ res = __ieee754_atan2f (rx + r1,
 169                                                  __copysignf (ix + r2,
 170                                                               __imag__ x));
 171               else
 172                 __imag__ res = __ieee754_atan2f (ix + r2, rx + r1);
 173             }
 174         }
 175       else
 176         {
 177           float s = __ieee754_hypotf (1.0f, rx);
 178
 179           __real__ res = __log1pf (2.0f * rx * (rx + s)) / 2.0f;
 180           if (adj)
 181             __imag__ res = __ieee754_atan2f (s, __imag__ x);
 182           else
 183             __imag__ res = __ieee754_atan2f (ix, s);
 184         }
 185       if (__real__ res < FLT_MIN)
 186         {
 187           volatile float force_underflow = __real__ res * __real__ res;
 188           (void) force_underflow;
 189         }
 190     }
 191   else
 192     {
 193       __real__ y = (rx - ix) * (rx + ix) + 1.0f;
 194       __imag__ y = 2.0f * rx * ix;
 195
 196       y = __csqrtf (y);
 197
 198       __real__ y += rx;
 199       __imag__ y += ix;
 200
 201       if (adj)
 202         {
 203           float t = __real__ y;
 204           __real__ y = __copysignf (__imag__ y, __imag__ x);
 205           __imag__ y = t;
 206         }
 207
 208       res = __clogf (y);
 209     }
 210
 211   /* Give results the correct sign for the original argument.  */
 212   __real__ res = __copysignf (__real__ res, __real__ x);
 213   __imag__ res = __copysignf (__imag__ res, (adj ? 1.0f : __imag__ x));
 214
 215   return res;
 216 }