# Filename: public_html/files/FP28
# Compute operator L for FP28. Verify Riccati polynomial of transformed operator LL.
# Use Maple Release 5 (or higher).
#
# Mark van Hoeij, november 1998, MSRI Berkeley.

z:=RootOf( x^32+x^28-x^20-x^16-x^12+x^4+1 ); # zeta 120.
zz:=evalf(exp(2*Pi*I/120));

a:=matrix(4,4,[z^15,0,0,0, 0,z^15,0,0, 0,0,0,z^15, 0,0,z^15,0]);
b:=matrix(4,4,[z^30,0,0,0,
	0,-1/3*z^30,2/3*z^30,2/3*z^30,
	0,2/3*z^30,-1/3*z^30,2/3*z^30,
	0,2/3*z^30,2/3*z^30,-1/3*z^30]);
c:=matrix(4,4,[-1/4*z^30, 1/2*z^28 - 1/2*z^20 -z^16 -1/2*z^12 - 1/2*z^8 +z^4 +3/4, 0,0,
	1/2*z^28 - 1/2*z^20-z^16-1/2*z^12-1/2*z^8+z^4+3/4, 1/4*z^30 ,0,0,
	0,0,0,z^30,
	0,0,z^30,0]);
d:=linalg[diag](1,1,z^20-1,-z^20);

la:=proc(A)
	[seq([seq(evala(Normal(A[i,j])),j=1..linalg[coldim](A))],i=1..linalg[rowdim](A))]
end:

groep:=proc(gens)
	W:={};
do
	W_old:=nops(W);
	W:={seq(seq(la(linalg[multiply](i,j)),i=gens),j=W),op(W),op(map(la,gens))};
	if nops(W)=W_old then break fi
od;
	W
end:

# nops(groep({a,b,c})); # 48 elementen.
# nops(groep({a,b,c,d})); # 480 elementen.
# groep({a, evalm(b &* d)}); # 96 elementen.
# groep({evalm( a &* c), evalm(b &* d)}); # 96 elementen.
# nops( groep({a, evalm(b &* c &* d )}) ); 
# groep({a, evalm(b &* d &* c)}); # 480 elementen.
# nops( groep({a, evalm(b &* c &* d )}) ); # 480 elementen.


A:=evalm(a);
B:=map(evala, evalm(b &* c &* d ));
C:=map(evala, evalm( A^(-1) &* B^(-1)));

CA=evala(Factor(convert(linalg[charpoly](  A ,x),RootOf)));
CB=evala(Sqrfree(convert(linalg[charpoly]( B ,x),RootOf))); # -> determine exponents
CC=evala(Sqrfree(convert(linalg[charpoly]( C ,x),RootOf))); #    from these charpoly's.

# This implements the formula for the multiplicity that Marius van der Put
# gave me.
multipliciteit:=proc(A,B)
	global geconjugeerde,z,zz;
	local cd,som,a,cp,n,x,i,k;
	cd:=linalg[coldim](A);
	som:=-cd;
	for a in [A,B,map(evala,evalm( A^(-1) &* B^(-1)))] do
		cp:=evala(Expand( convert(linalg[charpoly](a,x),RootOf) ));
		n:=1;
		while evala(Rem( (x^n-1)^cd , cp,x))<>0 do n:=n+1 od;
		lprint(n);
		som:=som+
		add(i/n * 1/n * add(exp(2.0*Pi*I*(-i*k)/n)*
		 subs(`if`(geconjugeerde=true,{RootOf(_Z^2+1)=I},{I=-I,RootOf(_Z^2+1)=-I})
		  ,subs(z=zz,linalg[trace](`if`(k=0,linalg[diag](1$cd),a^k))) ),k=0..n-1),i=0..n-1)
	od;
	evalf(som)
end:

# multipliciteit(A,B); # multiplicity = 3, so I can add 2 to the exponents.



# This is the generic diff eq with 3 regular singularities of order 4.
equa:=DF^4+(a1*x+a0)/x/(x-1)*DF^3+(b2*x^2+b1*x+b0)/x^2/(x-1)^2*DF^2
+(c3*x^3+c2*x^2+c1*x+c0)/x^3/(x-1)^3*DF
+(d4*x^4+d3*x^3+d2*x^2+d1*x+d0)/x^4/(x-1)^4;


with(DEtools);
_Envdiffopdomain:=[DF,x];

singularities:=[0,1,infty];
# Want the following exponents:
# e[p,i], p in singularities, i in {1,2,3,4}

fuchs_relation:=0; # -1 * sum of exponents e[p,i]
eqns:={};
for p in singularities do
indeq[p]:=subs(_Z=x,op(gen_exp(equa,T,x=subs(infty=infinity,p))[1][1]));
indeq[p]:=collect(indeq[p]/lcoeff(indeq[p],x),x,normal);
fuchs_relation:=normal(fuchs_relation + coeff(indeq[p],x,3));
eqns:=eqns union {coeffs(collect(mul(x-e[p,i],i=1..4)-indeq[p],x),x)}
od;

vars:=indets(equa) minus indets([_Envdiffopdomain,singularities]);
eqns:=map(numer@normal,subs(e[infty,1]=solve(
	fuchs_relation+add(add(e[p,i],i=1..4),p=singularities)
	,e[infty,1]),eqns));
sol:=readlib(`solve/linear`)(eqns,vars);
L:=collect(subs(sol,equa),DF,normal);

L:=
subs(e[0,1]=1/8-1,e[0,2]=1/8,e[0,3]=1/8+1 ,e[0,4]=5/8+1,
     e[1,1]=1/10 ,  e[1,2]=3/10, e[1,3]=7/10,e[1,4]=9/10,
     e[infty,1]=1/8, e[infty,2]=3/8,e[infty,3]=5/8,e[infty,4]=7/8,
	L);

v:=seq(`if`(has(i,ln(T)),i,NULL),i=formal_sol(L,T,x=0,terms=3));
SS:={ seq(subs(T=1,coeff( convert(numer(v[1]),polynom),ln(T),i)),i=1..2),
	   seq(subs(T=1,coeff( convert(numer(v[2]),polynom),ln(T),i)),i=1..2)
};
S:={seq(`if`(degree(i,indets(S))=1,i,NULL),i=SS)};
S:=solve(S);
L:=collect(subs(S,L),DF,normal);

v:=seq(`if`(has(i,ln(T)),i,NULL),i=formal_sol(L,T,x=0,terms=3));
S:=SS union { subs(T=5,coeff( convert(numer(v[1]),polynom),ln(T))) };
S:=solve(S);
L:=collect(subs(S,L),DF,normal);
v:=seq(`if`(has(i,ln(T)),i,NULL),i=formal_sol(L,T,x=0,terms=3));
# OK, logs are gone.


LL:=symmetric_product(L,DF + (1/8)/x + (1/10)/(x-1) );
 # should have invariant of degree 5, can determine that from A,B,C.

if false then # skip computation of semi-invariants of the group.
read Matrices; # load code for Sympower's and Hom's of representations.

N:=5; # search semi-invariant of L, corresponding to invariant of LL.

aa:=matrix(1,1,[ -z^15 ]); # search semi-invariant, if you make all
bb:=matrix(1,1,[ z^30 ]);  # these matrices [1] then you're looking
cc:=matrix(1,1,[ z^30 ]);  # for invariants, and you'll find a 5-dim
dd:=matrix(1,1,[ 1 ]);     # space of invariants of degree 8, and
                           # nothing of lower degree.


aN:=map(evala,Sympower(a,N));
bN:=map(evala,Sympower(b,N));
cN:=map(evala,Sympower(c,N));
dN:=map(evala,Sympower(d,N));

S:=map(vect2matrix,linalg[nullspace](
 linalg[blockmatrix](4,1,[ Hom(aa,aN), Hom(bb,bN) ,Hom(cc,cN), Hom(dd,dN) ])
),linalg[rowdim](aa),linalg[rowdim](aN) );
lprint(N,nops(S));
semi_invariants[N]:=map(op,map(columns2poly,S,4,N,X));

# For N=1..10 find only degree 8 invariants, 5 of them.
# Do they factor???
# 2 semi-invariants of degree 5. Take the one that has X[1] as factor.
# That should factor into linear factors.

fi;

# Here is LL again:
LL1:=
DF^4+1/10*(89*x-45)/x/(x-1)*DF^3+1/400*(7479*x^2-6895*x+600)/x^2/(x-1)^2*DF^2+
1/4000/x^2*(32781*x^2-40825*x+8300)/(x-1)^3*DF+1/20000*(3927*x-2725)/x^2/(x-1)
^3;
if normal(LL1-LL)<>0 then ERROR() fi;



# Now with the code of Jacques-Arthur and me I computed the following
# Riccati polynomial for LL (not for L, then you have to apply
# a substitution on P).

P:=
Z^5+(245366784*x^3+739444563*x^2+737492600*x-2787596875)/(-557519375+
184373150*x+246481521*x^2+122683392*x^3)/x*Z^4+2/5*(483065856*x^4+1875325881*x^
3+1373204105*x^2-16703581625*x+13937984375)/(-557519375+184373150*x+246481521*x
^2+122683392*x^3)/x^2/(x-1)*Z^3+2/25*(939294720*x^5+5380106667*x^4+5809988990*x
^3-91889491800*x^2+148598501250*x-69689921875)/(-557519375+184373150*x+
246481521*x^2+122683392*x^3)/x^3/(x-1)^2*Z^2+1/125*(1808142336*x^6+15408360513*
x^5+43238642915*x^4-522670373950*x^3+1182921498750*x^2-1068395471875*x+
348449609375)/(-557519375+184373150*x+246481521*x^2+122683392*x^3)/x^4/(x-1)^3*
Z+1/3125*(3451908096*x^7+41967774147*x^6+302211588068*x^5-3042266659625*x^4+
8490891950000*x^3-11023761109375*x^2+6968992187500*x-1742248046875)/(-557519375
+184373150*x+246481521*x^2+122683392*x^3)/x^5/(x-1)^4 ;

# Verifying:
#
# rightdivision(LL,DF-RootOf(P,Z)); # takes too long.
#
# A[0]:=u(x);
# for i from 1 to 4 do
# lprint(i);
# A[i]:=evala(subs(diff(u(x),x)=u(x)*RootOf(P,Z),  diff(A[i-1],x) )); # takes too long.
# od;

# Faster verification:
N:=10; # check up to order 10.
with(algcurves);
p:=puiseux(numer(P),x=2,Z,N)[1];
Normalizer:=evala;
p:=series(p,x=2,N);
p:=series(exp(int(p,x)),x=2,N);
series( add( coeff(LL,DF,i)*diff(p,x$i),i=1..4)+ coeff(LL,DF,0)*p ,x=2,N);
evala(Expand( subs(x=x+2,convert(%,polynom)) ));
if % <> 0 then ERROR(NOT_OK) else lprint(OK) fi;


# We had LL = symmetric_product(L,DF + (1/8)/x + (1/10)/(x-1) );
# so the Riccati polynomial for L should be:

Q:=collect( subs(Z=Z - ((1/8)/x + (1/10)/(x-1)) ,P) ,X,factor);
# Verify that this is a Riccati polynomial for the operator L.
N:=10; # check up to order 10.
p:=puiseux(numer(Q),x=2,Z,N)[1];
Normalizer:=evala;
p:=series(p,x=2,N);
p:=series(exp(int(p,x)),x=2,N);
series( add( coeff(L,DF,i)*diff(p,x$i),i=1..4)+ coeff(L,DF,0)*p ,x=2,N);
evala(Expand( subs(x=x+2,convert(%,polynom)) ));
if % <> 0 then ERROR(NOT_OK) else lprint(OK) fi;

# Summary:
" The differential equation ", diffop2de(L,y(x)),"has Galois group FP28 and has",
exp(Int(RootOf(Q,Z),x)),"as a Liouvillian (in fact: algebraic) solution.";