Integral $\int\limits_0^\infty \prod\limits_{k=0}^\infty\frac{1+\frac{x^2}{(b+1+k)^2}}{1+\frac{x^2}{(a+k)^2}} \ dx$

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.

Responding to the observation/comment by Random Variable from Thu Aug 8, what is missing above to turn this into a rigorous proof is an evaluation of the quantity $$R(a,b) = \sum_{m=a}^b \operatorname{Res}\left(\prod_{k=a}^b \frac{1}{k^2+x^2};x=im\right).$$

We put $$g_1(x) = \prod_{k=a}^b \frac{1}{k+ix} \quad \text{and} \quad g_2(x) = \prod_{k=a}^b \frac{1}{k-ix}.$$

The partial fraction decompositions of $g_1(x)$ and $g_2(x)$ are $$ g_1(x) = \frac{1}{(b-a)!} \sum_{k=a}^b (-1)^{k+1-a} {b-a \choose k-a} \frac{i}{x-ki}$$ and $$ g_2(x) = \frac{1}{(b-a)!} \sum_{k=a}^b (-1)^{k-a} {b-a \choose k-a} \frac{i}{x+ki}.$$ This implies that $$ R(a,b) = \sum_{m=a}^b \frac{i}{(b-a)!} (-1)^{m+1-a} {b-a \choose m-a} \left(\frac{1}{(b-a)!} \sum_{k=a}^b (-1)^{k-a} {b-a \choose k-a} \frac{1}{m+k} \right)$$ which is $$ \frac{i}{((b-a)!)^2} \sum_{m=a}^b \sum_{k=a}^b (-1)^{m+1+k-2a} {b-a \choose m-a} {b-a \choose k-a} \frac{1}{m+k}.$$ Now observe that $$ \sum_{m=a}^b \sum_{k=a}^b (-1)^{m+1+k-2a} {b-a \choose m-a} {b-a \choose k-a} x^{m+k-1} \\ = - x^{2a-1} \sum_{m=a}^b {b-a \choose m-a} (-1)^{m-a} x^{m-a} \sum_{k=a}^b {b-a \choose k-a} (-1)^{k-a} x^{k-a} \\ = - x^{2a-1} (1-x)^{2(b-a)}.$$ Integrating we find that $$ R(a,b) = - i \frac{1}{((b-a)!)^2} \operatorname{B}(2a, 2(b-a)+1).$$ Returning to $I(a,b)$ from the other post, we get that $$I(a,b) = \frac{1}{2} \frac{\Gamma(b+1)^2}{\Gamma(a)^2} \times 2\pi i \times \left(- i \frac{1}{((b-a)!)^2} \operatorname{B}(2a, 2(b-a)+1) \right) \\= \pi \frac{\Gamma(b+1)^2}{\Gamma(a)^2} \frac{1}{((b-a)!)^2} \operatorname{B}(2a, 2(b-a)+1).$$ Switching to gamma functions, this becomes $$ I(a,b) = \pi \frac{\Gamma(b+1)^2}{\Gamma(a)^2} \frac{1}{\Gamma(b-a+1)^2} \frac{\Gamma(2a)\Gamma(2(b-a)+1)}{\Gamma(2b+1)}.$$

To conclude we apply the duplication formula several times, getting $$ \pi \frac{\Gamma(b+1)^2}{\Gamma(a)^2 \Gamma(b-a+1)^2} \frac{\frac{2^{2a-1}}{\sqrt{\pi}} \Gamma(a) \Gamma(a+1/2) \frac{2^{2b-2a}}{\sqrt{\pi}} \Gamma(b-a+1/2) \Gamma(b-a+1)} {\frac{2^{2b}}{\sqrt{\pi}} \Gamma(b+1/2) \Gamma(b+1)},$$ which is $$ \sqrt{\pi} 2^{2a-1+2b-2a-2b} \frac{\Gamma(b+1)}{\Gamma(a) \Gamma(b-a+1)}\Gamma(a+1/2) \frac{ \Gamma(b-a+1/2)}{ \Gamma(b+1/2)},$$ which is indeed $$ \frac{\sqrt{\pi}}{2} \frac{\Gamma(b+1)\Gamma(a+1/2)\Gamma(b-a+1/2)} {\Gamma(a) \Gamma(b-a+1) \Gamma(b+1/2)},$$ as claimed.

To conclude this we actually do the simplification in terms of the $\Gamma$ function. We have $$\prod_{k=0}^{b-a-1}\frac{1}{2a+2k+1} = \frac{1}{2^{b-a}} \frac{\Gamma(a + 1/2)}{\Gamma(b+1/2)}.$$ Furthermore $$\prod_{k=0}^{b-a} (a+k) = \frac{\Gamma(b+1)}{\Gamma(a)}.$$ Finally, $$ \frac{1}{2^{b-a}} {2(b-a) - 1 \choose b-a} = \frac{1}{2^{b-a}} \frac{\Gamma(2b-2a)}{\Gamma(b-a+1)\Gamma(b-a)}.$$ Putting these all together, we have $$I(a, b) = \frac{\pi}{4^{b-a}} \frac{\Gamma(a + 1/2)}{\Gamma(b+1/2)} \frac{\Gamma(b+1)}{\Gamma(a)} \frac{\Gamma(2b-2a)}{\Gamma(b-a+1)\Gamma(b-a)}$$ Now by the duplication formula, $$ \Gamma(2b-2a) = \frac{2^{2(b-a)-1}}{\sqrt\pi} \Gamma(b-a) \Gamma(b-a+1/2)$$ so that finally $$I(a, b) = \frac{\pi}{4^{b-a}} \frac{\Gamma(a + 1/2)}{\Gamma(b+1/2)} \frac{\Gamma(b+1)}{\Gamma(a)} \frac{\Gamma(b-a+1/2)}{\Gamma(b-a+1)} \frac{2^{2(b-a)-1}}{\sqrt\pi} \\ = \frac{\sqrt\pi}{2} \frac{\Gamma(a + 1/2)}{\Gamma(b+1/2)} \frac{\Gamma(b+1)}{\Gamma(a)} \frac{\Gamma(b-a+1/2)}{\Gamma(b-a+1)}.$$

We have verfied the formula for $a$ and $b$ positive integers with $b>a.$ It should now follow by a continuity argument that we can extend it to real $a,b$ with $b>a>1.$