1 /* sys/invhyp.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include "gsl__config.h"
21 #include <math.h>
22 #include "gsl_math.h"
23 #include "gsl_machine.h"
24
25 double
gsl_acosh(const double x)26 gsl_acosh (const double x)
27 {
28 if (x > 1.0 / GSL_SQRT_DBL_EPSILON)
29 {
30 return log (x) + M_LN2;
31 }
32 else if (x > 2)
33 {
34 return log (2 * x - 1 / (sqrt (x * x - 1) + x));
35 }
36 else if (x > 1)
37 {
38 double t = x - 1;
39 return log1p (t + sqrt (2 * t + t * t));
40 }
41 else if (x == 1)
42 {
43 return 0;
44 }
45 else
46 {
47 return GSL_NAN;
48 }
49 }
50
51 double
gsl_asinh(const double x)52 gsl_asinh (const double x)
53 {
54 double a = fabs (x);
55 double s = (x < 0) ? -1 : 1;
56
57 if (a > 1 / GSL_SQRT_DBL_EPSILON)
58 {
59 return s * (log (a) + M_LN2);
60 }
61 else if (a > 2)
62 {
63 return s * log (2 * a + 1 / (a + sqrt (a * a + 1)));
64 }
65 else if (a > GSL_SQRT_DBL_EPSILON)
66 {
67 double a2 = a * a;
68 return s * log1p (a + a2 / (1 + sqrt (1 + a2)));
69 }
70 else
71 {
72 return x;
73 }
74 }
75
76 double
gsl_atanh(const double x)77 gsl_atanh (const double x)
78 {
79 double a = fabs (x);
80 double s = (x < 0) ? -1 : 1;
81
82 if (a > 1)
83 {
84 return GSL_NAN;
85 }
86 else if (a == 1)
87 {
88 return (x < 0) ? GSL_NEGINF : GSL_POSINF;
89 }
90 else if (a >= 0.5)
91 {
92 return s * 0.5 * log1p (2 * a / (1 - a));
93 }
94 else if (a > GSL_DBL_EPSILON)
95 {
96 return s * 0.5 * log1p (2 * a + 2 * a * a / (1 - a));
97 }
98 else
99 {
100 return x;
101 }
102 }
103