/* Various mathematical functions missing from C */
#include <math.h>

#include "support.h"

double sign(r)
double r;
{
   if (r < 0.0) return(-1.0);
   else if (r == 0.0) return (0.0);
   else return(1.0);
}

void Rect(r, phi, x, y)
double r, phi, *x, *y;
{
   *x = r * cos(phi);
   *y = r * sin(phi);
}

void Polar(x, y, r, phi)
double x, y, *r, *phi;
{
   *r = sqrt(x * x + y * y);
   *phi = atan2(y, x);
}

double stirling(n)
double n;
{
   double y = 1 / (12 * n);

   return (pow(n / E, n) * sqrt(2 * PI * n) * (1 + y * (1 + y * (0.5 - y * (4.6333333333333333 + y * 4.7583333333333333)))));
}

double gamma(x)
double x;
{
   double fx, tx, res, i;

   if (x >= 15.0) return(stirling(x - 1));
   else {
      if ((fx = modf(x, &tx)) < 0) { tx -= 1.0; fx += 1.0; } /* give real int & frac */

      if (fx == 0 && tx < 0) return(-HUGE);
      if (tx < -200) return(0.0); /* Underflow */

      res = stirling(fx + 14.0);
      for (i = 14.0; i >= tx; i -= 1.0) res /= i + fx;

      return(res);
   }
}

double factorial(x)
int x;
{
   double r = 1.0;

   if (x > 250) r = HUGE; /* Certainly too big */
   else for (; x > 0; x--) r *= x;

   return(r);
}

double Perm(x, y)
int x, y;
{
   double i, res = 1.0, lim = x - y;

   for (i = x; i > lim; i -= 1.0) res *= i;

   return(res);
}

double Comb(x, y)
int x,y;
{
   double i, lim = y, res = Perm(x, y);

   for (i = 1; i <= lim; i += 1.0) res /= i;

   return(res);
}

double hr(x)
double x;
{
   double h, m, s;

   /* f = modf(x, &i) returns the frcational part of x in f and the integral part in i (all double) */
   m = 100.0 * modf(x, &h);
   s = 100.0 * modf(m, &m);

   return(h + m / 60.0 + s / 3600.0);
}

double hms(x)
double x;
{
   double h, m, s;

   m = 60.0 * modf(x, &h);
   s = 60.0 * modf(m, &m);

   return(h + m / 100.0 + s / 10000.0);
}

double trunc(x)
double x;
{
   modf(x, &x);
   return(x);
}

double frac(x)
double x;
{
   return(modf(x, &x));
}

double asinh(x)
double x;
{
   return(log(x + sqrt(x * x + 1)));
}

double acosh(x)
double x;
{
   if (x < 1.0) return(0.0);
   else return(log(x + sqrt(x * x -1)));
}

double atanh(x)
double x;
{
   if (x > 1.0) return(0.0);
   else return(log((1.0 + x) / (1.0 - x)) / 2.0);
}

double neg(x)
double x;
{
    return(-x);
}


double sqr(x)
double x;
{
    return(x * x);
}

double exp10(x)
double x;
{
    return(pow(10.0, x));
}