/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:46 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "sgamma.h"
#include <float.h>
#include <stdlib.h>
/* PARAMETER translations */
#define C1 0.9189385332046727417803297e0
#define HALF 0.5e0
#define ONE 1.0e0
#define PI 3.1415926535897932384626434e0
#define TWELVE 12.0e0
#define TWO 2.0e0
#define ZERO 0.0e0
/* end of PARAMETER translations */
float /*FUNCTION*/ sgamma(
float x)
{
LOGICAL32 parity;
long int _l0, _l1, i, j, n;
float const_, del, f, fact, fp, res, sgamma_v, sum,
temp, x1, x2, xden, xnum, y, y1, ysq, z;
static float eps, xgbig, xminin;
static float xinf = 0.0e0;
static float p[8]={-1.71618513886549492533811e0,2.47656508055759199108314e1,
-3.79804256470945635097577e2,6.29331155312818442661052e2,8.66966202790413211295064e2,
-3.14512729688483675254357e4,-3.61444134186911729807069e4,6.64561438202405440627855e4};
static float q[8]={-3.08402300119738975254353e1,3.15350626979604161529144e2,
-1.01515636749021914166146e3,-3.10777167157231109440444e3,2.25381184209801510330112e4,
4.75584627752788110767815e3,-1.34659959864969306392456e5,-1.15132259675553483497211e5};
static float c[7]={-1.910444077728e-03,8.4171387781295e-04,-5.952379913043012e-04,
7.93650793500350248e-04,-2.777777777777681622553e-03,8.333333333333333331554247e-02,
5.7083835261e-03};
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const C = &c[0] - 1;
float *const P = &p[0] - 1;
float *const Q = &q[0] - 1;
/* end of OFFSET VECTORS */
/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
* ALL RIGHTS RESERVED.
* Based on Government Sponsored Research NAS7-03001.
*>> 2001-05-25 SGAMMA Krogh Minor change for making .f90 version.
*>> 1996-03-30 SGAMMA Krogh Added external statement.
*>> 1994-10-20 SGAMMA Krogh Changes to use M77CON
*>> 1991-10-21 SGAMMA CLL Eliminated DGAM1 as a separate subroutine.
*>> 1991-01-16 SGAMMA Lawson Replaced D2MACH/R2MACH with DGAM1.
*>> 1985-08-02 SGAMMA Lawson Initial code.
*--S replaces "?": ?GAMMA, ?ERM1, ?ERV1
*
* ----------------------------------------------------------------------
*
* THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A double precision
* ARGUMENT X. PERMITS NEGATIVE AS WELL AS POSITIVE X. NOTE
* THAT THE GAMMA FUNCTION HAS POLES AT ZERO AND AT NEGATIVE
* ARGUMENTS. COMPUTATION IS BASED ON AN ALGORITHN OUTLINED IN
* W.J.CODY, 'AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
* FUNCTIONS', LECTURE NOTES IN MATHEMATICS, 506, NUMERICAL ANALYSIS
* DUNDEE, 1975, G. A. WATSON (ED.),SPRINGER VERLAG, BERLIN, 1976.
* THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
* FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
* FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
* THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM HART, ET. AL.,
* COMPUTER APPROXIMATIONS, WILEY AND SONS, NEW YORK, 1968.
* LOWER ORDER APPROXIMATIONS CAN BE SUBSTITUTED FOR THESE ON
* MACHINES WITH LESS PRECISE ARITHMETIC.
*
* Designed & programmed by W.J.CODY, Argonne National Lab.,1982.
* Minor changes for the JPL library by C.L.LAWSON & S.CHAN,JPL,1983.
*
************************************************************************
*
*-- Begin mask code changes
* EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
*
* EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
* 1.0 + EPS .GT. 1.0
* XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER.
* XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
* both XMININ and 1/XMININ are representable.
* XGBIG - A value such that Gamma(XGBIG) = 0.875 * XINF.
* (Computed and used in [D/S]GAMMA.)
* XLBIG - A value such that LogGamma(XLBIG) = 0.875 * XINF.
* (Computed and used in [D/S]LGAMA.)
*
* Values of XINF, XGBIG, and XLBIG for some machines:
*
* XINF XGBIG XLBIG Machines
*
* 2**127 = 0.170e39 34.81 0.180e37 Vax SP & DP; Unisys SP
* 2**128 = 0.340e39 35.00 0.358e37 IEEE SP
* 2**252 = 0.723e76 57.54 0.376e74 IBM30xx DP
* 2**1023 = 0.899e308 171.46 0.112e306 Unisys DP
* 2**1024 = 0.180e309 171.60 0.2216e306 IEEE DP
* 2**1070 = 0.126e323 177.78 0.1501e320 CDC/7600 SP
* 2**8191 = 0.550e2466 966.94 0.8464e2462 Cray SP & DP
*-- End mask code changes
*
************************************************************************
*
* ERROR RETURNS
*
* THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
* WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
* TO BE FREE OF UNDERFLOW AND OVERFLOW.
*
* AUTHOR:W. J. CODY
* APPLIED MATHMATICS DIVISION
* ARGONNE NATIONAL LABORATORY
* ARGONNE, IL 60439
*
* LATEST MODIFICATION by Cody: MAY 18, 1982
*
* ------------------------------------------------------------------ */
/* C1 = LOG base e of SQRT(2*PI)
* */
/* ----------------------------------------------------------------------
* NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
* APPROXIMATION OVER (1,2).
* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
* COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
* */
if (xinf == ZERO)
{
eps = FLT_EPSILON;
xinf = FLT_MAX;
if (FLT_MIN*FLT_MAX >= ONE)
{
xminin = FLT_MIN;
}
else
{
xminin = ONE/FLT_MAX;
}
/* Compute XGBIG
*
* XGBIG will satisfy Gamma(XGBIG) = 0.875 * XINF.
* Use a Newton iteration and the following approximation for
* the gamma function:
* log(gamma(x)) ~ (x - .5)*log(x) - x + 0.5 * log(2 * PI)
* */
temp = logf( 0.875e0*xinf );
const_ = HALF*logf( TWO*PI ) - temp;
x1 = temp*0.34e0;
for (j = 1; j <= 7; j++)
{
f = (x1 - HALF)*logf( x1 ) - x1 + const_;
fp = ((x1 - HALF)/x1) + logf( x1 ) - ONE;
del = -f/fp;
x2 = x1 + del;
if (fabsf( del ) < 0.5e-5*x2)
goto L_45;
x1 = x2;
}
L_45:
;
xgbig = x2;
}
parity = FALSE;
fact = ONE;
n = 0;
y = x;
if (y > ZERO)
goto L_200;
/* ----------------------------------------------------------------------
* ARGUMENT IS NEGATIVE OR ZERO
* ---------------------------------------------------------------------- */
y = -x;
j = (long)( y );
res = y - (float)( j );
if (res == ZERO)
goto L_700;
if (j != (j/2)*2)
parity = TRUE;
fact = -PI/sinf( PI*res );
y += ONE;
/* ----------------------------------------------------------------------
* ARGUMENT IS POSITIVE
* ---------------------------------------------------------------------- */
L_200:
if (y < eps)
goto L_650;
if (y >= TWELVE)
goto L_300;
y1 = y;
if (y >= ONE)
goto L_210;
/* ----------------------------------------------------------------------
* 0.0 .LT. ARGUMENT .LT. 1.0
* ---------------------------------------------------------------------- */
z = y;
y += ONE;
goto L_250;
/* ----------------------------------------------------------------------
* 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
* ---------------------------------------------------------------------- */
L_210:
n = (long)( y ) - 1;
y -= (float)( n );
z = y - ONE;
/* ----------------------------------------------------------------------
* EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
* ---------------------------------------------------------------------- */
L_250:
xnum = ZERO;
xden = ONE;
for (i = 1; i <= 8; i++)
{
xnum = (xnum + P[i])*z;
xden = xden*z + Q[i];
}
res = (xnum/xden + HALF) + HALF;
if (y == y1)
goto L_900;
if (y1 > y)
goto L_280;
/* ----------------------------------------------------------------------
* ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
* ---------------------------------------------------------------------- */
res /= y1;
goto L_900;
/* ----------------------------------------------------------------------
* ADJUST RESULT FOR CASE 2.0 .LT. 12.0
* ---------------------------------------------------------------------- */
L_280:
for (i = 1; i <= n; i++)
{
res *= y;
y += ONE;
}
goto L_900;
/* ----------------------------------------------------------------------
* EVALUATE FOR ARGUMENT .GE. 12.0
* ---------------------------------------------------------------------- */
L_300:
if (y > xgbig)
goto L_720;
ysq = y*y;
sum = C[7];
for (i = 1; i <= 6; i++)
{
sum = sum/ysq + C[i];
}
sum = ((sum/y + C1) - y) + (y - HALF)*logf( y );
res = expf( sum );
goto L_900;
/* ----------------------------------------------------------------------
* ARGUMENT .LT. EPS
* ---------------------------------------------------------------------- */
L_650:
if (y < xminin)
goto L_740;
res = ONE/y;
goto L_900;
/* ----------------------------------------------------------------------
* RETURN FOR SINGULARITIES,EXTREME ARGUMENTS, ETC.
* ---------------------------------------------------------------------- */
L_700:
serm1( "SGAMMA", 1, 0, "POLE AT 0 AND NEG INTEGERS", "X", x, '.' );
goto L_780;
L_720:
serm1( "SGAMMA", 2, 0, "X SO LARGE VALUE WOULD OVERFLOW", "X",
x, ',' );
serv1( "LIMIT", xgbig, '.' );
goto L_780;
L_740:
serm1( "SGAMMA", 3, 0, "X TOO NEAR TO A SINGULARITY.VALUE WOULD OVERFLOW."
, "X", x, '.' );
L_780:
sgamma_v = xinf;
goto L_950;
/* ----------------------------------------------------------------------
* FINAL ADJUSTMENTS AND RETURN
* ---------------------------------------------------------------------- */
L_900:
if (parity)
res = -res;
if (fact != ONE)
res = fact/res;
sgamma_v = res;
L_950:
return( sgamma_v );
} /* end of function */