  Motivated by our previous work on hypergeometric functions and the parbelos
constant, we perform a deeper investigation on the interplay among generalized
complete elliptic integrals, Fourier-Legendre (FL) series expansions, and ${}_p
F_q$ series. We produce new hypergeometric transformations and closed-form
evaluations for new series involving harmonic numbers, through the use of the
integration method outlined as follows: Letting $K$ denote the complete
elliptic integral of the first kind, for a suitable function $g$ we evaluate
integrals such as $$ \int_{0}^{1} K\left( \sqrt{x} \right) g(x) \, dx $$ in two
different ways: (1) by expanding $K$ as a Maclaurin series, perhaps after a
transformation or a change of variable, and then integrating term-by-term; and
(2) by expanding $g$ as a shifted FL series, and then integrating term-by-term.
Equating the expressions produced by these two approaches often gives us new
closed-form evaluations, as in the formulas involving Catalan's constant $G$ $$
\sum _{n = 0}^{\infty } \binom{2 n}{n}^2 \frac{H_{n + \frac{1}{4}} -
  H_{n-\frac{1}{4}}}{16^{n} }
  = \frac{\Gamma^4 \left(\frac{1}{4}\right)}{8 \pi^2}-\frac{4 G}{\pi},$$
  $$ \sum _{m, n \geq 0} \frac{\binom{2 m}{m}^2 \binom{2 n}{n}^2 }{ 16^{m + n}
(m+n+1) (2 m+3) } =
  \frac{7 \zeta (3) - 4 G}{\pi^2}.$$
