I repleace the gauss quadrature with an adaptive Romberg in the 3D quadrature (page 169. in NR in C++) like this:

#include <cstdio>
#include <cmath>
#include "nr.h"

using namespace std;

DP yy1(const DP x)
{
return 0;
}

DP yy2(const DP x)
{
return 2.0;
}

DP z1(const DP x, const DP y)
{
return 0;
}

DP z2(const DP x, const DP y)
{
return 2.0;
}

DP myfunc(const DP x, const DP y, const DP z)
{
return 1.0;
}

namespace NRquad3d {
DP xsav,ysav;
DP (*nrfunc)(const DP, const DP, const DP);

DP f3(const DP z)
{
return nrfunc(xsav,ysav,z);
}

DP f2(const DP y)
{
ysav=y;
return NR::qromb(f3,z1(xsav,y),z2(xsav,y));
}

DP f1(const DP x)
{
xsav=x;
return NR::qromb(f2,yy1(x),yy2(x));
}
}

DP NR::quad3d(DP func(const DP, const DP, const DP), const DP x1, const DP x2)
{
NRquad3d::nrfunc=func;
return qromb(NRquad3d::f1,x1,x2);
}

void main()
{
DP result;

result = NR::quad3d(myfunc, 0, 2.0);

cout << result << endl;
}

Because I'm integrating the function f(x,y,z) = 1, it should give the volume of the integration domain, which is 8.0 (0 < z < 2 ;
0 < y < 2 ; 0 < x < 2) . But it gives 4.5!!!! I just couldn't understand it!