// 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 + 1] = {
  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[MAXAGE + 1] = {
  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 vn_precomputed[MAXAGE];
typedef struct {
    double abw;
    double ebw;
    double rbw;
} struct_barwert_result;
                                            
typedef struct { 
    double irate;
    /* indizes: gender (0/1), age (0..130), duration (0..50) */
    struct_barwert_result results[2][131][51];    
} struct_barwert_addrisk0_cache;
           
static struct_barwert_addrisk0_cache barwert_caches[2];
static int isBarwertPrecomputed=0;

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) {
    double q = mortality[gender][x] + addrisk; 
    return q > 1 ? 1 : q;
}
static double qx_addrisk0(int x, int gender) {
    return mortality[gender][x]; 
}

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;

    int over100 = (MAXAGE + 1) - x;
    int top = over100 < n ? over100 : n;

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

    // skip loop until `n', because in the next steps `p' would be zero anyways
    if (i < n) {
        *abw += p * vn(i+1);
        *rbw += p * vn(i);
        p = 0;
    }

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

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

    int over100 = (MAXAGE + 1) - x;
    int top = over100 < n ? over100 : n;

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

    // skip loop until `n', because in the next steps `p' would be zero anyways
    if (i < n) {
        *abw += p * vn(i+1);
        *rbw += p * vn(i);
        p = 0;
    }

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

static void barwert_addrisk0_fromCache(int x, int gender, double irate, int n, double *abw, double *ebw, double *rbw) {
    int ind_cache;
    struct_barwert_addrisk0_cache *pcache;               
    struct_barwert_result         *presult;
 
    for (ind_cache=0; ind_cache<2; ind_cache++) {
        pcache = &barwert_caches[ind_cache];
        if (pcache->irate==irate) {                
            presult = &(pcache->results[gender][x][n]);
            *abw = presult->abw;
            *ebw = presult->ebw;
            *rbw = presult->rbw;
            return;
        }
    }
    precompute_vn_irate(irate, n);
    barwert_addrisk0(x, gender, irate, n, abw, ebw, rbw);
}


static void precompute_barwert_addrisk0() {
    int ind_cache;
    double irate_arr[2] = { 0.025, 0.035 };
    double irate;                            
    struct_barwert_addrisk0_cache *pcache;
    struct_barwert_result *presults;
    int gender, age, n;
    
    for (ind_cache=0; ind_cache<2; ind_cache++) {
        pcache = &(barwert_caches[ind_cache]);
        irate = irate_arr[ind_cache];
        pcache->irate = irate;
        precompute_vn_irate(irate, MAXAGE-1);

        for (gender=0; gender<=1; gender++) {
            for (age=0; age<=130; age++) {
                for (n=0; n<=50; n++) {         
                    presults = &(pcache->results[gender][age][n]);
                    barwert_addrisk0(age, gender, irate, n, &(presults->abw), &(presults->ebw), &(presults->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;
                       
    precompute_vn_irate(irate, n);
    
    barwert(age, gender, 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++) {                
        barwert(age+t, gender, 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(int gender, int age, int n, double irate, 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;
                       
    if (!isBarwertPrecomputed) {
        precompute_barwert_addrisk0();
        isBarwertPrecomputed = 1;
    }
    
    barwert_addrisk0_fromCache(age, gender, irate, 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++) {                
        barwert_addrisk0_fromCache(age+t, gender, irate, 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(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(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;    
}