]> git.ipfire.org Git - thirdparty/glibc.git/blobdiff - math/k_casinhf.c
Update copyright dates with scripts/update-copyrights.
[thirdparty/glibc.git] / math / k_casinhf.c
index 4d1a92f8138a9cf663940a8321607362121735db..04c1a21f253a64e10c46d31e5f8dfbcd50660f2c 100644 (file)
@@ -1,7 +1,7 @@
 /* Return arc hyperbole sine for float value, with the imaginary part
    of the result possibly adjusted for use in computing other
    functions.
-   Copyright (C) 1997-2013 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
@@ -77,6 +77,39 @@ __kernel_casinhf (__complex__ float x, int adj)
       else
        __imag__ res = __ieee754_atan2f (s, rx);
     }
+  else if (ix > 1.0f && ix < 1.5f && rx < 0.5f)
+    {
+      if (rx < FLT_EPSILON * FLT_EPSILON)
+       {
+         float ix2m1 = (ix + 1.0f) * (ix - 1.0f);
+         float s = __ieee754_sqrtf (ix2m1);
+
+         __real__ res = __log1pf (2.0f * (ix2m1 + ix * s)) / 2.0f;
+         if (adj)
+           __imag__ res = __ieee754_atan2f (rx, __copysignf (s, __imag__ x));
+         else
+           __imag__ res = __ieee754_atan2f (s, rx);
+       }
+      else
+       {
+         float ix2m1 = (ix + 1.0f) * (ix - 1.0f);
+         float rx2 = rx * rx;
+         float f = rx2 * (2.0f + rx2 + 2.0f * ix * ix);
+         float d = __ieee754_sqrtf (ix2m1 * ix2m1 + f);
+         float dp = d + ix2m1;
+         float dm = f / dp;
+         float r1 = __ieee754_sqrtf ((dm + rx2) / 2.0f);
+         float r2 = rx * ix / r1;
+
+         __real__ res
+           = __log1pf (rx2 + dp + 2.0f * (rx * r1 + ix * r2)) / 2.0f;
+         if (adj)
+           __imag__ res = __ieee754_atan2f (rx + r1, __copysignf (ix + r2,
+                                                                  __imag__ x));
+         else
+           __imag__ res = __ieee754_atan2f (ix + r2, rx + r1);
+       }
+    }
   else if (ix == 1.0f && rx < 0.5f)
     {
       if (rx < FLT_EPSILON / 8.0f)
@@ -103,6 +136,58 @@ __kernel_casinhf (__complex__ float x, int adj)
            __imag__ res = __ieee754_atan2f (1.0f + s2, rx + s1);
        }
     }
+  else if (ix < 1.0f && rx < 0.5f)
+    {
+      if (ix >= FLT_EPSILON)
+       {
+         if (rx < FLT_EPSILON * FLT_EPSILON)
+           {
+             float onemix2 = (1.0f + ix) * (1.0f - ix);
+             float s = __ieee754_sqrtf (onemix2);
+
+             __real__ res = __log1pf (2.0f * rx / s) / 2.0f;
+             if (adj)
+               __imag__ res = __ieee754_atan2f (s, __imag__ x);
+             else
+               __imag__ res = __ieee754_atan2f (ix, s);
+           }
+         else
+           {
+             float onemix2 = (1.0f + ix) * (1.0f - ix);
+             float rx2 = rx * rx;
+             float f = rx2 * (2.0f + rx2 + 2.0f * ix * ix);
+             float d = __ieee754_sqrtf (onemix2 * onemix2 + f);
+             float dp = d + onemix2;
+             float dm = f / dp;
+             float r1 = __ieee754_sqrtf ((dp + rx2) / 2.0f);
+             float r2 = rx * ix / r1;
+
+             __real__ res
+               = __log1pf (rx2 + dm + 2.0f * (rx * r1 + ix * r2)) / 2.0f;
+             if (adj)
+               __imag__ res = __ieee754_atan2f (rx + r1,
+                                                __copysignf (ix + r2,
+                                                             __imag__ x));
+             else
+               __imag__ res = __ieee754_atan2f (ix + r2, rx + r1);
+           }
+       }
+      else
+       {
+         float s = __ieee754_hypotf (1.0f, rx);
+
+         __real__ res = __log1pf (2.0f * rx * (rx + s)) / 2.0f;
+         if (adj)
+           __imag__ res = __ieee754_atan2f (s, __imag__ x);
+         else
+           __imag__ res = __ieee754_atan2f (ix, s);
+       }
+      if (__real__ res < FLT_MIN)
+       {
+         volatile float force_underflow = __real__ res * __real__ res;
+         (void) force_underflow;
+       }
+    }
   else
     {
       __real__ y = (rx - ix) * (rx + ix) + 1.0f;