The R Project SVN R

Rev

Rev 77685 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 77685 Rev 79203
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
- 
 
 3 
 *  Copyright (C) 2006-2019 The R Core Team
3 
 *  Copyright (C) 2005-6 Morten Welinder <terra@gnome.org>
 4 
 *  Copyright (C) 2005-6 Morten Welinder <terra@gnome.org>
4 
 *  Copyright (C) 2005-10 The R Foundation
 5 
 *  Copyright (C) 2005-10 The R Foundation
5 
 *  Copyright (C) 2006-2015 The R Core Team
 - 
 
6 
 *
 6 
 *
7 
 *  This program is free software; you can redistribute it and/or modify
 7 
 *  This program is free software; you can redistribute it and/or modify
8 
 *  it under the terms of the GNU General Public License as published by
 8 
 *  it under the terms of the GNU General Public License as published by
9 
 *  the Free Software Foundation; either version 2 of the License, or
 9 
 *  the Free Software Foundation; either version 2 of the License, or
10 
 *  (at your option) any later version.
 10 
 *  (at your option) any later version.
Line 65... Line 65...
65 
static const double M_cutoff = M_LN2 * DBL_MAX_EXP / DBL_EPSILON;/*=3.196577e18*/
 65 
static const double M_cutoff = M_LN2 * DBL_MAX_EXP / DBL_EPSILON;/*=3.196577e18*/
66 
 
 66 
 
67 
/* Continued fraction for calculation of
 67 
/* Continued fraction for calculation of
68 
 *    1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ... = sum_{k=0}^Inf x^k/(i+k*d)
 68 
 *    1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ... = sum_{k=0}^Inf x^k/(i+k*d)
69 
 *
 69 
 *
70 
 * auxilary in log1pmx() and lgamma1p()
 70 
 * auxiliary in log1pmx() and lgamma1p()
71 
 */
 71 
 */
72 
static double
 72 
static double
73 
logcf (double x, double i, double d,
 73 
logcf (double x, double i, double d,
74 
       double eps /* ~ relative tolerance */)
 74 
       double eps /* ~ relative tolerance */)
75 
{
 75 
{
Line 143... Line 143...
143 
 
 143 
 
144 
 
 144 
 
145 
/* Compute  log(gamma(a+1))  accurately also for small a (0 < a < 0.5). */
 145 
/* Compute  log(gamma(a+1))  accurately also for small a (0 < a < 0.5). */
146 
double lgamma1p (double a)
 146 
double lgamma1p (double a)
147 
{
 147 
{
- 
 
 148 
    if (fabs (a) >= 0.5)
- 
 
 149 
	return lgammafn (a + 1);
- 
 
 150 
 
148 
    const double eulers_const =	 0.5772156649015328606065120900824024;
 151 
    const double eulers_const =	 0.5772156649015328606065120900824024;
149 
 
 152 
 
150 
    /* coeffs[i] holds (zeta(i+2)-1)/(i+2) , i = 0:(N-1), N = 40 : */
 153 
    /* coeffs[i] holds (zeta(i+2)-1)/(i+2) , i = 0:(N-1), N = 40 : */
151 
    const int N = 40;
 154 
    const int N = 40;
152 
    static const double coeffs[40] = {
 155 
    static const double coeffs[40] = {
Line 192... Line 195...
192 
	0.1109139947083452201658320007192334e-13/* = (zeta(40+1)-1)/(40+1) */
 195 
	0.1109139947083452201658320007192334e-13/* = (zeta(40+1)-1)/(40+1) */
193 
    };
 196 
    };
194 
 
 197 
 
195 
    const double c = 0.2273736845824652515226821577978691e-12;/* zeta(N+2)-1 */
 198 
    const double c = 0.2273736845824652515226821577978691e-12;/* zeta(N+2)-1 */
196 
    const double tol_logcf = 1e-14;
 199 
    const double tol_logcf = 1e-14;
197 
    double lgam;
 - 
 
198 
    int i;
 - 
 
199 
 
 - 
 
200 
    if (fabs (a) >= 0.5)
 - 
 
201 
	return lgammafn (a + 1);
 - 
 
202 
 
 200 
 
203 
    /* Abramowitz & Stegun 6.1.33 : for |x| < 2,
 201 
    /* Abramowitz & Stegun 6.1.33 : for |x| < 2,
204 
     * <==> log(gamma(1+x)) = -(log(1+x) - x) - gamma*x + x^2 * \sum_{n=0}^\infty c_n (-x)^n
 202 
     * <==> log(gamma(1+x)) = -(log(1+x) - x) - gamma*x + x^2 * \sum_{n=0}^\infty c_n (-x)^n
205 
     * where c_n := (Zeta(n+2) - 1)/(n+2)  = coeffs[n]
 203 
     * where c_n := (Zeta(n+2) - 1)/(n+2)  = coeffs[n]
206 
     *
 204 
     *
207 
     * Here, another convergence acceleration trick is used to compute
 205 
     * Here, another convergence acceleration trick is used to compute
208 
     * lgam(x) :=  sum_{n=0..Inf} c_n (-x)^n
 206 
     * lgam(x) :=  sum_{n=0..Inf} c_n (-x)^n
209 
     */
 207 
     */
210 
    lgam = c * logcf(-a / 2, N + 2, 1, tol_logcf);
 208 
    double lgam = c * logcf(-a / 2, N + 2, 1, tol_logcf);
211 
    for (i = N - 1; i >= 0; i--)
 209 
    for (int i = N - 1; i >= 0; i--)
212 
	lgam = coeffs[i] - a * lgam;
 210 
	lgam = coeffs[i] - a * lgam;
213 
 
 211 
 
214 
    return (a * lgam - eulers_const) * a - log1pmx (a);
 212 
    return (a * lgam - eulers_const) * a - log1pmx (a);
215 
} /* lgamma1p */
 213 
} /* lgamma1p */
216 
 
 214