Integral of a power and exponential of a sum of powers: generalizing Gradshteyn and Ryzhik's 3.478(4) Gradshteyn and Ryzhik in 3.478(4) give
$$\int_0^\infty x^{\nu-1} \exp\left( -\beta x^p -\gamma x^{-p} \right) dx = \frac{2}{p} \left(\frac{\gamma}{\beta}\right)^{\nu/2p} K_{\nu/p}\left( 2 \sqrt{\beta\gamma} \right) $$
for $\mathrm{Re}(\beta) > 0$ and $\mathrm{Re}(\gamma) > 0$.
My question is, is there a solution to the integral where we replace the $-p$ power with $-1$? I.e.,
$$\int_0^\infty x^{\nu-1} \exp\left( -\beta x^p -\gamma x^{-1}\right) dx = ?$$
(Or more generally any power, so replacing $\gamma x^{-p}$ in the original with $\gamma x^{-q}$?)

Some background—this integral arises when applying the distribution of the product of two random variables to

*

*a Gamma random variable multiplied by

*a Generalized Gamma random variable (specifically, a Gamma random variable to a power)

as described in this thread on the Stan forum. Specifically, when the second Gamma is raised to some integer power $n$, $p=1/n$ in my question while $q=1$.
Given that we can find closed-form expressions of the densities when $p=q=1$, i.e., when we have the product of two Gamma random variables, via G&R as well as right here on Math Stack Exchange, I'm somewhat hopeful we can find an analytical density for when $p\neq q$.
 A: Thanks to James' note in the comments, I retried computer algebra systems and, for the specific case where $p=1/n$ for integer $n$ and $q=1$, Sympy delivered: the following Python code,
import sympy as s

x, b, g, v, n = s.symbols('x beta gamma nu n', real=True, positive=True)
for power in range(1, 6):
  f = (x**(v - 1) * s.exp(-b * x**(1 / n) - g / x)).subs(n, power)
  res = s.integrate(f, (x, 0, s.oo)).simplify()

  replace = {v: 1.1, b: 0.9, g: .4}

  actual = float(res.subs(replace))
  expected = float(s.Integral(f.subs(replace), (x, 0, s.oo)).evalf())
  print(
      f'$p=1/{power}$: Quadrature: {expected:0.5g}, analytic: {actual:0.5g}, relative error: {(actual - expected) / expected:0.5g}'
  )
  print(f'\n$${s.latex(res)}$$\n')

prints the following: for a specific $\nu$, $\beta$ and $\gamma$, the analytical result as well as that found with quadrature integration:
$p=1/1$: Quadrature: 0.59315, analytic: 0.59315, relative error: 9.3587e-16
$$\frac{\pi \beta^{- \frac{\nu}{2}} \gamma^{\frac{\nu}{2}} \left(I_{- \nu}\left(2 \sqrt{\beta} \sqrt{\gamma}\right) - I_{\nu}\left(2 \sqrt{\beta} \sqrt{\gamma}\right)\right)}{\sin{\left(\pi \nu \right)}}$$
$p=1/2$: Quadrature: 2.2497, analytic: 2.2497, relative error: 0
$$\frac{\left(\frac{4}{\beta^{2}}\right)^{\nu} {G_{3, 0}^{0, 3}\left(\begin{matrix} 1 - \nu, \frac{1}{2} - \nu, 1 &  \\ &  \end{matrix} \middle| {\frac{4}{\beta^{2} \gamma}} \right)}}{\sqrt{\pi}}$$
$p=1/3$: Quadrature: 10.768, analytic: 10.768, relative error: 1.6496e-16
$$\frac{3^{3 \nu + \frac{1}{2}} \beta^{- 3 \nu} {G_{4, 0}^{0, 4}\left(\begin{matrix} 1 - \nu, \frac{2}{3} - \nu, \frac{1}{3} - \nu, 1 &  \\ &  \end{matrix} \middle| {\frac{27}{\beta^{3} \gamma}} \right)}}{2 \pi}$$
$p=1/4$: Quadrature: 63.698, analytic: 63.698, relative error: 0
$$\frac{2^{8 \nu - \frac{1}{2}} \beta^{- 4 \nu} {G_{5, 0}^{0, 5}\left(\begin{matrix} 1 - \nu, \frac{3}{4} - \nu, \frac{1}{2} - \nu, \frac{1}{4} - \nu, 1 &  \\ &  \end{matrix} \middle| {\frac{256}{\beta^{4} \gamma}} \right)}}{\pi^{\frac{3}{2}}}$$
$p=1/5$: Quadrature: 466.28, analytic: 466.28, relative error: -7.3145e-16
$$\frac{5^{5 \nu + \frac{1}{2}} \beta^{- 5 \nu} {G_{6, 0}^{0, 6}\left(\begin{matrix} 1 - \nu, \frac{4}{5} - \nu, \frac{3}{5} - \nu, \frac{2}{5} - \nu, \frac{1}{5} - \nu, 1 &  \\ &  \end{matrix} \middle| {\frac{3125}{\beta^{5} \gamma}} \right)}}{4 \pi^{2}}$$
The resulting expression seems to well-match numerical integration for the particular values of $\nu$, $\beta$ and $\gamma$ I picked here (I'm assuming that Sympy's evaluation of the Meijer G function isn't itself using numerical integration).
A: Let us take $x >0 $ and $y>0$ and $\alpha >0$, $\gamma_1>0$ and $\gamma_2>0$.
Then we define:
\begin{equation}
{\mathfrak J}_{\alpha,\gamma_1,\gamma_2}(x,y):= \int\limits_0^\infty t^{\alpha-1} \exp(-y t^{\gamma_1}) \exp(-\frac{x}{t^{\gamma_2}}) dt
\end{equation}
Then we have:
\begin{eqnarray}
&&{\mathfrak J}_{\alpha,\gamma_1,\gamma_2}(x,y) = \tag{1} \\
&&\frac{y^{-\frac{\alpha}{\gamma_1}}}{\gamma_1} \cdot \frac{1}{2\pi \imath} \int\limits_{-\imath \infty+ \zeta}^{\imath \infty +\zeta} \Gamma\left( \frac{\alpha+\gamma_2 s}{ \gamma_1} \right) \cdot \Gamma\left(s\right)  \cdot \left[ x y^{\frac{\gamma_2}{\gamma_1}} \right]^{-s} ds = \tag{2} \\
&& \frac{y^{-\frac{\alpha}{\gamma_1}}}{\gamma_1} \cdot \sum\limits_{n=0}^{\infty} \frac{(-1)^n}{n!} \left[
\frac{\gamma_1}{\gamma_2} \cdot \Gamma[-\frac{n \gamma_1+\alpha}{\gamma_2}] \cdot (x y^{\frac{\gamma_2}{\gamma_1}})^{\frac{n \gamma_1+\alpha}{\gamma_2}}
+
\Gamma[\frac{\alpha-n \gamma_2}{\gamma_1}] \cdot (x y^{\frac{\gamma_2}{\gamma_1}})^n
\right] \tag{3}
\end{eqnarray}
In $(2)$ we expressed the integrand through the inverse Mellin transform, which boils down to using the approach from here the first answer to this question. Here $\zeta > 0 $.Then in $(3)$ we used Cauchy theorem to compute the integral in question. The integral is then equal to the sum of the residues of the integrand at poles of order one being located at firstly $(\alpha + \gamma_2 s_n)/\gamma_1 = - n $ and secondly at $s_n = - n $ where $n \in {\mathbb N} $.
In[632]:= {alpha, g1, g2} = 
 RandomReal[{1, 2}, 3]; M = 20;(*g1=1;g2=1;*)
{x, y} = RandomReal[{1, 2}, 2];
NIntegrate[t^(alpha - 1) Exp[-y t^g1] Exp[-x/t^g2], {t, 0, Infinity}]
ofst = -alpha/2;
y^(-alpha/g1) /
  g1 NIntegrate[ 
   Gamma[(alpha + g2 (I s - ofst))/g1] Gamma[
     I s - ofst] (x y^(g2/g1))^(-I s + ofst), {s, -Infinity, 
    Infinity}]/(2 Pi)
Take[Accumulate[
  y^(-alpha/g1) /
    g1 Table[(-1)^n/n! ( 
      g1/g2 Gamma[-(n g1 + alpha)/
          g2] (x y^(g2/g1))^((n g1 + alpha)/g2) + 
       Gamma[(alpha - n g2)/g1] (x y^(g2/g1))^n), {n, 0, M}]], -5]


Out[633]= 0.137716

Out[635]= 0.137716 + 0. I

Out[636]= {0.137716, 0.137716, 0.137716, 0.137716, 0.137716}

