source: sasmodels/sasmodels/models/lib/sas_gamma.c @ e07b470

/*
The wrapper for gamma function from OpenCL and standard libraries
The OpenCL gamma function fails misserably 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)
*/

#ifdef tgamma
inline double sas_gamma( double x) { return tgamma(x+1)/x; }
#else
static double cephes_stirf(double x)
{
#define MAXSTIR 143.01608
#define 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);
}

double sas_gamma(double x) {
#define MAXGAM 171.624376956302725
#define LOGPI 1.14472988584940017414
        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