59 lines
1.7 KiB
C++
59 lines
1.7 KiB
C++
#include <math.h>
|
|
#include <numbers>
|
|
#include <stdio.h>
|
|
#include <vector>
|
|
|
|
// Calculate the coefficients used by Spouge's approximation (based on the C
|
|
// implemetation)
|
|
std::vector<double> CalculateCoefficients(int numCoeff)
|
|
{
|
|
std::vector<double> c(numCoeff);
|
|
double k1_factrl = 1.0;
|
|
c[0] = sqrt(2.0 * std::numbers::pi);
|
|
for(size_t k=1; k < numCoeff; k++)
|
|
{
|
|
c[k] = exp(numCoeff-k) * pow(numCoeff-k, k-0.5) / k1_factrl;
|
|
k1_factrl *= -(double)k;
|
|
}
|
|
return c;
|
|
}
|
|
|
|
// The Spouge approximation
|
|
double Gamma(const std::vector<double>& coeffs, double x)
|
|
{
|
|
const size_t numCoeff = coeffs.size();
|
|
double accm = coeffs[0];
|
|
for(size_t k=1; k < numCoeff; k++)
|
|
{
|
|
accm += coeffs[k] / ( x + k );
|
|
}
|
|
accm *= exp(-(x+numCoeff)) * pow(x+numCoeff, x+0.5);
|
|
return accm/x;
|
|
}
|
|
|
|
int main()
|
|
{
|
|
// estimate the gamma function with 1, 4, and 10 coefficients
|
|
const auto coeff1 = CalculateCoefficients(1);
|
|
const auto coeff4 = CalculateCoefficients(4);
|
|
const auto coeff10 = CalculateCoefficients(10);
|
|
|
|
const auto inputs = std::vector<double>{
|
|
0.001, 0.01, 0.1, 0.5, 1.0,
|
|
1.461632145, // minimum of the gamma function
|
|
2, 2.5, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100,
|
|
150 // causes overflow for this implemetation
|
|
};
|
|
|
|
printf("%16s%16s%16s%16s%16s\n", "gamma( x ) =", "Spouge 1", "Spouge 4", "Spouge 10", "built-in");
|
|
for(auto x : inputs)
|
|
{
|
|
printf("gamma(%7.3f) = %16.10g %16.10g %16.10g %16.10g\n",
|
|
x,
|
|
Gamma(coeff1, x),
|
|
Gamma(coeff4, x),
|
|
Gamma(coeff10, x),
|
|
std::tgamma(x)); // built-in gamma function
|
|
}
|
|
}
|