read(diffop);

# The following is the Hurwitz equation:

deq:=x^2*(x-1)^2*diff(y(x),x$3)+(7*x-4)*x*(x-1)*diff(y(x),x$2)+
(72/7*(x^2-x)-20/9*(x-1)+3/4*x)*diff(y(x),x)+(792/7^3*(x-1)+5/8+2/63)*y(x)=0;

f:=convert(deq,diffop);
sy4:=symmetric_power(f,4);

# The order of sy4 is 14. According to the theory (see: Singer and Ulmer,
# "Galois groups for second and third order linear differential equations",
# J. Symb. Comput. 16 (1993), pp. 9-36)
# this sy4 must have a factor of order 6 and a factor of order 8.

v:=factor_op(convert(sy4,diffop,x=infinity),`all right factors`,[]);

# There are 7 different exponential parts. This means that from these 7,
# at least one of them is not going to appear in the factor of order 6.
# It is a matter of trial and error to figure out which one is the good one.
# So we have to do 7 different computations, one computation for each of these
# local factors in the list v. To save some computation time, I'll take only
# one element from v, for which I know the factorization works.

for i in v do if nt(i,1)=xDF-6 then w:=i fi od;

`diffop/modulus2`:=nextprime(10000): # This is a global variable that needs
# to be set. Normally this is done by the procedure DFactor, but in this
# example we don't use DFactor because that takes too much time.

infolevel[diffop]:=10;

# Now the following command tries to search for a global factor of order 8
# for which the local factor w is a right hand factor. The first 8 means:
# search for order <= 8, and the last 8 (which is optional) means: search for
# order >= 8.
# The singularity we choose is infinity, these numbers 48 and 27 are bounds on
# the degrees of the global factor. The [] means that the ground field is Q.

# What will happen is the following: via modular arithmetic and the
# Pade algorithm (see the papers of Beckermann Labahn, and/or Derksen)
# a global factor of sy4 will be constructed modulo some power of a prime
# number. This prime power is being increased until the algorithm is able to
# construct a global factor in characteristic 0 (yes, I know, I should use
# chinese remaindering instead of just increasing the modulus.)

t:=try_factorization(w,8,[48$14],infinity,sy4,27,[],8);

# Now test if it is a factorization:
normal(sy4 - mult(t[1],t[2]));
# OK