dc8c814ca3
Fixed coding standard error in several files. Use of while( is incorrect; a space is required between while and (. Also ran tools/nxstyle and fix thoses complaints as well in most files.
Changes to comply with coding standard. Mostly focused on files with missing space after keyword in if(, switch(, and for(. Offending files also got changes to comply with tools nxstyle. If there were logs of nxstyle complaints, the file also got a taste of tools/indent.sh. Still need to fix occurrences of while( with missing space. There are a lot of them.
324 lines
8.7 KiB
C
324 lines
8.7 KiB
C
/****************************************************************************
|
|
* 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 <nuttx/config.h>
|
|
#include <nuttx/compiler.h>
|
|
|
|
#include <sys/types.h>
|
|
#include <math.h>
|
|
|
|
#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
|