A related function to the the gamma function is the beta function
and could be implemented in Mathematica as
betafn[n_, m_] := NIntegrate[x^(n-1) (1-x)^(m-1), {x,0,1}]
and can be compared with the built-in beta function
betafn[1,1] 1. Beta[1,1] 1 betafn[6,3] 0.00595238 Beta[6,3] 1 --- 168 5! 2! / 8! 1 --- 168 N[5! 2! / 8!] 0.00595238
One situation when the beta function makes its appearance is when calculating the total density of states by performing convolutions of classical densities of states which have a power law form. For example, two subsystems composed of different sets of harmonic oscillators
which leads to the general form above for 
harmonic oscillators.