/*
 * Calculate the integral I_n = \int_0^1 x^n e^(-x) dx
 * using the (stable) reverse recurrence relation
 */

#include <stdio.h>
#include <math.h>

#include <gsl/gsl_sf_gamma.h>

#define NTOP  23

int main(void)
{
    int i;
    double tst;
    double idir[NTOP+1];
    double iinv[NTOP+1];
    double iinv2[NTOP+1];

    iinv[NTOP] = 0.;
    iinv2[NTOP] = 1./(NTOP + 1.);
    idir[0] = 1. - 1./M_E;

    for(i = 1; i <= NTOP; i++)
    {
        idir[i] = -1./M_E + i*idir[i-1];
    }

    for(i = NTOP-1; i >= 0; i--)
    {
        iinv[i]  = (1./M_E +  iinv[i+1])/(i+1);
        iinv2[i] = (1./M_E + iinv2[i+1])/(i+1);
    }

    for(i = 1; i <= NTOP; i++)
    {
        tst = gsl_sf_gamma ((double) (i+1)) * gsl_sf_gamma_inc_P ((double) (i+1), 1.);
        printf("%3i   % 24.14f   % 24.14f   % 24.14f   % 24.14f\n",
               i, idir[i], iinv[i], iinv2[i], tst);
    }

    return 0;
}