M0 := Matrix(2,2, [ [1,1], [0,1]] );
M1 := Matrix(2,2, [ [1,0], [ alpha-2,1]] );
M2 := Matrix(2,2, [ [1,0], [-alpha-2,1]] );
M3 := (M0.M1.M2)^(-1);

# Input: V = [M1,M2,M3,M4]
# Output is simultaneously conjugated to V, and if two inputs are simultaneously
# conjugated, then their outputs are the same.
NormalizeV := proc(V) local M1,M2,v1,v2,a,P,i;
	M1 := Matrix(2,2,V[1]);
	v2 := Vector([1,0]);
	if convert(M1.v2, listlist) = [1,0] then
		v2 := Vector([0,1])
	fi;
	v1 := evala(M1.v2 - v2);
	for i from 2 to 4 do
		M2 := Matrix(2,2,V[i]);
		if evala(convert(M2.v1 - v1,listlist)) <> [0,0] then
			P := v2 + a * v1;
			P := solve({op(convert(M2 . P - P, listlist))}, {a});
			break
		fi
	od;
	P := <v1 | v2 + eval(a,P) * v1>;
	[seq( evala(convert(P^(-1) . Matrix(2,2,a) . P, listlist)), a=V)]
end:


# Braid entries i and i+1
Braid := proc(V, i) local M1,M2,M1new,s;
	M1 := Matrix(2,2,V[i]);
	M2 := Matrix(2,2,V[i+1]);
	s := 1;
	if i=3 then		# ConjugacyClass(M3) = -1 * ConjucacyClass(M4)
		s := -1		# so if we braid 3,4 then multiply by -1
	fi;
	M1new := convert(evala(M1 . M2 . M1^(-1)), listlist);
	NormalizeV(subsop(i=s*M1new, i+1 = s*V[i], V))
end:

BraidOrbit := proc(V) local W,S;
	W := {NormalizeV(V)};
	do
		S := map(Braid,W,1) union map(Braid,W,2) union map(Braid,W,3);
		if S minus W = {} then return [op(W)] fi;
		W := W union S
	od
end:

PureBraid := proc(V, i)
	Braid(Braid(V,i),i)
end:

PureBraidOrbit := proc(V) local W,S;
	W := {NormalizeV(V)};
	do
		S := map(PureBraid,W,1) union map(PureBraid,W,2);
		if S minus W = {} then return [op(W)] fi;
		W := W union S
	od
end:
	
O_zeta := proc(n) global alpha,M0,M1,M2,M3; local x, R, V;
	options remember;
	R := RootOf(NumberTheory:-CyclotomicPolynomial(n, x));
	V := subs(alpha = RootOf(sqrfree(evala(Norm(x - (R + 1/R))))[2,1,1]),
	  map(convert,[M0,M1,M2,M3],listlist));
	PureBraidOrbit(V)
end:
	# Note: the PureBraidOrbit is smaller than the BraidOrbit, but for the
	# differential equation the difference becomes irrelevant if we also
	# allow Mobius transformations to move the 4 singularities.

MatV := proc(V)
	map(Matrix,V)
end:


for V in O_zeta(8) do print(map(Matrix,convert(V,radical))) od;

# To get the full BraidOrbit, pick one of those V's and apply BraidOrbit(V).
	


FindIt := proc(V,O) local i;
	for i do if O[i]=V then return i fi od
end:

# Compute the "dessin d'enfant" (given by two permutations g0, g1) for a PureBraidOrbit
OrbitToDessin := proc(O) local V,i;
	seq(convert([seq(FindIt(PureBraid(V,i),O),V=O)],disjcyc), i=1..2)
end:

OrbitGenus := proc(O) local g0,g1,ginf,i;
	g0, g1 := OrbitToDessin(O);
	ginf := group[mulperms](g0, g1);
	1-nops(O) + add(nops(i)-1, i = map(op,[g0,g1,ginf]))/2
end:

seq([n,OrbitGenus(O_zeta(n))], n = [5, 8, 12]);

# By the "genus" of a differential equation, we mean the smallest genus of a field k
# such that we can Mobius-transform the differential equation to an element of k(x)[d/dx]
# in such a way that all singularities are in k.

# For n = 8, the genus of the orbit is 0.  So there must exist a genus-0 field k
# such that the differential equation can be reduced to k(x)[d/dx], and with all singularities
# in k.  And that is indeed the case, here with k = Q(p)

X8_rationalized_singularities := Dx^2+(16*p^4/(16*p^4*x-(p^2-1)^4)-4*p^3/(4*p^3*x+(p+1)*(p-1)^3)+1/x+1/(x-1))*Dx
+1/64*(-1024*p^7*x^2-4*p^3*(p^8-4*p^6+144*p^5-258*p^4-32*p^3+332*p^2-176*p+9)*x+(p+1)*(7*p^8+16*p^7-12*p^6-32*p^5+
226*p^4-48*p^3-60*p^2+15)*(p-1)^3)/(-16*p^4*x+(p^2-1)^4)/(x-1)/(4*p^3*x+(p+1)*(p-1)^3)/x ;

# The logarithmic singularities of this equation are located at 0, 1, infinity, and
x4 := (p^2-1)^4/(16*p^4);
# which is a Belyi map whose monodromy-matrices g0,g1 are:

g0,g1 := OrbitToDessin(O_zeta(8));

# Indeed, the 2 4-cycles in g0 correspond to the 2 roots of x4 of order 4,
# while the 2 2-cycles in g1 correspond to the 2 roots of 1-x4 of order 2.
factor(1-x4);


# For n = 5 the genus is 1.  This means that if we want all singularities to be k-rational, with genus(k) minimal,
# then genus(k) will be 1.

L_H5 := Dx^2+(2*x^3+(u^2+6*u+1)*x^2-16*u^2)/(x^3+(u^2+6*u+1)*x^2+8*u*(u+1)*x+16*u^2)/x*Dx
        +1/100*(25*x^2-(u+3)^2*x-40*u*(u+3))/(x^3+(u^2+6*u+1)*x^2+8*u*(u+1)*x+16*u^2)/x;

# Substituting  u = (p+1)/(p^3-p^2) makes one more singularity rational, which we then scale to x = 1
# (another logarithmic singularity is at x = infinity, and the apparant singularity is at x = 0)

dn := p^4*x^2+(p+1)*(p^3-p^2+3*p+1)*x+4*(p+1)^2;
L_H5_p := Dx^2+((2*p^4*x+(p+1)*(p^3-p^2+3*p+1))/dn+(-1+p)^2/((-1+p)^2*x+4)-1/x)*Dx
    +1/100*(25*(-1+p)^2*p^4*x^2-(3*p^3-3*p^2+p+1)^2*x-40*(p+1)*(3*p^3-3*p^2+p+1))/x/((-1+p)^2*x+4)/dn;

# The remaining singularity, x = RootOf(dn, x) requires an extension of Q(p) of genus:
algcurves[genus](dn, x, p);
# as expected.


# OrbitGenus(O_zeta(12)) is also 1, so we expect the same for its (not yet computed) differential equation.