- for (double r2 = 0.0; r2 < 3000.0; r2 += step_size) {
- double x = r2 + step_size*0.5;
- double probscore = prob_score_real(a, prodai_precompute, (r1 - x)*winfac);
- double z = (x - mu2)/sigma2;
- double gaussian = exp(-(z*z/2.0));
- sum += step_size * probscore * gaussian;
+
+ int n = int(3000.0 / int_step_size + 0.5);
+ double h = 3000.0 / double(n);
+ double sum = evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, 0.0);
+
+ for (int i = 1; i < n; i += 2) {
+ sum += 4.0 * evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, i * h);
+ }
+ for (int i = 2; i < n; i += 2) {
+ sum += 2.0 * evaluate_int_point(a, prodai_precompute, r1, mu2, sigma2, winfac, i * h);