gamma problems

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

gamma problems

L. Jonas Olsson
 I'm trying to get it to compile on 386BSD, its standard libm.a
does not have any gamma functions at all. But I also have the
GNU mathlib-2.1 installed as /usr/local/lib/libm.a, this library
comes with both gamma() and lgamma(), the problem is that this
lgamma() doesn't use the external variable signgam as it should.
 As this library is used by manu 386 Unices it's perhaps a good
idea to support despite these problems. Perhaps the best idea is
to take the gamma.c code from this library and add it to the
Octave distribution? I don't know if the use of gamma instead
of lgamma causes problems with overflow, does it?

        Jonas Olsson
        [hidden email]
   
(As there are few people on the list and the source isn't that long
I'll include it here. The mathlib distribution comes with GPL
copyright)

#ifndef lint
static char sccsid[] = "@(#)gamma.c 1.6 (ucb.beef) 10/3/89";
#endif /* !defined(lint) */
/*
 *  This routine calculates the GAMMA function for a "double" argument X.
 *    Computation is based on an algorithm outlined in reference 1.
 *    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 reference 2.
 *    The accuracy achieved depends on the arithmetic system, the
 *    compiler, the intrinsic functions, and proper selection of the
 *    machine-dependent constants.
 *
 *
 *********************************************************************
 *********************************************************************
 *
 *  Explanation of machine-dependent constants
 *
 *  beta   - radix for the floating-point representation
 *  maxexp - the smallest positive power of beta that overflows
 *  XBIG   - the largest argument for which GAMMA(X) is representable
 *           in the machine, i.e., the solution to the equation
 *                   GAMMA(XBIG) = beta**maxexp
 *  XINF   - the largest machine representable floating-point number;
 *           approximately beta**maxexp
 *  EPS    - the smallest positive floating-point number such that
 *           1.0+EPS .GT. 1.0
 *  XMININ - the smallest positive floating-point number such that
 *           1/XMININ is machine representable
 *
 *      Approximate values for some important machines are:
 *
 *                             beta       maxexp        XBIG
 *
 *  CRAY-1         (S.P.)        2         8191        967.095
 *  Cyber 180/855
 *    under NOS    (S.P.)        2         1070        177.980
 *  IEEE (IBM/XT,
 *    SUN, etc.)   (S.P.)        2          128        35.299
 *  IEEE (IBM/XT,
 *    SUN, etc.)   (D.P.)        2         1024        171.624
 *  IBM 3033       (D.P.)       16           63        57.801
 *  VAX D-Format   (D.P.)        2          127        34.844
 *  VAX G-Format   (D.P.)        2         1023        171.668
 *
 *                             XINF         EPS        XMININ
 *
 *  CRAY-1         (S.P.)   5.45E+2465   7.11E-15    1.84E-2466
 *  Cyber 180/855
 *    under NOS    (S.P.)   1.26E+322    3.55E-15    3.14E-294
 *  IEEE (IBM/XT,
 *    SUN, etc.)   (S.P.)   3.40E+38     1.19E-7     1.18E-38
 *  IEEE (IBM/XT,
 *    SUN, etc.)   (D.P.)   1.79D+308    2.22D-16    2.23D-308
 *  IBM 3033       (D.P.)   7.23D+75     2.22D-16    1.39D-76
 *  VAX D-Format   (D.P.)   1.70D+38     1.39D-17    5.88D-39
 *  VAX G-Format   (D.P.)   8.98D+307    1.11D-16    1.12D-308
 *
 *********************************************************************
 *********************************************************************
 *
 *  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.
 *
 *
 *   Intrinsic functions required are:
 *
 *      exp, floor, log, sin
 *
 *
 *   References: "An Overview of Software Development for Special
 *               Functions", W. J. Cody, Lecture Notes in Mathematics,
 *               506, Numerical Analysis Dundee, 1975, G. A. Watson
 *               (ed.), Springer Verlag, Berlin, 1976.
 *
 *               Computer Approximations, Hart, Et. Al., Wiley and
 *               sons, New York, 1968.
 *
 *   Latest modification: May 30, 1989
 *
 *   Authors: W. J. Cody and L. Stoltz
 *            Applied Mathematics Division
 *            Argonne National Laboratory
 *            Argonne, IL 60439
 */
#if defined(__STDC__) || defined(__GNUC__)
extern double floor(double),exp(double),log(double),sin(double);
#else
extern double floor(),exp(),log(),sin();
#endif
                                        /* Machine dependent parameters */
#if defined(vax) || defined(tahoe)
#define XBIG (double)34.844e0
#define XMININ (double)5.88e-39
#define EPS (double)1.39e-17
#define XINF (double)1.70e38
#else /* defined(vax) || defined(tahoe) */
#define XBIG (double)171.624e0
#define XMININ (double)2.23e-308
#define EPS (double)2.22e-16
#define XINF (double)1.79e308
#endif /* defined(vax) || defined(tahoe) */
                                        /* Mathematical constants */
#define ONE (double)1
#define HALF (double)0.5
#define TWELVE (double)12
#define ZERO (double)0
#define SQRTPI (double)0.9189385332046727417803297e0
#define PI (double)3.1415926535897932384626434e0

typedef int logical;
#define FALSE (logical)0
#define TRUE (logical)1

/*
 *   Numerator and denominator coefficients for rational minimax
 *      approximation over (1,2).
 */
static double P[] = {
        -1.71618513886549492533811e0,
         2.47656508055759199108314e1,
        -3.79804256470945635097577e2,
         6.29331155312818442661052e2,
         8.66966202790413211295064e2,
        -3.14512729688483675254357e4,
        -3.61444134186911729807069e4,
         6.64561438202405440627855e4,
};
static double Q[] = {
        -3.08402300119738975254353e1,
         3.15350626979604161529144e2,
        -1.01515636749021914166146e3,
        -3.10777167157231109440444e3,
         2.25381184209801510330112e4,
         4.75584627752788110767815e3,
        -1.34659959864969306392456e5,
        -1.15132259675553483497211e5,
};

/*
 *   Coefficients for minimax approximation over (12, INF).
 */
static double C[] = {
        -1.910444077728e-03,
         8.4171387781295e-04,
        -5.952379913043012e-04,
         7.93650793500350248e-04,
        -2.777777777777681622553e-03,
         8.333333333333333331554247e-02,
         5.7083835261e-03,
};

double
#if defined(__STDC__) || defined(__GNUC__)
gamma(double x)
#else
gamma(x)
double x;
#endif
{
        register i,n;
        logical parity;
        double fact,res,xden,xnum,dtmp1,dtmp2;

        parity = FALSE;
        fact = ONE;
        n = 0;
#define y x
        if (y <= ZERO) { /* Arg. is negative */
                y = -x;
#define y1 dtmp1
                y1 = floor(y);
                res = y-y1;
                if (res != ZERO) {
                        parity = floor(y1/2)*2 != y1;
                        fact = -PI/sin(PI*res);
                        y += ONE;
                }
                else
                        return XINF;
        }
                                                /* Arg. is positive */
        if (y < EPS) { /* Arg. < EPS */
                if (y >= XMININ)
                        res = ONE/y;
                else
                  return XINF;
        }
        else if (y < TWELVE) {
#define y1 dtmp1
#define z dtmp2
                y1 = y;
                if (y < ONE) { /* 0.0 < arg. < 1.0 */
                        z = y;
                        y += ONE;
                }
                else { /* 1.0 < arg. < 12.0 */

                        n = (int)y; n--;
                        y -= (double)n; /* reduce arg. if needed */
                        z = y-ONE;
                }
                        /* Evaluate approximation for 1.0 < arg. < 2.0 */
                xnum = ZERO;
                xden = ONE;
                for (i = 0; i <= 7; i++) {
                        xnum = (xnum+P[i])*z;
                        xden = xden*z+Q[i];
                }
                res = xnum/xden+ONE;
                if (y1 < y) /* Adjust res for 0.0 < arg. < 1.0 */
                        res /= y1;
                else if (y1 > y) { /* Adjust res for 2.0 < arg. < 12.0 */
                        for (i = 0; i <= n-1; i++) {
                                res *= y;
                                y += ONE;
                        }
                }
#undef z
#undef y1
        }
        else { /* Evaluate for arg. >= 12.0 */
                if (y <= XBIG) {
#define ysq dtmp1
#define sum dtmp2
                        ysq = y*y;
                        sum = C[6];
                        for (i = 0; i <= 5; i++)
                                sum = sum/ysq+C[i];
                        sum = sum/y-y; sum += SQRTPI;
                        sum += (y-HALF)*log(y);
                        res = exp(sum);
#undef sum
#undef ysq
                }
                else
                        return XINF;
        }
                                            /* Final adjustments and return */
        if (parity == TRUE)
                res = -res;
        if (fact != ONE)
                res = fact/res;
        return res;
#undef y
}

Loading...