Skip to content

Instantly share code, notes, and snippets.

@certik

certik/fm.c

Created Sep 21, 2012
Embed
What would you like to do?
Incomplete gamma function
/* This function calculates the integral Fm(t), it is based on [1], all
equations are references from there.
The idea is to use series expansion for F_maxm(t) and then the recursive
relation (24) downwards to calculate F_m(t) for m < maxm. For t >= maxt,
the series would take too many iterations to converge (for example with
maxt=20 and eps=1e-17, it could take longer than 82 iterations), so
we calculate F_0(t) directly and use (24) upwards.
[1] I. Shavitt: Methods in Computational Physics (Academic Press Inc., New
York, 1963), vol. 2
*/
double *Fm(int maxm, double t) {
double s, term, *F;
double maxt=20, eps=1e-17;
int m;
F = (double *)malloc((maxm+1)*sizeof(double));
if (t < maxt) {
// Series expansion for F_m(t), between equations (24) and (25)
// The worst case is with maxm=0, then after "m" iterations the "term"
// is:
// term = (2t)^m / (2m+1)!!
// With t=20, it takes m=82 to get term smaller than 1e-17:
// term = (2*20)^82 / (2*82+1)!! = 9.9104e-18 < 1e-17
// So in the worse case we will have 82 iterations.
term = 1.0 / (2*maxm + 1);
s = term;
m = 1;
while (fabs(term) > eps) {
term *= (2*t) / (2*maxm + 2 * m + 1);
s += term;
m++;
}
F[maxm] = s * exp(-t);
// Eq. (24) downwards:
for (m = maxm - 1; m >= 0; m--)
F[m] = (2*t*F[m + 1] + exp(-t)) / (2*m + 1);
} else {
// Eq. for F_0(t) on page 7:
F[0] = 0.5 * sqrt(M_PI/t) * erf(sqrt(t));
// Eq. (24) upwards, for t >= maxt=20, this converges well:
for (m = 0; m <= maxm - 1; m++)
F[m + 1] = ((2*m + 1)*F[m] - exp(-t)) / (2*t);
}
return F;
}
subroutine Fm(maxm, t, F)
! Calculates F_m(t) for m=0,1,..,maxm, where
!
! F_m(t) = \int_0^1 u^(2m) e^(-tu^2) du
!
! This routine is tested for accuracy at least 1e-15. Tests reveal that for any
! "t" one can go up to maxm=50. For t<20 or t>200, one can go up to maxm=500.
! This range should be sufficient for most applications. Besides this range,
! the accuracy might or might not be 1e-15. This test was done with maxt=20 and
! eps=1e-17 below.
!
! This function calculates the integral Fm(t), it is based on [1], all
! equations are references from there.
!
! The idea is to use series expansion for F_maxm(t) and then the recursive
! relation (24) downwards to calculate F_m(t) for m < maxm. For t >= maxt, the
! series would take too many iterations to converge (for example with maxt=20
! and eps=1e-17, it could take longer than 82 iterations), so we calculate
! F_0(t) directly and use (24) upwards.
!
! [1] I. Shavitt: Methods in Computational Physics (Academic Press Inc., New
! York, 1963), vol. 2
integer, intent(in) :: maxm
real(dp), intent(in) :: t
! The array F(0:maxm) will be equal to F(m) = F_m(t) for m = 0..maxm
real(dp), intent(out) :: F(0:)
real(dp), parameter :: maxt=20, eps=1e-17_dp
real(dp) :: s, term
integer :: m
if (ubound(F, 1) /= maxm) call stop_error("Fm: invalid bounds on F")
if (t < maxt) then
! Series expansion for F_m(t), between equations (24) and (25). The worst
! case is with maxm=0, then after "m" iterations the "term" is:
! term = (2t)^m / (2m+1)!!
! With t=20, it takes m=82 to get term smaller than 1e-17:
! term = (2*20)^82 / (2*82+1)!! = 9.9104e-18 < 1e-17
! So in the worse case we will have 82 iterations.
term = 1._dp / (2*maxm + 1)
s = term
m = 1
do while (abs(term) > eps)
term = term * (2*t) / (2*maxm + 2 * m + 1)
s = s + term
m = m + 1
end do
F(maxm) = s * exp(-t)
! Eq. (24) downwards:
do m = maxm-1, 0, -1
F(m) = (2*t*F(m + 1) + exp(-t)) / (2*m + 1)
end do
else
! Eq. for F_0(t) on page 7:
F(0) = 0.5_dp * sqrt(pi/t) * erf(sqrt(t))
! Eq. (24) upwards, for t >= maxt=20, this converges well:
do m = 0, maxm-1
F(m + 1) = ((2*m + 1)*F(m) - exp(-t)) / (2*t)
end do
endif
end subroutine
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment