dgamma.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/* dgamma.f -- translated by f2c (version 20041007).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

/** This routine has been editted to be thread safe **/

#define V3P_NETLIB_SRC
#include "v3p_netlib.h"

/* Table of constant values */

static integer c__3 = 3;
static integer c__42 = 42;
static integer c__4 = 4;
static integer c__2 = 2;
static integer c__1 = 1;

/* DECK DGAMMA */
doublereal dgamma_(doublereal *x)
{
    /* Initialized data */

    static doublereal gamcs[42] = { .008571195590989331421920062399942,
            .004415381324841006757191315771652,
            .05685043681599363378632664588789,
            -.004219835396418560501012500186624,
            .001326808181212460220584006796352,
            -1.893024529798880432523947023886e-4,
            3.606925327441245256578082217225e-5,
            -6.056761904460864218485548290365e-6,
            1.055829546302283344731823509093e-6,
            -1.811967365542384048291855891166e-7,
            3.117724964715322277790254593169e-8,
            -5.354219639019687140874081024347e-9,
            9.19327551985958894688778682594e-10,
            -1.577941280288339761767423273953e-10,
            2.707980622934954543266540433089e-11,
            -4.646818653825730144081661058933e-12,
            7.973350192007419656460767175359e-13,
            -1.368078209830916025799499172309e-13,
            2.347319486563800657233471771688e-14,
            -4.027432614949066932766570534699e-15,
            6.910051747372100912138336975257e-16,
            -1.185584500221992907052387126192e-16,
            2.034148542496373955201026051932e-17,
            -3.490054341717405849274012949108e-18,
            5.987993856485305567135051066026e-19,
            -1.027378057872228074490069778431e-19,
            1.762702816060529824942759660748e-20,
            -3.024320653735306260958772112042e-21,
            5.188914660218397839717833550506e-22,
            -8.902770842456576692449251601066e-23,
            1.527474068493342602274596891306e-23,
            -2.620731256187362900257328332799e-24,
            4.496464047830538670331046570666e-25,
            -7.714712731336877911703901525333e-26,
            1.323635453126044036486572714666e-26,
            -2.270999412942928816702313813333e-27,
            3.896418998003991449320816639999e-28,
            -6.685198115125953327792127999999e-29,
            1.146998663140024384347613866666e-29,
            -1.967938586345134677295103999999e-30,
            3.376448816585338090334890666666e-31,
            -5.793070335782135784625493333333e-32 };
    static doublereal pi = 3.1415926535897932384626433832795;
    static doublereal sq2pil = .91893853320467274178032973640562;
    // static logical first = TRUE_;

    /* System generated locals */
    integer i__1;
    real r__1;
    doublereal ret_val, d__1, d__2;

    /* Builtin functions */
    double sqrt(doublereal), d_int(doublereal *), log(doublereal), exp(
            doublereal), sin(doublereal);

    /* Local variables */
    integer i__, n;
    doublereal y;
    /* static */ integer ngam;
    /* static */ doublereal xmin, xmax, dxrel;
    extern doublereal d1mach_(integer *), d9lgmc_(doublereal *);
    extern /* Subroutine */ int dgamlm_(doublereal *, doublereal *);
    extern doublereal dcsevl_(doublereal *, doublereal *, integer *);
    extern integer initds_(doublereal *, integer *, real *);
    extern /* Subroutine */ int xermsg_(char *, char *, char *, integer *, 
            integer *, ftnlen, ftnlen, ftnlen);
    doublereal sinpiy;

/* ***BEGIN PROLOGUE  DGAMMA */
/* ***PURPOSE  Compute the complete Gamma function. */
/* ***LIBRARY   SLATEC (FNLIB) */
/* ***CATEGORY  C7A */
/* ***TYPE      DOUBLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C) */
/* ***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS */
/* ***AUTHOR  Fullerton, W., (LANL) */
/* ***DESCRIPTION */

/* DGAMMA(X) calculates the double precision complete Gamma function */
/* for double precision argument X. */

/* Series for GAM        on the interval  0.          to  1.00000E+00 */
/*                                        with weighted error   5.79E-32 */
/*                                         log weighted error  31.24 */
/*                               significant figures required  30.00 */
/*                                    decimal places required  32.05 */

/* ***REFERENCES  (NONE) */
/* ***ROUTINES CALLED  D1MACH, D9LGMC, DCSEVL, DGAMLM, INITDS, XERMSG */
/* ***REVISION HISTORY  (YYMMDD) */
/*   770601  DATE WRITTEN */
/*   890531  Changed all specific intrinsics to generic.  (WRB) */
/*   890911  Removed unnecessary intrinsics.  (WRB) */
/*   890911  REVISION DATE from Version 3.2 */
/*   891214  Prologue converted to Version 4.0 format.  (BAB) */
/*   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ) */
/*   920618  Removed space from variable name.  (RWC, WRB) */
/* ***END PROLOGUE  DGAMMA */

/* ***FIRST EXECUTABLE STATEMENT  DGAMMA */

    // d1mach has been made thread safe, so there is no need for the
    // statics in determining these values
//     if (first) {
//      r__1 = (real) d1mach_(&c__3) * .1f;
//      ngam = initds_(gamcs, &c__42, &r__1);

//      dgamlm_(&xmin, &xmax);
//      dxrel = sqrt(d1mach_(&c__4));
//     }
//     first = FALSE_;
    r__1 = (real) d1mach_(&c__3) * .1f;
    ngam = initds_(gamcs, &c__42, &r__1);
    dgamlm_(&xmin, &xmax);
    dxrel = sqrt(d1mach_(&c__4));

    y = abs(*x);
    if (y > 10.) {
        goto L50;
    }

/* COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND.  REDUCE INTERVAL AND FIND */
/* GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL. */

    n = (integer) (*x);
    if (*x < 0.) {
        --n;
    }
    y = *x - n;
    --n;
    d__1 = y * 2. - 1.;
    ret_val = dcsevl_(&d__1, gamcs, &ngam) + .9375;
    if (n == 0) {
        return ret_val;
    }

    if (n > 0) {
        goto L30;
    }

/* COMPUTE GAMMA(X) FOR X .LT. 1.0 */

    n = -n;
    if (*x == 0.) {
        xermsg_("SLATEC", "DGAMMA", "X IS 0", &c__4, &c__2, (ftnlen)6, (
                ftnlen)6, (ftnlen)6);
    }
    if (*x < 0.f && *x + n - 2 == 0.) {
        xermsg_("SLATEC", "DGAMMA", "X IS A NEGATIVE INTEGER", &c__4, &c__2, (
                ftnlen)6, (ftnlen)6, (ftnlen)23);
    }
    d__2 = *x - .5;
    if (*x < -.5 && (d__1 = (*x - d_int(&d__2)) / *x, abs(d__1)) < dxrel) {
        xermsg_("SLATEC", "DGAMMA", "ANSWER LT HALF PRECISION BECAUSE X TOO "
                "NEAR NEGATIVE INTEGER", &c__1, &c__1, (ftnlen)6, (ftnlen)6, (
                ftnlen)60);
    }

    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        ret_val /= *x + i__ - 1;
/* L20: */
    }
    return ret_val;

/* GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0 */

L30:
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        ret_val = (y + i__) * ret_val;
/* L40: */
    }
    return ret_val;

/* GAMMA(X) FOR ABS(X) .GT. 10.0.  RECALL Y = ABS(X). */

L50:
    if (*x > xmax) {
        xermsg_("SLATEC", "DGAMMA", "X SO BIG GAMMA OVERFLOWS", &c__3, &c__2, 
                (ftnlen)6, (ftnlen)6, (ftnlen)24);
    }

    ret_val = 0.;
    if (*x < xmin) {
        xermsg_("SLATEC", "DGAMMA", "X SO SMALL GAMMA UNDERFLOWS", &c__2, &
                c__1, (ftnlen)6, (ftnlen)6, (ftnlen)27);
    }
    if (*x < xmin) {
        return ret_val;
    }

    ret_val = exp((y - .5) * log(y) - y + sq2pil + d9lgmc_(&y));
    if (*x > 0.) {
        return ret_val;
    }

    d__2 = *x - .5;
    if ((d__1 = (*x - d_int(&d__2)) / *x, abs(d__1)) < dxrel) {
        xermsg_("SLATEC", "DGAMMA", "ANSWER LT HALF PRECISION, X TOO NEAR NE"
                "GATIVE INTEGER", &c__1, &c__1, (ftnlen)6, (ftnlen)6, (ftnlen)
                53);
    }

    sinpiy = sin(pi * y);
    if (sinpiy == 0.) {
        xermsg_("SLATEC", "DGAMMA", "X IS A NEGATIVE INTEGER", &c__4, &c__2, (
                ftnlen)6, (ftnlen)6, (ftnlen)23);
    }

    ret_val = -pi / (y * sinpiy * ret_val);

    return ret_val;
} /* dgamma_ */