diff options
Diffstat (limited to 'missing/acosh.c')
-rw-r--r-- | missing/acosh.c | 88 |
1 files changed, 88 insertions, 0 deletions
diff --git a/missing/acosh.c b/missing/acosh.c new file mode 100644 index 000000000..6a9daffa5 --- /dev/null +++ b/missing/acosh.c @@ -0,0 +1,88 @@ +/********************************************************************** + + acosh.c - + + $Author$ + $Date$ + created at: Fri Apr 12 00:34:17 JST 2002 + + public domain rewrite of acosh(3), asinh(3) and atanh(3) + +**********************************************************************/ + +#include <errno.h> +#include <float.h> +#include <math.h> + +/* DBL_MANT_DIG must be less than 4 times of bits of int */ +#ifdef DBL_MANT_DIG +#define DBL_MANT_DIG 53 /* in this case, at least 12 digit precision */ +#endif +#define BIG_CRITERIA_BIT (1<<DBL_MANT_DIG/2) +#if BIG_CRITERIA_BIT > 0 +#define BIG_CRITERIA (1.0*BIG_CRITERIA_BIT) +#else +#define BIG_CRITERIA (1.0*(1<<DBL_MANT_DIG/4)*(1<<(DBL_MANT_DIG/2+1-DBL_MANT_DIG/4))) +#endif +#define SMALL_CRITERIA_BIT (1<<(DBL_MANT_DIG/3)) +#if SMALL_CRITERIA_BIT > 0 +#define SMALL_CRITERIA (1.0/SMALL_CRITERIA_BIT) +#else +#define SMALL_CRITERIA (1.0*(1<<DBL_MANT_DIG/4)*(1<<(DBL_MANT_DIG/3+1-DBL_MANT_DIG/4))) +#endif + +#ifndef HAVE_ACOSH +double +acosh(x) + double x; +{ + if (x < 1) + x = -1; /* NaN */ + else if (x == 1) + return 0; + else if (x > BIG_CRITERIA) + x += x; + else + x += sqrt((x + 1) * (x - 1)); + return log(x); +} +#endif + +#ifndef HAVE_ASINH +double +asinh(x) + double x; +{ + int neg = x < 0; + double z = fabs(x); + + if (z < SMALL_CRITERIA) return x; + if (z < (1.0/(1<<DBL_MANT_DIG/5))) { + double x2 = z * z; + z *= 1 + x2 * (-1.0/6.0 + x2 * 3.0/40.0); + } + else if (z > BIG_CRITERIA) { + z = log(z + z); + } + else { + z = log(z + sqrt(z * z + 1)); + } + if (neg) z = -z; + return z; +} +#endif + +#ifndef HAVE_ATANH +double +atanh(x) + double x; +{ + int neg = x < 0; + double z = fabs(x); + + if (z < SMALL_CRITERIA) return x; + z = log(z > 1 ? -1 : (1 + z) / (1 - z)) / 2; + if (neg) z = -z; + return z; +} +#endif |