Неожиданное отклонение от абсолютного значения в N-мерной интеграции Монте-Карло в переполнении стека

Я написал код на С ++ для Монте-Карло в процессе интеграции, который отлично работал с моими другими функциями, когда я использовал двухмерную интеграцию. Я обобщил код для N-мерной интеграции, в данном конкретном случае я беру n = 10.
Я пытаюсь интегрировать простую функцию f = x1 + x2 + x3 + x4 + x5 + …. + x10, где x1 …. x10 попадает в пределы [-5.0, 5.0]. Я вижу большое отклонение в результате, когда я знаю, что абсолютный результат должен быть 0. Я буду очень признателен, если кто-нибудь любезно посмотрит на мой код и выяснит, где мой код разбивается. Я прилагаю код следующим образом:

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <cstdlib>
#include <ctime>

using namespace std;//define multivariate function F(x1, x2, ...xk)

double f(double x[], int n)
{
double y;
int j;
y = 0.0;
for (j = 0; j < n; j = j+1)
{
y = y + x[j];
}
y = y;

return y;
}
/*
*  Function f(x1, x2, ... xk)
*/

//define function for Monte Carlo Multidimensional integration

double int_mcnd(double(*fn)(double[],int),double a[], double b[], int n, int m)

{
double r, x[n], p;
int i, j;

// initial seed value (use system time)
srand(time(NULL));

r = 0.0;
p = 1.0;

// step 1: calculate the common factor p
for (j = 0; j < n; j = j+1)
{
p = p*(b[j]-a[j]);
}

// step 2: integration
for (i = 1; i <= m; i=i+1)
{
// calculate random x[] points
for (j = 0; j < n; j = j+1)
{
x[j] = a[j] + static_cast <double> (rand()) /( static_cast <double> (RAND_MAX/(b[j]-a[j])));
}
r = r + fn(x,n);
}
r = r*p/m;
return r;
}double f(double[], int);
double int_mcnd(double(*)(double[],int), double[], double[], int, int);

int main(int argc, char **argv)
{/* define how many integrals */
const int n = 10;

double b[n] = {5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0,5.0};
double a[n] = {-5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0,-5.0};

double result;
int i, m;
int N = 20;

cout.precision(6);
cout.setf(ios::fixed | ios::showpoint);

// current time in seconds (begin calculations)
time_t seconds_i;
seconds_i = time (NULL);

m = 2;                // initial number of intervals

// convert command-line input to N = number of points
//N = atoi( argv[1] );

for (i=0; i <=N; i=i+1)
{
result = int_mcnd(f, a, b, n, m);
cout << setw(30)  << m << setw(30) << result <<endl;
m = m*2;
}// current time in seconds (end of calculations)
time_t seconds_f;
seconds_f = time (NULL);
cout << endl << "total elapsed time = " << seconds_f - seconds_i << " seconds" << endl << endl;

return 0;
}

Вывод, который я получаю, выглядит так:

                            4           -41046426010.339691
8           -14557222958.913620
16            25601187040.145161
32            29498213233.367203
64            -2422980618.248888
128           -13400105151.286720
256           -11237568021.855265
512            -5950177645.396674
1024            -4726707072.013641
2048            -1240029475.829825
4096             1890210492.995555
8192              573067706.448856
16384              356227781.143659
32768             -343198855.224271
65536              171823353.999405
131072             -143383711.461758
262144             -197599063.607231
524288              -59641584.846697
1048576               10130826.266767
2097152              100880200.681037

total elapsed time = 1 seconds

что ничто близко к моему ожидаемому выходному нулю. Пожалуйста, помогите мне исправить код и спасибо заранее.