# I constructed the following operator with Galois group is H_{72}^{SL_3}
# with the methods of van der Put / Ulmer in:
# http://www.msri.org/publications/preprints/online/1998-058.html
# I selected the generators g0, g1, g_infty such that g0 and g1 have the
# same characteristic polynomial (same eigenvalues), and g_infty has
# "almost" the same characteristic polynomial (same up to a shift).
# This way all the coefficients can be obtained by exploiting the
# symmetry (first interchange g0 and g_infty, and *after* that interchange
# g0 and g1) and by exploiting the condition that there should be no
# logarithms in the local solutions at the apparent singularity. Then I
# moved the singularities a little in order to get rational coefficients;
# if g0 and g1 are located at 0, 1  then the apparent singularity is
# not rational and one gets algebraic coefficients.
#
# Mark van Hoeij, November 1998.

L:=
DF^3+(-1+21*x^2-24*x)/(x-1)/(3*x^2+1)*DF^2+1/48*(4437*x^3-5973*x^2+171*x-683)/
(x-1)/(3*x^2+1)^2*DF+1/216*(1983*x^2-7297*x+13338*x^4-22647*x^3-737)/(x-1)/(3*
x^2+1)^3;


# P is a Riccati polynomial of L. How such a polynomial P can be
# computed is theorem 2.1 in "Liouvillian solutions of linear differential
# equations of order three and higher".

P:=
X^9+30*x/(3*x^2+1)*X^8+1/4*(1599*x^2+1)/(3*x^2+1)^2*X^7+1/72*(223587*x^3+15*x
^2+441*x+5)/(3*x^2+1)^3*X^6+1/1152*(17855019*x^4+4800*x^3+73962*x^2+1600*x-21)/
(3*x^2+1)^4*X^5+1/3456*(178134957*x^5+119907*x^4+1290114*x^3+40062*x^2-1071*x+
31)/(3*x^2+1)^5*X^4+1/248832*(28420622685*x^6+38312352*x^5+323581833*x^4+
12864960*x^3-524073*x^2+31392*x+67)/(3*x^2+1)^6*X^3+1/497664*(80931106503*x^7+
191139507*x^6+1350862263*x^5+64704555*x^4-3559875*x^3+330129*x^2+1349*x-111)/(3
*x^2+1)^7*X^2+1/47775744*(6449857313913*x^8+24395769216*x^7+150200784540*x^6+
8354150784*x^5-580059306*x^4+73858176*x^3+437148*x^2-72576*x+505)/(3*x^2+1)^8*X
+1/644972544*(32112055562886*x^9+182329225533*x^8+1005364927404*x^7+63398867400
*x^6-5314028328*x^5+869144634*x^4+6682932*x^3-1669824*x^2+23490*x-127)/(3*x^2+1
)^9;

# P has genus 1, and j-invariant 0.

# Now L has a solution exp(int( RootOf(P,X) ,x)). This is an algebraic function,
# which (after multiplying by a suitable constant) is a root of the following
# polynomial:

A:=
Y^27-168/(3*x^2+1)^5*Y^18+405*(1+x)/(3*x^2+1)^7*Y^15-636/(3*x^2+1)^10*Y^9-324*
(1+x)/(3*x^2+1)^12*Y^6-243/4*(1+x)^2/(3*x^2+1)^14*Y^3-8/(3*x^2+1)^15 ;

# So RootOf(A,Y) is a solution of L. The Picard-Vessiot field of L is the
# splitting field of A over C(x), where C is the field of complex numbers.
# The Galois group of L and of A is H_{72}^{SL_3}, which has 72 * 3 = 216
# elements.

# Note that the Galois group of P over C(x) is H_72, it has 72 elements, and
# can be represented as a permutation group as follows:
with(group):
H72 := permgroup(9,{[], [[2, 8, 3, 9], [4, 6, 7, 5]], [[1, 9, 3, 7], [4, 8, 6, 5]]});
grouporder(H72);
# 
# The same holds for the polynomial
collect(rem(A,Y^3-Z,Y),Z,factor); # = A with Y^3 replaced by Z.

# The polynomial A has the following Galois group:

perm1:=[[1, 27, 2, 15, 19, 4, 24, 5, 25, 7, 17, 22],
        [3, 23, 14, 13, 26, 8, 21, 12, 11, 10, 9, 18], [6, 16, 20]];
perm2:=[[1, 27, 3, 17, 19, 4, 26, 2, 25, 7, 11, 24],
        [5, 20, 14, 8, 22, 6, 21, 10, 15, 16, 9, 23], [12, 18, 13]];

H72_SL3 := permgroup(27,{perm1,perm2});
grouporder(H72_SL3); # 3*72=216 elements.