STK_Funct_gamma.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*  Copyright (C) 2004-2016  Serge Iovleff, Université Lille 1, Inria
 
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.
 
    You should have received a copy of the GNU Lesser General Public
    License along with this program; if not, write to the
    Free Software Foundation, Inc.,
    59 Temple Place,
    Suite 330,
    Boston, MA 02111-1307
    USA
 
    Contact : S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
*/
 
/*
 * Project:  Analysis
 * Purpose:  Declaration of functions around the gamma function
 * Author:   Serge Iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 **/
 
#ifndef STK_FUNCT_GAMMA_H
#define STK_FUNCT_GAMMA_H
 
#include <vector>
 
#include "STK_Const_Math.h"
#include "STK_Const_Sequences.h"
 
namespace STK
{
 
namespace Funct
{
// forward declaration
Real factorial( int);
Real lfactorial( int);
 
Real factorial( Real const&);
Real lfactorial( Real const&);
 
Real gamma( Real const&);
Real lgamma( Real const&);
 
// raw versions
Real factorial_raw( int);
Real lfactorial_raw( int);
 
Real factorial_raw( Real const&);
Real lfactorial_raw( Real const&);
 
Real gamma_raw( Real const&);
Real lgamma_raw( Real const&);
 
inline Real lanczosSerie(Real const& z)
{
  Real sum = 0.0;
  for (int k=20; k>=0; k--)
    sum += Const::lanczosCoefArray[k]/(z+k);
  return 2.0240434640140357514731512432760e-10 + sum;
}
 
inline Real gammaLanczos(Real const& z)
{
  // 2 * sqrt(e/pi) = 1.86038273...
  return 1.8603827342052657173362492472666631120594218414085774528900013
       * exp((z-0.5)*(log(z+22.118910)-1.))
       * lanczosSerie(z);
}
 
inline double stirlingSerie(Real const& z)
{
  const Real z2 = z * z;
  return (z <= 50) ? ( Const::stirlingCoefArray[0]
                       + ( Const::stirlingCoefArray[1]
                         + ( Const::stirlingCoefArray[2]
                           + Const::stirlingCoefArray[3]/z2
                           )/z2
                         )/z2
                     )/z
                   : ( Const::stirlingCoefArray[0]
                       + ( Const::stirlingCoefArray[1]
                         + Const::stirlingCoefArray[2]/z2
                         )/z2
                     )/z
                   ;
}
 
inline Real gammaStirling( Real const& z)
{ return Const::_SQRT2PI_ * exp((z-0.5)*(log(z)-1.)+stirlingSerie(z)-0.5);}
 
inline Real lgammaStirling( Real const& z)
{ return ( Const::_LNSQRT2PI_ + (z-0.5)*log(z) - z + stirlingSerie(z) );}
 
inline Real lgammaStirlingError(Real const& z)
{
  int n = (int)std::floor(z);
  // z is a discrete value
  if (z == n)
  { return (n<100) ? Const::lgammaStirlingErrorArray[n] : stirlingSerie(z);}
  // z is a discrete value halves
  if (z== n + 0.5)
  { return (n<100) ? Const::lgammaStirlingErrorHalvesArray[n] : stirlingSerie(z);}
  // other cases
  return (z > 16) ? stirlingSerie(z) : lgamma(z) - (z-0.5)*log(z) + z - Const::_LNSQRT2PI_;
}
 
inline Real lgammaStirlingError(int n)
{ return (n<100.0) ? Const::lgammaStirlingErrorArray[n] : stirlingSerie(n);}
 
 
inline Real factorial( int n)
{
 // Check if z is Available and finite
 if (!Arithmetic<int>::isFinite(n)) return Arithmetic<Real>::NA();
 if (n < 0) { STKDOMAIN_ERROR_1ARG(factorial,n,"Negative argument");}
 return factorial_raw(n);
}
inline Real factorial_raw( int n)
{ return (n < 51) ? Const::factorialArray[n] : Funct::gamma_raw(n + 1.);}
 
inline Real factorial( Real const& z)
{
  // Check if z is Available and finite
  if (!Arithmetic<Real>::isFinite(z)) return z;
  // Negative integers or reals arguments not allowed
  if (z < 0 ||(z != std::floor(z))) { STKDOMAIN_ERROR_1ARG(factorial,z,"Negative or not discrete argument");}
  return factorial_raw(z);
}
inline Real factorial_raw( Real const& z)
{
  const int n = std::floor(z);
  return (n < 51) ? Const::factorialArray[n] : Funct::gamma_raw(z + 1.);
}
 
inline Real lfactorial(int n)
{
  // Check if z is Available and finite
  if (!Arithmetic<int>::isFinite(n)) return n;
  // Negative integers or reals arguments not allowed
  if (n < 0) { STKDOMAIN_ERROR_1ARG(lfactorial,n,"Negative integer argument");}
  // if n is a small number return a constant
  return lfactorial_raw(n);
}
inline Real lfactorial_raw(int n)
{ return (n < 51) ? Const::factorialLnArray[n] : lgamma( n + 1.);}
 
inline Real lfactorial(Real const& z)
{
  // Check if z is Available and finite
  if (!Arithmetic<Real>::isFinite(z)) return z;
  // Negative integers or real arguments not allowed
  if ((z < 0)||(z != std::floor(z)))
  { STKDOMAIN_ERROR_1ARG(lfactorial,z,"Negative integer or decimal argument");}
  // if n is a small number return a constant
  return lfactorial_raw(z);
}
inline Real lfactorial_raw(Real const& z)
{
  const int n = (int)std::floor(z);
  return (n < 51) ? Const::factorialLnArray[n] : lgamma(z + 1.);
}
 
inline Real gamma(Real const& z)
{
  // Check if z is Available and finite
  if (!Arithmetic<Real>::isFinite(z)) return z;
  // Negative integer argument not allowed
  if ( z<=0 && z == std::floor(z))
  { STKDOMAIN_ERROR_1ARG(Funct::gamma,z,"Negative integer or zero argument");}
  return gamma_raw(z);
}
 
inline Real gamma_raw(Real const& z)
{
  int n = (int)std::floor(z);
  // If z is an integer (z integer and z<0 is not possible)
  if (z == n)
  { return (z < 51) ? Const::factorialArray[(int)z-1] : gammaStirling(z);}
  // compute the absolute value of x -> y and compute the sign of the gamma function
  Real y = std::abs(z);
  int ny = std::floor(y);
  Real signgam = (z<0 && isEven(ny)) ? -1 : 1, value;
  // if y is an integer and a half, use reflection formula
  if (y == ny+0.5)
  { value = (ny<50) ?  Const::factorialHalvesArray[ny] : gammaStirling(y);}
  else // normal case
  {
    if (y <= 8)
    {
      Real r = y-ny;            // r in [0,1]
      value = gammaLanczos(r);  // use Lanzcos approximation
      // scale result
      for (int i=0; i<ny; i++) value *= (r+i);
    }
    else // shift value
    {
      if (y <=16)
      {
        value = gammaStirling(y+6.);
        for (int i=5; i>=0; --i) value /= (y+i);
      }
      else
        value = gammaStirling(y);
    }
  }
  // z >0 terminate
  if (z>0) return value;
  // z <0  -> use reflection formula and check overflow
  Real den = y*sin(Const::_PI_ *y)*value;
  return (den == 0.0) ?  signgam * Arithmetic<Real>::infinity() : -Const::_PI_/den;
}
 
inline Real lgamma(Real const& z)
{
  // Check if z is Available and finite
  if (!Arithmetic<Real>::isFinite(z)) return z;
  // Negative integer argument not allowed
  if ( z<=0 && z == std::floor(z))
  { STKDOMAIN_ERROR_1ARG(Funct::lgamma,z,"Negative integer or zero argument");}
  return lgamma_raw(z);
}
inline Real lgamma_raw(Real const& z)
{
  // compute the absolute value of x -> y
  Real y = std::abs(z), value;
  int ny = std::floor(y);
  // If x is an integer
  if (y == ny)
  { value = (ny<=50) ? Const::factorialLnArray[ny-1] : lgammaStirling(y);}
  else
  {
    // If x is an integer and a half
    if (y == ny+0.5)
    { value = (ny<50) ? Const::factorialLnHalvesArray[ny] : lgammaStirling(y);}
    else
    {
      // small values -> use gamma function
      if (y <= 16)  return log(std::abs(Funct::gamma_raw(z)));
      else          value = lgammaStirling(y);
    }
  }
  // z >0 terminate
  if (z>0) return value;
  // z <0  -> use reflectiong formula
  Real sinpiy = std::abs(sin(Const::_PI_ *y));
  // overflow
  return (sinpiy == 0.0) ? -Arithmetic<Real>::infinity()
                         : Const::_LNSQRTPI_2_+(z-0.5)*log(y)-z-log(sinpiy)+stirlingSerie(z);
}
 
} // namespace Funct
 
} // namespace STK
 
#endif // STK_FUNCT_GAMMA_H