lib/rubystats/modules.rb
# Gamma function.
# Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz<BR>
# Applied Mathematics Division<BR>
# Argonne National Laboratory<BR>
# Argonne, IL 60439<BR>
# <P>
# References:
# <OL>
# <LI>"An Overview of Software Development for Special Functions", W. J. Cody, Lecture Notes in Mathematics, 506, Numerical Analysis Dundee, 1975, G. A. Watson (ed.), Springer Verlag, Berlin, 1976.
# <LI>Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.
# </OL></P><P>
# From the original documentation:
# </P><P>
# This routine calculates the Gamma function for a real argument X.
# Computation is based on an algorithm outlined in reference 1.
# The program uses rational functions that approximate the Gamma
# function to at least 20 significant decimal digits.  Coefficients
# for the approximation over the interval (1,2) are unpublished.
# Those for the approximation for X .GE. 12 are from reference 2.
# The accuracy achieved depends on the arithmetic system, the
# compiler, the intrinsic functions, and proper selection of the
# machine-dependent constants.
# </P><P>
# Error returns:<BR>
# The program returns the value XINF for singularities or when overflow would occur.
# The computation is believed to be free of underflow and overflow.
# </P>
# Author:: Jaco van Kooten
Gamma function.
Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz
Applied Mathematics Division
Argonne National Laboratory
Argonne, IL 60439

References:

"An Overview of Software Development for Special Functions", W. J. Cody, Lecture Notes in Mathematics, 506, Numerical Analysis Dundee, 1975, G. A. Watson (ed.), Springer Verlag, Berlin, 1976.
Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.

From the original documentation:

This routine calculates the Gamma function for a real argument X.
Computation is based on an algorithm outlined in reference 1.
The program uses rational functions that approximate the Gamma
function to at least 20 significant decimal digits.  Coefficients
for the approximation over the interval (1,2) are unpublished.
Those for the approximation for X .GE. 12 are from reference 2.
The accuracy achieved depends on the arithmetic system, the
compiler, the intrinsic functions, and proper selection of the
machine-dependent constants.

Error returns:
The program returns the value XINF for singularities or when overflow would occur.
The computation is believed to be free of underflow and overflow.

Author:: Jaco van Kooten
Gamma function . Based on public domain NETLIB ( Fortran ) code by W . J . Cody and L . Stoltz Applied Mathematics Division Argonne National Laboratory Argonne IL 60439 References : " An Overview of Software Development for Special Functions " W . J . Cody Lecture Notes in Mathematics 506 Numerical Analysis Dundee 1975 G . A . Watson ( ed . ) Springer Verlag Berlin 1976 . Computer Approximations Hart Et . Al . Wiley and sons New York 1968 . From the original documentation : This routine calculates the Gamma function for a real argument X . Computation is based on an algorithm outlined in reference 1 . The program uses rational functions that approximate the Gamma function to at least 20 significant decimal digits . Coefficients for the approximation over the interval ( 1 2 ) are unpublished . Those for the approximation for X . GE . 12 are from reference 2 . The accuracy achieved depends on the arithmetic system the compiler the intrinsic functions and proper selection of the machine - dependent constants . Error returns : The program returns the value XINF for singularities or when overflow would occur . The computation is believed to be free of underflow and overflow . Author :: Jaco van Kooten
def orig_gamma(x) 
      g_p = [ -1.71618513886549492533811, 24.7656508055759199108314,
        -379.804256470945635097577, 629.331155312818442661052,
        866.966202790413211295064, -31451.2729688483675254357,
        -36144.4134186911729807069, 66456.1438202405440627855 ]
      g_q = [-30.8402300119738975254353, 315.350626979604161529144,
        -1015.15636749021914166146, -3107.77167157231109440444,
        22538.1184209801510330112, 4755.84627752788110767815,
        -134659.959864969306392456, -115132.259675553483497211 ] 
      g_c = [-0.001910444077728, 8.4171387781295e-4, -5.952379913043012e-4,
        7.93650793500350248e-4, -0.002777777777777681622553,
        0.08333333333333333331554247, 0.0057083835261 ]
      fact = 1.0
      i = 0
      n = 0
      y = x
      parity = false
      if y <= 0.0 
        y = -(x)
        y1 = y.to_i
        res = y - y1
        if res != 0.0 
          if y1 != (((y1*0.5).to_i) * 2.0)
            parity = true
            fact = -M_pi/sin(M_pi * res)
            y += 1
          end
        else
          return MAX_VALUE
        end
      end
      if y < EPS
        if y >= XMININ
          res = 1.0 / y
        else
          return MAX_VALUE
        end
      elsif y < 12.0
        y1 = y
        if y < 1.0
          z = y
          y += 1
        else 
          n = y.to_i - 1
          y -= n.to_f
          z = y - 1.0
        end
        xnum = 0.0
        xden = 1.0
        for i in (0...8) 
          xnum = (xnum + g_p[i]) * z
          xden = xden * z + g_q[i]
        end
        res = xnum / xden + 1.0
        if y1 < y
          res /= y1
        elsif y1 > y 
          for i in (0...n)
            res *= y
            y += 1
          end
        end
      else 
        if y <= GAMMA_X_MAX_VALUE
          ysq = y * y
          sum = g_c[6]
          for i in(0...6) 
            sum = sum / ysq + g_c[i]
            sum = sum / y - y + log(SQRT2PI)
            sum += (y - 0.5) * log(y)
            res = Math.exp(sum)
          end
        else
          return MAX_VALUE
        end
        if parity
          res = -res
          if fact != 1.0
            res = fact / res
            return res
          end
        end
      end
    end
def orig_gamma ( x ) g_p = [ - 1.71618513886549492533811 , 24.7656508055759199108314 , - 379.804256470945635097577 , 629.331155312818442661052 , 866.966202790413211295064 , - 31451.2729688483675254357 , - 36144.4134186911729807069 , 66456.1438202405440627855 ] g_q = [ - 30.8402300119738975254353 , 315.350626979604161529144 , - 1015.15636749021914166146 , - 3107.77167157231109440444 , 22538.1184209801510330112 , 4755.84627752788110767815 , - 134659.959864969306392456 , - 115132.259675553483497211 ] g_c = [ - 0.001910444077728 , 8.4171387781295e-4 , - 5.952379913043012e-4 , 7.93650793500350248e-4 , - 0.002777777777777681622553 , 0.08333333333333333331554247 , 0.0057083835261 ] fact = 1.0 i = 0 n = 0 y = x parity = false if y <= 0.0 y = - ( x ) y1 = y . to_i res = y - y1 if res != 0.0 if y1 != ( ( ( y1 * 0.5 ) . to_i ) * 2.0 ) parity = true fact = - M_pi / sin ( M_pi * res ) y += 1 end else return MAX_VALUE end end if y < EPS if y >= XMININ res = 1.0 / y else return MAX_VALUE end elsif y < 12.0 y1 = y if y < 1.0 z = y y += 1 else n = y . to_i - 1 y -= n . to_f z = y - 1.0 end xnum = 0.0 xden = 1.0 for i in ( 0 ... 8 ) xnum = ( xnum + g_p [ i ] ) * z xden = xden * z + g_q [ i ] end res = xnum / xden + 1.0 if y1 < y res /= y1 elsif y1 > y for i in ( 0 ... n ) res *= y y += 1 end end else if y <= GAMMA_X_MAX_VALUE ysq = y * y sum = g_c [ 6 ] for i in ( 0 ... 6 ) sum = sum / ysq + g_c [ i ] sum = sum / y - y + log ( SQRT2PI ) sum += ( y - 0.5 ) * log ( y ) res = Math . exp ( sum ) end else return MAX_VALUE end if parity res = - res if fact != 1.0 res = fact / res return res end end end end
Gamma function.
Gamma function .
# Gamma related constants # ---------------------------------------------------------------------- #  Argument is negative # ---------------------------------------------------------------------- # ---------------------------------------------------------------------- #  Argument is positive # ---------------------------------------------------------------------- # ---------------------------------------------------------------------- #  Argument .LT. EPS # ---------------------------------------------------------------------- #end # ---------------------------------------------------------------------- #  0.0 .LT. argument .LT. 1.0 # ---------------------------------------------------------------------- # ---------------------------------------------------------------------- #  1.0 .LT. argument .LT. 12.0, reduce argument if necessary # ---------------------------------------------------------------------- # ---------------------------------------------------------------------- #  Evaluate approximation for 1.0 .LT. argument .LT. 2.0 # ---------------------------------------------------------------------- # ---------------------------------------------------------------------- #  Adjust result for case  0.0 .LT. argument .LT. 1.0 # ---------------------------------------------------------------------- # ---------------------------------------------------------------------- #  Adjust result for case  2.0 .LT. argument .LT. 12.0 # ---------------------------------------------------------------------- # ---------------------------------------------------------------------- #  Evaluate for argument .GE. 12.0 # ---------------------------------------------------------------------- # ---------------------------------------------------------------------- #  Final adjustments and return # ----------------------------------------------------------------------
x
<not_specific>