Kai Tietz | 815a664 | 2009-02-19 10:54:09 +0000 | [diff] [blame] | 1 | /** |
| 2 | * This file has no copyright assigned and is placed in the Public Domain. |
| 3 | * This file is part of the w64 mingw-runtime package. |
Kai Tietz | fa0cfe3 | 2010-01-15 20:02:21 +0000 | [diff] [blame] | 4 | * No warranty is given; refer to the file DISCLAIMER.PD within this package. |
Kai Tietz | 815a664 | 2009-02-19 10:54:09 +0000 | [diff] [blame] | 5 | */ |
| 6 | #include <math.h> |
| 7 | #include <errno.h> |
| 8 | #include "fastmath.h" |
| 9 | |
| 10 | /* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */ |
| 11 | double asinh(double x) |
| 12 | { |
| 13 | double z; |
| 14 | if (!isfinite (x)) |
| 15 | return x; |
| 16 | z = fabs (x); |
| 17 | |
| 18 | /* Avoid setting FPU underflow exception flag in x * x. */ |
| 19 | #if 0 |
| 20 | if ( z < 0x1p-32) |
| 21 | return x; |
| 22 | #endif |
| 23 | |
| 24 | /* Use log1p to avoid cancellation with small x. Put |
| 25 | x * x in denom, so overflow is harmless. |
| 26 | asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0) |
| 27 | = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0)) */ |
| 28 | |
| 29 | z = __fast_log1p (z + z * z / (__fast_sqrt (z * z + 1.0) + 1.0)); |
| 30 | |
| 31 | return ( x > 0.0 ? z : -z); |
| 32 | } |
| 33 | |