/*							md_tanh.c
 *
 *	Hyperbolic tangent
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, md_tanh();
 *
 * y = md_tanh( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns hyperbolic tangent of argument in the range MINLOG to
 * MAXLOG.
 *
 * A rational function is used for |x| < 0.625.  The form
 * x + x**3 P(x)/Q(x) of Cody _& Waite is employed.
 * Otherwise,
 *    md_tanh(x) = md_sinh(x)/md_cosh(x) = 1  -  2/(md_exp(2x) + 1).
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC       -2,2        50000       3.3e-17     6.4e-18
 *    IEEE      -2,2        30000       2.5e-16     5.8e-17
 *
 */

/*
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/

#include "mconf.h"

#ifdef UNK
static double P[] = {
-9.64399179425052238628E-1,
-9.92877231001918586564E1,
-1.61468768441708447952E3
};
static double Q[] = {
/* 1.00000000000000000000E0,*/
 1.12811678491632931402E2,
 2.23548839060100448583E3,
 4.84406305325125486048E3
};
#endif
#ifdef DEC
static unsigned short P[] = {
0140166,0161335,0053753,0075126,
0141706,0111520,0070463,0040552,
0142711,0153001,0101300,0025430
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041741,0117624,0051300,0156060,
0043013,0133720,0071251,0127717,
0043227,0060201,0021020,0020136
};
#endif

#ifdef IBMPC
static unsigned short P[] = {
0x6f4b,0xaafd,0xdc5b,0xbfee,
0x682d,0x0e26,0xd26a,0xc058,
0x0563,0x3058,0x3ac0,0xc099
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x1b86,0x8a58,0x33f2,0x405c,
0x35fa,0x0e55,0x76fa,0x40a1,
0x040c,0x2442,0xec10,0x40b2
};
#endif

#ifdef MIEEE
static unsigned short P[] = {
0xbfee,0xdc5b,0xaafd,0x6f4b,
0xc058,0xd26a,0x0e26,0x682d,
0xc099,0x3ac0,0x3058,0x0563
};
static unsigned short Q[] = {
0x405c,0x33f2,0x8a58,0x1b86,
0x40a1,0x76fa,0x0e55,0x35fa,
0x40b2,0xec10,0x2442,0x040c
};
#endif

#ifdef ANSIPROT
extern double md_fabs ( double );
extern double md_exp ( double );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
#else
double md_fabs(), md_exp(), polevl(), p1evl();
#endif
extern double MAXLOG;

double md_tanh(x)
double x;
{
double s, z;

#ifdef MINUSZERO
if( x == 0.0 )
	return(x);
#endif
z = md_fabs(x);
if( z > 0.5 * MAXLOG )
	{
	if( x > 0 )
		return( 1.0 );
	else
		return( -1.0 );
	}
if( z >= 0.625 )
	{
	s = md_exp(2.0*z);
	z =  1.0  - 2.0/(s + 1.0);
	if( x < 0 )
		z = -z;
	}
else
	{
	if( x == 0.0 )
	  return(x);
	s = x * x;
	z = polevl( s, P, 2 )/p1evl(s, Q, 3);
	z = x * s * z;
	z = x + z;
	}
return( z );
}