gamma.c 2.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139
  1. /*
  2. * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
  3. * See the copyright notice in the ACK home directory, in the file "Copyright".
  4. *
  5. * Author: Ceriel J.H. Jacobs
  6. */
  7. /* $Id$ */
  8. #include <math.h>
  9. #include <errno.h>
  10. extern int errno;
  11. static double
  12. smallpos_gamma(x)
  13. double x;
  14. {
  15. /* Approximation of gamma function using
  16. gamma(x) = P(x-2) / Q(x-2) for x in [2,3]
  17. */
  18. /* Hart & Cheney # 5251 */
  19. static double p[11] = {
  20. -0.2983543278574342138830437659e+06,
  21. -0.2384953970018198872468734423e+06,
  22. -0.1170494760121780688403854445e+06,
  23. -0.3949445048301571936421824091e+05,
  24. -0.1046699423827521405330650531e+05,
  25. -0.2188218110071816359394795998e+04,
  26. -0.3805112208641734657584922631e+03,
  27. -0.5283123755635845383718978382e+02,
  28. -0.6128571763704498306889428212e+01,
  29. -0.5028018054416812467364198750e+00,
  30. -0.3343060322330595274515660112e-01
  31. };
  32. static double q[9] = {
  33. -0.2983543278574342138830438524e+06,
  34. -0.1123558608748644911342306408e+06,
  35. 0.5332716689118142157485686311e+05,
  36. 0.8571160498907043851961147763e+04,
  37. -0.4734865977028211706556819770e+04,
  38. 0.1960497612885585838997039621e+03,
  39. 0.1257733367869888645966647426e+03,
  40. -0.2053126153100672764513929067e+02,
  41. 0.1000000000000000000000000000e+01
  42. };
  43. double result = 1.0;
  44. while (x > 3) {
  45. x -= 1.0;
  46. result *= x;
  47. }
  48. while (x < 2) {
  49. result /= x;
  50. x += 1.0;
  51. }
  52. x -= 2.0;
  53. return result * POLYNOM10(x, p) / POLYNOM8(x, q);
  54. }
  55. #define log_sqrt_2pi 0.91893853320467274178032973640561763
  56. int signgam;
  57. static double
  58. bigpos_loggamma(x)
  59. double x;
  60. {
  61. /* computes the log(gamma(x)) function for big arguments
  62. using the Stirling form
  63. log(gamma(x)) = (x - 0.5)log(x) - x + log(sqrt(2*pi)) + fi(x)
  64. where fi(x) = (1/x)*P(1/(x*x))/Q(1/(x*x)) for x in [12,1000]
  65. */
  66. /* Hart & Cheney # 5468 */
  67. static double p[4] = {
  68. 0.12398282342474941538685913e+00,
  69. 0.67082783834332134961461700e+00,
  70. 0.64507302912892202513890000e+00,
  71. 0.66662907040200752600000000e-01
  72. };
  73. static double q[4] = {
  74. 0.14877938810969929846815600e+01,
  75. 0.80995271894897557472821400e+01,
  76. 0.79966911236636441947720000e+01,
  77. 0.10000000000000000000000000e+01
  78. };
  79. double rsq = 1.0/(x*x);
  80. extern double log();
  81. return (x-0.5)*log(x)-x+log_sqrt_2pi+POLYNOM3(rsq, p)/(x*POLYNOM3(rsq, q));
  82. }
  83. static double
  84. neg_loggamma(x)
  85. double x;
  86. {
  87. /* compute the log(gamma(x)) function for negative values of x,
  88. using the rule:
  89. -x*gamma(x)*gamma(-x) = pi/sin(z*pi)
  90. */
  91. extern double sin(), log();
  92. double sinpix;
  93. x = -x;
  94. sinpix = sin(M_PI * x);
  95. if (sinpix == 0.0) {
  96. errno = EDOM;
  97. return HUGE;
  98. }
  99. if (sinpix < 0) sinpix = -sinpix;
  100. else signgam = -1;
  101. return log(M_PI/(x * smallpos_gamma(x) * sinpix));
  102. }
  103. double
  104. gamma(x)
  105. double x;
  106. {
  107. /* Wrong name; Actually computes log(gamma(x))
  108. */
  109. extern double log();
  110. signgam = 1;
  111. if (x <= 0) {
  112. return neg_loggamma(x);
  113. }
  114. if (x > 12.0) {
  115. return bigpos_loggamma(x);
  116. }
  117. return log(smallpos_gamma(x));
  118. }