Sum(LegendreP(n,y) * z^n, n=0..infinity) = (1 - 2*y*z + z^2)^(-1/2);
# Check to order 10:
ShouldBeZero := simplify(series(lhs(%)-rhs(%), z=0, 10));

LegendreP(n,y) = hypergeom( [-n, n+1], [1], (1-y)/2 );
# Check for random n, say n = 11:
ShouldBeZero := simplify(eval(lhs(%)-rhs(%), n=11   ));

f := 2*(1-y^2)*z / (1-2*y^2*z + z^2):
Sum(LegendreP(n,y)^2 * z^n, n = 0..infinity)
  = (1 - 2*y^2*z + z^2)^(-1/2) * hypergeom( [1/4, 3/4], [1], f^2);
# Check to order 10:
ShouldBeZero := simplify( series(lhs(%)-rhs(%), z=0, 10) );

# Variation on previous formula:
f := normal(2*f/(f+1)):
Sum(LegendreP(n,y)^2 * z^n, n = 0..infinity)
  = (1 + (2-4*y^2)*z + z^2)^(-1/2) * hypergeom([1/2, 1/2],[1], f);
ShouldBeZero := simplify( series(lhs(%)-rhs(%), z=0, 10) );

########### Radius of convergence from Remark 1.

R_y := 1/4 * min( abs( y + (y^2-1)^(1/2) ), abs( y-(y^2-1)^(1/2) ) )^2;

#       R_y > 0 for any y in C
#       R_y = 1/4 if y is a real number in the interval [-1, 1]
#       R_y < 1/4 for any y not in that interval
#
# The series F(y,z) converges for any z in C with abs(z) < R_y.


########### Theorem 1, floating point check

F_numeric := proc(y0::complex(numeric), z0::complex(numeric))
	local R, S,u,n;
	global R_y, y, z;
	R := evalf(subs(y = y0, R_y),5);
	lprint("Radius of convergence |z| < ", R);
	lprint("|z|/Radius = ", evalf(abs(z0)/R, 3));
	S, n, u := 0, 0, 1;
	while abs(u) > 10^(-12) do
		# n'th term of series F(y,z)
		u := binomial(2*n,n) * LegendreP(n, y0)^2 * z0^n;
		S := S + u;
		n := n+1;
		if n > 1000 then error "Giving up" fi
	od;
	lprint(n, "terms used");
	S
end:

I_pm := proc(u, x, teken) local Boven, Onder, f, v;
	Digits := 40;
	Boven := 1 - u*v + teken * v * (2*u^2 - 2*u)^(1/2);
	Onder := v * (1-v) * ( (1-v)*(1-u^2*v)*(1+u*v)^2 + x*v*(1-u*v)^2 );
	# f := Boven/sqrt(Onder);
	# Same f but easier on the floating point integrator:
	f := evalf(Boven) / sqrt( v - v^2 ) / sqrt( evalf(collect( (1-v)*(1-u^2*v)*(1+u*v)^2 + x*v*(1-u*v)^2, v)) );
	1/Pi * evalf(Int(f, v=0..1), 10);
end:

w := sqrt( (1 + 4*z)^2 - 16 * y^2 * z) + 4 * y * sqrt(-z) ;

# Note 1: The expression inside the first square root in w is positive for
#         any z in (-1/4, R_y) and vanishes at R_y if y is a real number.
# Note 2: Since (2*u^2 - 2*u)^(1/2) in I_pm is multiplied by "teken" = +/- 1
#         the sign of this square root is always OK.
# Note 3: If z > 0 then the computation will involve complex numbers even
#         when y and z are both real.


# Now check Theorem 1 for a few values of y, z with abs(z) < R_y

y := 0.71; z := -0.12;
F_numeric(y, z)  =  w * I_pm(4*z, w^2, 1) * I_pm(4*z, w^2, -1);
diff1 := abs(lhs(%) - rhs(%)):

y := 0.55; z := 0.11;
F_numeric(y, z)  =  w * I_pm(4*z, w^2, 1) * I_pm(4*z, w^2, -1);
diff2 := abs(lhs(%) - rhs(%)):

y := 0.0001 * I - 0.4;  z := y^2 * (0.3 - 0.0002 * I);
F_numeric(y, z)  =  w * I_pm(4*z, w^2, 1) * I_pm(4*z, w^2, -1);
diff3 := abs(lhs(%) - rhs(%)):
lprint("Largest abs(difference) is ", max(diff1, diff2, diff3));