123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139 |
- /*
- * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
- * See the copyright notice in the ACK home directory, in the file "Copyright".
- *
- * Author: Ceriel J.H. Jacobs
- */
- /* $Id$ */
- #include <math.h>
- #include <errno.h>
- extern int errno;
- static double
- smallpos_gamma(x)
- double x;
- {
- /* Approximation of gamma function using
- gamma(x) = P(x-2) / Q(x-2) for x in [2,3]
- */
- /* Hart & Cheney # 5251 */
- static double p[11] = {
- -0.2983543278574342138830437659e+06,
- -0.2384953970018198872468734423e+06,
- -0.1170494760121780688403854445e+06,
- -0.3949445048301571936421824091e+05,
- -0.1046699423827521405330650531e+05,
- -0.2188218110071816359394795998e+04,
- -0.3805112208641734657584922631e+03,
- -0.5283123755635845383718978382e+02,
- -0.6128571763704498306889428212e+01,
- -0.5028018054416812467364198750e+00,
- -0.3343060322330595274515660112e-01
- };
- static double q[9] = {
- -0.2983543278574342138830438524e+06,
- -0.1123558608748644911342306408e+06,
- 0.5332716689118142157485686311e+05,
- 0.8571160498907043851961147763e+04,
- -0.4734865977028211706556819770e+04,
- 0.1960497612885585838997039621e+03,
- 0.1257733367869888645966647426e+03,
- -0.2053126153100672764513929067e+02,
- 0.1000000000000000000000000000e+01
- };
- double result = 1.0;
- while (x > 3) {
- x -= 1.0;
- result *= x;
- }
- while (x < 2) {
- result /= x;
- x += 1.0;
- }
- x -= 2.0;
- return result * POLYNOM10(x, p) / POLYNOM8(x, q);
- }
- #define log_sqrt_2pi 0.91893853320467274178032973640561763
- int signgam;
- static double
- bigpos_loggamma(x)
- double x;
- {
- /* computes the log(gamma(x)) function for big arguments
- using the Stirling form
- log(gamma(x)) = (x - 0.5)log(x) - x + log(sqrt(2*pi)) + fi(x)
- where fi(x) = (1/x)*P(1/(x*x))/Q(1/(x*x)) for x in [12,1000]
- */
- /* Hart & Cheney # 5468 */
- static double p[4] = {
- 0.12398282342474941538685913e+00,
- 0.67082783834332134961461700e+00,
- 0.64507302912892202513890000e+00,
- 0.66662907040200752600000000e-01
- };
- static double q[4] = {
- 0.14877938810969929846815600e+01,
- 0.80995271894897557472821400e+01,
- 0.79966911236636441947720000e+01,
- 0.10000000000000000000000000e+01
- };
- double rsq = 1.0/(x*x);
- extern double log();
- return (x-0.5)*log(x)-x+log_sqrt_2pi+POLYNOM3(rsq, p)/(x*POLYNOM3(rsq, q));
- }
- static double
- neg_loggamma(x)
- double x;
- {
- /* compute the log(gamma(x)) function for negative values of x,
- using the rule:
- -x*gamma(x)*gamma(-x) = pi/sin(z*pi)
- */
- extern double sin(), log();
- double sinpix;
- x = -x;
- sinpix = sin(M_PI * x);
- if (sinpix == 0.0) {
- errno = EDOM;
- return HUGE;
- }
- if (sinpix < 0) sinpix = -sinpix;
- else signgam = -1;
- return log(M_PI/(x * smallpos_gamma(x) * sinpix));
- }
- double
- gamma(x)
- double x;
- {
- /* Wrong name; Actually computes log(gamma(x))
- */
- extern double log();
- signgam = 1;
- if (x <= 0) {
- return neg_loggamma(x);
- }
- if (x > 12.0) {
- return bigpos_loggamma(x);
- }
- return log(smallpos_gamma(x));
- }
|