2020-12-31 07:32:16 -08:00
|
|
|
// Copyright © 2008-2021 Pioneer Developers. See AUTHORS.txt for details
|
2012-09-15 17:59:15 -07:00
|
|
|
// Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
|
2012-09-12 04:38:30 -07:00
|
|
|
|
2011-10-31 15:27:17 -07:00
|
|
|
// Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it.
|
|
|
|
|
2012-07-30 02:06:04 -07:00
|
|
|
#include "Win32Setup.h"
|
|
|
|
|
2011-10-31 15:27:17 -07:00
|
|
|
#include "Pi.h"
|
2019-01-02 08:59:07 -08:00
|
|
|
#include "libs.h"
|
2011-10-31 15:27:17 -07:00
|
|
|
#include <cmath>
|
|
|
|
#include <iostream>
|
2019-01-02 08:59:07 -08:00
|
|
|
#include <sstream>
|
2011-10-31 15:27:17 -07:00
|
|
|
#include <stdexcept>
|
|
|
|
|
2012-07-28 15:06:49 -07:00
|
|
|
#include "WinMath.h"
|
2011-10-31 15:27:17 -07:00
|
|
|
|
|
|
|
// Note that the functions Gamma and LogGamma are mutually dependent.
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
double Gamma(
|
|
|
|
double x // We require x > 0
|
2011-10-31 15:27:17 -07:00
|
|
|
)
|
|
|
|
{
|
2019-01-02 08:59:07 -08:00
|
|
|
if (x <= 0.0) {
|
2011-10-31 15:27:17 -07:00
|
|
|
std::stringstream os;
|
2019-01-02 08:59:07 -08:00
|
|
|
os << "Invalid input argument " << x << ". Argument must be positive.";
|
|
|
|
throw std::invalid_argument(os.str());
|
2011-10-31 15:27:17 -07:00
|
|
|
}
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
// Split the function domain into three intervals:
|
|
|
|
// (0, 0.001), [0.001, 12), and (12, infinity)
|
2011-10-31 15:27:17 -07:00
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
///////////////////////////////////////////////////////////////////////////
|
|
|
|
// First interval: (0, 0.001)
|
2011-10-31 15:27:17 -07:00
|
|
|
//
|
|
|
|
// For small x, 1/Gamma(x) has power series x + gamma x^2 - ...
|
|
|
|
// So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
|
|
|
|
// The relative error over this interval is less than 6e-7.
|
|
|
|
|
|
|
|
const double gamma = 0.577215664901532860606512090; // Euler's gamma constant
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
if (x < 0.001)
|
|
|
|
return 1.0 / (x * (1.0 + gamma * x));
|
2011-10-31 15:27:17 -07:00
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
///////////////////////////////////////////////////////////////////////////
|
|
|
|
// Second interval: [0.001, 12)
|
2012-07-04 12:16:48 -07:00
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
if (x < 12.0) {
|
|
|
|
// The algorithm directly approximates gamma over (1,2) and uses
|
|
|
|
// reduction identities to reduce other arguments to this interval.
|
2012-07-04 12:16:48 -07:00
|
|
|
|
2011-10-31 15:27:17 -07:00
|
|
|
double y = x;
|
2019-01-02 08:59:07 -08:00
|
|
|
int n = 0;
|
|
|
|
bool arg_was_less_than_one = (y < 1.0);
|
|
|
|
|
|
|
|
// Add or subtract integers as necessary to bring y into (1,2)
|
|
|
|
// Will correct for this below
|
|
|
|
if (arg_was_less_than_one) {
|
|
|
|
y += 1.0;
|
|
|
|
} else {
|
|
|
|
n = static_cast<int>(floor(y)) - 1; // will use n later
|
|
|
|
y -= n;
|
|
|
|
}
|
|
|
|
|
|
|
|
// numerator coefficients for approximation over the interval (1,2)
|
2019-01-21 01:19:07 -08:00
|
|
|
static const double p[] = {
|
|
|
|
-1.71618513886549492533811E+0,
|
|
|
|
2.47656508055759199108314E+1,
|
|
|
|
-3.79804256470945635097577E+2,
|
|
|
|
6.29331155312818442661052E+2,
|
|
|
|
8.66966202790413211295064E+2,
|
|
|
|
-3.14512729688483675254357E+4,
|
|
|
|
-3.61444134186911729807069E+4,
|
|
|
|
6.64561438202405440627855E+4
|
|
|
|
};
|
2019-01-02 08:59:07 -08:00
|
|
|
|
|
|
|
// denominator coefficients for approximation over the interval (1,2)
|
2019-01-21 01:19:07 -08:00
|
|
|
static const double q[] = {
|
|
|
|
-3.08402300119738975254353E+1,
|
|
|
|
3.15350626979604161529144E+2,
|
|
|
|
-1.01515636749021914166146E+3,
|
|
|
|
-3.10777167157231109440444E+3,
|
|
|
|
2.25381184209801510330112E+4,
|
|
|
|
4.75584627752788110767815E+3,
|
|
|
|
-1.34659959864969306392456E+5,
|
|
|
|
-1.15132259675553483497211E+5
|
|
|
|
};
|
2019-01-02 08:59:07 -08:00
|
|
|
|
|
|
|
double num = 0.0;
|
|
|
|
double den = 1.0;
|
|
|
|
int i;
|
|
|
|
|
|
|
|
double z = y - 1;
|
|
|
|
for (i = 0; i < 8; i++) {
|
|
|
|
num = (num + p[i]) * z;
|
|
|
|
den = den * z + q[i];
|
|
|
|
}
|
|
|
|
double result = num / den + 1.0;
|
|
|
|
|
|
|
|
// Apply correction if argument was not initially in (1,2)
|
|
|
|
if (arg_was_less_than_one) {
|
|
|
|
// Use identity gamma(z) = gamma(z+1)/z
|
|
|
|
// The variable "result" now holds gamma of the original y + 1
|
|
|
|
// Thus we use y-1 to get back the orginal y.
|
|
|
|
result /= (y - 1.0);
|
|
|
|
} else {
|
|
|
|
// Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
|
|
|
|
for (i = 0; i < n; i++)
|
|
|
|
result *= y++;
|
|
|
|
}
|
2011-10-31 15:27:17 -07:00
|
|
|
|
|
|
|
return result;
|
2019-01-02 08:59:07 -08:00
|
|
|
}
|
2011-10-31 15:27:17 -07:00
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
///////////////////////////////////////////////////////////////////////////
|
|
|
|
// Third interval: [12, infinity)
|
2011-10-31 15:27:17 -07:00
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
if (x > 171.624) {
|
2011-10-31 15:27:17 -07:00
|
|
|
// Correct answer too large to display. Force +infinity.
|
|
|
|
double temp = DBL_MAX;
|
2019-01-02 08:59:07 -08:00
|
|
|
return temp * 2.0;
|
|
|
|
}
|
2011-10-31 15:27:17 -07:00
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
return exp(LogGamma(x));
|
2011-10-31 15:27:17 -07:00
|
|
|
}
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
double LogGamma(
|
|
|
|
double x // x must be positive
|
2011-10-31 15:27:17 -07:00
|
|
|
)
|
|
|
|
{
|
2019-01-02 08:59:07 -08:00
|
|
|
if (x <= 0.0) {
|
2011-10-31 15:27:17 -07:00
|
|
|
std::stringstream os;
|
2019-01-02 08:59:07 -08:00
|
|
|
os << "Invalid input argument " << x << ". Argument must be positive.";
|
|
|
|
throw std::invalid_argument(os.str());
|
2011-10-31 15:27:17 -07:00
|
|
|
}
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
if (x < 12.0) {
|
|
|
|
return log(fabs(Gamma(x)));
|
|
|
|
}
|
2011-10-31 15:27:17 -07:00
|
|
|
|
|
|
|
// Abramowitz and Stegun 6.1.41
|
2019-01-02 08:59:07 -08:00
|
|
|
// Asymptotic series should be good to at least 11 or 12 figures
|
|
|
|
// For error analysis, see Whittiker and Watson
|
|
|
|
// A Course in Modern Analysis (1927), page 252
|
|
|
|
|
2019-01-21 01:19:07 -08:00
|
|
|
static const double c[8] = {
|
|
|
|
1.0 / 12.0,
|
|
|
|
-1.0 / 360.0,
|
|
|
|
1.0 / 1260.0,
|
|
|
|
-1.0 / 1680.0,
|
|
|
|
1.0 / 1188.0,
|
|
|
|
-691.0 / 360360.0,
|
|
|
|
1.0 / 156.0,
|
|
|
|
-3617.0 / 122400.0
|
|
|
|
};
|
2019-01-02 08:59:07 -08:00
|
|
|
double z = 1.0 / (x * x);
|
|
|
|
double sum = c[7];
|
|
|
|
for (int i = 6; i >= 0; i--) {
|
|
|
|
sum *= z;
|
|
|
|
sum += c[i];
|
|
|
|
}
|
|
|
|
double series = sum / x;
|
|
|
|
|
|
|
|
static const double halfLogTwoPi = 0.91893853320467274178032973640562;
|
|
|
|
double logGamma = (x - 0.5) * log(x) - x + halfLogTwoPi + series;
|
2011-10-31 15:27:17 -07:00
|
|
|
return logGamma;
|
|
|
|
}
|
|
|
|
|
|
|
|
// Can delete these functions and the #include for <iostream> if not testing.
|
|
|
|
|
|
|
|
void TestGamma()
|
|
|
|
{
|
2019-01-02 08:59:07 -08:00
|
|
|
struct TestCase {
|
2011-10-31 15:27:17 -07:00
|
|
|
double input;
|
|
|
|
double expected;
|
|
|
|
};
|
|
|
|
|
2019-01-21 01:19:07 -08:00
|
|
|
TestCase test[] = {
|
|
|
|
// Test near branches in code for (0, 0.001), [0.001, 12), (12, infinity)
|
|
|
|
{ 1e-20, 1e+20 },
|
|
|
|
{ 2.19824158876e-16, 4.5490905327e+15 }, // 0.99*DBL_EPSILON
|
|
|
|
{ 2.24265050974e-16, 4.45900953205e+15 }, // 1.01*DBL_EPSILON
|
|
|
|
{ 0.00099, 1009.52477271 },
|
|
|
|
{ 0.00100, 999.423772485 },
|
|
|
|
{ 0.00101, 989.522792258 },
|
|
|
|
{ 6.1, 142.451944066 },
|
|
|
|
{ 11.999, 39819417.4793 },
|
|
|
|
{ 12, 39916800.0 },
|
|
|
|
{ 12.001, 40014424.1571 },
|
|
|
|
{ 15.2, 149037380723.0 }
|
|
|
|
};
|
2011-10-31 15:27:17 -07:00
|
|
|
|
|
|
|
size_t numTests = sizeof(test) / sizeof(TestCase);
|
|
|
|
|
|
|
|
double worst_absolute_error = 0.0;
|
|
|
|
double worst_relative_error = 0.0;
|
|
|
|
size_t worst_absolute_error_case = 0;
|
|
|
|
size_t worst_relative_error_case = 0;
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
for (size_t t = 0; t < numTests; t++) {
|
|
|
|
double computed = Gamma(test[t].input);
|
2011-10-31 15:27:17 -07:00
|
|
|
double absolute_error = fabs(computed - test[t].expected);
|
|
|
|
double relative_error = absolute_error / test[t].expected;
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
if (absolute_error > worst_absolute_error) {
|
2011-10-31 15:27:17 -07:00
|
|
|
worst_absolute_error = absolute_error;
|
|
|
|
worst_absolute_error_case = t;
|
|
|
|
}
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
if (relative_error > worst_relative_error) {
|
2011-10-31 15:27:17 -07:00
|
|
|
worst_relative_error = absolute_error;
|
|
|
|
worst_relative_error_case = t;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
size_t t = worst_absolute_error_case;
|
|
|
|
double x = test[t].input;
|
|
|
|
double y = test[t].expected;
|
|
|
|
std::cout << "Worst absolute error: "
|
|
|
|
<< fabs(Gamma(x) - y)
|
|
|
|
<< "\nGamma( "
|
2019-01-02 08:59:07 -08:00
|
|
|
<< x
|
2012-07-04 12:16:48 -07:00
|
|
|
<< ") computed as "
|
|
|
|
<< Gamma(x)
|
2011-10-31 15:27:17 -07:00
|
|
|
<< " but exact value is "
|
|
|
|
<< y
|
|
|
|
<< "\n";
|
|
|
|
|
|
|
|
t = worst_relative_error_case;
|
|
|
|
x = test[t].input;
|
|
|
|
y = test[t].expected;
|
|
|
|
std::cout << "Worst relative error: "
|
|
|
|
<< (Gamma(x) - y) / y
|
2019-01-02 08:59:07 -08:00
|
|
|
<< "\nGamma( "
|
|
|
|
<< x
|
2012-07-04 12:16:48 -07:00
|
|
|
<< ") computed as "
|
|
|
|
<< Gamma(x)
|
2011-10-31 15:27:17 -07:00
|
|
|
<< " but exact value is "
|
|
|
|
<< y
|
|
|
|
<< "\n";
|
|
|
|
}
|
|
|
|
|
|
|
|
void TestLogGamma()
|
|
|
|
{
|
2019-01-02 08:59:07 -08:00
|
|
|
struct TestCase {
|
2011-10-31 15:27:17 -07:00
|
|
|
double input;
|
|
|
|
double expected;
|
|
|
|
};
|
|
|
|
|
2019-01-21 01:19:07 -08:00
|
|
|
TestCase test[] = {
|
|
|
|
{ 1e-12, 27.6310211159 },
|
|
|
|
{ 0.9999, 5.77297915613e-05 },
|
|
|
|
{ 1.0001, -5.77133422205e-05 },
|
|
|
|
{ 3.1, 0.787375083274 },
|
|
|
|
{ 6.3, 5.30734288962 },
|
|
|
|
{ 11.9999, 17.5020635801 },
|
|
|
|
{ 12, 17.5023078459 },
|
|
|
|
{ 12.0001, 17.5025521125 },
|
|
|
|
{ 27.4, 62.5755868211 }
|
|
|
|
};
|
2011-10-31 15:27:17 -07:00
|
|
|
|
|
|
|
size_t numTests = sizeof(test) / sizeof(TestCase);
|
|
|
|
|
|
|
|
double worst_absolute_error = 0.0;
|
|
|
|
double worst_relative_error = 0.0;
|
|
|
|
size_t worst_absolute_error_case = 0;
|
|
|
|
size_t worst_relative_error_case = 0;
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
for (size_t t = 0; t < numTests; t++) {
|
|
|
|
double computed = LogGamma(test[t].input);
|
2011-10-31 15:27:17 -07:00
|
|
|
double absolute_error = fabs(computed - test[t].expected);
|
|
|
|
double relative_error = absolute_error / test[t].expected;
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
if (absolute_error > worst_absolute_error) {
|
2011-10-31 15:27:17 -07:00
|
|
|
worst_absolute_error = absolute_error;
|
|
|
|
worst_absolute_error_case = t;
|
|
|
|
}
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
if (relative_error > worst_relative_error) {
|
2011-10-31 15:27:17 -07:00
|
|
|
worst_relative_error = absolute_error;
|
|
|
|
worst_relative_error_case = t;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
size_t t = worst_absolute_error_case;
|
|
|
|
double x = test[t].input;
|
|
|
|
double y = test[t].expected;
|
|
|
|
std::cout << "Worst absolute error: "
|
|
|
|
<< fabs(LogGamma(x) - y)
|
|
|
|
<< "\nGamma( "
|
2019-01-02 08:59:07 -08:00
|
|
|
<< x
|
2012-07-04 12:16:48 -07:00
|
|
|
<< ") computed as "
|
|
|
|
<< LogGamma(x)
|
2011-10-31 15:27:17 -07:00
|
|
|
<< " but exact value is "
|
|
|
|
<< y
|
|
|
|
<< "\n";
|
|
|
|
|
|
|
|
t = worst_relative_error_case;
|
|
|
|
x = test[t].input;
|
|
|
|
y = test[t].expected;
|
|
|
|
std::cout << "Worst relative error: "
|
|
|
|
<< (LogGamma(x) - y) / y
|
2019-01-02 08:59:07 -08:00
|
|
|
<< "\nGamma( "
|
|
|
|
<< x
|
2012-07-04 12:16:48 -07:00
|
|
|
<< ") computed as "
|
|
|
|
<< LogGamma(x)
|
2011-10-31 15:27:17 -07:00
|
|
|
<< " but exact value is "
|
|
|
|
<< y
|
|
|
|
<< "\n";
|
|
|
|
}
|
|
|
|
|
2019-01-02 08:59:07 -08:00
|
|
|
#if defined(_MSC_VER) && (_MSC_VER <= 1700)
|
2013-03-30 06:54:20 -07:00
|
|
|
// http://social.msdn.microsoft.com/Forums/en-US/Vsexpressvc/thread/25c923af-a824-40f8-8fd4-e5574bc147af/
|
2019-01-02 08:59:07 -08:00
|
|
|
double asinh(double value)
|
|
|
|
{
|
2013-03-30 06:54:20 -07:00
|
|
|
double returned;
|
2019-01-02 08:59:07 -08:00
|
|
|
if (value > 0.0)
|
2013-03-30 06:54:20 -07:00
|
|
|
returned = log(value + sqrt(value * value + 1.0));
|
|
|
|
else
|
|
|
|
returned = -log(-value + sqrt(value * value + 1.0));
|
|
|
|
return returned;
|
|
|
|
}
|
|
|
|
|
|
|
|
// http://stackoverflow.com/questions/15539116/atanh-arc-hyperbolic-tangent-function-missing-in-ms-visual-c
|
2019-01-02 08:59:07 -08:00
|
|
|
double atanh(double x) //implements: return (log(1+x) - log(1-x))/2
|
2013-03-30 06:54:20 -07:00
|
|
|
{
|
2019-01-02 08:59:07 -08:00
|
|
|
return (LogOnePlusX(x) - LogOnePlusX(-x)) / 2.0;
|
2013-03-30 06:54:20 -07:00
|
|
|
}
|
|
|
|
#endif
|