Method to evaluate an infinite sum of ratio of Gamma functions (how does Mathematica do it?) This question arose from Amdeberhan's question, the evaluation of a double integral, which can be reduced to the evaluation of this series:
$$\sum _{n=0}^{\infty } \frac{\Gamma \left(n+\frac{1}{2}\right)^2 \Gamma \left(n+\frac{s}{2}\right)}{\Gamma (n+1)^2 \Gamma (n+s)}=\frac{\pi ^2 2^{1-s} \Gamma \left(\frac{s}{2}\right)}{\left[\Gamma \left(\frac{3}{4}\right) \Gamma \left(\frac{s}{2} +\frac{1}{4}\right)\right]^2},\;\;{\rm Re}\,s>0.$$
The evaluation of the sum is Mathematica output.

Can someone enlighten me as to how this calculation proceeds?

I went so far as to pay for Wolfram Alpha Pro, hoping that it would disclose the steps, but to no avail. What is even more frustrating is that for $s=1$ the right-hand-side is the square of a complete elliptic integral, which is also recognized immediately by Mathematica and was the original question in the cited post, so far without a conclusive answer.
 A: This is a special case of Watson's Theorem
$$\def\h{\frac{1}{2}}
\def\g#1{\Gamma(#1)\,}
{}_3F_2\left({a,\ b,\ c\atop\h a+\h b+\h, 2c }\biggm| 1 \right)
  =\frac{\g\h\g{c+\h}\g{\h a+\h b +\h}\g{c+\h -\h a -\h b}}
    {\g{\h a +\h} \g{\h b+\h} \g{c+\h -\h a} \g{c+\h -\h b}}.
$$
A: FWIW, Maple (which gets the same result) says this comes from "definite summation using hypergeometric functions".
Hmmm: it looks like this comes from
$$ {}_{3}^{}F_{2}^{} \left(\frac{1}{2},\frac{1}{2},\frac{s}{2};1,s ;1\right)
 = 
\frac{\pi  2^{-s +1} \Gamma \! \left(s \right)}{\Gamma \! \left(\frac{3}{4}\right)^{2} \Gamma \! \left(\frac{1}{4}+\frac{s}{2}\right)^{2}}
$$
A: Let's rewrite the given problem
$$\sum _{n=0}^{\infty } \frac{\Gamma \left(n+\frac{1}{2}\right)^2 \Gamma \left(n+\frac{s}{2}\right)}{\Gamma (n+1)^2 \Gamma (n+s)}=\frac{\pi ^2 2^{1-s} \Gamma \left(\frac{s}{2}\right)}{\left[\Gamma \left(\frac{3}{4}\right) \Gamma \left(\frac{s}{2} +\frac{1}{4}\right)\right]^2}.\tag1$$
Change $s\rightarrow 2s$, then multiply equation (1) through by $\Gamma(s+1)$ to convert the problem to
$$\sum_{n=0}^{\infty}16^{-n}\binom{2n}n^2\binom{2s+n-1}{s+n-1}^{-1}
=\frac{3\pi}{2^{2s+1}}\binom{s}{s-\frac34}\binom{s-1}{s-\frac34}.\tag2$$
Define the two functions
\begin{align*}
F(s,n)&=16^{-n}\binom{2n}n^2\binom{2s+n-1}{s+n-1}^{-1}\frac{2^{2s+1}}{3\pi}\binom{s}{s-\frac34}^{-1}\binom{s-1}{s-\frac34}^{-1} \qquad \text{and} \\
G(s,n)&=\frac{n^2}{s(2s + n)}F(s,n).\end{align*}
It is (mechanically) routine to verify that $F(s+1,n)-F(s,n)=G(s,n+1)-G(s,n)$. Now, sum both sides over all $n\in\mathbb{Z}_{\geq0}$ resulting in the vanishing of the RHS. That means $f(s+1)=f(s)$ where $f(s)=\sum_{n\geq0}F(s,n)$ or $f(s)$ is a constant. To determine the constant, check that $f(1)=1$ because
$$\sum_{n=0}^{\infty}\frac1{16^n(n+1)}\binom{2n}n^2=\frac4{\pi}.$$
The latter evaluation can be justified in different ways.
Remark. This technique is called the Wilf-Zeilberger methodology. It could be what Mathematica or Maple is "hiding" from view. :-)
