418 lines
8.3 KiB
C
418 lines
8.3 KiB
C
#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 __x86_64__
|
|
#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*/
|
|
#if defined(__arm__) || defined(_ARM_)
|
|
#define MAXNUML 1.7976931348623158E308
|
|
#define MAXLOGL 7.09782712893383996843E2
|
|
#define MINLOGL -7.08396418532264106224E2
|
|
#define LOGE2L 6.93147180559945309417E-1
|
|
#define LOG2EL 1.44269504088896340736
|
|
#define PIL 3.14159265358979323846
|
|
#define PIO2L 1.57079632679489661923
|
|
#define PIO4L 7.85398163397448309616E-1
|
|
#else
|
|
#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
|
|
#endif /* defined(__arm__) || defined(_ARM_) */
|
|
|
|
#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);
|
|
}
|
|
|