mathplus.h

Go to the documentation of this file.

// $Id: mathplus.h,v 1.16.2.5 2019/05/28 11:13:02 garloff Exp $

#ifndef TBCI_MATHPLUS_H
#define TBCI_MATHPLUS_H  

/* These do NOT include basics.h and do NOT use namespace TBCI ! */
/* It should be ANSI C (not: C++) compliant */

#include <math.h>
#include "tbci/constants.h"

#ifndef MATH__
# define MATH__
#endif

#ifndef LONG_LONG
# ifdef __GNUC__
#  define LONG_LONG long long
# else
#  define LONG_LONG long
# endif
#endif

#if defined(__cplusplus) && defined(NAMESPACE_TBCI)
NAMESPACE_TBCI
#endif


inline double bernoulli(const double x)
{
  if (MATH__ fabs(x)<0.01) return ((-1.0/720.0*x*x+1.0/12.0)*x-0.5)*x +1.0;
  return (x/(MATH__ exp(x)-1));
}


#ifndef __USE_MISC
inline double asinh(const double x) 
{
        if(-x == MATH__ sqrt(1+x*x)) return MATH__ log (1/(2 * MATH__ fabs(x)));
        else return MATH__ log(x + MATH__ sqrt(1+x*x) );
}
#endif

inline unsigned long binom (const unsigned long x, unsigned long y)
{
   unsigned LONG_LONG e, d;
   unsigned long i;
   if (y > x) return 0;
   if (y == x || y == 0) return 1;
   if (y > x/2) y = x - y;
   e = x; d = y;
   { for (i = (x-1); i > (x-y); i--) e *= i; }
   { for (i = (y-1); i > 1; i--) d *= i; }
   return (e / d);
}

inline unsigned long binomial (const unsigned long y, const unsigned long z)
{
   unsigned LONG_LONG e, d;
   unsigned long i, yy, x = y + z;
   if (y == x || y == 0) return 1;
   if (y > x/2) yy = x - y; else yy = y;
   e = x; d = yy;
   { for (i = (x-1); i > (x-yy); i--) e *= i; }
   { for (i = (yy-1); i > 1; i--) d *= i; }
   return (e / d);
}


inline unsigned LONG_LONG fac (const unsigned char x)
{
   unsigned char i;
   unsigned LONG_LONG f = 1;
   for (i = 2; i <= x; i++) f *= i;
   return f;
}


inline unsigned long trinomial (const unsigned long x, const unsigned long y, const unsigned long z)
{
  return fac((unsigned char)(x+y+z)) / ( fac((unsigned char)x)
                                        *fac((unsigned char)y)
                                        *fac((unsigned char)z) );
} 


inline long double fact (const double x)
{
  long double corr, r = 1.0/x;
  //Abramov/Stegun (51840/139 = 372.94964..)
  //corr = 1.0 + r/12 + r/(x*288) - 139.0*r/(x*x*51840);
  //optimized for [42;43]
  corr = 1.0 + r/12.00000034 + r/(x*287.97452028) - r/(x*x*371.045958768);
  return (MATH__ pow(x/E1,x) * MATH__ sqrt(x*(pi*2.0)) * corr);
}

inline long double ldgamma (const double x)
{
  long double res; double x0 = x - 1.0;
  if (x0 >= 42.0 && x0 <= 43.0) return fact (x0);
  x0 += 42.0 - MATH__ floor(x0); // ]42.0;43.0[
  res = fact (x0);
  if ((x-x0) < 0.5) while ((x-x0) < 0.5) { res /= x0; x0 -= 1.0; }
  else while ((x-x0) > 1.5) { x0 += 1.0; res *= x0; }
  return res;
}

inline double poisson (const double x, const double la)
{
        return (MATH__ pow(la,x) * MATH__ exp(-la) / ldgamma(x+1.0));
}

inline double chi_s (const double x, const double n)
{
   double k = 1.0 / (MATH__ pow(2.0,n/2.0)*ldgamma(n/2.0));
   return k * MATH__ pow(x,(n-2.0)/2.0) * MATH__ exp (-x/2.0);
}

inline double erfc3 (const double x, const double c, const double s)
{
        double d = (x-c)/s;
        return MATH__ exp (-0.5*d*d) / (s * MATH__ sqrt(2.0*pi));
}

#if defined(__cplusplus) && defined(NAMESPACE_TBCI)
NAMESPACE_END
#endif
#endif  /* TBCI_MATHPLUS_H */