| #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, | 
 | #ifdef _WIN64 | 
 | #define XPD_SHORT 0, 0, | 
 | #define XPD_LONG 0, | 
 | #else | 
 | #define XPD_SHORT | 
 | #define XPD_LONG | 
 | #endif | 
 |  | 
 | #if UNK | 
 | typedef union uLD { long double ld; unsigned short sh[8]; long lo[4]; } uLD; | 
 | typedef union uD { double d; unsigned short sh[4]; } uD; | 
 | #elif IBMPC | 
 | typedef union uLD { unsigned short sh[8]; long double ld; long lo[4]; } uLD; | 
 | typedef union uD { unsigned short sh[4]; double d; } uD; | 
 | #elif MIEEE | 
 | typedef union uLD { long lo[4]; long double ld; unsigned short sh[8]; } uLD; | 
 | typedef union uD { unsigned short sh[4]; double d; } uD; | 
 | #else | 
 | #error Unknown uLD/uD type definition | 
 | #endif | 
 |  | 
 | #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(double x, const uD *p, int n) | 
 | { | 
 | 	register double y; | 
 |  | 
 | 	y = p->d; | 
 | 	p++; | 
 | 	do | 
 | 	{ | 
 | 		y = y * x + p->d; | 
 | 		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(double x, const uD *p, int n) | 
 | { | 
 | 	register double y; | 
 |  | 
 | 	n -= 1; | 
 | 	y = x + p->d; p++; | 
 | 	do | 
 | 	{ | 
 | 		y = y * x + p->d; 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(long double x, const uLD *p, int n) | 
 | { | 
 | 	register long double y; | 
 |  | 
 | 	y = p->ld; | 
 | 	p++; | 
 | 	do | 
 | 	{ | 
 | 		y = y * x + p->ld; | 
 | 		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(long double x, const uLD *p, int n) | 
 | { | 
 | 	register long double y; | 
 |  | 
 | 	n -= 1; | 
 | 	y = x + p->ld; | 
 | 	p++; | 
 |  | 
 | 	do | 
 | 	{ | 
 | 		y = y * x + p->ld; | 
 | 		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); | 
 | } | 
 |  |