#include <math.h>
#include <errno.h>
#define IBMPC 1
#define ANSIPROT 1
#define MINUSZERO 1
#define INFINITIES 1
#define NANS 1
#define DENORMAL 1
#define VOLATILE
#define mtherr(fname, code)
#define XPD 0,
//#define _CEPHES_USE_ERRNO
#ifdef _CEPHES_USE_ERRNO
#define _SET_ERRNO(x) errno = (x)
#else
#define _SET_ERRNO(x)
#endif
/* constants used by cephes functions */
/* double */
#define MAXNUM 1.7976931348623158E308
#define MAXLOG 7.09782712893383996843E2
#define MINLOG -7.08396418532264106224E2
#define LOGE2 6.93147180559945309417E-1
#define LOG2E 1.44269504088896340736
#define PI 3.14159265358979323846
#define PIO2 1.57079632679489661923
#define PIO4 7.85398163397448309616E-1
#define NEGZERO (-0.0)
#undef NAN
#undef INFINITY
#if (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ > 2))
#define INFINITY __builtin_huge_val()
#define NAN __builtin_nan("")
#else
extern double __INF;
#define INFINITY (__INF)
extern double __QNAN;
#define NAN (__QNAN)
#endif
/*long double*/
#define MAXNUML 1.189731495357231765021263853E4932L
#define MAXLOGL 1.1356523406294143949492E4L
#define MINLOGL -1.13994985314888605586758E4L
#define LOGE2L 6.9314718055994530941723E-1L
#define LOG2EL 1.4426950408889634073599E0L
#define PIL 3.1415926535897932384626L
#define PIO2L 1.5707963267948966192313L
#define PIO4L 7.8539816339744830961566E-1L
#define isfinitel isfinite
#define isinfl isinf
#define isnanl isnan
#define signbitl signbit
#define NEGZEROL (-0.0L)
#undef NANL
#undef INFINITYL
#if (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ > 2))
#define INFINITYL __builtin_huge_vall()
#define NANL __builtin_nanl("")
#else
extern long double __INFL;
#define INFINITYL (__INFL)
extern long double __QNANL;
#define NANL (__QNANL)
#endif
/* float */
#define MAXNUMF 3.4028234663852885981170418348451692544e38F
#define MAXLOGF 88.72283905206835F
#define MINLOGF -103.278929903431851103F /* log(2^-149) */
#define LOG2EF 1.44269504088896341F
#define LOGE2F 0.693147180559945309F
#define PIF 3.141592653589793238F
#define PIO2F 1.5707963267948966192F
#define PIO4F 0.7853981633974483096F
#define isfinitef isfinite
#define isinff isinf
#define isnanf isnan
#define signbitf signbit
#define NEGZEROF (-0.0F)
#undef NANF
#undef INFINITYF
#if (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ > 2))
#define INFINITYF __builtin_huge_valf()
#define NANF __builtin_nanf("")
#else
extern float __INFF;
#define INFINITYF (__INFF)
extern float __QNANF;
#define NANF (__QNANF)
#endif
/* double */
/*
Cephes Math Library Release 2.2: July, 1992
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* polevl.c
* p1evl.c
*
* Evaluate polynomial
*
*
*
* SYNOPSIS:
*
* int N;
* double x, y, coef[N+1], polevl[];
*
* y = polevl( x, coef, N );
*
*
*
* DESCRIPTION:
*
* Evaluates polynomial of degree N:
*
* 2 N
* y = C + C x + C x +...+ C x
* 0 1 2 N
*
* Coefficients are stored in reverse order:
*
* coef[0] = C , ..., coef[N] = C .
* N 0
*
* The function p1evl() assumes that coef[N] = 1.0 and is
* omitted from the array. Its calling arguments are
* otherwise the same as polevl().
*
*
* SPEED:
*
* In the interest of speed, there are no checks for out
* of bounds arithmetic. This routine is used by most of
* the functions in the library. Depending on available
* equipment features, the user may wish to rewrite the
* program in microcode or assembly language.
*
*/
/* Polynomial evaluator:
* P[0] x^n + P[1] x^(n-1) + ... + P[n]
*/
static __inline__ double polevl( x, p, n )
double x;
const void *p;
int n;
{
register double y;
register double *P = (double *)p;
y = *P++;
do
{
y = y * x + *P++;
}
while( --n );
return(y);
}
/* Polynomial evaluator:
* x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n]
*/
static __inline__ double p1evl( x, p, n )
double x;
const void *p;
int n;
{
register double y;
register double *P = (double *)p;
n -= 1;
y = x + *P++;
do
{
y = y * x + *P++;
}
while( --n );
return( y );
}
/* long double */
/*
Cephes Math Library Release 2.2: July, 1992
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* polevll.c
* p1evll.c
*
* Evaluate polynomial
*
*
*
* SYNOPSIS:
*
* int N;
* long double x, y, coef[N+1], polevl[];
*
* y = polevll( x, coef, N );
*
*
*
* DESCRIPTION:
*
* Evaluates polynomial of degree N:
*
* 2 N
* y = C + C x + C x +...+ C x
* 0 1 2 N
*
* Coefficients are stored in reverse order:
*
* coef[0] = C , ..., coef[N] = C .
* N 0
*
* The function p1evll() assumes that coef[N] = 1.0 and is
* omitted from the array. Its calling arguments are
* otherwise the same as polevll().
*
*
* SPEED:
*
* In the interest of speed, there are no checks for out
* of bounds arithmetic. This routine is used by most of
* the functions in the library. Depending on available
* equipment features, the user may wish to rewrite the
* program in microcode or assembly language.
*
*/
/* Polynomial evaluator:
* P[0] x^n + P[1] x^(n-1) + ... + P[n]
*/
static __inline__ long double polevll( x, p, n )
long double x;
const void *p;
int n;
{
register long double y;
register long double *P = (long double *)p;
y = *P++;
do
{
y = y * x + *P++;
}
while( --n );
return(y);
}
/* Polynomial evaluator:
* x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n]
*/
static __inline__ long double p1evll( x, p, n )
long double x;
const void *p;
int n;
{
register long double y;
register long double *P = (long double *)p;
n -= 1;
y = x + *P++;
do
{
y = y * x + *P++;
}
while( --n );
return( y );
}
/* Float version */
/* polevlf.c
* p1evlf.c
*
* Evaluate polynomial
*
*
*
* SYNOPSIS:
*
* int N;
* float x, y, coef[N+1], polevlf[];
*
* y = polevlf( x, coef, N );
*
*
*
* DESCRIPTION:
*
* Evaluates polynomial of degree N:
*
* 2 N
* y = C + C x + C x +...+ C x
* 0 1 2 N
*
* Coefficients are stored in reverse order:
*
* coef[0] = C , ..., coef[N] = C .
* N 0
*
* The function p1evl() assumes that coef[N] = 1.0 and is
* omitted from the array. Its calling arguments are
* otherwise the same as polevl().
*
*
* SPEED:
*
* In the interest of speed, there are no checks for out
* of bounds arithmetic. This routine is used by most of
* the functions in the library. Depending on available
* equipment features, the user may wish to rewrite the
* program in microcode or assembly language.
*
*/
/*
Cephes Math Library Release 2.1: December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
static __inline__ float polevlf(float x, const float* coef, int N )
{
float ans;
float *p;
int i;
p = (float*)coef;
ans = *p++;
/*
for( i=0; i<N; i++ )
ans = ans * x + *p++;
*/
i = N;
do
ans = ans * x + *p++;
while( --i );
return( ans );
}
/* p1evl() */
/* N
* Evaluate polynomial when coefficient of x is 1.0.
* Otherwise same as polevl.
*/
static __inline__ float p1evlf( float x, const float *coef, int N )
{
float ans;
float *p;
int i;
p = (float*)coef;
ans = x + *p++;
i = N-1;
do
ans = ans * x + *p++;
while( --i );
return( ans );
}