Evaluate $ \int^\infty_0\int^\infty_0 x^a y^{1-a} (1+x)^{-b-1}(1+y)^{-b-1} \exp(-c\frac{x}{y})dxdy $ 
Evaluate
  $$ \int^\infty_0\int^\infty_0 x^a y^{1-a} (1+x)^{-b-1}(1+y)^{-b-1} \exp(-c\frac{x}{y})dxdy $$
  under the condition $a>1$, $b>0$, $c>0$. Note that none of $a$, $b$ and $c$ is integer.

Mathematica found the following form, but I prefer more compact expression for the faster numerical evaluation.
$$
 \bigg(2 c^2 \Gamma (a-2) \Gamma (b-2) \Gamma (b-1) \Gamma (b+1) \Gamma
 (-a+b+1) \, _2F_2(3,b+1;3-a,3-b;-c)
 \nonumber\\ ~~~
 -\pi  \csc (\pi  a)
   \Big(\pi  c^b \Gamma (b+1) \Gamma (2 b-1) (\cot (\pi  (a-b))+\cot
 (\pi  b)) \, _2F_2(b+1,2 b-1;b-1,-a+b+1;-c)
 \nonumber\\ ~~~~~ ~~~~~ 
 +(a-1) a c^a \Gamma (b-1) \Gamma(b-a) \Gamma (-a+b+1) \Gamma (a+b-1) \,
 _2F_2(a+1,a+b-1;a-1,a-b+1;-c)\Big)\Bigg)
 \nonumber\\ ~~~~~ ~~~~~
\bigg/\Big(c^a{\Gamma (b-1) \Gamma (b+1)^2 \Gamma (-a+b+1)}\Big)
$$
 A: $$
\dfrac{d^\alpha}{dc^{\alpha}}\mathrm{e}^{-c\frac{x}{y}} = (-1)^{\alpha}x^\alpha y^{-\alpha}\mathrm{e}^{-c\frac{x}{y}} 
$$
thus
$$
x^\alpha y^{-\alpha}\mathrm{e}^{-c\frac{x}{y}} = (-1)^{\alpha}\dfrac{d^\alpha}{dc^{\alpha}}\mathrm{e}^{-c\frac{x}{y}} 
$$
thus your integral looks like this now
$$
(-1)^{\alpha}\dfrac{d^\alpha}{dc^{\alpha}}\int \int y(1-x)^{-(b+1)}(1-y)^{-(b+1)}\mathrm{e}^{-c\frac{x}{y}} dxdy
$$
so I suggest trying to find the reduced integral (may be simpler to look at to the eye at least) and remember to take the derivative $\alpha$ times with respect to $c$.
A: $\int_0^\infty\int_0^\infty x^ay^{1-a}(1+x)^{-b-1}(1+y)^{-b-1}e^{-c\frac{x}{y}}~dx~dy$
$=\Gamma(a+1)\int_0^\infty y^{1-a}(1+y)^{-b-1}U\left(a+1,a-b+1,\dfrac{c}{y}\right)dy$ (according to http://en.wikipedia.org/wiki/Confluent_hypergeometric_function#Integral_representations)
$=\dfrac{\Gamma(a+1)\Gamma(b-a)}{\Gamma(b+1)}\int_0^\infty y^{1-a}(1+y)^{-b-1}{_1F_1}\left(a+1,a-b+1,\dfrac{c}{y}\right)dy+c^{b-a}\Gamma(a-b)\int_0^\infty y^{1-b}(1+y)^{-b-1}{_1F_1}\left(b+1,b-a+1,\dfrac{c}{y}\right)dy$ (according to http://en.wikipedia.org/wiki/Confluent_hypergeometric_function#Kummer.27s_equation)
$=\dfrac{\Gamma(a+1)\Gamma(b-a)}{\Gamma(b+1)}\int_\infty^0\left(\dfrac{1}{y}\right)^{1-a}\left(1+\dfrac{1}{y}\right)^{-b-1}{_1F_1}(a+1,a-b+1,cy)~d\left(\dfrac{1}{y}\right)+c^{b-a}\Gamma(a-b)\int_\infty^0\left(\dfrac{1}{y}\right)^{1-b}\left(1+\dfrac{1}{y}\right)^{-b-1}{_1F_1}(b+1,b-a+1,cy)~d\left(\dfrac{1}{y}\right)$
$=\dfrac{\Gamma(a+1)\Gamma(b-a)}{\Gamma(b+1)}\int_0^\infty y^{a+b-2}(y+1)^{-b-1}{_1F_1}(a+1,a-b+1,cy)~dy+c^{b-a}\Gamma(a-b)\int_0^\infty y^{2b-2}(y+1)^{-b-1}{_1F_1}(b+1,b-a+1,cy)~dy$
$=\dfrac{\Gamma(a+1)\Gamma(b-a)}{\Gamma(b+1)}\int_0^\infty\sum\limits_{n=0}^\infty\dfrac{(a+1)_nc^ny^{a+b+n-2}(y+1)^{-b-1}}{(a-b+1)_nn!}dy+c^{b-a}\Gamma(a-b)\int_0^\infty\sum\limits_{n=0}^\infty\dfrac{(b+1)_nc^ny^{2b+n-2}(y+1)^{-b-1}}{(b-a+1)_nn!}dy$
$=\dfrac{\Gamma(a+1)\Gamma(b-a)}{\Gamma(b+1)}\sum\limits_{n=0}^\infty\dfrac{(a+1)_nc^nB(a+b+n-1,2-a-n)}{(a-b+1)_nn!}+c^{b-a}\Gamma(a-b)\sum\limits_{n=0}^\infty\dfrac{(b+1)_nc^nB(2b+n-1,2-b-n)}{(b-a+1)_nn!}$ (according to http://en.wikipedia.org/wiki/Beta_function#Properties)
$=\sum\limits_{n=0}^\infty\dfrac{\Gamma(a+1)\Gamma(b-a)(a+1)_n\Gamma(a+b+n-1)\Gamma(2-a-n)c^n}{\Gamma(b+1)\Gamma(b+1)(a-b+1)_nn!}+\sum\limits_{n=0}^\infty\dfrac{\Gamma(a-b)(b+1)_n\Gamma(2b+n-1)\Gamma(2-b-n)c^{b-a+n}}{\Gamma(b+1)(b-a+1)_nn!}$
$=\sum\limits_{n=0}^\infty\dfrac{\Gamma(a+1)\Gamma(2-a)\Gamma(b-a)\Gamma(a+b-1)(a+1)_n(a+b-1)_n(-1)^nc^n}{(\Gamma(b+1))^2(a-1)_n(a-b+1)_nn!}+\sum\limits_{n=0}^\infty\dfrac{\Gamma(a-b)\Gamma(2-b)\Gamma(2b-1)(b+1)_n(2b-1)_n(-1)^nc^{b-a+n}}{\Gamma(b+1)(b-1)_n(b-a+1)_nn!}$ (according to http://en.wikipedia.org/wiki/Pochhammer_symbol#Properties)
$=\dfrac{\Gamma(a+1)\Gamma(2-a)\Gamma(b-a)\Gamma(a+b-1)}{(\Gamma(b+1))^2}{_2F_2}(a+1,a+b-1;a-1,a-b+1;-c)+\dfrac{\Gamma(a-b)\Gamma(2-b)\Gamma(2b-1)c^{b-a}}{\Gamma(b+1)}{_2F_2}(b+1,2b-1;b-1,b-a+1;-c)$
