| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133 |
- /* asinh.c
- *
- * Inverse hyperbolic sine
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, asinh();
- *
- * y = asinh( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns inverse hyperbolic sine of argument.
- *
- * If |x| < 0.5, the function is approximated by a rational
- * form x + x**3 P(x)/Q(x). Otherwise,
- *
- * asinh(x) = log( x + sqrt(1 + x*x) ).
- *
- *
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * DEC -3,3 75000 4.6e-17 1.1e-17
- * IEEE -1,1 30000 3.7e-16 7.8e-17
- * IEEE 1,3 30000 2.5e-16 6.7e-17
- *
- */
- /* asinh.c */
- /*
- Cephes Math Library Release 2.8: June, 2000
- Copyright 1984, 1995, 2000 by Stephen L. Moshier
- */
- #include "mconf.h"
- #ifdef UNK
- static double P[] = {-4.33231683752342103572E-3, -5.91750212056387121207E-1,
- -4.37390226194356683570E0, -9.09030533308377316566E0,
- -5.56682227230859640450E0};
- static double Q[] = {
- /* 1.00000000000000000000E0,*/
- 1.28757002067426453537E1, 4.86042483805291788324E1,
- 6.95722521337257608734E1, 3.34009336338516356383E1};
- #endif
- #ifdef DEC
- static unsigned short P[] = {0136215, 0173033, 0110410, 0105475, 0140027,
- 0076361, 0020056, 0164520, 0140613, 0173401,
- 0160136, 0053142, 0141021, 0070744, 0000503,
- 0176261, 0140662, 0021550, 0073106, 0133351};
- static unsigned short Q[] = {
- /* 0040200,0000000,0000000,0000000,*/
- 0041116, 0001336, 0034120, 0173054, 0041502, 0065300, 0013144, 0021231,
- 0041613, 0022376, 0035516, 0153063, 0041405, 0115216, 0054265, 0004557};
- #endif
- #ifdef IBMPC
- static unsigned short P[] = {0x1168, 0x7221, 0xbec3, 0xbf71, 0xdd2a,
- 0x2405, 0xef9e, 0xbfe2, 0xcacc, 0x3c0b,
- 0x7ee0, 0xc011, 0x7f96, 0x8028, 0x2e3c,
- 0xc022, 0xd6dd, 0x0ec8, 0x446d, 0xc016};
- static unsigned short Q[] = {
- /* 0x0000,0x0000,0x0000,0x3ff0,*/
- 0x1ec5, 0xc70a, 0xc05b, 0x4029, 0x8453, 0x02cc, 0x4d58, 0x4048,
- 0xdac6, 0xc769, 0x649f, 0x4051, 0xa12e, 0xcb16, 0xb351, 0x4040};
- #endif
- #ifdef MIEEE
- static unsigned short P[] = {0xbf71, 0xbec3, 0x7221, 0x1168, 0xbfe2,
- 0xef9e, 0x2405, 0xdd2a, 0xc011, 0x7ee0,
- 0x3c0b, 0xcacc, 0xc022, 0x2e3c, 0x8028,
- 0x7f96, 0xc016, 0x446d, 0x0ec8, 0xd6dd};
- static unsigned short Q[] = {0x4029, 0xc05b, 0xc70a, 0x1ec5, 0x4048, 0x4d58,
- 0x02cc, 0x8453, 0x4051, 0x649f, 0xc769, 0xdac6,
- 0x4040, 0xb351, 0xcb16, 0xa12e};
- #endif
- #ifdef ANSIPROT
- extern double polevl(double, void *, int);
- extern double p1evl(double, void *, int);
- extern double sqrt(double);
- extern double log(double);
- #else
- double log(), sqrt(), polevl(), p1evl();
- #endif
- extern double LOGE2, INFINITY;
- double asinh(xx) double xx;
- {
- double a, z, x;
- int sign;
- #ifdef MINUSZERO
- if (xx == 0.0)
- return (xx);
- #endif
- if (xx < 0.0) {
- sign = -1;
- x = -xx;
- } else {
- sign = 1;
- x = xx;
- }
- if (x > 1.0e8) {
- #ifdef INFINITIES
- if (x == INFINITY)
- return (xx);
- #endif
- return (sign * (log(x) + LOGE2));
- }
- z = x * x;
- if (x < 0.5) {
- a = (polevl(z, P, 4) / p1evl(z, Q, 4)) * z;
- a = a * x + x;
- if (sign < 0)
- a = -a;
- return (a);
- }
- a = sqrt(z + 1.0);
- return (sign * log(x + a));
- }
|