Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Fetching contributors…

Cannot retrieve contributors at this time

file 68 lines (60 sloc) 2.103 kb
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
/* lgamma_r.c - public domain implementation of error function lgamma_r(3m)

lgamma_r() is based on gamma(). modified by Tanaka Akira.

reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
(New Algorithm handbook in C language) (Gijyutsu hyouron
sha, Tokyo, 1991) [in Japanese]
http://oku.edu.mie-u.ac.jp/~okumura/algo/
*/

/***********************************************************
gamma.c -- Gamma function
***********************************************************/
#include <math.h>
#include <errno.h>
#define PI 3.14159265358979324 /* $\pi$ */
#define LOG_2PI 1.83787706640934548 /* $\log 2\pi$ */
#define LOG_PI 1.14472988584940017 /* $\log_e \pi$ */
#define N 8

#define B0 1 /* Bernoulli numbers */
#define B1 (-1.0 / 2.0)
#define B2 ( 1.0 / 6.0)
#define B4 (-1.0 / 30.0)
#define B6 ( 1.0 / 42.0)
#define B8 (-1.0 / 30.0)
#define B10 ( 5.0 / 66.0)
#define B12 (-691.0 / 2730.0)
#define B14 ( 7.0 / 6.0)
#define B16 (-3617.0 / 510.0)

static double
loggamma(double x) /* the natural logarithm of the Gamma function. */
{
    double v, w;

    if (x == 1.0 || x == 2.0) return 0.0;

    v = 1;
    while (x < N) { v *= x; x++; }
    w = 1 / (x * x);
    return ((((((((B16 / (16 * 15)) * w + (B14 / (14 * 13))) * w
                + (B12 / (12 * 11))) * w + (B10 / (10 * 9))) * w
                + (B8 / ( 8 * 7))) * w + (B6 / ( 6 * 5))) * w
                + (B4 / ( 4 * 3))) * w + (B2 / ( 2 * 1))) / x
                + 0.5 * LOG_2PI - log(v) - x + (x - 0.5) * log(x);
}

/* the natural logarithm of the absolute value of the Gamma function */
double
lgamma_r(double x, int *signp)
{
    if (x <= 0) {
        double i, f, s;
        f = modf(-x, &i);
        if (f == 0.0) { /* pole error */
            *signp = 1;
            errno = ERANGE;
            return HUGE_VAL;
        }
        *signp = (fmod(i, 2.0) != 0.0) ? 1 : -1;
        s = sin(PI * f);
        if (s < 0) s = -s;
        return LOG_PI - log(s) - loggamma(1 - x);
    }
    *signp = 1;
    return loggamma(x);
}
Something went wrong with that request. Please try again.