Let us first assume $(a,b,c)\in \Bbb {(R^+)}^3$. Then
$$
\int_0^1\int_0^1
\frac{x^a(1-y)^b}{(1-xy)^c}dxdy
{=
\int_0^1\int_0^1
\sum_{n=0}^\infty\frac{\Gamma(n+c)}{n!\Gamma(c)}x^{n+a}(1-y)^by^ndxdy
\\=
\sum_{n=0}^\infty\frac{1}{n+a+1}\frac{\Gamma(n+c)}{n!\Gamma(c)}\int_0^1(1-y)^by^ndy
\\=
\sum_{n=0}^\infty\frac{1}{n+a+1}\frac{\Gamma(n+c)}{n!\Gamma(c)}\frac{\Gamma(b+1)n!}{\Gamma(b+n+2)}
\\=
\frac{\Gamma(b+1)}{\Gamma(c)}
\sum_{n=0}^\infty\frac{1}{n+a+1}\frac{\Gamma(n+c)}{\Gamma(b+n+2)},
}
$$
where we used $(1-u)^{-c}=\sum_{n=0}^\infty\frac{\Gamma(n+c)}{n!\Gamma(c)}u^n$ for $|u|<1$. At this point, due to the non-negativity of the arguments of the summands and integrals, we can deduce the conditions of convergence, which is $c<b+2$. From now on, we assume $(a,b,c)\in (\Bbb Z^{\ge 0})^3$ and $c\ge 1$. Therefore,
$$
\frac{\Gamma(n+c)}{\Gamma(b+n+2)}
{
=
\frac{(n+c-1)!}{(b+n+1)!}
=\sum_{i=0}^\infty\frac{A_i}{b+n+1-i},
}
$$
where
$$
A_i
{=
\lim_{n\to i-b-1}\frac{(b+n+1-i)(n+c-1)!}{(b+n+1)!}
\\=
\lim_{n\to i-b-1}\frac{(b+n+1-i)}{\prod_{j=0}^{b-c+1}(b+n+1-j)}
\\=
\frac{(-1)^{b-c-i+1}}{i!(b-c-i+1)!}.
}
$$
Replacement results in
$$
\int_0^1\int_0^1
\frac{x^a(1-y)^b}{(1-xy)^c}dxdy
{
=
\frac{\Gamma(b+1)}{\Gamma(c)}
\sum_{n=0}^\infty\frac{1}{n+a+1}\frac{\Gamma(n+c)}{\Gamma(b+n+2)}
\\=
\frac{b!}{(c-1)!}
\sum_{i=0}^{b-c+1}\sum_{n=1}^\infty\frac{1}{n+a}\frac{A_i}{b+n-i}
\\=
\frac{b!}{(c-1)!}
\sum_{i=0,i\ne b-a}^{b-c+1}A_i\frac{H_a-H_{b-i}}{a-b+i}
\\+
\frac{b!}{(c-1)!}
A_{b-a}\left(\frac{\pi^2}{6}-H^{(2)}_a\right)\Bbb{I}_{c-1\le a\le b}
\\=
\binom{b}{c-1}
\sum_{i=0,i\ne a-c+1}^{b-c+1}\binom{b-c+1}{i}(-1)^{i}\frac{H_a-H_{c+i-1}}{a-c-i+1}
\\+
\binom{b}{a}\binom{a}{c-1}
(-1)^{a-c+1}\left(\frac{\pi^2}{6}-H^{(2)}_a\right)\Bbb{I}_{c-1\le a\le b}.
}
$$
Therefore, finally
$$
\int_0^1\int_0^1
\frac{x^a(1-y)^b}{(1-xy)^c}dxdy
{
=
\binom{b}{c-1}
\sum_{i=0,i\ne a-c+1}^{b-c+1}\binom{b-c+1}{i}(-1)^{i}\frac{H_a-H_{c+i-1}}{a-c-i+1}
\\+
\binom{b}{a}\binom{a}{c-1}
(-1)^{a-c+1}\left(\frac{\pi^2}{6}-H^{(2)}_a\right)\Bbb{I}_{c-1\le a\le b}.
}
$$
Some special cases
$$
{
\int_0^1\int_0^1
\frac{x^a(1-y)^b}{(1-xy)^{b+1}}dxdy
=
\frac{H_a-H_{b}}{a-b}\Bbb{I}_{a\ne b}
+
\left(\frac{\pi^2}{6}-H^{(2)}_a\right)\Bbb{I}_{a=b}
\\
\int_0^1\int_0^1
\frac{x^a(1-y)^b}{(1-xy)^{a+1}}dxdy
=
\binom{b}{a}\left(\frac{\pi^2}{6}-H^{(2)}_a\right)
+
\binom{b}{a}
\sum_{i=1}^{b-a}\binom{b-a}{i}(-1)^{i}\frac{H_{a+i}-H_a}{i}
}
$$
Python code for checking the integrity of the results
import numpy as np
from math import factorial
from numpy import pi
import matplotlib.pyplot as plt
from tqdm import tqdm
def H(n):
return sum(1/np.arange(1,n+1))
def H2(n):
return sum(1/np.arange(1,n+1)**2)
def binom(n,m):
return factorial(n)/factorial(m)/factorial(n-m)
def GammaFuncValue(a,b,c,N=10):
return sum([1/(n+a+1)/np.prod(np.arange(n+c,n+b+2)) for n in range(0,N)])*factorial(b)/factorial(c-1)
def MonteCarloValue(a,b,c,N=10000):
x=np.random.rand(N)
y=np.random.rand(N)
z=xa*(1-y)b/(1-xy)*c
return sum(z)/N
def ClosedForm(a,b,c):
values=[]
for i in range(b-c+2):
if i==a-c+1:
values.append(
binom(b,a)*binom(a,c-1)*np.cos(np.pi*(a-c+1))*(pi**2/6-H2(a))
)
else:
values.append(
binom(b,c-1)*binom(b-c+1,i)*np.cos(pi*i)*(H(a)-H(c+i-1))/(a-c-i+1)
)
return sum(values)
if name=="main":
n=100
progbar=tqdm(total=n,position=0)
monte_carlo_values=[]
gamma_func_values=[]
closed_form_values=[]
for i in range(n):
a=np.random.randint(1,10)
b=np.random.randint(1,10)
c=np.random.randint(1,b+2)
# print("===========================")
# print(f"a = {a} , b = {b} , c = {c}")
# print(f"MonteCarloValue = {MonteCarloValue(a,b,c,100000)}")
# print(f"GammaFuncValue = {GammaFuncValue(a,b,c,10000)}")
# print(f"ClosedForm = {ClosedForm(a,b,c)}")
monte_carlo_values.append(MonteCarloValue(a,b,c,100000))
gamma_func_values.append(GammaFuncValue(a,b,c,10000))
closed_form_values.append(ClosedForm(a,b,c))
progbar.update(1)
plt.plot(monte_carlo_values)
plt.plot(gamma_func_values)
plt.plot(closed_form_values)
