113618394Sjtc /* @(#)w_gamma.c 5.1 93/09/24 */
213618394Sjtc /*
313618394Sjtc * ====================================================
413618394Sjtc * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
513618394Sjtc *
613618394Sjtc * Developed at SunPro, a Sun Microsystems, Inc. business.
713618394Sjtc * Permission to use, copy, modify, and distribute this
813618394Sjtc * software is freely granted, provided that this notice
913618394Sjtc * is preserved.
1013618394Sjtc * ====================================================
1113618394Sjtc */
1213618394Sjtc
13118c71a6Slukem #include <sys/cdefs.h>
14d1f06e0bSjtc #if defined(LIBM_SCCS) && !defined(lint)
15*aa30599eSwiz __RCSID("$NetBSD: w_gamma.c,v 1.11 2002/05/26 22:02:00 wiz Exp $");
16bc3f7bf6Sjtc #endif
17bc3f7bf6Sjtc
1813618394Sjtc /* double gamma(double x)
1913618394Sjtc * Return the logarithm of the Gamma function of x.
2013618394Sjtc *
2113618394Sjtc * Method: call gamma_r
2213618394Sjtc */
2313618394Sjtc
248346e333Sjtc #include "math.h"
258346e333Sjtc #include "math_private.h"
2613618394Sjtc
27*aa30599eSwiz double
gamma(double x)28*aa30599eSwiz gamma(double x)
2913618394Sjtc {
3013618394Sjtc #ifdef _IEEE_LIBM
3197740060Sjtc return __ieee754_lgamma_r(x,&signgam);
3213618394Sjtc #else
3313618394Sjtc double y;
3497740060Sjtc y = __ieee754_lgamma_r(x,&signgam);
3513618394Sjtc if(_LIB_VERSION == _IEEE_) return y;
3613618394Sjtc if(!finite(y)&&finite(x)) {
3713618394Sjtc if(floor(x)==x&&x<=0.0)
3813618394Sjtc return __kernel_standard(x,x,41); /* gamma pole */
3913618394Sjtc else
4013618394Sjtc return __kernel_standard(x,x,40); /* gamma overflow */
4113618394Sjtc } else
4213618394Sjtc return y;
4313618394Sjtc #endif
4413618394Sjtc }
45