// vim: sw=4 ts=2 sta et
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
/* life insurance calculation */

#define MAXAGE 100
static double mortality_male[MAXAGE + 2] = {
  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, 1.0
};

static double mortality_female[MAXAGE + 2] = {
  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, 1.0
};
              
static double mortality_max[2] = { 0.37116, 0.38934 };

#define GENDER_MALE   0
#define GENDER_FEMALE 1

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

static double irate_cache[2] = { 0.025, 0.035 };
/* indizes: irate-ind, gender (0/1), age (0..80) */
/* results: [prem_n10 res_n10(0)..res_n10(10) .... prem_n50 res_n50 resn50(0) ... res_n50(50)] */
/*          duration 10..50, duration+2 values per duration */
/*           -> 12 + 13 + ... + 52 = 52*53/2 - 11*12/2 = 1312 values */
double results_cache[2][2][81][1312] =
#include "pre.txt"
;
/* indizes: irate-ind, x (0..50) */
static double vn_cache[2][51] =
#include "prev.txt"
;
static double vn_precomputed[51];

static double _pow(double v, int p) {
    switch(p) {
        case 0: return 1;
	case 1: return v;
	case 2: return v*v;
        case 3: return v*v*v;
        default: return exp(p * log(v));
    }
}

static double *precompute_vn_irate(double irate, int n) {
    int i, ind_cache;
    for (ind_cache=0; ind_cache<2; ind_cache++) {
        if (irate_cache[ind_cache]==irate) {
            return vn_cache[ind_cache];
        }
    }    

    double t = 1 + irate;
    vn_precomputed[0] = 1;
    for (i=1; i<=n; i++) {
        vn_precomputed[i] = 1 / _pow(t, (double) i);
    }                  
    return vn_precomputed;
}

static double qx(double *qx_values, int x, double addrisk) {
    double q = qx_values[x] + addrisk; 
    return q > 1 ? 1 : q;
}                           

static void barwert(int x, int gender, double* vn_irate, double addrisk, int n, double *abw, double *ebw, double *rbw) {
    int i;
    int age;

    double p;
    double q;
    double vn0, vn1;
    
    double ret_abw;
    double ret_rbw; 
    
    double *qx_values;
                   
    if (n==0) {
        *abw=0;
        *rbw=0;
        *ebw=1;
        return;
    } else if (x>100) {
        *abw=vn_irate[1];
        *rbw=1;
        *ebw=0;
        return;
    }                 
    
    qx_values = mortality[gender];
    p=1;                         
                   
    if (x+n>102) {
        n = 102-x;
    }           
    if (mortality_max[gender] + addrisk <= 1) {
        qx_values[MAXAGE+1] = 1.0 - addrisk;
        q = qx_values[x] + addrisk;
        vn0 = vn_irate[1];
        ret_abw = q * vn0;
        ret_rbw = 1;
        p = (1-q);
        q = qx_values[x+1] + addrisk;
        for (i=1; i < n; i++) {
            vn1 = vn_irate[i+1];
            ret_abw += p * q * vn1;
            ret_rbw += p * vn0;
            p *= (1-q); 
            vn0 = vn1;
            q = qx_values[x+i+1] + addrisk;
        }
    } else {
        ret_abw = 0;
        ret_rbw = 0; 
        vn0 = 1;
        qx_values[MAXAGE+1] = 1.0;
        for (i=0, age=x; i < n; i++, age++) {
            q = qx(qx_values, age, addrisk);
            vn1 = vn_irate[i+1];
            ret_abw += p * q * vn1;
            ret_rbw += p * vn0;
            p *= (1-q); 
            vn0 = vn1;
        }
    }
    *ebw = p * vn_irate[n];
    *abw = ret_abw;
    *rbw = ret_rbw;            
}

static double doit_common(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;
    
    double *vn_irate = precompute_vn_irate(irate, n);
    
    barwert(age, gender, vn_irate, addrisk, n, &abw, &ebw, &rbw);
    
    /* 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++) {                
        if (t>0) {
            barwert(age+t, gender, vn_irate, addrisk, n-t, &abw, &ebw, &rbw);
        }
        res = abw + ebw;
        if (t<n) {
            res -= premnorm * rbw; 
        }          
        res *= vsum;
        if (res<0) {
            res = 0;
        }
        res_out[t] = res;
        totalsum += res;
    }
    return totalsum;
}

static double doit_addrisk0_postcache(double *pre, int n, double vsum, double* res_out) {
    int t;                     
    double res;
    double premnorm;
    double prem;
    double totalsum;
    
    /* compute premium */
    premnorm = pre[0];
    prem     = vsum * premnorm;
    totalsum = prem;
                
    /* compute reserve values over time */
    for (t=0; t<=n; t++) {                
        res = pre[t+1] * vsum;
        res_out[t] = res;
        totalsum += res;
    }
    return totalsum;
}

static double doit_addrisk0_withcache(int gender, int age, int n, double irate, double vsum, double* res_out) {
    int ind_cache;
    double *pre;                                       
    int start;

    for (ind_cache=0; ind_cache<2; ind_cache++) {
        if (irate_cache[ind_cache]==irate) {
            pre = results_cache[ind_cache][gender][age];
            /* startposition of data: 12 + 13 + ... + (n+2-1) = (n+2-1)*(n+2) - (11*12)/2 */
            start = (n+1)*(n+2)/2 - 66;
            return doit_addrisk0_postcache(pre+start, n, vsum, res_out);
        }
    }
    return doit_common (gender, age, n, irate, 0, vsum, res_out);
}
    
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) {
    return addrisk==0.0 
           ? doit_addrisk0_withcache(gender, age, n, irate, vsum, res_out)
           : doit_common  (gender, age, n, irate, addrisk, vsum, res_out);
}
                                          
                                          
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[MAXAGE]; /* 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;    
}