- return simpson_integrate(ProbScoreEvaluator(k, a, binomial_precompute, r1, mu2, sigma2, winfac), 0.0, 6000.0, int_step_size);
+ int sz = (6000.0 - 0.0) / int_step_size;
+ double h = (6000.0 - 0.0) / sz;
+
+ fftw_plan f1, f2, b;
+ complex<double> *func1, *func2, *res;
+
+ func1 = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
+ func2 = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
+ res = reinterpret_cast<complex<double> *>(fftw_malloc(sz*2*sizeof(complex<double>)));
+ f1 = fftw_plan_dft_1d(sz*2,
+ reinterpret_cast<fftw_complex*>(func1),
+ reinterpret_cast<fftw_complex*>(func1),
+ FFTW_FORWARD,
+ FFTW_MEASURE);
+ f2 = fftw_plan_dft_1d(sz*2,
+ reinterpret_cast<fftw_complex*>(func2),
+ reinterpret_cast<fftw_complex*>(func2),
+ FFTW_FORWARD,
+ FFTW_MEASURE);
+ b = fftw_plan_dft_1d(sz*2,
+ reinterpret_cast<fftw_complex*>(res),
+ reinterpret_cast<fftw_complex*>(res),
+ FFTW_BACKWARD,
+ FFTW_MEASURE);
+
+ // start off by zero
+ for (int i = 0; i < sz*2; ++i) {
+ func1[i].real() = func1[i].imag() = func2[i].real() = func2[i].imag() = 0.0;
+ }
+
+ for (int i = 0; i < sz; ++i) {
+ double x1 = 0.0 + h*i;
+ double z = (x1 - mu2)/sigma2;
+ func1[i].real() = exp(-(z*z/2.0));
+
+ double x2 = -3000.0 + h*i;
+ func2[(i - sz/2 + sz*2)%(sz*2)].real() = prob_score_real(k, a, binomial_precompute, x2*winfac);
+ }
+
+ result.reserve(sz*2);
+
+ // convolve
+ fftw_execute(f1);
+ fftw_execute(f2);
+ for (int i = 0; i < sz*2; ++i) {
+ res[i] = func1[i] * func2[i];
+ }
+ fftw_execute(b);
+ for (int i = 0; i < sz; ++i) {
+ double r1 = i*h;
+ result.push_back(make_pair(r1, abs(res[i])));
+ }