################### Routines to prove/check the formulas from
## our ISSAC'05 paper
## M. van Hoeij and J.A Weil
###################
# initialization
restart:  
with(DEtools):
_Envdiffopdomain:=[Dx,x]:
 with(diffalg):
pullback := proc(L, a)
 local i, f;
 global x, Dx;
     f := add(
         mult(subs(x = a, coeff(L, Dx, i)), Dx/diff(a, x) $ i)
         , i = 0 .. degree(L, Dx));
     sort(collect(f/lcoeff(f, Dx), Dx, X->factor(simplify(X,symbolic))), Dx)
end;
singularites:=proc(L)
local co, liste;
    co := lcoeff(collect(numer(normal(L)), Dx), Dx);
    liste := select(X -> has(X, x), factors(co)[2]);
    [seq(RootOf(i[1], x), i = liste), infinity]
end proc:
localex:=proc(L) seq(gen_exp(L,_t,x=p),p=singularites(L)) end proc:
################################################################################################
################### Primitive groups  ###################
## former standard equations with lowest semi-invariant with value 1 (section 4.3)
s:=(6*lambda-1)*(6*lambda+1)/144;
St:=pullback( Dx^2+(7*x-4)/6/x/(x-1)*Dx-s/x/(x-1), 1/(x+1) );
St_A4:=subs(lambda=1/3,St);
St_S4:=subs(lambda=1/4,St);
St_A5:=subs(lambda=1/5,St);
##### former standard equations:
A4S:=St_A4:
S4S:=St_S4:
A5S:=St_A5:
###################################
## standard equations with lowest invariant with value 1 (section 5)
b6:=1/2/((x+1)*x):
L:=symmetric_product(St_A4,Dx+b6/6); 
## better choice 1 (gives nicer formula).
A4I:=pullback(L,-2*x/(x+1));
verif:=coeff(symmetric_power(A4I,6),Dx,0);
## other choice 2: normalise the invariant of degree 8 to 1 so that
## the standard equation is a simple pullback of S4I. Gives nicer formula
b8:=-1/3/(x+1);
L:=symmetric_product(St_A4,Dx+b8/8); 
## better choice (makes it a pullback of S4I)
A4I8:=pullback(L,-2/(x+1));
verif:=coeff(symmetric_power(A4I8,8),Dx,0);
###
b8:=-1/24/(x+1);
L:=symmetric_product(St_S4,Dx+b8);
## better choice
S4I:=pullback(L,x/(-x+1));
    ## Dx^2+1/4*(5*x-2)/(-1+x)/x*Dx-7/576/(-1+x)^2/x

#verif:=normal(A4I8-pullback(S4I,x^2));

###### pullback formula for rational S4 (section 5.2.2)
L:=S4I; 
PL:=pullback(L,F(x)):
A0,A1:=coeff(PL,Dx,0),coeff(PL,Dx,1);
sys:=map(numer,{A1-a1(x),A0-a0(x)});
R := differential_ring(ranking=[F,a1,a0], derivations=[x],  notation=diff):
GE:= Rosenfeld_Groebner(sys, R);
EQ:=rewrite_rules(GE[1]);
## simplification: may set G=2a1+diff(a0)/a0 to simplify the formulas.
the_f:=normal(subs(diff(a0(x),x)=G(x)*a0(x)-2*a1(x)*a0(x) , rhs(rewrite_rules(GE[1])[1]) )):
the_f:=collect(the_f,G(x));
## result: -1-144/7*a0(x)^3/G(x)^2

#### same process for rational A4 (section 5.2.3)
L:=A4I; 
PL:=pullback(L,F(x)):
A0,A1:=coeff(PL,Dx,0),coeff(PL,Dx,1);
sys:=map(numer,{A1-a1(x),A0-a0(x)});
R := differential_ring(ranking=[F,a1,a0], derivations=[x],  notation=diff):
GE:= Rosenfeld_Groebner(sys, R);
rewrite_rules(GE[1]);
## simplification: may set G=2a1+diff(a0)/a0 to simplify the formulas.
the_pol:=collect(normal(subs(diff(a0(x),x)=G(x)*a0(x)-2*a1(x)*a0(x) ,    
		rhs(rewrite_rules(GE[1])[1]) )),F(x),X->collect(X,G(x)) );
PF:=subs(F(x)=F,F(x)^2-the_pol);
a0(x)^3=-5/64*G(x)^2+1/64*Delta(x)^2;
PF2:=subs(a0(x)^3=-5/64*G(x)^2+1/64*Delta(x)^2,PF); 
the_f2:=simplify(solve(PF2,F)[1],symbolic);
## result is 1/2*(-Delta(x)+5^(1/2)*G(x))/Delta(x)

## Alternative approach for rational A4:
## compute the invariant of degree 8, normalise to 1
## Use the S4 formula to compute the pullback to S4I and then to A4I8: sqre root.

#### From A4 to D2 (section 5.2.4)
Dn:=Dx^2+x*Dx/((x-1)*(x+1))-1/4/((x-1)*n^2*(x+1)):
Q:=eval(Dn,n=2);
localex(A4S); ## check local exponents to see which transform to apply
b1:=1/12/(x+1)-1/12/x;
L1:=symmetric_product(A4S,Dx-b1);
localex(L1); ## check local exponents
 ## arrange exponents are 0,1/3 at zero and infinity so they will become ordinary:
L2:=pullback(L1,x-1);
 ##
L3:=pullback(L2,x^3);
localex(L3); ## shows that this is a standard equations with group D_2.
alias(J=RootOf(_Z^2+_Z+1)): 
            ## rearrange the singular points to get closer to desired equation.
p3:=(x+1+2*J)/(x-1-2*J) # sends to 1,1,inf ; inverse: (2*J+1)*(1+y)/(-1+y)
L4:=pullback(L3,p3): L4:=collect(L4,Dx,X->factor(evala(expand(X))));;
localex(L4);
b4:=1/12/(x-1)+1/12/(x+1);
L5:=symmetric_product(L4,Dx-b4);
 ## this is the standard equation for D2.
 ## but the inversion of the above steps gives an ugly formula so let's not use it


## one could maybe find an explicit formula as in the A4 case but this might be
## useless.

##############################################################################
## Formulas for the primitive case, using semi-invariants: section 4.3
#St:=Dx^2+1/6*(8*x+3)/(x+1)/x*Dx+_s/(x+1)^2/x ;
L:=pullback(Dx^2+(7*x-4)/6/x/(x-1)*Dx-_s/x/(x-1), 1/(x+1) );
PSt:=pullback(L,F(x)):
A1:=collect(coeff(PSt,Dx,1),diff(F(x),x$2),factor);
A0:=coeff(PSt,Dx,0);
K := field_extension(transcendental_elements=[_s]);
R := differential_ring(ranking=[F,a1,a0], derivations=[x],field_of_constants=K, notation=diff):
sys:=map(numer,{A1-a1(x),A0-a0(x)}):
GE:= Rosenfeld_Groebner(sys, R);
rewrite_rules(GE[1]);
## result letting G=A0'/A0+2*A1;
pullback_formula:=F(x)=factor(normal( subs(a1(x)=1/2*(G(x)-diff(a0(x),x)/a0(x) ) , rhs(rewrite_rules(GE[1])[1]))));
## F(x) = 9*S*G(x)^2/a0(x);

################################################################################################
######## FORMULAS FOR IMPRIMITIVE GROUPS
Dn:=Dx^2+x*Dx/((x-1)*(x+1))-1/4/((x-1)*n^2*(x+1)):
Q:=eval(Dn,n=2);
SQ:= symmetric_product(Q,Dx-coeff(Q,Dx,1)/2):
r1:=4*x/(x-1)/(x+1);
D2I:=symmetric_product(SQ,Dx+r1/6);
aa0,aa1:=coeff(D2I,Dx,0),coeff(D2I,Dx,1):
gd2I:=normal(2*aa1+diff(aa0,x)/aa0);
G2:=eval(gd2I,x=f(x))*diff(f(x),x);
aa0:=coeff(pullback(D2I,f(x)),Dx,0);
 rel1:=solve(G2=G,diff(f(x),x));
 rel2:=solve(aa0=_a0,diff(f(x),x)^2);
P:=(primpart( numer(normal(rel1^2-rel2)), G)): P:=collect(P,f(x),factor);
PP:=subs(f(x)=F^(1/2),P);