/**************************************************************************** * libs/libc/math/lib_gamma.c * * Ported to NuttX from FreeBSD by Alan Carvalho de Assis: * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in * the documentation and/or other materials provided with the * distribution. * 3. Neither the name NuttX nor the names of its contributors may be * used to endorse or promote products derived from this software * without specific prior written permission. * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. * ****************************************************************************/ /* "A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964) * "Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001) * "An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004) * * Approximation method: * * (x - 0.5) S(x) * Gamma(x) = (x + g - 0.5) * ---------------- * exp(x + g - 0.5) * * with * a1 a2 a3 aN * S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ] * x + 1 x + 2 x + 3 x + N * * with a0, a1, a2, a3,.. aN constants which depend on g. * * for x < 0 the following reflection formula is used: * * Gamma(x)*Gamma(-x) = -pi/(x sin(pi x)) * * most ideas and constants are from boost and python */ /**************************************************************************** * Included Files ****************************************************************************/ #include #include #include #include #ifdef CONFIG_HAVE_DOUBLE /**************************************************************************** * Pre-processor Definitions ****************************************************************************/ #define FORCE_EVAL(x) \ do \ { \ if (sizeof(x) == sizeof(float)) \ { \ volatile float __x; \ UNUSED(__x); \ __x = (x); \ } \ else if (sizeof(x) == sizeof(double)) \ { \ volatile double __x; \ UNUSED(__x); \ __x = (x); \ } \ else \ { \ volatile long double __x; \ UNUSED(__x); \ __x = (x); \ } \ } \ while (0) #define N 12 /**************************************************************************** * Private Data ****************************************************************************/ static const double pi = 3.141592653589793238462643383279502884; static const double g_gmhalf = 5.524680040776729583740234375; static const double g_snum[N + 1] = { 23531376880.410759688572007674451636754734846804940, 42919803642.649098768957899047001988850926355848959, 35711959237.355668049440185451547166705960488635843, 17921034426.037209699919755754458931112671403265390, 6039542586.3520280050642916443072979210699388420708, 1439720407.3117216736632230727949123939715485786772, 248874557.86205415651146038641322942321632125127801, 31426415.585400194380614231628318205362874684987640, 2876370.6289353724412254090516208496135991145378768, 186056.26539522349504029498971604569928220784236328, 8071.6720023658162106380029022722506138218516325024, 210.82427775157934587250973392071336271166969580291, 2.5066282746310002701649081771338373386264310793408, }; static const double g_sden[N + 1] = { 0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535, 2637558, 357423, 32670, 1925, 66, 1, }; /* n! for small integer n */ static const double g_fact[] = { 1, 1, 2, 6, 24, 120, 720, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0, 20922789888000.0, 355687428096000.0, 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0, }; /* S(x) rational function for positive x */ /**************************************************************************** * Private Functions ****************************************************************************/ /* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */ static double sinpi(double x) { int n; /* argument reduction: x = |x| mod 2 */ /* spurious inexact when x is odd int */ x = x * 0.5; x = 2 * (x - floor(x)); /* reduce x into [-.25,.25] */ n = 4 * x; n = (n + 1) / 2; x -= n * 0.5; x *= pi; switch (n) { default: /* case 4 */ case 0: return __sin(x, 0, 0); case 1: return __cos(x, 0); case 2: return __sin(-x, 0, 0); case 3: return -__cos(x, 0); } } static double s(double x) { double num = 0; double den = 0; int i; /* to avoid overflow handle large x differently */ if (x < 8) { for (i = N; i >= 0; i--) { num = num * x + g_snum[i]; den = den * x + g_sden[i]; } } else { for (i = 0; i <= N; i++) { num = num / x + g_snum[i]; den = den / x + g_sden[i]; } } return num/den; } /**************************************************************************** * Public Functions ****************************************************************************/ double tgamma(double x) { union { double f; uint64_t i; } u; u.f = x; double absx; double y; double dy; double z; double r; uint32_t ix = u.i >> 32 & 0x7fffffff; int sign = u.i >> 63; /* special cases */ if (ix >= 0x7ff00000) { /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */ return x + INFINITY; } if (ix < (0x3ff - 54) << 20) { /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */ return 1 / x; } /* integer arguments */ /* raise inexact when non-integer */ if (x == floor(x)) { if (sign) { return 0 / 0.0; } if (x <= sizeof g_fact / sizeof *g_fact) { return g_fact[(int)x - 1]; } } /* x >= 172: tgamma(x)=inf with overflow */ /* x =< -184: tgamma(x)=+-0 with underflow */ if (ix >= 0x40670000) { /* |x| >= 184 */ if (sign) { FORCE_EVAL((float)(0x1p-126 / x)); if (floor(x) * 0.5 == floor(x * 0.5)) { return 0; } return -0.0; } x *= 0x1p1023; return x; } absx = sign ? -x : x; /* handle the error of x + g - 0.5 */ y = absx + g_gmhalf; if (absx > g_gmhalf) { dy = y - absx; dy -= g_gmhalf; } else { dy = y - g_gmhalf; dy -= absx; } z = absx - 0.5; r = s(absx) * exp(-y); if (x < 0) { /* reflection formula for negative x */ /* sinpi(absx) is not 0, integers are already handled */ r = -pi / (sinpi(absx) * absx * r); dy = -dy; z = -z; } r += dy * (g_gmhalf + 0.5) * r / y; z = pow(y, 0.5 * z); y = r * z * z; return y; } double gamma(double x) { return tgamma(x); } #endif