#include  <math.h>

extern int errno;

/* This math function file is to make up for deficiencies in some
 * library math functions.  The following are missing for some vendors:
 *
 *   (double) asinh( double )   inverse hyperbolic sin
 *   (double) acosh( double )   inverse hyperbolic cos
 *   (double) atanh( double )   inverse hyperbolic tan
 *   (double) cbrt(  double )   cubic root
 *   (double) rint(  double )   nearest rounded integer
 *   (double) trunc( double )   truncated integer
 *
 * For Linux, it's rint, trunc, cbrt, gamma.
 *
 * Some versions of SunOS lack trunc().
 */

#ifdef AMIGA
double gamma(double x) {
	return exp(lgamma(x));
}
#endif

#ifdef NOGAMMA
double gamma(double x)
{
    double const cof[] = { 76.18009173E0, -86.50532033E0, 24.01409822E0,
			  -1.231739516E0, 0.120858003E-2, -0,536382E-5 } ;
    double const stp = 2.50662827465, fpf = 5.5, pi=3.1415926535897932;

    int i;
    double tmp,ser;

printf("gamma(%f) = ", x);
    if (x < 1.)
	return log(pi/sin(pi*(1.-x))) - gamma(1.-x);
    else
    {
	x = x - 1.0;
	tmp = x + fpf;
	tmp = (x+.5)*log(tmp)-tmp;
	ser = 1.0;
	for (i=0; i<6; i++)
	{
	    x = x + 1.0;
	    ser += cof[i]/x;
	}
printf("%f\n", tmp+log(stp*ser));
	return tmp+log(stp*ser);
    }
}
#endif

#ifdef NORINT
double rint(double x)
{
     double diff;

     diff = x - floor(x);
     if (diff >= 0.5) {
	 	return(ceil(x));
	 }
     else {
	   	return(floor(x));
	 }
}
#endif

#ifdef NOTRUNC
double trunc(double x)
{
    if (x < 0)
		return(ceil(x));
    else
		return(floor(x));
}
#endif

#ifdef NOCBRT
#define  ONE_THIRD 0.3333333333333333333333333333
double cbrt(double x)
{
    double y, val, pow();

    y = fabs(x);
    val = pow(y, ONE_THIRD) * (x/y);
    return(val);
}
#endif

#ifdef NOHYP_TRIGO
double acosh(double x)
{
	double val;

	val = log(x) + log(1.0 + sqrt(1.0 - 1.0/(x*x)));
	return(val);
}

double asinh(double x)
{
	double val;

    if(x == 0.0)
		return(0.0);
	val = log(fabs(x)) + log(1.0 + sqrt(1.0 + 1.0/(x*x)));
	if(x < 0)
		val = -val;
	return(val);
}

double atanh(double x)
{
	double val;

	val = 0.5 * log((1.0 + x)/(1.0 - x));
	return(val);
}

#endif