slatec-bessel-cpp/f2c/dgamma.cpp

/* dgamma.f -- translated by f2c (version 20100827).
   This file no longer depends on f2c.
*/

#include "slatec-internal.hpp"
#include <tuple>

/* Table of constant values */

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

    /* Initialized data */

static double const 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 double const pi = 3.1415926535897932384626433832795;
static double const sq2pil = .91893853320467274178032973640562;

static float const r__1 = (float) d1mach_(3) * (float).1;
static integer const ngam = initds_(gamcs, &c__42, &r__1);
static double const dxrel = sqrt(d1mach_(4));
static std::tuple<double, double> const xboth = ([]() { double a, b; dgamlm_(&a, &b); return std::make_tuple(a, b); })();
static double const xmin = std::get<0>(xboth);
static double const xmax = std::get<1>(xboth);

double dgamma_(double const *x)
{
    /* System generated locals */
    integer i__1;
    double ret_val, d__1, d__2;

    /* Local variables */
    integer i__, n;
    double y;
    double 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 */
    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 < (float)0. && *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 - f2c::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 - f2c::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_ */