/* 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.3: March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/
#include "mconf.h"
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
};
#ifndef ANSIPROT
double log(), sqrt(), polevl(), p1evl();
#endif
extern double LOGE2;
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 )
{
if(!finite(x))
return(xx);
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) );
}