|
0
|
1 #include "gamma.hpp"
|
|
|
2
|
|
|
3 /** The digamma function in long double precision.
|
|
|
4 * @param x the real value of the argument
|
|
|
5 * @return the value of the digamma (psi) function at that point
|
|
|
6 * @author Richard J. Mathar
|
|
|
7 * @since 2005-11-24
|
|
|
8 */
|
|
|
9 long double digammal(long double x)
|
|
|
10 {
|
|
|
11 /* force into the interval 1..3 */
|
|
|
12 if( x < 0.0L )
|
|
|
13 return digammal(1.0L-x)+M_PIl/tanl(M_PIl*(1.0L-x)) ; /* reflection formula */
|
|
|
14 else if( x < 1.0L )
|
|
|
15 return digammal(1.0L+x)-1.0L/x ;
|
|
|
16 else if ( x == 1.0L)
|
|
|
17 return -M_GAMMAl ;
|
|
|
18 else if ( x == 2.0L)
|
|
|
19 return 1.0L-M_GAMMAl ;
|
|
|
20 else if ( x == 3.0L)
|
|
|
21 return 1.5L-M_GAMMAl ;
|
|
|
22 else if ( x > 3.0L)
|
|
|
23 /* duplication formula */
|
|
|
24 return 0.5L*(digammal(x/2.0L)+digammal((x+1.0L)/2.0L))+M_LN2l ;
|
|
|
25 else
|
|
|
26 {
|
|
|
27 /* Just for your information, the following lines contain
|
|
|
28 * the Maple source code to re-generate the table that is
|
|
|
29 * eventually becoming the Kncoe[] array below
|
|
|
30 * interface(prettyprint=0) :
|
|
|
31 * Digits := 63 :
|
|
|
32 * r := 0 :
|
|
|
33 *
|
|
|
34 * for l from 1 to 60 do
|
|
|
35 * d := binomial(-1/2,l) :
|
|
|
36 * r := r+d*(-1)^l*(Zeta(2*l+1) -1) ;
|
|
|
37 * evalf(r) ;
|
|
|
38 * print(%,evalf(1+Psi(1)-r)) ;
|
|
|
39 *o d :
|
|
|
40 *
|
|
|
41 * for N from 1 to 28 do
|
|
|
42 * r := 0 :
|
|
|
43 * n := N-1 :
|
|
|
44 *
|
|
|
45 * for l from iquo(n+3,2) to 70 do
|
|
|
46 * d := 0 :
|
|
|
47 * for s from 0 to n+1 do
|
|
|
48 * d := d+(-1)^s*binomial(n+1,s)*binomial((s-1)/2,l) :
|
|
|
49 * od :
|
|
|
50 * if 2*l-n > 1 then
|
|
|
51 * r := r+d*(-1)^l*(Zeta(2*l-n) -1) :
|
|
|
52 * fi :
|
|
|
53 * od :
|
|
|
54 * print(evalf((-1)^n*2*r)) ;
|
|
|
55 *od :
|
|
|
56 *quit :
|
|
|
57 */
|
|
|
58 static long double Kncoe[] = { .30459198558715155634315638246624251L,
|
|
|
59 .72037977439182833573548891941219706L, -.12454959243861367729528855995001087L,
|
|
|
60 .27769457331927827002810119567456810e-1L, -.67762371439822456447373550186163070e-2L,
|
|
|
61 .17238755142247705209823876688592170e-2L, -.44817699064252933515310345718960928e-3L,
|
|
|
62 .11793660000155572716272710617753373e-3L, -.31253894280980134452125172274246963e-4L,
|
|
|
63 .83173997012173283398932708991137488e-5L, -.22191427643780045431149221890172210e-5L,
|
|
|
64 .59302266729329346291029599913617915e-6L, -.15863051191470655433559920279603632e-6L,
|
|
|
65 .42459203983193603241777510648681429e-7L, -.11369129616951114238848106591780146e-7L,
|
|
|
66 .304502217295931698401459168423403510e-8L, -.81568455080753152802915013641723686e-9L,
|
|
|
67 .21852324749975455125936715817306383e-9L, -.58546491441689515680751900276454407e-10L,
|
|
|
68 .15686348450871204869813586459513648e-10L, -.42029496273143231373796179302482033e-11L,
|
|
|
69 .11261435719264907097227520956710754e-11L, -.30174353636860279765375177200637590e-12L,
|
|
|
70 .80850955256389526647406571868193768e-13L, -.21663779809421233144009565199997351e-13L,
|
|
|
71 .58047634271339391495076374966835526e-14L, -.15553767189204733561108869588173845e-14L,
|
|
|
72 .41676108598040807753707828039353330e-15L, -.11167065064221317094734023242188463e-15L } ;
|
|
|
73
|
|
|
74 register long double Tn_1 = 1.0L ; /* T_{n-1}(x), started at n=1 */
|
|
|
75 register long double Tn = x-2.0L ; /* T_{n}(x) , started at n=1 */
|
|
|
76 register long double resul = Kncoe[0] + Kncoe[1]*Tn ;
|
|
|
77
|
|
|
78 x -= 2.0L ;
|
|
|
79 int n ;
|
|
|
80
|
|
|
81 for( n = 2 ; n < sizeof(Kncoe)/sizeof(long double) ;n++)
|
|
|
82 {
|
|
|
83 const long double Tn1 = 2.0L * x * Tn - Tn_1 ; /* Chebyshev recursion, Eq. 22.7.4 Abramowitz-Stegun */
|
|
|
84 resul += Kncoe[n]*Tn1 ;
|
|
|
85 Tn_1 = Tn ;
|
|
|
86 Tn = Tn1 ;
|
|
|
87 }
|
|
|
88 return resul ;
|
|
|
89 }
|
|
|
90 }
|
|
|
91
|
|
|
92
|
|
|
93
|
|
|
94 double trigamma ( double x, int *ifault )
|
|
|
95
|
|
|
96 //****************************************************************************
|
|
|
97 // purpose:
|
|
|
98 //
|
|
|
99 // trigamma calculates trigamma(x) = d**2 log(gamma(x)) / dx**2
|
|
|
100 //
|
|
|
101 // licensing:
|
|
|
102 //
|
|
|
103 // this code is distributed under the gnu lgpl license.
|
|
|
104 //
|
|
|
105 // modified:
|
|
|
106 //
|
|
|
107 // 19 january 2008
|
|
|
108 //
|
|
|
109 // author:
|
|
|
110 //
|
|
|
111 // original fortran77 version by be schneider.
|
|
|
112 // c++ version by john burkardt.
|
|
|
113 //
|
|
|
114 // reference:
|
|
|
115 //
|
|
|
116 // be schneider,
|
|
|
117 // algorithm as 121:
|
|
|
118 // trigamma function,
|
|
|
119 // applied statistics,
|
|
|
120 // volume 27, number 1, pages 97-99, 1978.
|
|
|
121 //
|
|
|
122 // parameters:
|
|
|
123 //
|
|
|
124 // input, double x, the argument of the trigamma function.
|
|
|
125 // 0 < x.
|
|
|
126 //
|
|
|
127 // output, int *ifault, error flag.
|
|
|
128 // 0, no error.
|
|
|
129 // 1, x <= 0.
|
|
|
130 //
|
|
|
131 // output, double trigamma, the value of the trigamma function at x.
|
|
|
132 //
|
|
|
133 {
|
|
|
134 double a = 0.0001;
|
|
|
135 double b = 5.0;
|
|
|
136 double b2 = 0.1666666667;
|
|
|
137 double b4 = -0.03333333333;
|
|
|
138 double b6 = 0.02380952381;
|
|
|
139 double b8 = -0.03333333333;
|
|
|
140 double value;
|
|
|
141 double y;
|
|
|
142 double z;
|
|
|
143 //
|
|
|
144 // check the input.
|
|
|
145 //
|
|
|
146 if ( x <= 0.0 )
|
|
|
147 {
|
|
|
148 *ifault = 1;
|
|
|
149 value = 0.0;
|
|
|
150 return value;
|
|
|
151 }
|
|
|
152
|
|
|
153 *ifault = 0;
|
|
|
154 z = x;
|
|
|
155 //
|
|
|
156 // use small value approximation if x <= a.
|
|
|
157 //
|
|
|
158 if ( x <= a )
|
|
|
159 {
|
|
|
160 value = 1.0 / x / x;
|
|
|
161 return value;
|
|
|
162 }
|
|
|
163 //
|
|
|
164 // increase argument to ( x + i ) >= b.
|
|
|
165 //
|
|
|
166 value = 0.0;
|
|
|
167
|
|
|
168 while ( z < b )
|
|
|
169 {
|
|
|
170 value = value + 1.0 / z / z;
|
|
|
171 z = z + 1.0;
|
|
|
172 }
|
|
|
173 //
|
|
|
174 // apply asymptotic formula if argument is b or greater.
|
|
|
175 //
|
|
|
176 y = 1.0 / z / z;
|
|
|
177
|
|
|
178 value = value + 0.5 *
|
|
|
179 y + ( 1.0 + y * ( b2+ y * ( b4 + y * ( b6+ y * b8 )))) / z;
|
|
|
180
|
|
|
181 return value;
|
|
|
182 }
|
|
|
183
|
|
|
184
|
|
|
185 double LogGammaDensity( double x, double k, double theta )
|
|
|
186 {
|
|
|
187 return -k * log( theta ) + ( k - 1 ) * log( x ) - x / theta - lgamma( k ) ;
|
|
|
188 }
|
|
|
189
|
|
|
190 double MixtureGammaAssignment( double x, double pi, double* k, double *theta )
|
|
|
191 {
|
|
|
192 if ( pi == 1 )
|
|
|
193 return 0 ;
|
|
|
194 else if ( pi == 0 )
|
|
|
195 return 1 ;
|
|
|
196
|
|
|
197 double lf0 = LogGammaDensity( x, k[0], theta[0] ) ;
|
|
|
198 double lf1 = LogGammaDensity( x, k[1], theta[1] ) ;
|
|
|
199
|
|
|
200 return (double)1.0 / ( 1.0 + exp( lf1 + log( 1 - pi ) - lf0 - log( pi ) ) ) ;
|
|
|
201 }
|
