# From Möller_Zagier-CompMath2016.pdf page 2299.

A := t*(t-1)*(t-lambda)*(t-1/lambda);
B := 3/4 * (3*t^2-(beta+gamm)*t+gamm);
lambda := (31 - 7*sqrt(17))/2;
beta := lambda + 1/lambda + 1;
gamm := (27 - 5*sqrt(17))/4;

with(DEtools):
_Envdiffopdomain := [Dt,t];
L_MZ := mult(Dt, A, Dt) + B;
L_MZ := collect(L_MZ/A, Dt, factor);

# Next example on page 2299 has more singularities:
Atilde := A / (t^2 + (137-95*sqrt(17))/128 * t + 1);
Btilde := ( 1/4 * t^4 + (1113 - 399*sqrt(17))/512 * t^3 - (260375 - 69633*sqrt(17))/16384 * t^2
 - (1387 - 301*sqrt(17))/128 * t + (23-5*sqrt(17))/8 )  /  (t^2 + (137-95*sqrt(17))/128 * t + 1)^2;
Ltilde := mult(Dt, Atilde, Dt) + Btilde;
Ltilde := collect(Ltilde/Atilde, Dt, factor);


convert(formal_sol(L_MZ, T, t = 0, terms = 20)[1,2],polynom);
ifactor(denom(%)); # Arithmetic
convert(formal_sol(Ltilde, T, t = 0, terms = 20)[1,2],polynom);
ifactor(denom(%)); # Arithmetic


hypergeometricsols(L_MZ);
hypergeometricsols(Ltilde);
# Contains solid code for <= 5 true singularities, so we can conclude there are no 2F1 sol.

L_Heun := Dt^2+(434*t*17^(1/2)-768*t^2-217*17^(1/2)-1278*t+895)/(-256*t-895+217*17^(1/2))/t/(t-1)*Dt+60*(-4*t-3+17^(1/2))/(-256*t-895+217*17^(1/2))/t/(t-1) ;
# Local exponents: [0,0] at the 3 finite points, and [3/4, 5/4] at infinty.


# Verify that MZ is pullback of L_Heun
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:
symmetric_product(chvar(L_Heun, t, Dt, t = -1/128*(-64*t+7*17^(1/2)+31)*(2*t-31+7*17^(1/2))/(t+1)^2, {}), Dt + (3/2)/(t+1));
ShouldBeZero := evala(% - L_MZ);

# The singularities of L_Heun are:  0, 1, -895/256+217/256*17^(1/2), infinity
# Since -895/256+217/256*17^(1/2) = -0.001117497 is so close to 0, we replaced 17^(1/2) by
# its Galois conjugate -17^(1/2) before computing the floating point monodromy.
# We applied a floating-point version of NormalizeV from Section 4.2 to that, and then we
# used algebraic-number reconstruction to find the matrices in Section 5.1.

V := [[[1, 1], [0, 1]], [[1, 0], [-7/2-1/2*17^(1/2), 1]], [[-3/2-1/2*17^(1/2), 9/8+1/8*17^(1/2)], [-13/2-3/2*17^(1/2), 7/2+1/2*17^(1/2)]]];
map(Matrix,V);

# If one applies Braid(V,1) n times to those matrices, then the result is that
# the third matrix V[3] is conjugated by P^n for some matrix P.
# In particular, entries of that matrix satisfy a simple recurrence, with constant coefficients.
# That makes it easy to prove that the sequence is not periodic (if it was, then all roots
# of the polynomial belonging to the minimal constant recurrence should be roots of unity, and
# this is not the case).
# Therefore, PureBraidOrbit(V) is infinite, and indeed, if we call that code on this V, it runs forever.