Does anybody know how to prove this identity?
$$\int_0^\infty \prod_{k=0}^\infty\frac{1+\frac{x^2}{(b+1+k)^2}}{1+\frac{x^2}{(a+k)^2}} \ dx=\frac{\sqrt{\pi}}{2}\frac{\Gamma \left(a+\frac{1}{2}\right)\Gamma(b+1)\Gamma \left(b-a+\frac{1}{2}\right)}{\Gamma(a)\Gamma \left(b+\frac{1}{2}\right)\Gamma(b-a+1)}$$
I found this on Wikipedia.
While this is fairly advanced we can at least bring it into a manageable form which allows specific values to be calculated, assuming $a$ and $b$ are positive integers. Start by evaluating the inner product, which has two parts.
First, $$f_1(x) = \prod_{k=0}^\infty \frac{1}{1+ \frac{x^2}{(a+k)^2}} = \prod_{k=1}^{a-1}\left( 1+ \frac{x^2}{k^2} \right) \prod_{k=1}^\infty \frac{1}{1+ \frac{x^2}{k^2}} = \prod_{k=1}^{a-1} \frac{(k+ix)(k-ix)}{k^2} \frac{\pi x}{\sinh \pi x} \\ = \frac{1}{\Gamma(a)^2} \prod_{k=1}^{a-1} (k+ix)(k-ix) \frac{\pi x}{\sinh \pi x} = \frac{1}{\Gamma(a)^2} \frac{\Gamma(a+ix)\Gamma(a-ix)}{\Gamma(1+ix)\Gamma(1-ix)} \frac{\pi x}{\sinh \pi x} \\ = \frac{\Gamma(a+ix)\Gamma(a-ix)}{\Gamma(a)^2},$$ where we have used Euler's reflection formula in the last step. Next, by the same reasoning, $$f_2(x) = \prod_{k=0}^\infty \left( 1+ \frac{x^2}{(b+1+k)^2}\right) = \frac{\Gamma(b+1)^2}{\Gamma(b+1+ix)\Gamma(b+1-ix)},$$ so the integral becomes $$ I(a, b) = \frac{\Gamma(b+1)^2}{\Gamma(a)^2} \int_0^\infty \frac{\Gamma(a+ix)\Gamma(a-ix)}{\Gamma(b+1+ix)\Gamma(b+1-ix)} dx \\ = \frac{1}{2} \frac{\Gamma(b+1)^2}{\Gamma(a)^2} \int_{-\infty}^\infty \frac{\Gamma(a+ix)\Gamma(a-ix)}{\Gamma(b+1+ix)\Gamma(b+1-ix)} dx.$$ Now suppose that $b>a.$ We have $$ \frac{\Gamma(a+ix)}{\Gamma(b+1+ix)} = \prod_{k=a}^b \frac{1}{k+ix} \quad \text{and} \quad \frac{\Gamma(a-ix)}{\Gamma(b+1-ix)} = \prod_{k=a}^b \frac{1}{k-ix}.$$ We evaluate the inner integral $J(a,b)$ with the Cauchy Residue Theorem using a contour consisting in the limit of the x axis and a half circle in the upper half plane. The poles are simple and located at $x=im$ with $a\le m\le b$ and we get $$ J(a,b) = 2\pi i \sum_{m=a}^b \operatorname{Res} \left(\prod_{k=a}^b \frac{1}{k+ix} \prod_{k=a}^b \frac{1}{k-ix}; x=im \right) \\= 2\pi i \sum_{m=a}^b \operatorname{Res} \left(\prod_{k=a}^b \frac{1}{k^2+x^2}; x=im \right).$$ Using the fact that $$ I(a,b) = \frac{1}{2} \frac{\Gamma(b+1)^2}{\Gamma(a)^2} J(a,b)$$ we may now calculate specific values. Suppose that $b=a+1$ and obtain $$ I(a, a+1) = \frac{\pi}{2} \frac{a(a+1)}{2a+1}.$$ For $b=a+2$ we get $$ I(a, a+2) = \frac{3\pi}{4} \frac{a(a+1)(a+2)}{(2a+1)(2a+3)}.$$ Setting $b = a+3$ we find $$ I(a, a+3) = \frac{5\pi}{4} \frac{a(a+1)(a+2)(a+3)}{(2a+1)(2a+3)(2a+5)}.$$ Setting $b = a+4$ we find $$ I(a, a+4) = \frac{35\pi}{16} \frac{a(a+1)(a+2)(a+3)(a+4)}{(2a+1)(2a+3)(2a+5)(2a+7)}.$$ The pattern here is readily apparent, with the coefficient at the front of $I(a,a+n)$ being $$\frac{1}{2^n}{2n-1\choose n}.$$ These formulas including the factorials in the binomial coefficient may of course be rewritten in terms of $\Gamma$ functions. Edit. Note that the product in the original source becomes finite in this case which might make it more manageable. This latter representation should of course simplify to the same result as the long computation and indeed by the looks of it this transformation should be a trivial computation.