173 lines
3.3 KiB
C
173 lines
3.3 KiB
C
|
/* sinhl.c
|
|||
|
|
*
|
|||
|
|
* Hyperbolic sine, long double precision
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* SYNOPSIS:
|
|||
|
|
*
|
|||
|
|
* long double x, y, sinhl();
|
|||
|
|
*
|
|||
|
|
* y = sinhl( x );
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* DESCRIPTION:
|
|||
|
|
*
|
|||
|
|
* Returns hyperbolic sine of argument in the range MINLOGL to
|
|||
|
|
* MAXLOGL.
|
|||
|
|
*
|
|||
|
|
* The range is partitioned into two segments. If |x| <= 1, a
|
|||
|
|
* rational function of the form x + x**3 P(x)/Q(x) is employed.
|
|||
|
|
* Otherwise the calculation is sinh(x) = ( exp(x) - exp(-x) )/2.
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* ACCURACY:
|
|||
|
|
*
|
|||
|
|
* Relative error:
|
|||
|
|
* arithmetic domain # trials peak rms
|
|||
|
|
* IEEE -2,2 10000 1.5e-19 3.9e-20
|
|||
|
|
* IEEE +-10000 30000 1.1e-19 2.8e-20
|
|||
|
|
*
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
/*
|
|||
|
|
Cephes Math Library Release 2.7: January, 1998
|
|||
|
|
Copyright 1984, 1991, 1998 by Stephen L. Moshier
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
/*
|
|||
|
|
Modified for mingw
|
|||
|
|
2002-07-22 Danny Smith <dannysmith@users.sourceforge.net>
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
#ifdef __MINGW32__
|
|||
|
|
#include "cephes_mconf.h"
|
|||
|
|
#else
|
|||
|
|
#include "mconf.h"
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifndef _SET_ERRNO
|
|||
|
|
#define _SET_ERRNO(x)
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef UNK
|
|||
|
|
static long double P[] = {
|
|||
|
|
1.7550769032975377032681E-6L,
|
|||
|
|
4.1680702175874268714539E-4L,
|
|||
|
|
3.0993532520425419002409E-2L,
|
|||
|
|
9.9999999999999999998002E-1L,
|
|||
|
|
};
|
|||
|
|
static long double Q[] = {
|
|||
|
|
1.7453965448620151484660E-8L,
|
|||
|
|
-5.9116673682651952419571E-6L,
|
|||
|
|
1.0599252315677389339530E-3L,
|
|||
|
|
-1.1403880487744749056675E-1L,
|
|||
|
|
6.0000000000000000000200E0L,
|
|||
|
|
};
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef IBMPC
|
|||
|
|
static const unsigned short P[] = {
|
|||
|
|
0xec6a,0xd942,0xfbb3,0xeb8f,0x3feb, XPD
|
|||
|
|
0x365e,0xb30a,0xe437,0xda86,0x3ff3, XPD
|
|||
|
|
0x8890,0x01f6,0x2612,0xfde6,0x3ff9, XPD
|
|||
|
|
0x0000,0x0000,0x0000,0x8000,0x3fff, XPD
|
|||
|
|
};
|
|||
|
|
static const unsigned short Q[] = {
|
|||
|
|
0x4edd,0x4c21,0xad09,0x95ed,0x3fe5, XPD
|
|||
|
|
0x4376,0x9b70,0xd605,0xc65c,0xbfed, XPD
|
|||
|
|
0xc8ad,0x5d21,0x3069,0x8aed,0x3ff5, XPD
|
|||
|
|
0x9c32,0x6374,0x2d4b,0xe98d,0xbffb, XPD
|
|||
|
|
0x0000,0x0000,0x0000,0xc000,0x4001, XPD
|
|||
|
|
};
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef MIEEE
|
|||
|
|
static long P[] = {
|
|||
|
|
0x3feb0000,0xeb8ffbb3,0xd942ec6a,
|
|||
|
|
0x3ff30000,0xda86e437,0xb30a365e,
|
|||
|
|
0x3ff90000,0xfde62612,0x01f68890,
|
|||
|
|
0x3fff0000,0x80000000,0x00000000,
|
|||
|
|
};
|
|||
|
|
static long Q[] = {
|
|||
|
|
0x3fe50000,0x95edad09,0x4c214edd,
|
|||
|
|
0xbfed0000,0xc65cd605,0x9b704376,
|
|||
|
|
0x3ff50000,0x8aed3069,0x5d21c8ad,
|
|||
|
|
0xbffb0000,0xe98d2d4b,0x63749c32,
|
|||
|
|
0x40010000,0xc0000000,0x00000000,
|
|||
|
|
};
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifndef __MINGW32__
|
|||
|
|
extern long double MAXNUML, MAXLOGL, MINLOGL, LOGE2L;
|
|||
|
|
#ifdef ANSIPROT
|
|||
|
|
extern long double fabsl ( long double );
|
|||
|
|
extern long double expl ( long double );
|
|||
|
|
extern long double polevll ( long double, void *, int );
|
|||
|
|
extern long double p1evll ( long double, void *, int );
|
|||
|
|
#else
|
|||
|
|
long double fabsl(), expl(), polevll(), p1evll();
|
|||
|
|
#endif
|
|||
|
|
#ifdef INFINITIES
|
|||
|
|
extern long double INFINITYL;
|
|||
|
|
#endif
|
|||
|
|
#ifdef NANS
|
|||
|
|
extern long double NANL;
|
|||
|
|
#endif
|
|||
|
|
#endif /* __MINGW32__ */
|
|||
|
|
|
|||
|
|
long double sinhl(x)
|
|||
|
|
long double x;
|
|||
|
|
{
|
|||
|
|
long double a;
|
|||
|
|
|
|||
|
|
#ifdef MINUSZERO
|
|||
|
|
if( x == 0.0 )
|
|||
|
|
return(x);
|
|||
|
|
#endif
|
|||
|
|
#ifdef NANS
|
|||
|
|
if (isnanl(x))
|
|||
|
|
{
|
|||
|
|
_SET_ERRNO(EDOM);
|
|||
|
|
}
|
|||
|
|
#endif
|
|||
|
|
a = fabsl(x);
|
|||
|
|
if( (x > (MAXLOGL + LOGE2L)) || (x > -(MINLOGL-LOGE2L) ) )
|
|||
|
|
{
|
|||
|
|
mtherr( "sinhl", DOMAIN );
|
|||
|
|
_SET_ERRNO(ERANGE);
|
|||
|
|
#ifdef INFINITIES
|
|||
|
|
if( x > 0.0L )
|
|||
|
|
return( INFINITYL );
|
|||
|
|
else
|
|||
|
|
return( -INFINITYL );
|
|||
|
|
#else
|
|||
|
|
if( x > 0.0L )
|
|||
|
|
return( MAXNUML );
|
|||
|
|
else
|
|||
|
|
return( -MAXNUML );
|
|||
|
|
#endif
|
|||
|
|
}
|
|||
|
|
if( a > 1.0L )
|
|||
|
|
{
|
|||
|
|
if( a >= (MAXLOGL - LOGE2L) )
|
|||
|
|
{
|
|||
|
|
a = expl(0.5L*a);
|
|||
|
|
a = (0.5L * a) * a;
|
|||
|
|
if( x < 0.0L )
|
|||
|
|
a = -a;
|
|||
|
|
return(a);
|
|||
|
|
}
|
|||
|
|
a = expl(a);
|
|||
|
|
a = 0.5L*a - (0.5L/a);
|
|||
|
|
if( x < 0.0L )
|
|||
|
|
a = -a;
|
|||
|
|
return(a);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
a *= a;
|
|||
|
|
return( x + x * a * (polevll(a,P,3)/polevll(a,Q,4)) );
|
|||
|
|
}
|