Blame view

GAMMA.cpp 1.9 KB
da3d4e0e   dmayerich   Initial commit.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
  //  gamma.cpp -- computation of gamma function.
  //      Algorithms and coefficient values from "Computation of Special
  //      Functions", Zhang and Jin, John Wiley and Sons, 1996.
  //
  //  (C) 2003, C. Bond. All rights reserved.
  //
  // Returns gamma function of argument 'x'.
  //
  // NOTE: Returns 1e308 if argument is a negative integer or 0,
  //      or if argument exceeds 171.
  //
  #define _USE_MATH_DEFINES
  #include <math.h>
  double gamma(double x)
  {
      int i,k,m;
      double ga,gr,r,z;
  
      static double g[] = {
          1.0,
          0.5772156649015329,
         -0.6558780715202538,
         -0.420026350340952e-1,
          0.1665386113822915,
         -0.421977345555443e-1,
         -0.9621971527877e-2,
          0.7218943246663e-2,
         -0.11651675918591e-2,
         -0.2152416741149e-3,
          0.1280502823882e-3,
         -0.201348547807e-4,
         -0.12504934821e-5,
          0.1133027232e-5,
         -0.2056338417e-6,
          0.6116095e-8,
          0.50020075e-8,
         -0.11812746e-8,
          0.1043427e-9,
          0.77823e-11,
         -0.36968e-11,
          0.51e-12,
         -0.206e-13,
         -0.54e-14,
          0.14e-14};
  
      if (x > 171.0) return 1e308;    // This value is an overflow flag.
      if (x == (int)x) {
          if (x > 0.0) {
              ga = 1.0;               // use factorial
              for (i=2;i<x;i++) {
                 ga *= i;
              }
           }
           else
              ga = 1e308;
       }
       else {
          if (fabs(x) > 1.0) {
              z = fabs(x);
              m = (int)z;
              r = 1.0;
              for (k=1;k<=m;k++) {
                  r *= (z-k);
              }
              z -= m;
          }
          else
              z = x;
          gr = g[24];
          for (k=23;k>=0;k--) {
              gr = gr*z+g[k];
          }
          ga = 1.0/(gr*z);
          if (fabs(x) > 1.0) {
              ga *= r;
              if (x < 0.0) {
                  ga = -M_PI/(x*ga*sin(M_PI*x));
              }
          }
      }
      return ga;
  }