/*							mconf.h
 *
 *	Common include file for math routines
 *
 *
 *
 * SYNOPSIS:
 *
 * #include "mconf.h"
 *
 *
 *
 * DESCRIPTION:
 *
 * This file contains definitions for error codes that are
 * passed to the common error handling routine mtherr()
 * (which see).
 *
 * The file also includes a conditional assembly definition
 * for the type of computer arithmetic (IEEE, DEC, Motorola
 * IEEE, or UNKnown).
 *
 * For Digital Equipment PDP-11 and VAX computers, certain
 * IBM systems, and others that use numbers with a 56-bit
 * significand, the symbol DEC should be defined.  In this
 * mode, most floating point constants are given as arrays
 * of octal integers to eliminate decimal to binary conversion
 * errors that might be introduced by the compiler.
 *
 * For computers, such as IBM PC, that follow the IEEE 
 * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
 * Std 754-1985), the symbol IBMPC should be defined.  These
 * numbers have 53-bit significands.  In this mode, constants
 * are provided as arrays of hexadecimal 16 bit integers.
 *
 * To accommodate other types of computer arithmetic, all
 * constants are also provided in a normal decimal radix
 * which one can hope are correctly converted to a suitable
 * format by the available C language compiler.  To invoke
 * this mode, the symbol UNK is defined.
 *
 * An important difference among these modes is a predefined
 * set of machine arithmetic constants for each.  The numbers
 * MACHEP (the machine roundoff error), MAXNUM (largest number
 * represented), and several other parameters are preset by
 * the configuration symbol.  Check the file const.c to
 * ensure that these values are correct for your computer.
 *
 */

/*
Cephes Math Library Release 2.0:  April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


/* Constant definitions for math error conditions
 */

#define DOMAIN		1	/* argument domain error */
#define SING		2	/* argument singularity */
#define OVERFLOW	3	/* overflow range error */
#define UNDERFLOW	4	/* underflow range error */
#define TLOSS		5	/* total loss of precision */
#define PLOSS		6	/* partial loss of precision */

#define EDOM		33
#define ERANGE		34

typedef struct
	{
	double r;
	double i;
	}cmplx;

/* Type of computer arithmetic */

/* PDP-11, Pro350, VAX:
 */
/*define DEC 1*/

/* Intel IEEE, low order words come first:
 */
#define IBMPC 1

/* Motorola IEEE, high order words come first
 * (Sun workstation):
 */
/*define MIEEE 1*/

/* UNKnown arithmetic, invokes coefficients given in
 * normal decimal format.  Beware of range boundary
 * problems (MACHEP, MAXLOG, etc. in const.c) and
 * roundoff problems in pow.c:
 /*define UNK 1*/

/*							const.c
 *
 *	Globally declared constants
 *
 *
 *
 * SYNOPSIS:
 *
 * extern double nameofconstant;
 *
 *
 *
 *
 * DESCRIPTION:
 *
 * This file contains a number of mathematical constants and
 * also some needed size parameters of the computer arithmetic.
 * The values are supplied as arrays of hexadecimal integers
 * for IEEE arithmetic; arrays of octal constants for DEC
 * arithmetic; and in a normal decimal scientific notation for
 * other machines.  The particular notation used is determined
 * by a symbol (DEC, IBMPC, or UNK) defined in the include file
 * mconf.h.
 *
 * The default size parameters are as follows.
 *
 * For DEC and UNK modes:
 * MACHEP =  1.38777878078144567553E-17       2**-56
 * MAXLOG =  8.8029691931113054295988E1       log(2**127)
 * MINLOG = -8.872283911167299960540E1        log(2**-128)
 * MAXNUM =  1.701411834604692317316873e38    2**127
 *
 * For IEEE arithmetic (IBMPC):
 * MACHEP =  1.11022302462515654042E-16       2**-53
 * MAXLOG =  7.09782712893383996843E2         log(2**1024)
 * MINLOG = -7.08396418532264106224E2         log(2**-1022)
 * MAXNUM =  1.7976931348623158E308           2**1024
 *
 * The global symbols for mathematical constants are
 * PI     =  3.14159265358979323846           pi
 * PIO2   =  1.57079632679489661923           pi/2
 * PIO4   =  7.85398163397448309616E-1        pi/4
 * SQRT2  =  1.41421356237309504880           sqrt(2)
 * SQRTH  =  7.07106781186547524401E-1        sqrt(2)/2
 * LOG2E  =  1.4426950408889634073599         1/log(2)
 * SQ2OPI =  7.9788456080286535587989E-1      sqrt( 2/pi )
 * LOGE2  =  6.93147180559945309417E-1        log(2)
 * LOGSQ2 =  3.46573590279972654709E-1        log(2)/2
 * THPIO4 =  2.35619449019234492885           3*pi/4
 * TWOOPI =  6.36619772367581343075535E-1     2/pi
 *
 * These lists are subject to change.
 */

/*							const.c */

/*
Cephes Math Library Release 2.0:  April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#include "mconf.h"

#ifdef UNK
double MACHEP =  1.38777878078144567553E-17;   /* 2**-56 */
double MAXLOG =  8.8029691931113054295988E1;   /* log(2**127) */
double MINLOG = -8.872283911167299960540E1;    /* log(2**-128) */
double MAXNUM =  1.701411834604692317316873e38; /* 2**127 */
double PI     =  3.14159265358979323846;       /* pi */
double PIO2   =  1.57079632679489661923;       /* pi/2 */
double PIO4   =  7.85398163397448309616E-1;    /* pi/4 */
double SQRT2  =  1.41421356237309504880;       /* sqrt(2) */
double SQRTH  =  7.07106781186547524401E-1;    /* sqrt(2)/2 */
double LOG2E  =  1.4426950408889634073599;     /* 1/log(2) */
double SQ2OPI =  7.9788456080286535587989E-1;  /* sqrt( 2/pi ) */
double LOGE2  =  6.93147180559945309417E-1;    /* log(2) */
double LOGSQ2 =  3.46573590279972654709E-1;    /* log(2)/2 */
double THPIO4 =  2.35619449019234492885;       /* 3*pi/4 */
double TWOOPI =  6.36619772367581343075535E-1; /* 2/pi */
#endif

#ifdef IBMPC
			/* 2**-53 =  1.11022302462515654042E-16 */
short MACHEP[4] = {0x0000,0x0000,0x0000,0x3ca0};
			/* log(2**1024) =   7.09782712893383996843E2 */
short MAXLOG[4] = {0x39ef,0xfefa,0x2e42,0x4086};
			/* log(2**-1022) = - 7.08396418532264106224E2 */
short MINLOG[4] = {0xbcd2,0xdd7a,0x232b,0xc086};
			/* 2**1024*(1-MACHEP) =  1.7976931348623158E308 */
short MAXNUM[4] = {0xffff,0xffff,0xffff,0x7fef};
short PI[4]     = {0x2d18,0x5444,0x21fb,0x4009};
short PIO2[4]   = {0x2d18,0x5444,0x21fb,0x3ff9};
short PIO4[4]   = {0x2d18,0x5444,0x21fb,0x3fe9};
short SQRT2[4]  = {0x3bcd,0x667f,0xa09e,0x3ff6};
short SQRTH[4]  = {0x3bcd,0x667f,0xa09e,0x3fe6};
short LOG2E[4]  = {0x82fe,0x652b,0x1547,0x3ff7};
short SQ2OPI[4] = {0x3651,0x33d4,0x8845,0x3fe9};
short LOGE2[4]  = {0x39ef,0xfefa,0x2e42,0x3fe6};
short LOGSQ2[4] = {0x39ef,0xfefa,0x2e42,0x3fd6};
short THPIO4[4] = {0x21d2,0x7f33,0xd97c,0x4002};
short TWOOPI[4] = {0xc883,0x6dc9,0x5f30,0x3fe4};
#endif

#ifdef MIEEE
			/* 2**-53 =  1.11022302462515654042E-16 */
short MACHEP[4] = {
0x3ca0,0x0000,0x0000,0x0000
};
			/* log(2**1024) =   7.09782712893383996843E2 */
short MAXLOG[4] = {
0x4086,0x2e42,0xfefa,0x39ef
};
			/* log(2**-1022) = - 7.08396418532264106224E2 */
short MINLOG[4] = {
0xc086,0x232b,0xdd7a,0xbcd2
};
			/* 2**1024*(1-MACHEP) =  1.7976931348623158E308 */
short MAXNUM[4] = {
0x7fef,0xffff,0xffff,0xffff
};
short PI[4]     = {
0x4009,0x21fb,0x5444,0x2d18
};
short PIO2[4]   = {
0x3ff9,0x21fb,0x5444,0x2d18
};
short PIO4[4]   = {
0x3fe9,0x21fb,0x5444,0x2d18
};
short SQRT2[4]  = {
0x3ff6,0xa09e,0x667f,0x3bcd
};
short SQRTH[4]  = {
0x3fe6,0xa09e,0x667f,0x3bcd
};
short LOG2E[4]  = {
0x3ff7,0x1547,0x652b,0x82fe
};
short SQ2OPI[4] = {
0x3fe9,0x8845,0x33d4,0x3651
};
short LOGE2[4]  = {
0x3fe6,0x2e42,0xfefa,0x39ef
};
short LOGSQ2[4] = {
0x3fd6,0x2e42,0xfefa,0x39ef
};
short THPIO4[4] = {
0x4002,0xd97c,0x7f33,0x21d2
};
short TWOOPI[4] = {
0x3fe4,0x5f30,0x6dc9,0xc883
};
#endif

#ifdef DEC
			/* 2**-56 =  1.38777878078144567553E-17 */
short MACHEP[4] = {0022200,0000000,0000000,0000000};
			/* log 2**127 = 88.029691931113054295988 */
short MAXLOG[4] = {041660,007463,0143742,025733,};
			/* log 2**-128 = -88.72283911167299960540 */
short MINLOG[4] = {0141661,071027,0173721,0147572,};
			/* 2**127 = 1.701411834604692317316873e38 */
short MAXNUM[4] = {077777,0177777,0177777,0177777,};
short PI[4]     = {040511,007732,0121041,064302,};
short PIO2[4]   = {040311,007732,0121041,064302,};
short PIO4[4]   = {040111,007732,0121041,064302,};
short SQRT2[4]  = {040265,002363,031771,0157145,};
short SQRTH[4]  = {040065,002363,031771,0157144,};
short LOG2E[4]  = {040270,0125073,024534,013761,};
short SQ2OPI[4] = {040114,041051,0117241,0131204,};
short LOGE2[4]  = {040061,071027,0173721,0147572,};
short LOGSQ2[4] = {037661,071027,0173721,0147572,};
short THPIO4[4] = {040426,0145743,0174631,007222,};
short TWOOPI[4] = {040042,0174603,067116,042025,};
#endif

extern short MACHEP[];
extern short MAXLOG[];
extern short MINLOG[];
extern short MAXNUM[];
extern short PI[];
extern short PIO2[];
extern short PIO4[];
extern short SQRT2[];
extern short SQRTH[];
extern short LOG2E[];
extern short SQ2OPI[];
extern short LOGE2[];
extern short LOGSQ2[];
extern short THPIO4[];
extern short TWOOPI[];
static char x[] =
"Cephes Math Library Copyright 1987 by Stephen L. Moshier";

/*							polevl.c
 *							p1evl.c
 *
 *	Evaluate polynomial
 *
 *
 *
 * SYNOPSIS:
 *
 * int N;
 * double x, y, coef[N+1], polevl[];
 *
 * y = polevl( C, x, 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.0:  April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


double polevl( x, coef, N )
double x;
double coef[];
int N;
{
double ans;
short i;
register double *p;

p = coef;
ans = *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.
 */

double p1evl( x, coef, N )
double x;
double coef[];
int N;
{
double ans;
register double *p;
short i;

p = coef;
ans = x + *p++;
i = N-1;

do
	ans = ans * x  + *p++;
while( --i );

return( ans );
}


/*							chbevl.c
 *
 *	Evaluate Chebyshev series
 *
 *
 *
 * SYNOPSIS:
 *
 * int N;
 * double x, y, coef[N], chebevl();
 *
 * y = chbevl( x, coef, N );
 *
 *
 *
 * DESCRIPTION:
 *
 * Evaluates the series
 *
 *        N-1
 *         - '
 *  y  =   >   coef[i] T (x/2)
 *         -            i
 *        i=0
 *
 * of Chebyshev polynomials Ti at argument x/2.
 *
 * Coefficients are stored in reverse order, i.e. the zero
 * order term is last in the array.  Note N is the number of
 * coefficients, not the order.
 *
 * If coefficients are for the interval a to b, x must
 * have been transformed to x -> 2(2x - b - a)/(b-a) before
 * entering the routine.  This maps x from (a, b) to (-1, 1),
 * over which the Chebyshev polynomials are defined.
 *
 * If the coefficients are for the inverted interval, in
 * which (a, b) is mapped to (1/b, 1/a), the transformation
 * required is x -> 2(2ab/x - b - a)/(b-a).  If b is infinity,
 * this becomes x -> 4a/x - 1.
 *
 *
 *
 * SPEED:
 *
 * Taking advantage of the recurrence properties of the
 * Chebyshev polynomials, the routine requires one more
 * addition per loop than evaluating a nested polynomial of
 * the same degree.
 *
 */
/*							chbevl.c	*/

/*
Cephes Math Library Release 2.0:  April, 1987
Copyright 1985, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

double chbevl( x, array, n )
double x;
double array[];
int n;
{
double b0, b1, b2, *p;
int i;

p = array;
b0 = *p++;
b1 = 0.0;
i = n - 1;

do
	{
	b2 = b1;
	b1 = b0;
	b0 = x * b1  -  b2  + *p++;
	}
while( --i );

return( 0.5*(b0-b2) );
}

/*							mtherr.c
 *
 *	Library common error handling routine
 *
 *
 *
 * SYNOPSIS:
 *
 * char *fctnam;
 * int code;
 * void mtherr();
 *
 * mtherr( fctnam, code );
 *
 *
 *
 * DESCRIPTION:
 *
 * This routine may be called to report one of the following
 * error conditions (in the include file mconf.h).
 *  
 *   Mnemonic        Value          Significance
 *
 *    DOMAIN            1       argument domain error
 *    SING              2       function singularity
 *    OVERFLOW          3       overflow range error
 *    UNDERFLOW         4       underflow range error
 *    TLOSS             5       total loss of precision
 *    PLOSS             6       partial loss of precision
 *    EDOM             33       Unix domain error code
 *    ERANGE           34       Unix range error code
 *
 * The default version of the file prints the function name,
 * passed to it by the pointer fctnam, followed by the
 * error condition.  The display is directed to the standard
 * output device.  The routine then returns to the calling
 * program.  Users may wish to modify the program to abort by
 * calling exit() under severe error conditions such as domain
 * errors.
 *
 * Since all error conditions pass control to this function,
 * the display may be easily changed, eliminated, or directed
 * to an error logging device.
 *
 * SEE ALSO:
 *
 * mconf.h
 *
 */

/*
Cephes Math Library Release 2.0:  April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#include "mconf.h"

/* Notice: the order of appearance of the following
 * messages is bound to the error codes defined
 * in mconf.h.
 */
static char *ermsg[7] = {
"unknown",      /* error code 0 */
"domain",       /* error code 1 */
"singularity",  /* et seq.      */
"overflow",
"underflow",
"total loss of precision",
"partial loss of precision"
};



mtherr( name, code )
char *name;
int code;
{

/* Display string passed by calling program,
 * which is supposed to be the name of the
 * function in which the error occurred:
 */
printf( "\n%s ", name );

/* Display error message defined
 * by the code argument.
 */
if( (code <= 0) || (code >= 6) )
	code = 0;
printf( "%s error\n", ermsg[code] );

/* Return to calling
 * program
 */
}