/*							md_asinh.c
 *
 *	Inverse hyperbolic sine
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, md_asinh();
 *
 * y = md_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,
 *
 *     md_asinh(x) = md_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
 *
 */

/*						md_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 md_log ( double );
#else
double md_log(), sqrt(), polevl(), p1evl();
#endif
extern double LOGE2, INFINITY;

double md_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 * (md_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 * md_log(x + a) );
}