2D численное интегрирование с бесконечным пределом (C ++)

Для того, чтобы интегрировать двумерную функцию вида

$$\int_{1}^\infty \int_{-\sqrt{x^2-1}}^{\sqrt{x^2-1}} e^{-x} \rm{d}y \rm{d}x,$$

Я пытался использовать следующий код (написанный на C ++), взятый в основном из книги «Числовые рецепты», которая вызывает гауссовскую квадратурную процедуру для интеграции:

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

float quad2d(float (*func)(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)
{
return (*nrfunc)(xsav,y);
}

Этот код прекрасно работает для двумерных интегралов с конечными пределами, но не работает, поскольку внешний предел переводится в бесконечность. Чтобы объяснить это, я попытался использовать изменение переменных:

#define FUNC(x) ((*funk)(-log(x))/(x))

float qgaus(float (*funk)(float), float aa, float bb)
{
int j;
float xr,xm,dx,s,a,b;

b=exp(-aa);
a=0.0;

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 f(float x, float y)
{
float a = exp(-x);
return a;
}

float yy1(float x)
{
float y = -sqrt(x*x-1);
return y;
}

float yy2(float x)
{
float y = sqrt(x*x-1);
return y;
}static float xsav;
static float (*nrfunc)(float, float);

float quad2d(float (*func)(float, float), float x1, float x2)
{
float qgaus(float (*func)(float), float aa, float bb);
float f1(float x);

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

float f1(float x)
{
float qgaus(float (*func)(float), float aa, float bb);
float f2(float y);
float yy1(float);
float yy2(float);

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

float f2(float y)
{
float k = (*nrfunc)(xsav,y);
return k;
}

int main ()
{
float z;
z = quad2d(f, 1.0, 20.0);
cout << z << endl;
}

но это все еще не дает правильный ответ. Так должно быть

$2 \times \rm{BesselK}[1,1] \approx 1.20381$

но вместо этого дает

2.15501

Будем очень благодарны за любые предложения о том, как я мог бы изменить этот код для учета бесконечного лимита!