source: sasmodels/sasmodels/models/lib/sas_gamma.c @ 1596de3

/*
The wrapper for gamma function from OpenCL and standard libraries
The OpenCL gamma function fails miserably on values lower than 1.0
while works fine on larger values.
We use gamma definition Gamma(t + 1) = t * Gamma(t) to compute
to function for values lower than 1.0. Namely Gamma(t) = 1/t * Gamma(t + 1)
*/

#if defined(NEED_TGAMMA)
static double cephes_stirf(double x)
{
        const double MAXSTIR=143.01608;
        const double SQTPI=2.50662827463100050242E0;
        double y, w, v;

        w = 1.0 / x;

        w = ((((
                7.87311395793093628397E-4*w
                - 2.29549961613378126380E-4)*w
                - 2.68132617805781232825E-3)*w
                + 3.47222221605458667310E-3)*w
                + 8.33333333333482257126E-2)*w
                + 1.0;
        y = exp(x);
        if (x > MAXSTIR)
        { /* Avoid overflow in pow() */
                v = pow(x, 0.5 * x - 0.25);
                y = v * (v / y);
        }
        else
        {
                y = pow(x, x - 0.5) / y;
        }
        y = SQTPI * y * w;
        return(y);
}

static double tgamma(double x) {
        double p, q, z;
        int sgngam;
        int i;

        sgngam = 1;
        if (isnan(x))
                return(x);
        q = fabs(x);

        if (q > 33.0)
        {
                if (x < 0.0)
                {
                        p = floor(q);
                        if (p == q)
                        {
                                return (NAN);
                        }
                        i = p;
                        if ((i & 1) == 0)
                                sgngam = -1;
                        z = q - p;
                        if (z > 0.5)
                        {
                                p += 1.0;
                                z = q - p;
                        }
                        z = q * sin(M_PI * z);
                        if (z == 0.0)
                        {
                                return(NAN);
                        }
                        z = fabs(z);
                        z = M_PI / (z * cephes_stirf(q));
                }
                else
                {
                        z = cephes_stirf(x);
                }
                return(sgngam * z);
        }

        z = 1.0;
        while (x >= 3.0)
        {
                x -= 1.0;
                z *= x;
        }

        while (x < 0.0)
        {
                if (x > -1.E-9)
                        goto small;
                z /= x;
                x += 1.0;
        }

        while (x < 2.0)
        {
                if (x < 1.e-9)
                        goto small;
                z /= x;
                x += 1.0;
        }

        if (x == 2.0)
                return(z);

        x -= 2.0;
        p = (((((
                +1.60119522476751861407E-4*x
                + 1.19135147006586384913E-3)*x
                + 1.04213797561761569935E-2)*x
                + 4.76367800457137231464E-2)*x
                + 2.07448227648435975150E-1)*x
                + 4.94214826801497100753E-1)*x
                + 9.99999999999999996796E-1;
        q = ((((((
                -2.31581873324120129819E-5*x
                + 5.39605580493303397842E-4)*x
                - 4.45641913851797240494E-3)*x
                + 1.18139785222060435552E-2)*x
                + 3.58236398605498653373E-2)*x
                - 2.34591795718243348568E-1)*x
                + 7.14304917030273074085E-2)*x
                + 1.00000000000000000320E0;
        return(z * p / q);

small:
        if (x == 0.0)
        {
                return (NAN);
        }
        else
                return(z / ((1.0 + 0.5772156649015329 * x) * x));
}
#endif // NEED_TGAMMA


inline double sas_gamma( double x) { return tgamma(x+1)/x; }