fix lossdistrib_truncated

1 files changed, 18 insertions, 20 deletions

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index 61f3c84..d28b1b2 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -145,7 +145,7 @@ void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaul
     int i, j;
     double* pshocked = malloc(sizeof(double) * (*np) * (*nZ));
 
-#pragma omp parallel for private(j)
+    #pragma omp parallel for private(j)
     for(i = 0; i < *nZ; i++){
         for(j = 0; j < *np; j++){
             pshocked[j + (*np) * i] = shockprob(p[j], rho[j], Z[i], 0);
@@ -172,34 +172,32 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
 
     int i, j, d1, d2, M;
     double lu, d, p1, p2;
-    double *q1, *q2;
+    double *qtemp;
+    int one = 1;
+    double pbar;
+    int bound;
 
     lu = 1./(*N-1);
-    q1 = calloc( *T, sizeof(double));
-    q2 = calloc( *T, sizeof(double));
+    qtemp = malloc( *T * sizeof(double));
     q[0] = 1;
     M = 1;
     for(i=0; i<(*np); i++){
         d = (*defaultflag)? w[i] / lu : S[i] * w[i] / lu;
-        d1 = floor(d);
-        d2 = ceil(d);
-        p1 = p[i] * (d2-d);
+        d1 = (int)floor(d);
+        d2 = d1 + 1;
+        p1 = p[i] * ((double)d2-d);
         p2 = p[i] - p1;
-        for(j=0; j < MIN(M, *T); j++){
-            q1[j] = p1 * q[j];
-            q2[j] = p2 * q[j];
-            q[j] = (1-p[i]) * q[j];
-        }
-        for(j=0; j < MIN(M, *T-d1); j++){
-            q[d1+j] += q1[j];
-        };
-        for(j=0; j < MIN(M, *T-d2); j++){
-            q[d2+j] += q2[j];
-        };
+        memcpy(qtemp, q, MIN(M, *T) * sizeof(double));
+        pbar = 1-p[i];
+        bound = MIN(M, *T);
+        dscal_(&bound, &pbar, q, &one);
+        bound = MIN(M, *T-d1);
+        daxpy_(&bound, &p1, qtemp, &one, q+d1, &one);
+        bound = MIN(M, *T-d2);
+        daxpy_(&bound, &p2, qtemp, &one, q+d2, &one);
         M += d2;
     }
-    free(q1);
-    free(q2);
+    free(qtemp);
 }
 
 void lossdistrib_joint(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {