# Same as in the previous file:
L_plus_x := Dx^2 + ( 1/x - 1/(x + (t+1)*(t-1)^3) + 2*(x+t^4+2*t^2-1)/(x^2 + (2*t^4+4*t^2-2)*x + (t^2-1)^4) )*Dx
  + ( 4*x^2 + (8*t^4-16*t^3-4*t^2+20*t-9)*x+(4*t^4-8*t^3+12*t^2+4*t-5)*(t+1)*(t-1)^3 )
  / ( 16 * x * (x^2+(2*t^4+4*t^2-2)*x+(t^2-1)^4) * (x+(t+1)*(t-1)^3) ) ;
L_minus_x := subs(t = -t, L_plus_x);



with(DEtools): _Envdiffopdomain := [Dx,x];
chvar := proc(L, y, Dy, sb1,sb2)
	local i, f, x, Dx, s;
	Dx, x := op(_Envdiffopdomain);
	f := add( subs(sb1,sb2,coeff(L, Dy, i)) * mult( Dx/diff(subs(sb1,sb2,y),x) $ i), i = 0 .. degree(L, Dy));
	s := proc(a) factor(simplify(normal(a,expanded),symbolic)) end:
	sort(collect(f/lcoeff(f, Dx), Dx, s), Dx);
end:

m := 1/(64*x) - (t+1)*(t-1)^3;
Lm := chvar(L_plus_x, x,Dx, x = m, {});


gfun[diffeqtorec](diffop2de(symmetric_product(Lm,Dx+1/2/x),y(x)),y(x),u(n));
subs(n=n-2,%):
collect(%,indets(%,function),factor);

# Typed from paper:
(n+1)^2 * u(n+1)
-2^2 * (16 * (t^4 - 6*t^3 - 4*t^2 + 6*t - 1)*(n^2+n) + 4*t^4 - 24*t^3 - 12*t^2 + 20*t - 3)*u(n)
-2^11 * t * (t-1)^3 * (t+1) * ( 8*(t^2 + 2*t - 1)*n^2 - 2*t^2 - 6*t + 3) * u(n-1)
+2^18 * t^2 * (t-1)^6 * (t+1)^2 * (2*n+1) * (2*n-3)*u(n-2);

# Check that the computed recurrence and the one in the paper are the same:
ShouldBeZero := normal(%% - %);


# Typed from paper:
Sol_L2 := (Rminus^(1/8) + Rplus^(1/8)) / (2 * sqrt(S));
RSA :=
	Rplus = 1/Rminus,
	Rminus  = A - sqrt(A^2-1),
	S = (1-16*(t+1)*(t-1)^3 * x) * (1 + 2^6 * (t^3+t^2-t) * x - 2^10*(t^3+t^2)*(t-1)^3 * x^2),
	A = 1 + 2^7*(2*t-1)^2 * x - 2^11*(t-1)^3*(2*t-1)*(2*t^2+5*t-1) * x^2
	      + 2^17 * t*(t-1)^6*(2*t^2 + 2*t - 1) * x^3 - 2^21 * (t^3+t^2)*(t-1)^9 * x^4;


sort(collect(convert(eval(series(subs(RSA,Sol_L2),x=0,6), csgn = proc() 1 end),polynom),x,factor),x,ascending);
sort(collect(convert(formal_sol(symmetric_product(Lm,Dx+1/2/x),T,x=0)[1,2],polynom),T,factor),T,ascending);
ShouldBeZero := seq(normal(coeff(%,T,i) - binomial(2*i,i) * coeff(%%,x,i)),i=0..5);

Relation_Zx := ( ((4*S*Z^2 - 2)^2 - 2)^2 - 2 - 2*A ) / (1-16*(t+1)*(t-1)^3 * x);
LongForm_Zx := normal(subs(RSA, %)):

# From computation, and compared with paper:
param := 
	zeta = 1/16/t*(t^2-v^2)/(t*v^4+(t+1)*(t-1)^2*(t^2-2*v^2-t)),
	Zsquare = 1/2/t*(t*v^4+(t+1)*(t-1)^2*(t^2-2*v^2-t))^2*(t-1+v)^2*(v-t)/v^2/(t^2+t*v-1)/(v^4-2*t^2*v^2+(t^2-1)^2)^2;

# Verify parametrization for some random parameter-values:
subs( t = 99, v = 76, subs(param, subs(Z = sqrt(Zsquare), x=zeta, LongForm_Zx)));

IntegrandSquare := factor( subs(param, Zsquare/(zeta * (4*x-zeta))) * diff(subs(param,zeta),v)^2 );
TypedFromPaper := (-2) * (v+t-1)^2 / ( (v+t)*(v+t-t^(-1))*(v^2-t^2 + 64*t*x*( t*v^4 + (t+1)*(t-1)^2*(t^2-2*v^2-t))) );

# Check that the computed IntegrandSquare is the same as the one from the paper:
ShouldBeZero := normal(%% - %);

###############  Transformations between this Integrand and the one in Theorem 1 ##############

m := 1/(64*x) - (t+1)*(t-1)^3;
m_inv := subs(X=x,solve(m=X,x));
IntegrandSquare := factor(subs(x = m_inv, IntegrandSquare));

# After rescaling x and a Mobius transformation w.r.t. we get
	IS := subs(x = -x/(4*u^2), u = 1/(1-2*t^2), IntegrandSquare);
	# Transformation in equation (23), start of Section 3.
	mv := 2*v*t/(2*t^2+v-1)-t;
	TransformedIntegrandSquare := factor(subs(v = mv, IS) * diff(mv,v)^2);

# Square of integrand in Theorem 1 is:
	Boven := 1 - u*v + v * 2 * t * u;   # From Theorem 1, replaced (2*u^2 - 2*u)^(1/2) = 2*t*u
	Onder := v * (1-v) * ( (1-v)*(1-u^2*v)*(1+u*v)^2 + x*v*(1-u*v)^2 );
	IntegrandSquareTheorem1 := factor(subs(u = 1/(1-2*t^2), Boven^2/Onder));

Ratio := factor(IntegrandSquareTheorem1 / TransformedIntegrandSquare);
indets(%);
# This Ratio is independent of v, and can thus be moved to the algebraic
# prefactor in front of the integral.