#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
/* life insurance calculation */

static int mortality_male[101] = {
    744,    46,    31,    22,    19,    20,    20,    18,    17,    17
,    17,    18,    19,    23,    29,    36,    46,    59,    75,    90
,   103,   111,   114,   113,   112,   115,   120,   126,   132,   139
,   147,   158,   171,   187,   205,   224,   244,   264,   287,   313
,   344,   381,   427,   478,   535,   596,   659,   725,   793,   867
,   946,  1030,  1119,  1211,  1306,  1402,  1503,  1613,  1733,  1865
,  2008,  2162,  2326,  2500,  2687,  2889,  3105,  3340,  3592,  3864
,  4159,  4479,  4824,  5197,  5597,  6034,  6513,  7049,  7668,  8373
,  9149,  9993, 10823, 11645, 12476, 13228, 14256, 15356, 16530, 17781
, 19112, 20526, 22024, 23608, 25279, 27037, 28883, 30816, 32834, 34935
, 37116 
};


static int mortality_female[101] = {
    615,    33,    24,    18,    16,    14,    12,    11,    11,    12
,    12,    12,    13,    15,    18,    23,    27,    28,    27,    28
,    29,    29,    27,    26,    26,    28,    30,    32,    35,    37
,    40,    43,    46,    50,    55,    60,    67,    75,    85,    97
,   112,   129,   148,   169,   191,   214,   237,   263,   290,   319
,   351,   384,   418,   453,   489,   527,   568,   611,   657,   706
,   758,   815,   880,   953,  1036,  1131,  1238,  1360,  1501,  1666
,  1856,  2076,  2323,  2603,  2921,  3282,  3695,  4171,  4724,  5357
,  6060,  6840,  7650,  8499,  9402, 10310, 11440, 12669, 14002, 15444
, 16997, 18666, 20452, 22357, 24383, 26527, 28789, 31166, 33653, 36245
, 38934
};

#define GENDER_MALE   0
#define GENDER_FEMALE 1

static int *mortality[] ={ mortality_male, mortality_female };
static double vn_precomputed[100];

static void precompute_vn_irate(double irate, int n) {
    int i;
    for (i=0; i<=n; i++) {
        vn_precomputed[i] = 1 / pow(1 + irate, (double) i);
    }
}
                                        
static double vn(int n) {
   return vn_precomputed[n];   
}
static double qx(int x, int gender, double addrisk) {
    if (x>100) {
        return 1;
    }
    double q = (double) mortality[gender][x] / 100000.0 + addrisk; 
    return q > 1 ? 1 : q;
}
#if 0
static double px(int x, int gender, double addrisk) {
    return (1 - qx(x, gender, addrisk));
}
#endif

#if 0
static double npx(int x, int gender, double addrisk, int n) {
    int t;
    double p=1;
    for (t=0; t<n; t++) {
        p *= px(x+t, gender, addrisk);
    }    
    return p;
}
#endif

#if 0
static double nAx(int x, int gender, double irate, double addrisk, int n) {
    int i;
    double p;
    double ret = 0;
    for (i=0; i<n; i++) {
        p = npx(x, gender, addrisk, i);
        ret += p * qx(x+i, gender, addrisk) * vn(irate, i+1);
     }
     return ret;
}
#endif

#if 0
static double nEx(int x, int gender, double irate, double addrisk, int n) {
    double p;
    double ret;
    
    p = npx(x, gender, addrisk, n);
    ret = p * vn(irate, n);
    return ret;
}
#endif

#if 0
static double aexn(int x, int gender, double irate, double addrisk, int n) {
    int i;
    double p;
    double ret;
    ret = 0;
    for (i=0; i<n; i++) {
        p = npx(x, gender, addrisk, i);
        ret += p * vn(irate, i);
    }
    return ret;
}
#endif

static void barwert(int x, int gender, double irate, double addrisk, int n, double *abw, double *ebw, double *rbw) {
    int i;
    double p=1;                         
    double q;
    
    *abw = 0;
    *rbw = 0;

    /* nAx, aexn */
    for (i=0; i<n; i++) {
        q = qx(x+i, gender, addrisk);
        *abw += p * q * vn(i+1);
        *rbw += p * vn(i);
        p *= (1-q); 
    }

    /* nEx */
    *ebw = p * vn(n);
}


static double doit(int gender, int age, int n, double irate, double addrisk, double vsum, double* res_out) __attribute__ ((noinline));
static double doit(int gender, int age, int n, double irate, double addrisk, double vsum, double* res_out) {
    int t;                     
    double gamma = 0.004;                                
    double alpha = 0.02;

    double res;
    double abw, ebw, rbw;
    
    double premnorm;
    double prem;
    double totalsum;
                       
    precompute_vn_irate(irate, n);
    
    barwert(age, gender, irate, addrisk, n, &abw, &ebw, &rbw);
    
#if 0
    fprintf(stdout, "abw=%.15f\n", abw);
    fprintf(stdout, "ebw=%.15f\n", ebw);
    fprintf(stdout, "rbw=%.15f\n", rbw);
#endif
    
    /* compute premium */
    premnorm = (abw + ebw + gamma*rbw) / (rbw - alpha);
    prem     = vsum * premnorm;
    totalsum = prem;
                
    /* compute reserve values over time */
    for (t=0; t<=n; t++) {
        barwert(age+t, gender, irate, addrisk, n-t, &abw, &ebw, &rbw);
#if 0
        fprintf(stdout, "abw %d=%.15f\n", t, abw);
        fprintf(stdout, "ebw %d=%.15f\n", t, ebw);
        fprintf(stdout, "rbw %d=%.15f\n", t, rbw);
#endif
        res = abw + ebw;
        if (t<n) {
            res -= premnorm * rbw; 
        }          
        res *= vsum;
        if (res<0) {
            res = 0;
        }
        res_out[t] = res;
        totalsum += res;
    }
    return totalsum;
}
                                          
                                          
int main (int argc, char **argv)  {
    /* input for one calculatoin */
    int gender;      /* 0:male, 1:female                   */
    int age;         /* age in years at begin of contract  */
    int n;           /* contract duration                  */
    double addrisk;  /* additional risk, between 0 and 1   */
    double irate;    /* interest rate, between 0 and 1     */

    /* output */
    double ret;    
    double res[100]; /* maximum duration is 100 years */
    double vsum;

    int i, r;
    long double totalsum=0;  /* sum all computed premiums to compare results */

    clock_t clockStart, clockEnd;
    clockStart = clock();
 
    srand(42);    

#if 0
    ret = doit(GENDER_MALE, 20, 40, 0.035, 0, 100000, res);
    fprintf(stdout, "prem 40,20,0.035,0,100000=%.15f\n", ret);
    for (i=0; i<=40; i++) {
        fprintf(stdout, "res %d=%.15f\n", i, res[i]);
    }
    exit(0);
#endif

    for (i=0; i<100000; i++) {
       /* generate random input for one case */
       gender  = rand() % 2;       /* 0, 1 */
       age     = rand() % 81;      /* 0 .. 80 */
       n       = rand() % 41 + 10; /* 10 .. 50 */
       vsum    = (double) rand() / RAND_MAX * 100000 + 90000;  /* 10000 ... 60000 */
       irate   = rand() % 2 ? 0.025 : 0.035;        /* 0.025 (50%), 0.035 (50%) */
       r = rand();
       addrisk = r%2 ? (double) r / RAND_MAX * 0.3 : 0.0;  /* 0 .. 0.30 (50%), 0 (50%) */
                                               
       /* action: compute premium and reserves */                                               
       ret = doit(gender, age, n, irate, addrisk, vsum, res);
                                                   
       if (totalsum>1e9) {
          fprintf(stdout, "reset, i=%d, totalsum=%.15Lf\n", i, totalsum);
          totalsum=0;  
       }
       totalsum += ret;
    }
    clockEnd = clock();

    fprintf(stdout, "totalsum = %.15Lf\n", totalsum);
    fprintf(stderr, "totalsum = %.15Lf\n", totalsum);
    fprintf(stderr, "took %ld ms\n", (clockEnd - clockStart) / (CLOCKS_PER_SEC / 1000));
    
    return 0;    
}