Please don't be frightened by the length of this, I just wanted to provide ample detail. If you want, you can skip the derivation and go straight to the result at the bottom.
I have
$ \begin{multline}\int_{-1}^{1}\left(\sum_{n_1-l_1=0}^{\infty}C_{n_1-l_1}^{l_1+1}(x)r^{n_1-l_1}\right)\left(\sum_{n_2-l_2=0}^{\infty}C_{n_2-l_2}^{l_2+1}(x)r^{n_1-l_1}\right)\times\\ \left(\sum_{n_3-l_3=0}^{\infty}C_{n_3-l_3}^{l_3+1}(x)r^{n_3-l_3}\right)(1-x^2)^{(l_1+l_2+l_3+1)/2}dx \end{multline} $ Applying the identity $\left(\sum_{n=0}^\infty a_nr^n\right)\left(\sum_{n=0}^\infty b_nr^n\right)=\sum_{n=0}^\infty\sum_{k=0}^n a_kb_{n-k}r^n$ twice to the above gives: $ =\sum_{i=0}^{\infty}\sum_{j=0}^{i}\sum_{k=0}^{i-j}r^{i}\int_{-1}^{1}C_{j}^{l_1+1}(x)C_{k}^{l_2+1}(x)C_{i-j-k}^{l_3+1}(x)(1-x^2)^{(l_1+l_2+l_3+1)/2}dx $
Using the generating function for the Gegenbauer Polynomials: $ (1-2xr+r^2)^{-\lambda}=\sum_{n=0}^{\infty}C_{n}^{\lambda}(x)r^n\ \ \ ,\left|r\right|<1 $ And writing $L=l_1+l_2+l_3$, we get: $ =\int_{-1}^{1}(1-2xr+r^2)^{-(L+3)}(1-x^2)^{(L+3)/2} $ $ =\sqrt{\pi}\frac{\Gamma(\frac{L+3}{2})}{\Gamma(\frac{L+4}{2})}(1-r^2)^{-(L+3)} $ Rewriting $(1-r^2)^{-(L+3)}=(1-r^2)^{-(l_1+1)}(1-r^2)^{-(l_2+1)}(1-r^2)^{-(l_3+1)}$ and applying the binomial theorem gives: $ \begin{multline} \sqrt{\pi}\frac{\Gamma\left(\frac{L+3}{2}\right)}{\Gamma\left(\frac{L+4}{2}\right)}\left(\sum_{n=0}^{\infty}\left(-1\right)^n\left(\begin{array}{c}-l_1-1\\n\end{array}\right)r^{2n}\right)\times\\ \left(\sum_{n=0}^{\infty}\left(-1\right)^n\left(\begin{array}{c}-l_2-1\\n\end{array}\right)r^{2n}\right)\left(\sum_{n=0}^{\infty}\left(-1\right)^n\left(\begin{array}{c}-l_3-1\\n\end{array}\right)r^{2n}\right) \end{multline} $ Applying the product of series identity twice again gives: $ =\sum_{n=0}^{\infty}\sum_{l=0}^{n}\sum_{m=0}^{n-l}r^{2n}(-1)^n\left(\begin{array}{c}-l_1-1\\l\end{array}\right)\left(\begin{array}{c}-l_2-1\\m\end{array}\right)\left(\begin{array}{c}-l_3-1\\n-l-m\end{array}\right)\sqrt{\pi}\frac{\Gamma\left(\frac{L+3}{2}\right)}{\Gamma\left(\frac{L+4}{2}\right)} $
Comparing the series in $i,j,k$ to this one in $n,l,m$ I find that $i=2n$ (so that the original series is zero unless i is an even integer). Then comparing the coefficients of $r$ I get: $ \begin{multline} \sum_{j=0}^{i}\sum_{k=0}^{i-j}\int_{-1}^{1}C_{j}^{l_1+1}(x)C_{k}^{l_2+1}(x)C_{i-j-k}^{l_3+1}(x)(1-x^2)^{(l_1+l_2+l_3+1)/2}dx\\=\sum_{l=0}^{n}\sum_{m=0}^{n-l}(-1)^n\left(\begin{array}{c}-l_1-1\\l\end{array}\right)\left(\begin{array}{c}-l_2-1\\m\end{array}\right)\left(\begin{array}{c}-l_3-1\\n-l-m\end{array}\right)\sqrt{\pi}\frac{\Gamma\left(\frac{L+3}{2}\right)}{\Gamma\left(\frac{L+4}{2}\right)} \end{multline} $ I then set $ \begin{align} i=2n\\j=l=n_1-l_1\\k=m=n_2-l_2\\i-j-k=n_3-l_3 \end{align} $ Whence I get the $\bf{FINAL\ RESULT}$: $ \int_{-1}^{1}C_{n_1-l_1}^{l_1+1}(x)C_{n_2-l_2}^{l_2+1}(x)C_{n_3-l_3}^{l_3+1}(x)(1-x^2)^{(l_1+l_2+l_3+1)/2}dx= $ $ (-1)^{(N-L)/2}\left(\begin{array}{c}-l_1-1\\ n_1-l_1\end{array}\right)\left(\begin{array}{c}-l_2-1\\n_2-l_2\end{array}\right)\left(\begin{array}{c}-l_3-1\\-\frac{1}{2}((n_1-l_1)+(n_2-l_2)-(n_3-l_3))\end{array}\right)\sqrt{\pi}\frac{\Gamma\left(\frac{L+3}{2}\right)}{\Gamma\left(\frac{L+4}{2}\right)} $
So, \bf{HERE's} my question:
Why is this incorrect? I've tested it and it is not true for all values of $n_1,l_1,n_2,l_2,n_3,l_3$.
I've now spent the better part of a week trying to figure out why this doesn't work, but I can't figure it out. It seems like maybe it has something to do with the way I equate the indices, or perhaps a theorem about the convergence of the product of power series that would make something I did illegal, but I can't find anything so far that would tell me that there's anything wrong with it. In case it's helpful the answer seems to be off by a multiplicative factor that is a function of the indices (a gamma function maybe?).