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

#define GENDER_MALE   0
#define GENDER_FEMALE 1

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) {
    double q;
    if (x>100) {
        return 1;
    }
    switch(x + gender*101) {
        case 0  : q=0.00744; break;
        case 1  : q=0.00046; break;
        case 2  : q=0.00031; break;
        case 3  : q=0.00022; break;
        case 4  : q=0.00019; break;
        case 5  : q=0.00020; break;
        case 6  : q=0.00020; break;
        case 7  : q=0.00018; break;
        case 8  : q=0.00017; break;
        case 9  : q=0.00017; break;
        case 10 : q=0.00017; break;
        case 11 : q=0.00018; break;
        case 12 : q=0.00019; break;
        case 13 : q=0.00023; break;
        case 14 : q=0.00029; break;
        case 15 : q=0.00036; break;
        case 16 : q=0.00046; break;
        case 17 : q=0.00059; break;
        case 18 : q=0.00075; break;
        case 19 : q=0.00090; break;
        case 20 : q=0.00103; break;
        case 21 : q=0.00111; break;
        case 22 : q=0.00114; break;
        case 23 : q=0.00113; break;
        case 24 : q=0.00112; break;
        case 25 : q=0.00115; break;
        case 26 : q=0.00120; break;
        case 27 : q=0.00126; break;
        case 28 : q=0.00132; break;
        case 29 : q=0.00139; break;
        case 30 : q=0.00147; break;
        case 31 : q=0.00158; break;
        case 32 : q=0.00171; break;
        case 33 : q=0.00187; break;
        case 34 : q=0.00205; break;
        case 35 : q=0.00224; break;
        case 36 : q=0.00244; break;
        case 37 : q=0.00264; break;
        case 38 : q=0.00287; break;
        case 39 : q=0.00313; break;
        case 40 : q=0.00344; break;
        case 41 : q=0.00381; break;
        case 42 : q=0.00427; break;
        case 43 : q=0.00478; break;
        case 44 : q=0.00535; break;
        case 45 : q=0.00596; break;
        case 46 : q=0.00659; break;
        case 47 : q=0.00725; break;
        case 48 : q=0.00793; break;
        case 49 : q=0.00867; break;
        case 50 : q=0.00946; break;
        case 51 : q=0.01030; break;
        case 52 : q=0.01119; break;
        case 53 : q=0.01211; break;
        case 54 : q=0.01306; break;
        case 55 : q=0.01402; break;
        case 56 : q=0.01503; break;
        case 57 : q=0.01613; break;
        case 58 : q=0.01733; break;
        case 59 : q=0.01865; break;
        case 60 : q=0.02008; break;
        case 61 : q=0.02162; break;
        case 62 : q=0.02326; break;
        case 63 : q=0.02500; break;
        case 64 : q=0.02687; break;
        case 65 : q=0.02889; break;
        case 66 : q=0.03105; break;
        case 67 : q=0.03340; break;
        case 68 : q=0.03592; break;
        case 69 : q=0.03864; break;
        case 70 : q=0.04159; break;
        case 71 : q=0.04479; break;
        case 72 : q=0.04824; break;
        case 73 : q=0.05197; break;
        case 74 : q=0.05597; break;
        case 75 : q=0.06034; break;
        case 76 : q=0.06513; break;
        case 77 : q=0.07049; break;
        case 78 : q=0.07668; break;
        case 79 : q=0.08373; break;
        case 80 : q=0.09149; break;
        case 81 : q=0.09993; break;
        case 82 : q=0.10823; break;
        case 83 : q=0.11645; break;
        case 84 : q=0.12476; break;
        case 85 : q=0.13228; break;
        case 86 : q=0.14256; break;
        case 87 : q=0.15356; break;
        case 88 : q=0.16530; break;
        case 89 : q=0.17781; break;
        case 90 : q=0.19112; break;
        case 91 : q=0.20526; break;
        case 92 : q=0.22024; break;
        case 93 : q=0.23608; break;
        case 94 : q=0.25279; break;
        case 95 : q=0.27037; break;
        case 96 : q=0.28883; break;
        case 97 : q=0.30816; break;
        case 98 : q=0.32834; break;
        case 99 : q=0.34935; break;
        case 100: q=0.37116; break; 
        case 101: q=0.00615; break;
        case 102: q=0.00033; break;
        case 103: q=0.00024; break;
        case 104: q=0.00018; break;
        case 105: q=0.00016; break;
        case 106: q=0.00014; break;
        case 107: q=0.00012; break;
        case 108: q=0.00011; break;
        case 109: q=0.00011; break;
        case 110: q=0.00012; break;
        case 111: q=0.00012; break;
        case 112: q=0.00012; break;
        case 113: q=0.00013; break;
        case 114: q=0.00015; break;
        case 115: q=0.00018; break;
        case 116: q=0.00023; break;
        case 117: q=0.00027; break;
        case 118: q=0.00028; break;
        case 119: q=0.00027; break;
        case 120: q=0.00028; break;
        case 121: q=0.00029; break;
        case 122: q=0.00029; break;
        case 123: q=0.00027; break;
        case 124: q=0.00026; break;
        case 125: q=0.00026; break;
        case 126: q=0.00028; break;
        case 127: q=0.00030; break;
        case 128: q=0.00032; break;
        case 129: q=0.00035; break;
        case 130: q=0.00037; break;
        case 131: q=0.00040; break;
        case 132: q=0.00043; break;
        case 133: q=0.00046; break;
        case 134: q=0.00050; break;
        case 135: q=0.00055; break;
        case 136: q=0.00060; break;
        case 137: q=0.00067; break;
        case 138: q=0.00075; break;
        case 139: q=0.00085; break;
        case 140: q=0.00097; break;
        case 141: q=0.00112; break;
        case 142: q=0.00129; break;
        case 143: q=0.00148; break;
        case 144: q=0.00169; break;
        case 145: q=0.00191; break;
        case 146: q=0.00214; break;
        case 147: q=0.00237; break;
        case 148: q=0.00263; break;
        case 149: q=0.00290; break;
        case 150: q=0.00319; break;
        case 151: q=0.00351; break;
        case 152: q=0.00384; break;
        case 153: q=0.00418; break;
        case 154: q=0.00453; break;
        case 155: q=0.00489; break;
        case 156: q=0.00527; break;
        case 157: q=0.00568; break;
        case 158: q=0.00611; break;
        case 159: q=0.00657; break;
        case 160: q=0.00706; break;
        case 161: q=0.00758; break;
        case 162: q=0.00815; break;
        case 163: q=0.00880; break;
        case 164: q=0.00953; break;
        case 165: q=0.01036; break;
        case 166: q=0.01131; break;
        case 167: q=0.01238; break;
        case 168: q=0.01360; break;
        case 169: q=0.01501; break;
        case 170: q=0.01666; break;
        case 171: q=0.01856; break;
        case 172: q=0.02076; break;
        case 173: q=0.02323; break;
        case 174: q=0.02603; break;
        case 175: q=0.02921; break;
        case 176: q=0.03282; break;
        case 177: q=0.03695; break;
        case 178: q=0.04171; break;
        case 179: q=0.04724; break;
        case 180: q=0.05357; break;
        case 181: q=0.06060; break;
        case 182: q=0.06840; break;
        case 183: q=0.07650; break;
        case 184: q=0.08499; break;
        case 185: q=0.09402; break;
        case 186: q=0.10310; break;
        case 187: q=0.11440; break;
        case 188: q=0.12669; break;
        case 189: q=0.14002; break;
        case 190: q=0.15444; break;
        case 191: q=0.16997; break;
        case 192: q=0.18666; break;
        case 193: q=0.20452; break;
        case 194: q=0.22357; break;
        case 195: q=0.24383; break;
        case 196: q=0.26527; break;
        case 197: q=0.28789; break;
        case 198: q=0.31166; break;
        case 199: q=0.33653; break;
        case 200: q=0.36245; break;
        case 201: q=0.38934; break;
    }
    q += addrisk; 
    return q > 1 ? 1 : q;
}

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(irate, i+1);
        *rbw += p * vn(irate, i);
        p *= (1-q); 
    }

    /* nEx */
    *ebw = p * vn(irate, 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;
                                                                
    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;
}
                                          
                                          
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;    
}