1 module rmathd.gamma;
2 
3 public import rmathd.common;
4 public import rmathd.normal;
5 //import rmathd.uniform;
6 
7 /*
8 ** dmd gamma.d uniform.d normal.d common.d && ./gamma
9 */
10 
11 T dgamma(T: double)(T x, T shape, T scale, int give_log)
12 {
13     T pr;
14     if (isNaN(x) || isNaN(shape) || isNaN(scale))
15         return x + shape + scale;
16 
17     mixin R_D__0!give_log;
18 
19     if (shape < 0 || scale <= 0)
20     	return T.nan;
21     if (x < 0)
22 	    return R_D__0;
23     if (shape == 0) /* point mass at 0 */
24 	    return (x == 0)? T.infinity : R_D__0;
25     if (x == 0) {
26 	    if (shape < 1)
27 	    	return T.infinity;
28 	    if (shape > 1)
29 	    	return R_D__0;
30 	    return give_log ? -log(scale) : 1 / scale;
31     }
32 
33     if (shape < 1) {
34 	    pr = dpois_raw!T(shape, x/scale, give_log);
35 	    return give_log ?  pr + log(shape/x) : pr*shape/x;
36     }
37     /* else  shape >= 1 */
38     pr = dpois_raw!T(shape - 1, x/scale, give_log);
39     return give_log ? pr - log(scale) : pr/scale;
40 }
41 
42 /* Continued fraction for calculation of
43  *    1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ... = sum_{k=0}^Inf x^k/(i+k*d)
44  *
45  * auxilary in log1pmx() and lgamma1p()
46  */
47 static T logcf(T)(T x, T i, T d, T eps /* ~ relative tolerance */)
48 {
49     T c1 = 2 * d;
50     T c2 = i + d;
51     T c4 = c2 + d;
52     T a1 = c2;
53     T b1 = i * (c2 - i * x);
54     T b2 = d * d * x;
55     T a2 = c4 * c2 - b2;
56 
57     //#if 0
58     //    assert(i > 0);
59     //    assert(d >= 0);
60     //#endif
61 
62     b2 = c4 * b1 - i * b2;
63 
64     while (fabs(a2 * b1 - a1 * b2) > fabs(eps * b1 * b2)) {
65 	T c3 = c2*c2*x;
66 	c2 += d;
67 	c4 += d;
68 	a1 = c4 * a2 - c3 * a1;
69 	b1 = c4 * b2 - c3 * b1;
70 
71 	c3 = c1 * c1 * x;
72 	c1 += d;
73 	c4 += d;
74 	a2 = c4 * a1 - c3 * a2;
75 	b2 = c4 * b1 - c3 * b2;
76 
77 	if (fabs (b2) > scalefactor) {
78 	    a1 /= scalefactor;
79 	    b1 /= scalefactor;
80 	    a2 /= scalefactor;
81 	    b2 /= scalefactor;
82 	} else if (fabs (b2) < 1 / scalefactor) {
83 	    a1 *= scalefactor;
84 	    b1 *= scalefactor;
85 	    a2 *= scalefactor;
86 	    b2 *= scalefactor;
87 	}
88     }
89 
90     return a2 / b2;
91 }
92 
93 /* Accurate calculation of log(1+x)-x, particularly for small x.  */
94 T log1pmx(T)(T x)
95 {
96     static const(T) minLog1Value = -0.79149064;
97 
98     if (x > 1 || x < minLog1Value)
99 	    return log1p(x) - x;
100     else { /* -.791 <=  x <= 1  -- expand in  [x/(2+x)]^2 =: y :
101 	    * log(1+x) - x =  x/(2+x) * [ 2 * y * S(y) - x],  with
102 	    * ---------------------------------------------
103 	    * S(y) = 1/3 + y/5 + y^2/7 + ... = \sum_{k=0}^\infty  y^k / (2k + 3)
104 	   */
105 	    T r = x / (2 + x), y = r * r;
106 	    if (fabs(x) < 1e-2) {
107 	        static const(T) two = 2;
108 	        return r * ((((two / 9 * y + two / 7) * y + two / 5) * y +
109 	    		    two / 3) * y - x);
110 	    } else {
111 	        static const(T) tol_logcf = 1e-14;
112 	        return r * (2 * y * logcf!T(y, 3, 2, tol_logcf) - x);
113 	    }
114     }
115 }
116 
117 /* Compute  log(gamma(a+1))  accurately also for small a (0 < a < 0.5). */
118 T lgamma1p(T)(T a)
119 {
120     const(T) eulers_const = 0.5772156649015328606065120900824024;
121 
122     /* coeffs[i] holds (zeta(i+2)-1)/(i+2) , i = 0:(N-1), N = 40 : */
123     const(int) N = 40;
124     static const(T)[40] coeffs = [
125 	0.3224670334241132182362075833230126e-0,/* = (zeta(2)-1)/2 */
126 	0.6735230105319809513324605383715000e-1,/* = (zeta(3)-1)/3 */
127 	0.2058080842778454787900092413529198e-1,
128 	0.7385551028673985266273097291406834e-2,
129 	0.2890510330741523285752988298486755e-2,
130 	0.1192753911703260977113935692828109e-2,
131 	0.5096695247430424223356548135815582e-3,
132 	0.2231547584535793797614188036013401e-3,
133 	0.9945751278180853371459589003190170e-4,
134 	0.4492623673813314170020750240635786e-4,
135 	0.2050721277567069155316650397830591e-4,
136 	0.9439488275268395903987425104415055e-5,
137 	0.4374866789907487804181793223952411e-5,
138 	0.2039215753801366236781900709670839e-5,
139 	0.9551412130407419832857179772951265e-6,
140 	0.4492469198764566043294290331193655e-6,
141 	0.2120718480555466586923135901077628e-6,
142 	0.1004322482396809960872083050053344e-6,
143 	0.4769810169363980565760193417246730e-7,
144 	0.2271109460894316491031998116062124e-7,
145 	0.1083865921489695409107491757968159e-7,
146 	0.5183475041970046655121248647057669e-8,
147 	0.2483674543802478317185008663991718e-8,
148 	0.1192140140586091207442548202774640e-8,
149 	0.5731367241678862013330194857961011e-9,
150 	0.2759522885124233145178149692816341e-9,
151 	0.1330476437424448948149715720858008e-9,
152 	0.6422964563838100022082448087644648e-10,
153 	0.3104424774732227276239215783404066e-10,
154 	0.1502138408075414217093301048780668e-10,
155 	0.7275974480239079662504549924814047e-11,
156 	0.3527742476575915083615072228655483e-11,
157 	0.1711991790559617908601084114443031e-11,
158 	0.8315385841420284819798357793954418e-12,
159 	0.4042200525289440065536008957032895e-12,
160 	0.1966475631096616490411045679010286e-12,
161 	0.9573630387838555763782200936508615e-13,
162 	0.4664076026428374224576492565974577e-13,
163 	0.2273736960065972320633279596737272e-13,
164 	0.1109139947083452201658320007192334e-13/* = (zeta(40+1)-1)/(40+1) */
165     ];
166 
167     const(T) c = 0.2273736845824652515226821577978691e-12;/* zeta(N+2)-1 */
168     const(T) tol_logcf = 1e-14;
169     T lgam;
170     int i;
171 
172     if (fabs (a) >= 0.5)
173 	    return lgammafn!T(a + 1);
174 
175     /* Abramowitz & Stegun 6.1.33 : for |x| < 2,
176      * <==> log(gamma(1+x)) = -(log(1+x) - x) - gamma*x + x^2 * \sum_{n=0}^\infty c_n (-x)^n
177      * where c_n := (Zeta(n+2) - 1)/(n+2)  = coeffs[n]
178      *
179      * Here, another convergence acceleration trick is used to compute
180      * lgam(x) :=  sum_{n=0..Inf} c_n (-x)^n
181      */
182     lgam = c * logcf!T(-a / 2, N + 2, 1, tol_logcf);
183     for (i = N - 1; i >= 0; i--)
184 	    lgam = coeffs[i] - a * lgam;
185 
186     return (a * lgam - eulers_const) * a - log1pmx!T(a);
187 } /* lgamma1p */
188 
189 
190 /*
191  * Compute the log of a sum from logs of terms, i.e.,
192  *
193  *     log (exp (logx) + exp (logy))
194  *
195  * without causing overflows and without throwing away large handfuls
196  * of accuracy.
197  */
198 T logspace_add(T)(T logx, T logy)
199 {
200     return fmax2!T(logx, logy) + log1p (exp (-fabs (logx - logy)));
201 }
202 
203 /*
204  * Compute the log of a difference from logs of terms, i.e.,
205  *
206  *     log (exp (logx) - exp (logy))
207  *
208  * without causing overflows and without throwing away large handfuls
209  * of accuracy.
210  */
211 T logspace_sub(T)(T logx, T logy)
212 {
213     return logx + R_Log1_Exp!T(logy - logx);
214 }
215 
216 /*
217  * Compute the log of a sum from logs of terms, i.e.,
218  *
219  *     log (sum_i  exp (logx[i]) ) =
220  *     log (e^M * sum_i  e^(logx[i] - M) ) =
221  *     M + log( sum_i  e^(logx[i] - M)
222  *
223  * without causing overflows or throwing much accuracy.
224  */
225 T logspace_sum(T)(const(T)* logx, int n)
226 {
227     if(n == 0)
228         return -T.infinity; // = log( sum(<empty>) )
229     if(n == 1)
230     	return logx[0];
231     if(n == 2)
232     	return logspace_add(logx[0], logx[1]);
233     // else (n >= 3) :
234     int i;
235     // Mx := max_i log(x_i)
236     T Mx = logx[0];
237     for(i = 1; i < n; i++)
238     	if(Mx < logx[i]) Mx = logx[i];
239     real s = 0.;
240     for(i = 0; i < n; i++)
241     	s += exp(logx[i] - Mx);
242     return Mx + cast(T)log(s);
243 }
244 
245 /* dpois_wrap (x__1, lambda) := dpois(x__1 - 1, lambda);  where
246  * dpois(k, L) := exp(-L) L^k / gamma(k+1)  {the usual Poisson probabilities}
247  *
248  * and  dpois*(.., give_log = TRUE) :=  log( dpois*(..) )
249 */
250 static T dpois_wrap(T)(T x_plus_1, T lambda, int give_log)
251 {
252 	mixin R_D__0!(give_log);
253 	mixin M_cutoff!T;
254     if (!isFinite(lambda))
255 	    return R_D__0;
256     if (x_plus_1 > 1)
257 	    return dpois_raw!T(x_plus_1 - 1, lambda, give_log);
258     if (lambda > fabs(x_plus_1 - 1) * M_cutoff)
259 	    return R_D_exp!T(-lambda - lgammafn!T(x_plus_1), give_log);
260     else {
261 	T d = dpois_raw!T(x_plus_1, lambda, give_log);
262 	return give_log
263 	    ? d + log(x_plus_1 / lambda)
264 	    : d * (x_plus_1 / lambda);
265     }
266 }
267 
268 /*
269  * Abramowitz and Stegun 6.5.29 [right]
270  */
271 static T pgamma_smallx(T)(T x, T alph, int lower_tail, int log_p)
272 {
273     T sum = 0, c = alph, n = 0, term;
274 
275     /*
276      * Relative to 6.5.29 all terms have been multiplied by alph
277      * and the first, thus being 1, is omitted.
278      */
279 
280     do {
281 	    n++;
282 	    c *= -x / n;
283 	    term = c / (alph + n);
284 	    sum += term;
285     } while (fabs(term) > EPSILON!T * fabs(sum));
286 
287     if (lower_tail) {
288 	    T f1 = log_p ? log1p(sum) : 1 + sum;
289 	    T f2;
290 	    if (alph > 1) {
291 	        f2 = dpois_raw!T(alph, x, log_p);
292 	        f2 = log_p ? f2 + x : f2 * exp (x);
293 	    } else if (log_p)
294 	        f2 = alph * log (x) - lgamma1p!T(alph);
295 	    else
296 	        f2 = pow(x, alph) / exp (lgamma1p!T(alph));
297         
298 	    return log_p ? f1 + f2 : f1 * f2;
299     } else {
300 	    T lf2 = alph * log (x) - lgamma1p!T(alph);
301         
302 	    if (log_p)
303 	        return R_Log1_Exp!T(log1p(sum) + lf2);
304 	    else {
305 	        T f1m1 = sum;
306 	        T f2m1 = expm1 (lf2);
307 	        return -(f1m1 + f2m1 + f1m1 * f2m1);
308 	    }
309     }
310 } /* pgamma_smallx() */
311 
312 
313 static T pd_upper_series(T)(T x, T y, int log_p)
314 {
315     T term = x / y;
316     T sum = term;
317 
318     do {
319 	    y++;
320 	    term *= x / y;
321 	    sum += term;
322     } while (term > sum * EPSILON!T);
323 
324     /* sum =  \sum_{n=1}^ oo  x^n     / (y*(y+1)*...*(y+n-1))
325      *	   =  \sum_{n=0}^ oo  x^(n+1) / (y*(y+1)*...*(y+n))
326      *	   =  x/y * (1 + \sum_{n=1}^oo	x^n / ((y+1)*...*(y+n)))
327      *	   ~  x/y +  o(x/y)   {which happens when alph -> Inf}
328      */
329     return log_p ? log(sum) : sum;
330 }
331 
332 
333 immutable(string) NEEDED_SCALE(string ctrl)()
334 {
335 	immutable(string) output = ctrl ~ `(b2 > scalefactor) {
336 	                         a1 /= scalefactor;
337 	                         b1 /= scalefactor;
338 	                         a2 /= scalefactor;
339 	                         b2 /= scalefactor;
340 	                     }`;
341 	return output;
342 }
343 
344 
345 static T pd_lower_cf (T)(T y, T d)
346 {
347     T f = 0.0 /* -Wall */, of, f0;
348     T i, c2, c3, c4,  a1, b1,  a2, b2;
349 
350     enum max_it = 200000;
351 
352     if (y == 0)
353     	return 0;
354 
355     f0 = y/d;
356     /* Needed, e.g. for  pgamma(10^c(100,295), shape= 1.1, log=TRUE): */
357     if(fabs(y - 1) < fabs(d) * EPSILON!T) { /* includes y < d = Inf */
358 	    return (f0);
359     }
360 
361     if(f0 > 1.) f0 = 1.;
362     c2 = y;
363     c4 = d; /* original (y,d), *not* potentially scaled ones!*/
364 
365     a1 = 0; b1 = 1;
366     a2 = y; b2 = d;
367 
368     mixin (NEEDED_SCALE!("while")());
369 
370     i = 0; of = -1.; /* far away */
371     while (i < max_it) {
372 	    i++;	c2--;	c3 = i * c2;	c4 += 2;
373 	    /* c2 = y - i,  c3 = i(y - i),  c4 = d + 2i,  for i odd */
374 	    a1 = c4 * a2 + c3 * a1;
375 	    b1 = c4 * b2 + c3 * b1;
376         
377 	    i++;	c2--;	c3 = i * c2;	c4 += 2;
378 	    /* c2 = y - i,  c3 = i(y - i),  c4 = d + 2i,  for i even */
379 	    a2 = c4 * a1 + c3 * a2;
380 	    b2 = c4 * b1 + c3 * b2;
381         
382 	    mixin (NEEDED_SCALE!("if")());
383         
384 	    if (b2 != 0) {
385 	        f = a2 / b2;
386 	        /* convergence check: relative; "absolute" for very small f : */
387 	        if (fabs (f - of) <= EPSILON!T * fmax2(f0, fabs(f))) {
388 	    	    return f;
389 	        }
390 	        of = f;
391 	    }
392     }
393 
394     return f;/* should not happen ... */
395 }
396 
397 
398 static T pd_lower_series(T)(T lambda, T y)
399 {
400     T term = 1, sum = 0;
401 
402     while (y >= 1 && term > sum * EPSILON!T) {
403 	    term *= y / lambda;
404 	    sum += term;
405 	    y--;
406     }
407     /* sum =  \sum_{n=0}^ oo  y*(y-1)*...*(y - n) / lambda^(n+1)
408      *	   =  y/lambda * (1 + \sum_{n=1}^Inf  (y-1)*...*(y-n) / lambda^n)
409      *	   ~  y/lambda + o(y/lambda)
410      */
411 
412     if (y != floor(y)) {
413 	    /*
414 	     * The series does not converge as the terms start getting
415 	     * bigger (besides flipping sign) for y < -lambda.
416 	     */
417 	    T f;
418 	    /* FIXME: in quite few cases, adding  term*f  has no effect (f too small)
419 	     *	  and is unnecessary e.g. for pgamma(4e12, 121.1) */
420 	    f = pd_lower_cf!T(y, lambda + 1 - y);
421 	    sum += term * f;
422     }
423 
424     return sum;
425 } /* pd_lower_series() */
426 
427 
428 /*
429  * Compute the following ratio with higher accuracy that would be had
430  * from doing it directly.
431  *
432  *		 dnorm (x, 0, 1, FALSE)
433  *	   ----------------------------------
434  *	   pnorm (x, 0, 1, lower_tail, FALSE)
435  *
436  * Abramowitz & Stegun 26.2.12
437  */
438 static T dpnorm(T)(T x, int lower_tail, T lp)
439 {
440     /*
441      * So as not to repeat a pnorm call, we expect
442      *
443      *	 lp == pnorm (x, 0, 1, lower_tail, TRUE)
444      *
445      * but use it only in the non-critical case where either x is small
446      * or p==exp(lp) is close to 1.
447      */
448 
449     if (x < 0) {
450 	    x = -x;
451 	    lower_tail = !lower_tail;
452     }
453 
454     if (x > 10 && !lower_tail) {
455 	    T term = 1 / x;
456 	    T sum_ = term;
457 	    T x2 = x * x;
458 	    T i = 1;
459         
460 	    do {
461 	        term *= -i / x2;
462 	        sum_ += term;
463 	        i += 2;
464 	    } while(fabs(term) > EPSILON!T*sum_);
465         
466 	    return 1/sum_;
467     } else {
468 	    T d = dnorm(x, 0., 1., 0);
469 	    return d/exp(lp);
470     }
471 }
472 
473 /*
474  * Asymptotic expansion to calculate the probability that Poisson variate
475  * has value <= x.
476  * Various assertions about this are made (without proof) at
477  * http://members.aol.com/iandjmsmith/PoissonApprox.htm
478  */
479 static T ppois_asymp(T)(T x, T lambda, int lower_tail, int log_p)
480 {
481     static const(T)[8] coefs_a = [
482 	-1e99, /* placeholder used for 1-indexing */
483 	2/3.,
484 	-4/135.,
485 	8/2835.,
486 	16/8505.,
487 	-8992/12629925.,
488 	-334144/492567075.,
489 	698752/1477701225.
490     ];
491 
492     static const(T)[8] coefs_b = [
493 	-1e99, /* placeholder */
494 	1/12.,
495 	1/288.,
496 	-139/51840.,
497 	-571/2488320.,
498 	163879/209018880.,
499 	5246819/75246796800.,
500 	-534703531/902961561600.
501     ];
502 
503     T elfb, elfb_term;
504     T res12, res1_term, res1_ig, res2_term, res2_ig;
505     T dfm, pt_, s2pt, f, np;
506     int i;
507 
508     dfm = lambda - x;
509     /* If lambda is large, the distribution is highly concentrated
510        about lambda.  So representation error in x or lambda can lead
511        to arbitrarily large values of pt_ and hence divergence of the
512        coefficients of this approximation.
513     */
514     pt_ = - log1pmx (dfm / x);
515     s2pt = sqrt (2 * x * pt_);
516     if (dfm < 0)
517     	s2pt = -s2pt;
518 
519     res12 = 0;
520     res1_ig = res1_term = sqrt (x);
521     res2_ig = res2_term = s2pt;
522     for (i = 1; i < 8; i++) {
523 	    res12 += res1_ig * coefs_a[i];
524 	    res12 += res2_ig * coefs_b[i];
525 	    res1_term *= pt_ / i ;
526 	    res2_term *= 2 * pt_ / (2 * i + 1);
527 	    res1_ig = res1_ig / x + res1_term;
528 	    res2_ig = res2_ig / x + res2_term;
529     }
530 
531     elfb = x;
532     elfb_term = 1;
533     for (i = 1; i < 8; i++) {
534 	    elfb += elfb_term * coefs_b[i];
535 	    elfb_term /= x;
536     }
537     if (!lower_tail)
538     	elfb = -elfb;
539 
540     f = res12 / elfb;
541 
542     np = pnorm(s2pt, 0.0, 1.0, !lower_tail, log_p);
543 
544     if (log_p) {
545 	    T n_d_over_p = dpnorm (s2pt, !lower_tail, np);
546 	    return np + log1p (f * n_d_over_p);
547     } else {
548 	    T nd = dnorm (s2pt, 0., 1., log_p);
549 	    return np + f * nd;
550     }
551 } /* ppois_asymp() */
552 
553 T pgamma_raw(T)(T x, T alph, int lower_tail, int log_p)
554 {
555 /* Here, assume that  (x,alph) are not NA  &  alph > 0 . */
556 
557     T res;
558     mixin R_D__1!log_p;
559 
560     mixin R_P_bounds_01!(x, 0., T.infinity);
561 
562     if (x < 1) {
563 	    res = pgamma_smallx!T(x, alph, lower_tail, log_p);
564     } else if (x <= alph - 1 && x < 0.8 * (alph + 50)) {
565 	    /* incl. large alph compared to x */
566 	    T sum = pd_upper_series!T(x, alph, log_p);/* = x/alph + o(x/alph) */
567 	    T d = dpois_wrap!T(alph, x, log_p);
568 	    if (!lower_tail)
569 	        res = log_p ? R_Log1_Exp!T(d + sum) : 1 - d * sum;
570 	    else
571 	        res = log_p ? sum + d : sum * d;
572     } else if (alph - 1 < x && alph < 0.8 * (x + 50)) {
573 	    /* incl. large x compared to alph */
574 	    T sum;
575 	    T d = dpois_wrap!T(alph, x, log_p);
576 	    if (alph < 1) {
577 	        if (x * EPSILON!T > 1 - alph)
578 	    	sum = R_D__1;
579 	        else {
580 	    	T f = pd_lower_cf!T(alph, x - (alph - 1)) * x / alph;
581 	    	/* = [alph/(x - alph+1) + o(alph/(x-alph+1))] * x/alph = 1 + o(1) */
582 	    	sum = log_p ? log (f) : f;
583 	        }
584 	    } else {
585 	        sum = pd_lower_series!T(x, alph - 1);/* = (alph-1)/x + o((alph-1)/x) */
586 	        sum = log_p ? log1p (sum) : 1 + sum;
587 	    }
588 	    if (!lower_tail)
589 	        res = log_p ? sum + d : sum * d;
590 	    else
591 	        res = log_p ? R_Log1_Exp!T(d + sum) : 1 - d * sum;
592     } else { /* x >= 1 and x fairly near alph. */
593 	    res = ppois_asymp!T(alph - 1, x, !lower_tail, log_p);
594     }
595 
596     /*
597      * We lose a fair amount of accuracy to underflow in the cases
598      * where the final result is very close to DBL_MIN.	 In those
599      * cases, simply redo via log space.
600      */
601     if (!log_p && res < MIN!T / EPSILON!T) {
602 	    /* with(.Machine, double.xmin / double.eps) #|-> 1.002084e-292 */
603 	    return exp(pgamma_raw!T(x, alph, lower_tail, 1));
604     } else
605 	    return res;
606 }
607 
608 
609 T pgamma(T: double)(T x, T alph, T scale, int lower_tail, int log_p)
610 {
611     if (isNaN(x) || isNaN(alph) || isNaN(scale))
612 	    return x + alph + scale;
613     if(alph < 0. || scale <= 0.)
614 	    return T.nan;
615     x /= scale;
616 
617     if (isNaN(x)) /* eg. original x = scale = +Inf */
618 	    return x;
619 
620     if(alph == 0.) /* limit case; useful e.g. in pnchisq() */
621 	    return (x <= 0) ? R_DT_0!T(lower_tail, log_p): R_DT_1!T(lower_tail, log_p); /* <= assert  pgamma(0,0) ==> 0 */
622     return pgamma_raw!T(x, alph, lower_tail, log_p);
623 }
624 
625 
626 T qchisq_appr(T)(T p, T nu, T g /* = log Gamma(nu/2) */,
627 		   int lower_tail, int log_p, T tol /* EPS1 */)
628 {
629     enum C7	 = 4.67;
630     enum C8	 = 6.66;
631     enum C9	 = 6.73;
632     enum C10 = 13.32;
633 
634     T alpha, a, c, ch, p1;
635     T p2, q, t, x;
636 
637     /* test arguments and initialise */
638 
639     if (isNaN(p) || isNaN(nu))
640 	   return p + nu;
641     mixin (R_Q_P01_check!(p)());
642     if (nu <= 0)
643     	return T.nan;
644 
645     alpha = 0.5 * nu;/* = [pq]gamma() shape */
646     c = alpha - 1;
647 
648     if(nu < (-1.24)*(p1 = R_DT_log!T(p, lower_tail))) {	/* for small chi-squared */
649 	    /* log(alpha) + g = log(alpha) + log(gamma(alpha)) =
650 	     *        = log(alpha*gamma(alpha)) = lgamma(alpha+1) suffers from
651 	     *  catastrophic cancellation when alpha << 1
652 	     */
653 	    T lgam1pa = (alpha < 0.5) ? lgamma1p!T(alpha) : (log(alpha) + g);
654 	    ch = exp((lgam1pa + p1)/alpha + M_LN2);
655     } else if(nu > 0.32) {	/*  using Wilson and Hilferty estimate */
656         
657 	    x = qnorm!T(p, 0, 1, lower_tail, log_p);
658 	    p1 = 2./(9*nu);
659 	    ch = nu*pow(x*sqrt(p1) + 1-p1, 3);
660         
661 	    /* approximation for p tending to 1: */
662 	    if( ch > 2.2*nu + 6 )
663 	        ch = -2*(R_DT_Clog!T(p, lower_tail, log_p) - c*log(0.5*ch) + g);
664         
665     } else { /* "small nu" : 1.24*(-log(p)) <= nu <= 0.32 */
666         
667 	    ch = 0.4;
668 	    a = R_DT_Clog!T(p, lower_tail, log_p) + g + c*M_LN2;
669         
670 	    do {
671 	        q = ch;
672 	        p1 = 1. / (1 + ch*(C7 + ch));
673 	        p2 = ch*(C9 + ch*(C8 + ch));
674 	        t = -0.5 +(C7 + 2*ch)*p1 - (C9 + ch*(C10 + 3*ch))/p2;
675 	        ch -= (1 - exp(a + 0.5*ch)*p2*p1)/t;
676 	    } while(fabs(q - ch) > tol * fabs(ch));
677     }
678 
679     return ch;
680 }
681 
682 T qgamma(T: double)(T p, T alpha, T scale, int lower_tail, int log_p)
683 /*			shape = alpha */
684 {
685     enum EPS1 = 1e-2;
686     enum EPS2 = 5e-7; /* final precision of AS 91 */
687     enum EPS_N = 1e-15; /* precision of Newton step / iterations */
688     enum LN_EPS = -36.043653389117156; /* = log(.Machine$double.eps) iff IEEE_754 */
689 
690     enum MAXIT = 1000; /* was 20 */
691 
692     enum pMIN = 1e-100; /* was 0.000002 = 2e-6 */
693     enum pMAX = (1-1e-14); /* was (1-1e-12) and 0.999998 = 1 - 2e-6 */
694 
695     mixin R_D__0!log_p;
696 
697     const static T
698 	i420  = 1./ 420.,
699 	i2520 = 1./ 2520.,
700 	i5040 = 1./ 5040;
701 
702     T p_, a, b, c, g, ch, ch0, p1;
703     T p2, q, s1, s2, s3, s4, s5, s6, t, x;
704     int i, max_it_Newton = 1;
705 
706     /* test arguments and initialise */
707 
708     if (isNaN(p) || isNaN(alpha) || isNaN(scale))
709 	    return p + alpha + scale;
710 
711     mixin R_Q_P01_boundaries!(p, 0., T.infinity);
712 
713     if (alpha < 0 || scale <= 0)
714     	return T.nan;
715 
716     if (alpha == 0) /* all mass at 0 : */ 
717     	return 0.;
718 
719     if (alpha < 1e-10) {
720     /* Warning seems unnecessary now: */
721 	    max_it_Newton = 7;/* may still be increased below */
722     }
723 
724     mixin R_DT_qIv!(p);
725     p_ = R_DT_qIv;/* lower_tail prob (in any case) */
726 
727     g = lgammafn!T(alpha);/* log Gamma(v/2) */
728 
729     /*----- Phase I : Starting Approximation */
730     ch = qchisq_appr!T(p, /* nu= 'df' =  */ 2*alpha, /* lgamma(nu/2)= */ g,
731 		     lower_tail, log_p, /* tol= */ EPS1);
732     if(!isFinite(ch)) {
733 	    /* forget about all iterations! */
734 	    max_it_Newton = 0; goto END;
735     }
736     if(ch < EPS2) {/* Corrected according to AS 91; MM, May 25, 1999 */
737 	    max_it_Newton = 20;
738 	    goto END;/* and do Newton steps */
739     }
740 
741     /* FIXME: This (cutoff to {0, +Inf}) is far from optimal
742      * -----  when log_p or !lower_tail, but NOT doing it can be even worse */
743     if(p_ > pMAX || p_ < pMIN) {
744 	    /* did return ML_POSINF or 0.;	much better: */
745 	    max_it_Newton = 20;
746 	    goto END;/* and do Newton steps */
747     }
748 
749 /*----- Phase II: Iteration
750  *	Call pgamma() [AS 239]	and calculate seven term taylor series
751  */
752     c = alpha - 1;
753     s6 = (120 + c*(346 + 127*c)) * i5040; /* used below, is "const" */
754 
755     ch0 = ch;/* save initial approx. */
756     for(i=1; i <= MAXIT; i++ ) {
757 	q = ch;
758 	p1 = 0.5*ch;
759 	p2 = p_ - pgamma_raw!T(p1, alpha, /*lower_tail*/1, /*log_p*/0);
760 
761 	if(!isFinite(p2) || ch <= 0)
762 	{
763 		ch = ch0; max_it_Newton = 27; goto END; 
764 	}/*was  return ML_NAN;*/
765 
766 	t = p2*exp(alpha*M_LN2 + g + p1 - c*log(ch));
767 	b = t/ch;
768 	a = 0.5*t - b*c;
769 	s1 = (210 + a*(140 + a*(105 + a*(84 + a*(70 + 60*a))))) * i420;
770 	s2 = (420 + a*(735 + a*(966 + a*(1141 + 1278*a)))) * i2520;
771 	s3 = (210 + a*(462 + a*(707 + 932*a))) * i2520;
772 	s4 = (252 + a*(672 + 1182*a) + c*(294 + a*(889 + 1740*a))) * i5040;
773 	s5 = (84 + 2264*a + c*(1175 + 606*a)) * i2520;
774 
775 	ch += t*(1 + 0.5*t*s1 - b*c*(s1 - b*(s2 - b*(s3 - b*(s4 - b*(s5 - b*s6))))));
776 	if(fabs(q - ch) < EPS2*ch)
777 	    goto END;
778 	if(fabs(q - ch) > 0.1*ch) {/* diverging? -- also forces ch > 0 */
779 	    if(ch < q) ch = 0.9 * q; else ch = 1.1 * q;
780 	}
781     }
782 /* no convergence in MAXIT iterations -- but we add Newton now... */
783 /* was
784  *    ML_ERROR(ME_PRECISION, "qgamma");
785  * does nothing in R !*/
786 
787 END:
788 /* PR# 2214 :	 From: Morten Welinder <terra@diku.dk>, Fri, 25 Oct 2002 16:50
789    --------	 To: R-bugs@biostat.ku.dk     Subject: qgamma precision
790 
791    * With a final Newton step, double accuracy, e.g. for (p= 7e-4; nu= 0.9)
792    *
793    * Improved (MM): - only if rel.Err > EPS_N (= 1e-15);
794    *		    - also for lower_tail = FALSE	 or log_p = TRUE
795    *		    - optionally *iterate* Newton
796    */
797     x = 0.5*scale*ch;
798     if(max_it_Newton) {
799 	/* always use log scale */
800 	if (!log_p) {
801 	    p = log(p);
802 	    log_p = 1;
803 	}
804 	if(x == 0) {
805 	    const T _1_p = 1. + 1e-7;
806 	    const T _1_m = 1. - 1e-7;
807 	    x = MIN!T;
808 	    p_ = pgamma!T(x, alpha, scale, lower_tail, log_p);
809 	    if(( lower_tail && p_ > p * _1_p) ||
810 	       (!lower_tail && p_ < p * _1_m))
811 		return(0.);
812 	    /* else:  continue, using x = DBL_MIN instead of  0  */
813 	}
814 	else
815 	    p_ = pgamma!T(x, alpha, scale, lower_tail, log_p);
816 	if(p_ == -T.infinity)
817 		return 0; /* PR#14710 */
818 	for(i = 1; i <= max_it_Newton; i++) {
819 	    p1 = p_ - p;
820 	    if(fabs(p1) < fabs(EPS_N * p))
821 		break;
822 	    /* else */
823 	    if((g = dgamma!T(x, alpha, scale, log_p)) == R_D__0) {
824 		    break;
825 	    }
826 	    /* else :
827 	     * delta x = f(x)/f'(x);
828 	     * if(log_p) f(x) := log P(x) - p; f'(x) = d/dx log P(x) = P' / P
829 	     * ==> f(x)/f'(x) = f*P / P' = f*exp(p_) / P' (since p_ = log P(x))
830 	     */
831 	    t = log_p ? p1*exp(p_ - g) : p1/g ;/* = "delta x" */
832 	    t = lower_tail ? x - t : x + t;
833 	    p_ = pgamma!T(t, alpha, scale, lower_tail, log_p);
834 	    if (fabs(p_ - p) > fabs(p1) ||
835 		(i > 1 && fabs(p_ - p) == fabs(p1)) /* <- against flip-flop */) {
836 		    /* no improvement */
837 		    break;
838 	    } /* else : */
839 //#ifdef Harmful_notably_if_max_it_Newton_is_1
840 	    /* control step length: this could have started at
841 	       the initial approximation */
842 	    //if(t > 1.1*x) t = 1.1*x;
843 	    //else if(t < 0.9*x) t = 0.9*x;
844 //#endif
845 	    x = t;
846 	}
847     }
848 
849     return x;
850 }
851 
852 
853 T rgamma(T: double)(T a, T scale)
854 {
855 /* Constants : */
856     const static T sqrt32 = 5.656854;
857     const static T exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
858 
859     /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k))
860      * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k)
861      * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k)
862      */
863     const static T q1 = 0.04166669;
864     const static T q2 = 0.02083148;
865     const static T q3 = 0.00801191;
866     const static T q4 = 0.00144121;
867     const static T q5 = -7.388e-5;
868     const static T q6 = 2.4511e-4;
869     const static T q7 = 2.424e-4;
870 
871     const static T a1 = 0.3333333;
872     const static T a2 = -0.250003;
873     const static T a3 = 0.2000062;
874     const static T a4 = -0.1662921;
875     const static T a5 = 0.1423657;
876     const static T a6 = -0.1367177;
877     const static T a7 = 0.1233795;
878 
879     /* State variables [FIXME for threading!] :*/
880     static T aa = 0.;
881     static T aaa = 0.;
882     static T s, s2, d;    /* no. 1 (step 1) */
883     static T q0, b, si, c;/* no. 2 (step 4) */
884 
885     T e, p, q, r, t, u, v, w, x, ret_val;
886 
887     if (isNaN(a) || isNaN(scale))
888 	    return T.nan;
889     if (a <= 0.0 || scale <= 0.0) {
890 	if(scale == 0. || a == 0.) return 0.;
891 	    return T.nan;
892     }
893     if(!isFinite(a) || !isFinite(scale))
894     	return T.infinity;
895 
896     if (a < 1.) { /* GS algorithm for parameters a < 1 */
897 	    e = 1.0 + exp_m1 * a;
898 	    for(;;) {
899 	        p = e * unif_rand!T();
900 	        if (p >= 1.0) {
901 	    	x = -log((e - p) / a);
902 	    	if (exp_rand!T() >= (1.0 - a) * log(x))
903 	    	    break;
904 	        } else {
905 	    	x = exp(log(p) / a);
906 	    	if (exp_rand!T() >= x)
907 	    	    break;
908 	        }
909 	    }
910 	    return scale * x;
911     }
912 
913     /* --- a >= 1 : GD algorithm --- */
914 
915     /* Step 1: Recalculations of s2, s, d if a has changed */
916     if (a != aa) {
917 	    aa = a;
918 	    s2 = a - 0.5;
919 	    s = sqrt(s2);
920 	    d = sqrt32 - s * 12.0;
921     }
922     /* Step 2: t = standard normal deviate,
923                x = (s,1/2) -normal deviate. */
924 
925     /* immediate acceptance (i) */
926     t = norm_rand!T();
927     x = s + 0.5 * t;
928     ret_val = x * x;
929     if (t >= 0.0)
930 	return scale * ret_val;
931 
932     /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
933     u = unif_rand!T();
934     if (d * u <= t * t * t)
935 	return scale * ret_val;
936 
937     /* Step 4: recalculations of q0, b, si, c if necessary */
938 
939     if (a != aaa) {
940 	    aaa = a;
941 	    r = 1.0 / a;
942 	    q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
943 	           + q2) * r + q1) * r;
944         
945 	    /* Approximation depending on size of parameter a */
946 	    /* The constants in the expressions for b, si and c */
947 	    /* were established by numerical experiments */
948         
949 	    if (a <= 3.686) {
950 	        b = 0.463 + s + 0.178 * s2;
951 	        si = 1.235;
952 	        c = 0.195 / s - 0.079 + 0.16 * s;
953 	    } else if (a <= 13.022) {
954 	        b = 1.654 + 0.0076 * s2;
955 	        si = 1.68 / s + 0.275;
956 	        c = 0.062 / s + 0.024;
957 	    } else {
958 	        b = 1.77;
959 	        si = 0.75;
960 	        c = 0.1515 / s;
961 	    }
962     }
963     /* Step 5: no quotient test if x not positive */
964 
965     if (x > 0.0) {
966 	/* Step 6: calculation of v and quotient q */
967 	v = t / (s + s);
968 	if (fabs(v) <= 0.25)
969 	    q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v
970 				      + a3) * v + a2) * v + a1) * v;
971 	else
972 	    q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
973 
974 
975 	/* Step 7: quotient acceptance (q) */
976 	if (log(1.0 - u) <= q)
977 	    return scale * ret_val;
978     }
979 
980     for(;;) {
981 	/* Step 8: e = standard exponential deviate
982 	 *	u =  0,1 -uniform deviate
983 	 *	t = (b,si)-double exponential (laplace) sample */
984 	e = exp_rand!T();
985 	u = unif_rand!T();
986 	u = u + u - 1.0;
987 	if (u < 0.0)
988 	    t = b - si * e;
989 	else
990 	    t = b + si * e;
991 	/* Step	 9:  rejection if t < tau(1) = -0.71874483771719 */
992 	if (t >= -0.71874483771719) {
993 	    /* Step 10:	 calculation of v and quotient q */
994 	    v = t / (s + s);
995 	    if (fabs(v) <= 0.25)
996 		q = q0 + 0.5 * t * t *
997 		    ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
998 		      + a2) * v + a1) * v;
999 	    else
1000 		q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
1001 	    /* Step 11:	 hat acceptance (h) */
1002 	    /* (if q not positive go to step 8) */
1003 	    if (q > 0.0) {
1004 		w = expm1(q);
1005 		/*  ^^^^^ original code had approximation with rel.err < 2e-7 */
1006 		/* if t is rejected sample again at step 8 */
1007 		if (c * fabs(u) <= w * exp(e - 0.5 * t * t))
1008 		    break;
1009 	    }
1010 	}
1011     } /* repeat .. until  `t' is accepted */
1012     x = s + 0.5 * t;
1013     return scale * x * x;
1014 }
1015 
1016 
1017 void test_gamma()
1018 {
1019 	import std.stdio: write, writeln;
1020 	writeln("dgamma: ", dgamma(1., 1., 1., 0));
1021 	writeln("pgamma: ", pgamma(1., 1., 1., 1, 0));
1022 	writeln("qgamma: ", qgamma(0.5, 1., 1., 1, 0));
1023 	writeln("rgamma: ", rgamma(1., 1.), ", rgamma: ", rgamma(1., 1.) , ", rgamma: ", rgamma(1., 1.));
1024 }