##
## Draft Implementation of Hoeij-Weil Issac'05
## Note: this is not at all optimized: in its current state, it is slower than the
## fast kovacicsols implemented by MvH in Maple; 
## for test purpose only, at this stage.
#restart:  
with(DEtools):
_Envdiffopdomain:=[Dx,x]:
##


## Here should be: front procedure with type check, normalization, and reducibility cases.
## L is assumed to be a polynomial in Dx=d/dx, independent variable is x.
pullback_solve := proc(L)
 local ex, v, NL;
    ex := expsols(diffop2de(L, y(x)), y(x));
    if ex <> [] then print(reducible); RETURN(ex) end if;
    v := -1/2*coeff(L, Dx, 1);
    userinfo(2,DEtools,"shift by ",-1/2*coeff(L, Dx, 1));
    NL:=symmetric_product(L, Dx + v);
    userinfo(3,DEtools,"shifted unimodular operator is ",NL);
    pullback_solve_irred(NL);
    map(X -> simplify( X*factor(expand(exp(int(v, x)))) , symbolic), %);
end proc:


## Input: L:=Dx^2+a1*Dx+a0 is a unimodular operator. 
## I assume the reducible case has already been treated
pullback_solve_irred := proc(L)
local a0, a1, inv,m; 
    a0, a1 := seq(coeff(L, Dx, i), i = [0, 1]);
    m:=4; print(m);
    inv := invar2(L, m);
    if nops([inv]) = 1 then userinfo(2,DEtools,"projective group is dihedral");RETURN(dihedral(a0, a1, inv)) end if;
    if nops([inv]) = 2 then userinfo(2,DEtools,"projective group is D2");
    RETURN(quaternion(a0, a1, inv)) end if;
    m:=6:print(m);
    inv := invar2(L, m);
    if nops([inv]) = 1 then userinfo(2,DEtools,"projective group is tetraedral");
        userinfo(3,DEtools,"found invariant of degree ",m," with value",inv[1]);
        RETURN(tetraedral(a0, a1, inv)) end if;
    m:=8:print(m);
    inv := invar2(L, m);
    if nops([inv]) = 1 then userinfo(2,DEtools,"projective group is octaedral");
        userinfo(3,DEtools,"found invariant of degree ",m," with value",inv[1]);
        RETURN(octaedral(a0, a1, inv)) end if;
    m:=12:print(m);
    inv := invar2(L, m);
    if nops([inv]) = 1 then userinfo(2,DEtools,"projective group is icosaedral");
        userinfo(3,DEtools,"found invariant of degree ",m," with value",inv[1]);
        RETURN(icosaedral(a0, a1, inv)) end if;
    RETURN([0])
end proc:
oursolver:=pullback_solve:

###############################################################################
#############################  Ancillaries: symmetric_powers and invariants.
# Symmetric- power via bmw - there is also a square free faster one
sympow := proc (L, m) local a0, a1, M, i; 
    a0, a1 := seq(coeff(L, Dx, i), i = [0, 1]); 
    M[0] := 1; M[1] := Dx; 
    for i to m do 
        M[i+1] := collect(mult(Dx, M[i])+i*a1*M[i]+i*(m-i+1)*a0*M[i-1], Dx, normal) 
    end do 
end proc:
## computation of invariant values; one can do much better for C(x)
invar2 := proc(L, m) op(ratsols(diffop2de(sympow(L, m), y(x)), y(x))) end proc:
semiinvar2 := proc(L, m) op(expsols(diffop2de(sympow(L, m), y(x)), y(x))) end proc:

## note: the next is homemade, should rather use PDEtools[dchange]
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->normal(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:
# local exponents at singular points
localex:=proc(L) seq(gen_exp(L,_t,x=p),p=singularites(L)) end proc:
###############################################################################
#############################  List of Standard Equations and their solutions

###############################################################################################
## general principles:
## 0) already known: pullback can be found from invariants but (1) need proper pullback for #that
## or (2) a formula that uses an invariant fraction which is homogeneous of degree 0 - but #then
#we have the problem that the relevant invariant are known only up to constants -> need to use
#logarithmic derivatives.
## 1) good normal form if proper pullback -> set invariant value to 1
###############################################################################################
################ DIHEDRAL CASE
## L=D^2+a1*D+a0, inv is the invariant of degree 4. Returns pullback & sol if Dn or sol otherwise.
dihedral := proc(a0, a1, inv)
local b, A, n, A0;
    b := 1/4*diff(inv, x)/inv;
    A0 := a0 + a1*b + diff(b, x) + b^2;
    n := find_n(A0);
    if n < infinity then 
        map(X -> X*inv^(1/4), algebraic_diedral(A0, -1/2*diff(A0, x)/A0, n))
    else map(X -> X*inv^(1/4), [exp(int(sqrt(-A0), x)), exp(int(-sqrt(-A0), x))])
    end if
end proc:



## decide if exp(int(-A,x)) is algebraic and measure/bound its degree n
## not yet implemented.
find_n:=proc(A)
infinity:
end:

algebraic_diedral := proc(a0, a1, n)
local L2n, inv, B, f;
    L2n := diff(y(x), x $ 2) + a1*diff(y(x), x) + 4*n^2*a0*y(x);
    inv := ratsols(L2n, y(x))[1];
    B := 1/4*diff(inv, x)^2/(inv^2*n^2*a0);
    f:=simplify(1/sqrt(B + 1),symbolic): ## this is the pullback function
    userinfo(2,DEtools,"found invariant of degree ",n," with value",inv);
    userinfo(2,DEtools,"pullback map is ",f);
    eval([(x+(x^2-1)^(1/2))^(1/2/n), (x+(x^2-1)^(1/2))^(-1/2/n)],x=f);
end proc:



################ QUATERNION CASE
quaternion := proc(a0, a1, inv)
local i6, b, NL, G, f, F; #option trace;
    i6 := invar2(Dx^2 + a1*dx + a0, 6);
    userinfo(3,DEtools,"found invariant of degree ",6," with value",i6);
    b := diff(i6, x)/i6;
    NL := symmetric_product(Dx^2 + a1*Dx + a0, Dx + 1/6*b);
    A0, A1 := seq(coeff(NL, Dx, i), i = [0, 1]);
    G := 2*A1 + diff(A0, x)/A0;
    f := simplify(1/5*(25 + 320*A0/G^2)^(1/2), symbolic);
    F := solve(3*sqrt(-3)*(X^2 - 1) - (X^3 - 9*X)*f, X)[1];
    userinfo(2,DEtools,"pullback map is ",F);
    eval([hypergeom([-1/4, 1/4], [1/2], x/2 + 1/2)/((x + 1)^(1/12)*(x - 1)^(1/12)),
        (x/2 + 1/2)^(1/2)*hypergeom([1/4, 3/4], [3/2], x/2 + 1/2)/((x + 1)^(1/12)*(x - 1)^(1/12))],
        x = F)
end proc:

## shorter version,using only tetraedral solution: easier but less minimal.
## can be used to prove the formula of approach 2.
quaternion1 := proc(a0, a1, inv) 
  local i6; i6 := invar2(Dx^2 + a1*dx + a0, 6); tetraedral(a0, a1, inv) 
end proc:



################ TETRAEDRAL CASE
tetraedral := proc(a0, a1, inv)
local A0, A1, G, f; #option trace;
    symmetric_product(Dx^2 + a1*Dx + a0, Dx + 1/6*diff(inv, x)/inv);
    A0, A1 := coeff(%, Dx, 0), coeff(%, Dx, 1);
    G := 2*A1 + diff(A0, x)/A0;
    f := simplify(1/5*(25 + 320*A0/G^2)^(1/2), symbolic);
    userinfo(2,DEtools,"pullback map is ",f);
    eval([hypergeom([-1/12, 1/4], [1/2], 2*x/(x - 1))/(x/(x - 1))^(1/12),
        2^(1/2)*(x/(x - 1))^(5/12)*hypergeom([5/12, 3/4], [3/2], 2*x/(x - 1))], x = f)
end proc:



################ OCTAEDRAL CASE
octaedral := proc(a0, a1, inv)
local A0, A1, G, f;
    symmetric_product(Dx^2 + a1*Dx + a0, Dx + 1/8*diff(inv, x)/inv);
    A0, A1 := coeff(%, Dx, 0), coeff(%, Dx, 1);
    G := 2*A1 + diff(A0, x)/A0;
    f := normal(-7/144*G^2/A0);
    userinfo(2,DEtools,"pullback map is ",f);
    eval([hypergeom([-1/24, 5/24], [1/2], x)/(x - 1)^(1/24),
        x^(1/2)*hypergeom([11/24, 17/24], [3/2], x)/(x - 1)^(1/24)], x = f)
end proc:

################ ICOSAEDRAL CASE
icosaedral := proc(a0, a1, inv)
local A0, A1, G, pullback_A5S;
    symmetric_product(Dx^2 + a1*Dx + a0, Dx + 1/12*diff(inv, x)/inv);
    A0, A1 := coeff(%, Dx, 0), coeff(%, Dx, 1);
    G := 2*A1 + diff(A0, x)/A0;
    f := normal(11/400*G^2/A0);
    userinfo(2,DEtools,"pullback map is ",f);
    eval([hypergeom([-1/60, 19/60], [1/2], -x)/(x + 1)^(1/60),
        (-x)^(1/2)*hypergeom([29/60, 49/60], [3/2], -x)/(x + 1)^(1/60)], x = f)
end proc:

###################################################
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); 
A4I:=pullback(L,-2*x/(x+1));
b8:=-1/24/(x+1);
L:=symmetric_product(St_S4,Dx+b8);
S4I:=pullback(L,x/(-x+1));
######
St_Dn:=Dx^2+x*Dx/((x-1)*(x+1))-1/4/((x-1)*n^2*(x+1)):
Dn:=St_Dn:
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);