Incomplete gamma function

fm.c

/* 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;
}

fm.f90

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