/* * Mathlib : A C Library of Special Functions * Copyright (C) 2005-6 Morten Welinder * Copyright (C) 2005-6 The R Foundation * Copyright (C) 2006 The R Core Development Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, a copy is available at * http://www.r-project.org/Licenses/ * * SYNOPSIS * * #include * * double pgamma (double x, double alph, double scale, * int lower_tail, int log_p) * * double log1pmx (double x) * double lgamma1p (double a) * * double logspace_add (double logx, double logy) * double logspace_sub (double logx, double logy) * * * DESCRIPTION * * This function computes the distribution function for the * gamma distribution with shape parameter alph and scale parameter * scale. This is also known as the incomplete gamma function. * See Abramowitz and Stegun (6.5.1) for example. * * NOTES * * Complete redesign by Morten Welinder, originally for Gnumeric. * Improvements (e.g. "while NEEDED_SCALE") by Martin Maechler * The old version can be activated by compiling with -DR_USE_OLD_PGAMMA * * REFERENCES * */ #include "nmath.h" #include "dpq.h" /*----------- DEBUGGING ------------- * make CFLAGS='-DDEBUG_p -g -I/usr/local/include -I../include' */ /* Scalefactor:= (2^32)^8 = 2^256 = 1.157921e+77 */ #define SQR(x) ((x)*(x)) static const double scalefactor = SQR(SQR(SQR(4294967296.0))); #undef SQR /* If |x| > |k| * M_cutoff, then log[ exp(-x) * k^x ] =~= -x */ static const double M_cutoff = M_LN2 * DBL_MAX_EXP / DBL_EPSILON;/*=3.196577e18*/ /* Continued fraction for calculation of * 1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ... = sum_{k=0}^Inf x^k/(i+k*d) * * auxilary in log1pmx() and lgamma1p() */ static double logcf (double x, double i, double d, double eps /* ~ relative tolerance */) { double c1 = 2 * d; double c2 = i + d; double c4 = c2 + d; double a1 = c2; double b1 = i * (c2 - i * x); double b2 = d * d * x; double a2 = c4 * c2 - b2; #if 0 assert (i > 0); assert (d >= 0); #endif b2 = c4 * b1 - i * b2; while (fabs(a2 * b1 - a1 * b2) > fabs(eps * b1 * b2)) { double c3 = c2*c2*x; c2 += d; c4 += d; a1 = c4 * a2 - c3 * a1; b1 = c4 * b2 - c3 * b1; c3 = c1 * c1 * x; c1 += d; c4 += d; a2 = c4 * a1 - c3 * a2; b2 = c4 * b1 - c3 * b2; if (fabs (b2) > scalefactor) { a1 /= scalefactor; b1 /= scalefactor; a2 /= scalefactor; b2 /= scalefactor; } else if (fabs (b2) < 1 / scalefactor) { a1 *= scalefactor; b1 *= scalefactor; a2 *= scalefactor; b2 *= scalefactor; } } return a2 / b2; } /* Accurate calculation of log(1+x)-x, particularly for small x. */ double log1pmx (double x) { static const double minLog1Value = -0.79149064; if (x > 1 || x < minLog1Value) return log1p(x) - x; else { /* expand in [x/(2+x)]^2 */ double term = x / (2 + x); double y = term * term; if (fabs(x) < 1e-2) { static const double two = 2; return term * ((((two / 9 * y + two / 7) * y + two / 5) * y + two / 3) * y - x); } else { static const double tol_logcf = 1e-14; return term * (2 * y * logcf (y, 3, 2, tol_logcf) - x); } } } /* Compute log(gamma(a+1)) accurately also for small a (0 < a < 0.5). */ double lgamma1p (double a) { const double eulers_const = 0.5772156649015328606065120900824024; /* coeffs[i] holds (zeta(i+2)-1)/(i+2) , i = 0:(N-1), N = 40 : */ const int N = 40; static const double coeffs[40] = { 0.3224670334241132182362075833230126e-0,/* = (zeta(2)-1)/2 */ 0.6735230105319809513324605383715000e-1,/* = (zeta(3)-1)/3 */ 0.2058080842778454787900092413529198e-1, 0.7385551028673985266273097291406834e-2, 0.2890510330741523285752988298486755e-2, 0.1192753911703260977113935692828109e-2, 0.5096695247430424223356548135815582e-3, 0.2231547584535793797614188036013401e-3, 0.9945751278180853371459589003190170e-4, 0.4492623673813314170020750240635786e-4, 0.2050721277567069155316650397830591e-4, 0.9439488275268395903987425104415055e-5, 0.4374866789907487804181793223952411e-5, 0.2039215753801366236781900709670839e-5, 0.9551412130407419832857179772951265e-6, 0.4492469198764566043294290331193655e-6, 0.2120718480555466586923135901077628e-6, 0.1004322482396809960872083050053344e-6, 0.4769810169363980565760193417246730e-7, 0.2271109460894316491031998116062124e-7, 0.1083865921489695409107491757968159e-7, 0.5183475041970046655121248647057669e-8, 0.2483674543802478317185008663991718e-8, 0.1192140140586091207442548202774640e-8, 0.5731367241678862013330194857961011e-9, 0.2759522885124233145178149692816341e-9, 0.1330476437424448948149715720858008e-9, 0.6422964563838100022082448087644648e-10, 0.3104424774732227276239215783404066e-10, 0.1502138408075414217093301048780668e-10, 0.7275974480239079662504549924814047e-11, 0.3527742476575915083615072228655483e-11, 0.1711991790559617908601084114443031e-11, 0.8315385841420284819798357793954418e-12, 0.4042200525289440065536008957032895e-12, 0.1966475631096616490411045679010286e-12, 0.9573630387838555763782200936508615e-13, 0.4664076026428374224576492565974577e-13, 0.2273736960065972320633279596737272e-13, 0.1109139947083452201658320007192334e-13/* = (zeta(40+1)-1)/(40+1) */ }; const double c = 0.2273736845824652515226821577978691e-12;/* zeta(N+2)-1 */ const double tol_logcf = 1e-14; double lgam; int i; if (fabs (a) >= 0.5) return lgammafn (a + 1); /* Abramowitz & Stegun 6.1.33, * also http://functions.wolfram.com/06.11.06.0008.01 */ lgam = c * logcf(-a / 2, N + 2, 1, tol_logcf); for (i = N - 1; i >= 0; i--) lgam = coeffs[i] - a * lgam; return (a * lgam - eulers_const) * a - log1pmx (a); } /* lgamma1p */ /* * Compute the log of a sum from logs of terms, i.e., * * log (exp (logx) + exp (logy)) * * without causing overflows and without throwing away large handfuls * of accuracy. */ double logspace_add (double logx, double logy) { return fmax2 (logx, logy) + log1p (exp (-fabs (logx - logy))); } /* * Compute the log of a difference from logs of terms, i.e., * * log (exp (logx) - exp (logy)) * * without causing overflows and without throwing away large handfuls * of accuracy. */ double logspace_sub (double logx, double logy) { return logx + log1p (-exp (logy - logx)); } #ifndef R_USE_OLD_PGAMMA /* dpois_wrap (x_P_1, lambda, g_log) == * dpois (x_P_1 - 1, lambda, g_log) */ static double dpois_wrap (double x_plus_1, double lambda, int give_log) { #ifdef DEBUG_p REprintf (" dpois_wrap(x+1=%.14g, lambda=%.14g, log=%d)\n", x_plus_1, lambda, give_log); #endif if (!R_FINITE(lambda)) return R_D__0; if (x_plus_1 > 1) return dpois_raw (x_plus_1 - 1, lambda, give_log); if (lambda > fabs(x_plus_1 - 1) * M_cutoff) return R_D_exp(-lambda - lgammafn(x_plus_1)); else { double d = dpois_raw (x_plus_1, lambda, give_log); #ifdef DEBUG_p REprintf (" -> d=dpois_raw(..)=%.14g\n", d); #endif return give_log ? d + log (x_plus_1 / lambda) : d * (x_plus_1 / lambda); } } /* * Abramowitz and Stegun 6.5.29 [right] */ static double pgamma_smallx (double x, double alph, int lower_tail, int log_p) { double sum = 0, c = alph, n = 0, term; #ifdef DEBUG_p REprintf (" pg_smallx(x=%.12g, alph=%.12g): ", x, alph); #endif /* * Relative to 6.5.29 all terms have been multiplied by alph * and the first, thus being 1, is omitted. */ do { n++; c *= -x / n; term = c / (alph + n); sum += term; } while (fabs (term) > DBL_EPSILON * fabs (sum)); #ifdef DEBUG_p REprintf (" conv.sum=%g;", sum); #endif if (lower_tail) { double f1 = log_p ? log1p (sum) : 1 + sum; double f2; if (alph > 1) { f2 = dpois_raw (alph, x, log_p); f2 = log_p ? f2 + x : f2 * exp (x); } else if (log_p) f2 = alph * log (x) - lgamma1p (alph); else f2 = pow (x, alph) / exp (lgamma1p (alph)); #ifdef DEBUG_p REprintf (" (f1,f2)= (%g,%g)\n", f1,f2); #endif return log_p ? f1 + f2 : f1 * f2; } else { double lf2 = alph * log (x) - lgamma1p (alph); #ifdef DEBUG_p REprintf (" 1:%.14g 2:%.14g\n", alph * log (x), lgamma1p (alph)); REprintf (" sum=%.14g log(1+sum)=%.14g lf2=%.14g\n", sum, log1p (sum), lf2); #endif if (log_p) return R_Log1_Exp (log1p (sum) + lf2); else { double f1m1 = sum; double f2m1 = expm1 (lf2); return -(f1m1 + f2m1 + f1m1 * f2m1); } } } /* pgamma_smallx() */ static double pd_upper_series (double x, double y, int log_p) { double term = x / y; double sum = term; do { y++; term *= x / y; sum += term; } while (term > sum * DBL_EPSILON); /* sum = \sum_{n=1}^ oo x^n / (y*(y+1)*...*(y+n-1)) * = \sum_{n=0}^ oo x^(n+1) / (y*(y+1)*...*(y+n)) * = x/y * (1 + \sum_{n=1}^oo x^n / ((y+1)*...*(y+n))) * ~ x/y + o(x/y) {which happens when alph -> Inf} */ return log_p ? log (sum) : sum; } /* Continued fraction for calculation of * ??? * = (i / d) + o(i/d) */ static double pd_lower_cf (double i, double d) { double f = 0, of; double c1, c2, c3, c4, a1, b1, a2, b2; #define NEEDED_SCALE \ (b2 > scalefactor) { \ a1 /= scalefactor; \ b1 /= scalefactor; \ a2 /= scalefactor; \ b2 /= scalefactor; \ } #define max_it 200000 #ifdef DEBUG_p REprintf("pd_lower_cf(i=%.14g, d=%.14g)\n", i, d); #endif if (i < d * 1e-20) /* includes d = Inf, or i = 0 < d */ return (i/d); a1 = 0; b1 = 1; a2 = i; b2 = d; while NEEDED_SCALE if(a2 == 0) return 0;/* just in case, e.g. d=i=0 */ c2 = a2; c4 = b2; c1 = 0; while (c1 < max_it) { c1++; c2--; c3 = c1 * c2; c4 += 2; a1 = c4 * a2 + c3 * a1; b1 = c4 * b2 + c3 * b1; c1++; c2--; c3 = c1 * c2; c4 += 2; a2 = c4 * a1 + c3 * a2; b2 = c4 * b1 + c3 * b2; if NEEDED_SCALE if (b2 != 0) { of = f; f = a2 / b2; /* convergence check: relative; absolute for small f : */ if (fabs (f - of) <= DBL_EPSILON * fmax2(1., fabs(f))) return f; } } MATHLIB_WARNING(" ** NON-convergence in pgamma()'s pd_lower_cf() f= %g.\n", f); return f;/* should not happen ... */ } /* pd_lower_cf() */ #undef NEEDED_SCALE static double pd_lower_series (double lambda, double y) { double term = 1, sum = 0; #ifdef DEBUG_p REprintf("pd_lower_series(lam=%.14g, y=%.14g) ...", lambda, y); #endif while (y >= 1 && term > sum * DBL_EPSILON) { term *= y / lambda; sum += term; y--; } /* sum = \sum_{n=0}^ oo y*(y-1)*...*(y - n) / lambda^(n+1) * = y/lambda * (1 + \sum_{n=1}^Inf (y-1)*...*(y-n) / lambda^n * ~ y/lambda + o(y/lambda) */ #ifdef DEBUG_p REprintf(" done: term=%g, sum=%g, y= %g\n", term, sum, y); #endif if (y != floor (y)) { /* * The series does not converge as the terms start getting * bigger (besides flipping sign) for y < -lambda. */ double f; #ifdef DEBUG_p REprintf(" y not int: add another term "); #endif /* FIXME: in quite few cases, adding term*f has no effect (f too small) * and unnecessary e.g. for pgamma(4e12, 121.1) */ f = pd_lower_cf (y, lambda + 1 - y); #ifdef DEBUG_p REprintf(" (= %.14g) * term = %.14g to sum %g\n", f, term * f, sum); #endif sum += term * f; } return sum; } /* pd_lower_series() */ /* * Compute the following ratio with higher accuracy that would be had * from doing it directly. * * dnorm (x, 0, 1, FALSE) * ---------------------------------- * pnorm (x, 0, 1, lower_tail, FALSE) * * Abramowitz & Stegun 26.2.12 */ static double dpnorm (double x, int lower_tail, double lp) { /* * So as not to repeat a pnorm call, we expect * * lp == pnorm (x, 0, 1, lower_tail, TRUE) * * but use it only in the non-critical case where either x is small * or p==exp(lp) is close to 1. */ if (x < 0) { x = -x; lower_tail = !lower_tail; } if (x > 10 && !lower_tail) { double term = 1 / x; double sum = term; double x2 = x * x; double i = 1; do { term *= -i / x2; sum += term; i += 2; } while (fabs (term) > DBL_EPSILON * sum); return 1 / sum; } else { double d = dnorm (x, 0, 1, FALSE); return d / exp (lp); } } /* * Asymptotic expansion to calculate the probability that Poisson variate * has value <= x. * Various assertions about this are made (without proof) at * http://members.aol.com/iandjmsmith/PoissonApprox.htm */ static double ppois_asymp (double x, double lambda, int lower_tail, int log_p) { static const double coefs_a[8] = { -1e99, /* placeholder used for 1-indexing */ 2/3., -4/135., 8/2835., 16/8505., -8992/12629925., -334144/492567075., 698752/1477701225. }; static const double coefs_b[8] = { -1e99, /* placeholder */ 1/12., 1/288., -139/51840., -571/2488320., 163879/209018880., 5246819/75246796800., -534703531/902961561600. }; double elfb, elfb_term; double res12, res1_term, res1_ig, res2_term, res2_ig; double dfm, pt_, s2pt, f, np; int i; dfm = lambda - x; /* If lambda is large, the distribution is highly concentrated about lambda. So representation error in x or lambda can lead to arbitrarily large values of pt_ and hence divergence of the coefficients of this approximation. */ pt_ = - log1pmx (dfm / x); s2pt = sqrt (2 * x * pt_); if (dfm < 0) s2pt = -s2pt; res12 = 0; res1_ig = res1_term = sqrt (x); res2_ig = res2_term = s2pt; for (i = 1; i < 8; i++) { res12 += res1_ig * coefs_a[i]; res12 += res2_ig * coefs_b[i]; res1_term *= pt_ / i ; res2_term *= 2 * pt_ / (2 * i + 1); res1_ig = res1_ig / x + res1_term; res2_ig = res2_ig / x + res2_term; } elfb = x; elfb_term = 1; for (i = 1; i < 8; i++) { elfb += elfb_term * coefs_b[i]; elfb_term /= x; } if (!lower_tail) elfb = -elfb; #ifdef DEBUG_p REprintf ("res12 = %.14g elfb=%.14g\n", elfb, res12); #endif f = res12 / elfb; np = pnorm (s2pt, 0.0, 1.0, !lower_tail, log_p); if (log_p) { double n_d_over_p = dpnorm (s2pt, !lower_tail, np); #ifdef DEBUG_p REprintf ("pp*_asymp(): f=%.14g np=e^%.14g nd/np=%.14g f*nd/np=%.14g\n", f, np, n_d_over_p, f * n_d_over_p); #endif return np + log1p (f * n_d_over_p); } else { double nd = dnorm (s2pt, 0., 1., log_p); #ifdef DEBUG_p REprintf ("pp*_asymp(): f=%.14g np=%.14g nd=%.14g f*nd=%.14g\n", f, np, nd, f * nd); #endif return np + f * nd; } } /* ppois_asymp() */ double pgamma_raw (double x, double alph, int lower_tail, int log_p) { /* Here, assume that (x,alph) are not NA & alph > 0 . */ double res; #ifdef DEBUG_p REprintf("pgamma_raw(x=%.14g, alph=%.14g, low=%d, log=%d)\n", x, alph, lower_tail, log_p); #endif R_P_bounds_01(x, 0., ML_POSINF); if (x < 1) { res = pgamma_smallx (x, alph, lower_tail, log_p); } else if (x <= alph - 1 && x < 0.8 * (alph + 50)) { /* incl. large alph compared to x */ double sum = pd_upper_series (x, alph, log_p);/* = x/alph + o(x/alph) */ double d = dpois_wrap (alph, x, log_p); #ifdef DEBUG_p REprintf(" alph `large': sum=pd_upper*()= %.12g, d=dpois_w(*)= %.12g ", sum, d); #endif if (!lower_tail) res = log_p ? R_Log1_Exp (d + sum) : 1 - d * sum; else res = log_p ? sum + d : sum * d; } else if (alph - 1 < x && alph < 0.8 * (x + 50)) { /* incl. large x compared to alph */ double sum; double d = dpois_wrap (alph, x, log_p); #ifdef DEBUG_p REprintf(" x `large': d=dpois_w(*)= %.14g ", d); #endif if (alph < 1) { if (x * DBL_EPSILON > 1 - alph) sum = R_D__1; else { double f = pd_lower_cf (alph, x - (alph - 1)) * x / alph; /* = [alph/(x - alph+1) + o(alph/(x-alph+1))] * x/alph = 1 + o(1) */ sum = log_p ? log (f) : f; } } else { sum = pd_lower_series (x, alph - 1);/* = (alph-1)/x + o((alph-1)/x) */ sum = log_p ? log1p (sum) : 1 + sum; } #ifdef DEBUG_p REprintf(", sum= %.14g\n", sum); #endif if (!lower_tail) res = log_p ? sum + d : sum * d; else res = log_p ? R_Log1_Exp (d + sum) : 1 - d * sum; } else { /* x > 1 and x fairly near alph. */ #ifdef DEBUG_p REprintf(" using ppois_asymp()\n"); #endif res = ppois_asymp (alph - 1, x, !lower_tail, log_p); } /* * We lose a fair amount of accuracy to underflow in the cases * where the final result is very close to DBL_MIN. In those * cases, simply redo via log space. */ if (!log_p && res < DBL_MIN / DBL_EPSILON) { /* with(.Machine, double.xmin / double.eps) #|-> 1.002084e-292 */ #ifdef DEBUG_p REprintf(" very small res=%.14g; -> recompute via log\n", res); #endif return exp (pgamma_raw (x, alph, lower_tail, 1)); } else return res; } double pgamma(double x, double alph, double scale, int lower_tail, int log_p) { #ifdef IEEE_754 if (ISNAN(x) || ISNAN(alph) || ISNAN(scale)) return x + alph + scale; #endif if(alph <= 0. || scale <= 0.) ML_ERR_return_NAN; x /= scale; #ifdef IEEE_754 if (ISNAN(x)) /* eg. original x = scale = +Inf */ return x; #endif return pgamma_raw (x, alph, lower_tail, log_p); } /* From: terra@gnome.org (Morten Welinder) * To: R-bugs@biostat.ku.dk * Cc: maechler@stat.math.ethz.ch * Subject: Re: [Rd] pgamma discontinuity (PR#7307) * Date: Tue, 11 Jan 2005 13:57:26 -0500 (EST) * this version of pgamma appears to be quite good and certainly a vast * improvement over current R code. (I last looked at 2.0.1) Apart from * type naming, this is what I have been using for Gnumeric 1.4.1. * This could be included into R as-is, but you might want to benefit from * making logcf, log1pmx, lgamma1p, and possibly logspace_add/logspace_sub * available to other parts of R. * MM: I've not (yet?) taken logcf(), but the other four */ #else /* R_USE_OLD_PGAMMA */ /* * Copyright (C) 1998 Ross Ihaka * Copyright (C) 1999-2000 The R Development Core Team * Copyright (C) 2003-2004 The R Foundation * based on AS 239 (C) 1988 Royal Statistical Society * * ................ * * NOTES * * This function is an adaptation of Algorithm 239 from the * Applied Statistics Series. The algorithm is faster than * those by W. Fullerton in the FNLIB library and also the * TOMS 542 alorithm of W. Gautschi. It provides comparable * accuracy to those algorithms and is considerably simpler. * * REFERENCES * * Algorithm AS 239, Incomplete Gamma Function * Applied Statistics 37, 1988. */ /* now would need this here: */ double attribute_hidden pgamma_raw(x, alph, lower_tail, log_p) { return pgamma(x, alph, 1, lower_tail, log_p); } double pgamma(double x, double alph, double scale, int lower_tail, int log_p) { const static double xbig = 1.0e+8, xlarge = 1.0e+37, /* normal approx. for alph > alphlimit */ alphlimit = 1e5;/* was 1000. till R.1.8.x */ double pn1, pn2, pn3, pn4, pn5, pn6, arg, a, b, c, an, osum, sum; long n; int pearson; /* check that we have valid values for x and alph */ #ifdef IEEE_754 if (ISNAN(x) || ISNAN(alph) || ISNAN(scale)) return x + alph + scale; #endif #ifdef DEBUG_p REprintf("pgamma(x=%4g, alph=%4g, scale=%4g): ",x,alph,scale); #endif if(alph <= 0. || scale <= 0.) ML_ERR_return_NAN; x /= scale; #ifdef DEBUG_p REprintf("-> x=%4g; ",x); #endif #ifdef IEEE_754 if (ISNAN(x)) /* eg. original x = scale = Inf */ return x; #endif if (x <= 0.) return R_DT_0; #define USE_PNORM \ pn1 = sqrt(alph) * 3. * (pow(x/alph, 1./3.) + 1. / (9. * alph) - 1.); \ return pnorm(pn1, 0., 1., lower_tail, log_p); if (alph > alphlimit) { /* use a normal approximation */ USE_PNORM; } if (x > xbig * alph) { if (x > DBL_MAX * alph) /* if x is extremely large __compared to alph__ then return 1 */ return R_DT_1; else { /* this only "helps" when log_p = TRUE */ USE_PNORM; } } if (x <= 1. || x < alph) { pearson = 1;/* use pearson's series expansion. */ arg = alph * log(x) - x - lgammafn(alph + 1.); #ifdef DEBUG_p REprintf("Pearson arg=%g ", arg); #endif c = 1.; sum = 1.; a = alph; do { a += 1.; c *= x / a; sum += c; } while (c > DBL_EPSILON * sum); } else { /* x >= max( 1, alph) */ pearson = 0;/* use a continued fraction expansion */ arg = alph * log(x) - x - lgammafn(alph); #ifdef DEBUG_p REprintf("Cont.Fract. arg=%g ", arg); #endif a = 1. - alph; b = a + x + 1.; pn1 = 1.; pn2 = x; pn3 = x + 1.; pn4 = x * b; sum = pn3 / pn4; for (n = 1; ; n++) { a += 1.;/* = n+1 -alph */ b += 2.;/* = 2(n+1)-alph+x */ an = a * n; pn5 = b * pn3 - an * pn1; pn6 = b * pn4 - an * pn2; if (fabs(pn6) > 0.) { osum = sum; sum = pn5 / pn6; if (fabs(osum - sum) <= DBL_EPSILON * fmin2(1., sum)) break; } pn1 = pn3; pn2 = pn4; pn3 = pn5; pn4 = pn6; if (fabs(pn5) >= xlarge) { /* re-scale the terms in continued fraction if they are large */ #ifdef DEBUG_p REprintf(" [r] "); #endif pn1 /= xlarge; pn2 /= xlarge; pn3 /= xlarge; pn4 /= xlarge; } } } arg += log(sum); lower_tail = (lower_tail == pearson); if (log_p && lower_tail) return(arg); /* else */ /* sum = exp(arg); and return if(lower_tail) sum else 1-sum : */ return (lower_tail) ? exp(arg) : (log_p ? R_Log1_Exp(arg) : -expm1(arg)); } #endif /* R_USE_OLD_PGAMMA */