Ошибка сегментации: ошибка числовой интеграции 11 — 4D (C ++)

Я написал трехмерный код интеграции, который вызывает гауссовскую квадратурную подпрограмму для выполнения интеграции, как показано ниже:

#include <iostream>
#include <cmath>
using namespace std;

float qgaus(float (*func)(float), float a, float b)
{
int j;
float xr,xm,dx,s;
static float x[]={0.0,0.1488743389,0.4333953941,
0.6794095682,0.8650633666,0.9739065285};
static float w[]={0.0,0.2955242247,0.2692667193,
0.2190863625,0.1494513491,0.0666713443};

xm = 0.5*(b+a);
xr = 0.5*(b-a);
s = 0;

for (j=1;j<=5;j++)
{
dx=xr*x[j];
s += w[j]*((*func)(xm+dx)+(*func)(xm-dx));
}
return s *= xr;
}

float func(float x,float y,float z)
{
float f = x*y*z;
return f;
}

float yy1(float x)
{
float y = x;
return y;
}

float yy2(float x)
{
float y = 2*x;
return y;
}

float z1(float x,float y)
{
float z = 5*x*y;
return z;
}

float z2(float x,float y)
{
float z = 10*x*x*y;
return z;
}

static float xsav,ysav;
static float (*nrfunc)(float,float,float);

float quad3d(float (*func)(float, float, float), float x1, float x2)
{
float qgaus(float (*func)(float), float a, float b);
float f1(float x);

nrfunc=func;
return qgaus(f1,x1,x2);
}

float f1(float x)
{
float qgaus(float (*func)(float), float a, float b);
float f2(float y);
float yy1(float),yy2(float);

xsav=x;
return qgaus(f2,yy1(x),yy2(x));
}

float f2(float y)
{
float qgaus(float (*func)(float), float a, float b);
float f3(float z);
float z1(float,float),z2(float,float);

ysav = y;
return qgaus(f3,z1(xsav,y),z2(xsav,y));
}

float f3(float z)
{
return (*nrfunc)(xsav,ysav,z);
}

int main ()
{
float R;
R = quad3d(func, 0, 1);
cout << R << endl;
}

Этот код отлично работает для любой трехмерной функции, с которой я его тестировал. Я попытался изменить его для вычисления четырехмерной функции, заменив 3d-процедуру на 4d:

static float wsav,xsav,ysav;
static float (*nrfunc)(float,float,float,float);

float quad4d(float (*func)(float, float, float, float), float w1, float w2)
{
float qtrap(float (*func)(float), float a, float b);
float f1(float w);

nrfunc=func;
return qtrap(f1,w1,w2);
}

float f1(float w)
{
float qtrap(float (*func)(float), float a, float b);
float f2(float x);
float x1(float),x2(float);

wsav = w;
return qtrap(f2,x1(w),x2(w));
}

float f2(float x)
{
float qtrap(float (*func)(float), float a, float b);
float f3(float y);
float yy1(float,float),yy2(float,float);

xsav = x;
return qtrap(f2,yy1(wsav,x),yy2(wsav,x));
}

float f3(float y)
{
float qtrap(float (*func)(float), float a, float b);
float f4(float z);
float z1(float,float,float),z2(float,float,float);

ysav = y;
return qtrap(f3,z1(wsav,xsav,y),z2(wsav,xsav,y));
}

float f4(float z)
{
float t = (*nrfunc)(wsav,xsav,ysav,z);
return t;
}

Этот код компилируется правильно, но при запуске выдает «Ошибка сегментации: 11». Из того, что я понимаю, это подразумевает, что есть либо какая-то проблема с массивами, либо ошибка выделения памяти, но ни одна из них не имеет смысла, так как в случае 3d не было никаких проблем. Любая помощь с этим будет принята с благодарностью.