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

static double mortality_male[101] = {
  0.00744, 0.00046, 0.00031, 0.00022, 0.00019, 0.00020, 0.00020, 0.00018, 0.00017, 0.00017
, 0.00017, 0.00018, 0.00019, 0.00023, 0.00029, 0.00036, 0.00046, 0.00059, 0.00075, 0.00090
, 0.00103, 0.00111, 0.00114, 0.00113, 0.00112, 0.00115, 0.00120, 0.00126, 0.00132, 0.00139
, 0.00147, 0.00158, 0.00171, 0.00187, 0.00205, 0.00224, 0.00244, 0.00264, 0.00287, 0.00313
, 0.00344, 0.00381, 0.00427, 0.00478, 0.00535, 0.00596, 0.00659, 0.00725, 0.00793, 0.00867
, 0.00946, 0.01030, 0.01119, 0.01211, 0.01306, 0.01402, 0.01503, 0.01613, 0.01733, 0.01865
, 0.02008, 0.02162, 0.02326, 0.02500, 0.02687, 0.02889, 0.03105, 0.03340, 0.03592, 0.03864
, 0.04159, 0.04479, 0.04824, 0.05197, 0.05597, 0.06034, 0.06513, 0.07049, 0.07668, 0.08373
, 0.09149, 0.09993, 0.10823, 0.11645, 0.12476, 0.13228, 0.14256, 0.15356, 0.16530, 0.17781
, 0.19112, 0.20526, 0.22024, 0.23608, 0.25279, 0.27037, 0.28883, 0.30816, 0.32834, 0.34935
, 0.37116 
};


static double mortality_female[101] = {
  0.00615, 0.00033, 0.00024, 0.00018, 0.00016, 0.00014, 0.00012, 0.00011, 0.00011, 0.00012
, 0.00012, 0.00012, 0.00013, 0.00015, 0.00018, 0.00023, 0.00027, 0.00028, 0.00027, 0.00028
, 0.00029, 0.00029, 0.00027, 0.00026, 0.00026, 0.00028, 0.00030, 0.00032, 0.00035, 0.00037
, 0.00040, 0.00043, 0.00046, 0.00050, 0.00055, 0.00060, 0.00067, 0.00075, 0.00085, 0.00097
, 0.00112, 0.00129, 0.00148, 0.00169, 0.00191, 0.00214, 0.00237, 0.00263, 0.00290, 0.00319
, 0.00351, 0.00384, 0.00418, 0.00453, 0.00489, 0.00527, 0.00568, 0.00611, 0.00657, 0.00706
, 0.00758, 0.00815, 0.00880, 0.00953, 0.01036, 0.01131, 0.01238, 0.01360, 0.01501, 0.01666
, 0.01856, 0.02076, 0.02323, 0.02603, 0.02921, 0.03282, 0.03695, 0.04171, 0.04724, 0.05357
, 0.06060, 0.06840, 0.07650, 0.08499, 0.09402, 0.10310, 0.11440, 0.12669, 0.14002, 0.15444
, 0.16997, 0.18666, 0.20452, 0.22357, 0.24383, 0.26527, 0.28789, 0.31166, 0.33653, 0.36245
, 0.38934
};

#define GENDER_MALE   0
#define GENDER_FEMALE 1

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

static double hoch(double a, int b) {
   return pow(a, (double) b);
}
static double vn(double irate, int n) {
    return 1 / hoch(1 + irate, n);
}
static double qx(int x, int gender, double addrisk) {
    if (x>100) {
        return 1;
    }
    double q = mortality[gender][x] + addrisk; 
    return q > 1 ? 1 : q;
}
static double px(int x, int gender, double addrisk) {
    return (1 - qx(x, gender, addrisk));
}

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

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

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

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

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;
                                                                
    abw = nAx (age, gender, irate, addrisk, n);
    ebw = nEx (age, gender, irate, addrisk, n);
    rbw = aexn(age, gender, irate, addrisk, n);
    
    /* 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++) {
        abw = nAx (age+t, gender, irate, addrisk, n-t);
        ebw = nEx (age+t, gender, irate, addrisk, n-t);
        rbw = aexn(age+t, gender, irate, addrisk, n-t);
        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);    

    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;
    }
    totalsum += doit(0, 50, 50, 0.04, 0.70, 100000, res);
    totalsum += doit(1, 60, 50, 0.02, 0.69, 100000, res);
    totalsum += doit(1, 0, 40, 0.03, 0.0, 100000, res);

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