/* Computation of the u_1,...,u_l-free coefficients in E_l */;

%/* VV is the formal parameter for a V-tuple of unknown length;

A0(l,VV):= block ( [],
    if (l=0) then 1 else (
         -(1/l)*sum(A0(l-k, VV)*u[0]*VV(k), k, 1, l)
    )
);

%/* U and V are the actual parameters on which A0 and A1, defined below, will be called */;

U(l) := u[l]; V(l):=v[l];

%/* demo calculations */;

A0(0,V);

A0(1,V);

expand(A0(2,V));

expand(A0(3,V));

/* Computation of the terms E_l, where u_{l-m} occurs linearly and no other u_{k} occurs */;

A1(l,m,VV):= block ( [],
    if (l=0) then 0 else (
      if (m>=1) then (
         -(1/l)*(sum(A0(l-k, VV)*binomial(k,l-m)*VV(k-(l-m)), k,l-m, l)
            +sum(A1(l-k, m-k, VV)*binomial(k,0)*u[0]*VV(k), k, 1, m))
      ) else (
            -(1/l)*sum(A0(l-k, VV)*binomial(k,l)*VV(k-l), k,l, l )
    )
)
);

;

/* demo calculations */;

A1(1, 0,V);

A1(l, 0,V);

A1(l, 1,V);

ratsimp(makefact(A1(l, 1,V)));

A1(l, 2,V);

ratsimp(makefact(A1(l, 2,V)));

A1(l, 3,V);

;

ratsimp(makefact(A1(l, 3,V)));

/* obtaining the coefficient for the u_{l-4}-linear term  of E_{l}*/;

ratsimp(makefact(A1(l, 4,V)));

%/* program for computing power sums */;

expon(x,y) := block([], 
if (y=0) then 1 else x^y
);

expon(x,0);

expon(0,0);

expon(x,1); expon(2,l);

 p(t1,t2,t3,l):=expon(t1,l)+expon(t2,l)+expon(t3,l);

%/* Calculations to study the solutions of ADE with the lhs below:;

lhs(x):=P(x)*P(x-2)*P(x-3)-2*P(x-1)*P(x-1)*P(x-3)+P(x-1)*P(x-2)*P(x-2)
+P(x)*P(x-1)*P(x-3)-2*P(x)*P(x-2)*P(x-2)+P(x-1)*P(x-1)*P(x-2);

;

%/* ab(i,j,k) is the coefficient of P(x-i)P(x-j)P(x-k) in the l.h.s. above;

ab(i,j,k) := block ( [],
    retVal:0,
    if (i=0 and j=2 and k=3) then ( retVal:1),
    if (i=1 and j=1 and k=3) then ( retVal:-2),
    if (i=1 and j=2 and k=2) then ( retVal:1),
    if (i=0 and j=1 and k=3) then ( retVal:1),
    if (i=0 and j=2 and k=2) then ( retVal:-2),
    if (i=1 and j=1 and k=2) then ( retVal:1),
    retVal 
);

%/* computing the power sums corresponding to the t=(i,j,k) of the given l.h.s., l is the power;

V023(l):=p(0,2,3,l);

V113(l):=p(1,1,3,l);
V122(l):=p(1,2,2,l);
V013(l):=p(0,1,3,l);
V022(l):=p(0,2,2,l);
V112(l):=p(1,1,2,l);

%/* demo calculations */;

A1(0, 0,V113);

A1(1, 0,V113);

A1(l, 0,V113);

%/* the coefficient of the linear occurence u_{l} in S_{l}(u0,u1,..,u_l), must be zero  */;

%/* m=0 */;

AS0(l):=ab(0,2,3)*A1(l, 0,V023)+ab(1,1,3)*A1(l, 0,V113)+ab(1,2,2)*A1(l, 0,V122)+
ab(0,1,3)*A1(l, 0,V013)+ab(0,2,2)*A1(l, 0,V022)+ab(1,1,2)*A1(l, 0,V112);

ratsimp(makefact(AS0(l)));

/* the coefficient of the linear occurence u_{l-1} in S_{l}(u0,u1,..,u_l), must be zero */;

%/* m=1 */;

AS1(l):=ab(0,2,3)*A1(l, 1,V023)+ab(1,1,3)*A1(l, 1,V113)+ab(1,2,2)*A1(l, 1,V122)+
ab(0,1,3)*A1(l, 1,V013)+ab(0,2,2)*A1(l, 1,V022)+ab(1,1,2)*A1(l, 1,V112);

ratsimp(makefact(AS1(l)));

/* the coefficient of the linear occurence u_{l-2}  in S_{l}(u0,u1,..,u_l), must be zero */;

%/* m=2 */;

AS2(l):=ab(0,2,3)*A1(l, 2,V023)+ab(1,1,3)*A1(l, 2,V113)+ab(1,2,2)*A1(l, 2,V122)+
ab(0,1,3)*A1(l, 2,V013)+ab(0,2,2)*A1(l, 2,V022)+ab(1,1,2)*A1(l, 2,V112);

ratsimp(makefact(AS2(l)));

/* the coefficient of the linear occurence u_{l-3} in S_{l}(u0,u1,..,u_l), must be zero */;

%/* m=3 */;

AS3(l):=ab(0,2,3)*A1(l, 3,V023)+ab(1,1,3)*A1(l, 3,V113)+ab(1,2,2)*A1(l, 3,V122)+
ab(0,1,3)*A1(l, 3,V013)+ab(0,2,2)*A1(l, 3,V022)+ab(1,1,2)*A1(l, 3,V112);

ratsimp(makefact(AS3(l))); ratsimp(makefact(AS3(6)));

/* the coefficient of the linear occurence u_{l-4} in E_{l} */;

%/* m=4 */;

AS4(l):=ab(0,2,3)*A1(l, 4,V023)+ab(1,1,3)*A1(l, 4,V113)+ab(1,2,2)*A1(l, 4,V122)+
ab(0,1,3)*A1(l, 4,V013)+ab(0,2,2)*A1(l, 4,V022)+ab(1,1,2)*A1(l, 4,V112);

ratsimp(makefact(AS4(l)));

%/* check-up calculations */

E(l, UU, VV) := block ( [],
if (l=0) then 
    1
    else ( 
        -(1/l)*sum(E(l-k,UU, VV)*sum(binomial(k,ll)*VV(k-ll)*UU(ll), ll, 0, k), k, 1, l)
    )
);

E(0, U, V);

E(1, U, V);

A1u1check:coeff(E(1, U, V), u[1]);

E3:expand(E(3, U, V));

E6:expand(E(6, U, V));

A6u3check:coeff(E6, u[3]);

A6u3:expand(ratsimp(A1(6,3, V)));

A6u3checkCorr:expand(ratsubst(0, u[2], ratsubst(0, u[1], A6u3check)));

ratsimp(A6u3checkCorr - A6u3);

A6u2check:coeff(E6, u[2]);

A6u2checkCorr:expand(ratsubst(0, u[4], (ratsubst(0, u[3], ratsubst(0, u[1], A6u2check)))));

A6u2:expand(ratsimp(A1(6,4, V)));

ratsimp(A6u2-A6u2checkCorr);

A60check1:expand(ratsubst(0, u[6], ratsubst(0, u[5], ratsubst(0, u[4], E(6, U, V) ))));

A60check:expand(ratsubst(0, u[3], ratsubst(0, u[2], ratsubst(0, u[1], A60check1 ))));

A60:expand(A0(6,V));

ratsimp(A60-A60check);

A21check:expand(coeff(expand(E(2, U, V)), u[1]));

A21:expand(A1(2,1,V));

A22check:expand(coeff(expand(E(2, U, V)), u[2]));

A22:expand(A1(2,0,V));

A31check:expand(coeff(expand(ratsubst(0, u[2],E(3, U, V))), u[1]));

A31:expand(A1(3,2,V));

A32check:expand(coeff(expand(ratsubst(0, u[1],E(3, U, V))), u[2]));

A32:expand(A1(3,1,V));

A33check:expand(coeff(expand(ratsubst(0, u[1],E(3, U, V))), u[3]));

A33:expand(A1(3,0,V));

%/* computing solutions of lhs(P(x),P(x-1),P(x-2),P(x-3)) = 1 */

%/*trying d=6 */

lhs(x):=P(x)*P(x-2)*P(x-3)-2*P(x-1)*P(x-1)*P(x-3)+P(x-1)*P(x-2)*P(x-2)
+P(x)*P(x-1)*P(x-3)-2*P(x)*P(x-2)*P(x-2)+P(x-1)*P(x-1)*P(x-2);

P(x):=x^6+a5*x^5+a4*x^4+a3*x^3+a2*x^2+a1*x+a0;

LHS:expand(lhs(x));

S18:ratcoeff(LHS,x^18);

S17:ratcoeff(LHS,x^17);S16:ratcoeff(LHS,x^16);S15:ratcoeff(LHS,x^15);S14:ratcoeff(LHS,x^14);

S13:ratcoeff(LHS,x^13);S12:ratcoeff(LHS,x^12);S11:ratcoeff(LHS,x^11);S10:ratcoeff(LHS,x^10);

solve([S10, S11, S12], [a2, a3,a4]);

S9:ratcoeff(LHS,x^9);S8:ratcoeff(LHS,x^8);S7:ratcoeff(LHS,x^7);S6:ratcoeff(LHS,x^6);

solve([S10, S11, S12, S9], [a2, a3,a4,a1]);

solve([S10, S11, S12, S8], [a2, a3,a4,a1]);

solve([S10, S11, S12, S9, S8], [a2, a3,a4,a1, a0]);

solve([S10, S11, S12, S9, S8, S7, S6], [a2, a3,a4,a1, a0]);

S5:ratcoeff(LHS,x^5);S4:ratcoeff(LHS,x^4);S3:ratcoeff(LHS,x^3);S2:ratcoeff(LHS,x^2);

solve([S10, S11, S12, S9, S8, S7, S6, S5], [a2, a3,a4,a1, a0]);

solve([S10, S11, S12, S9, S8, S7, S6, S5, S4], [a2, a3,a4,a1, a0]);

solve([S10, S11, S12, S9, S8, S7, S6, S5, S4, S3, S2], [a2, a3,a4,a1, a0]);

S1:ratcoeff(LHS,x^1);S0:ratcoeff(LHS,x, 0);

solve([S10, S11, S12, S9, S8, S7, S6, S5, S4, S3, S2, S1], [a2, a3,a4,a1, a0]);

/* solve([S10, S11, S12, S9, S8, S7, S6, S5, S4, S3, S2, S1, S0], [a2, a3,a4,a1, a0]);*/

solve([S0-1, a2=(5*a5^4-630*a5^2+6993)/432], [a0]); /* subsitute previous equations in this solution and solve w.r.t a5;*/;

T1:factor(subst((5*a5^4-630*a5^2+6993)/432, a2, S0-1));

T2:factor(subst((5*a5^3-315*a5)/54, a3, T1));

T3:factor(subst((5*a5^2-105)/12, a4, T2));

T4:factor(subst((a5^5-210*a5^3+6993*a5)/1296, a1, T3));

T5:factor(subst((a5^6-315*a5^4+20979*a5^2-164025)/46656, a0, T4));

%/* trying d=5 */

lhs(x):=P(x)*P(x-2)*P(x-3)-2*P(x-1)*P(x-1)*P(x-3)+P(x-1)*P(x-2)*P(x-2)
+P(x)*P(x-1)*P(x-3)-2*P(x)*P(x-2)*P(x-2)+P(x-1)*P(x-1)*P(x-2);

P(x):=x^5+a4*x^4+a3*x^3+a2*x^2+a1*x+a0;

LHS:expand(lhs(x));

S15:ratcoeff(LHS,x^15);S14:ratcoeff(LHS,x^14);

S13:ratcoeff(LHS,x^13);S12:ratcoeff(LHS,x^12);S11:ratcoeff(LHS,x^11);S10:ratcoeff(LHS,x^10);

;

S9:ratcoeff(LHS,x^9);S8:ratcoeff(LHS,x^8);S7:ratcoeff(LHS,x^7);S6:ratcoeff(LHS,x^6);

solve([S9, S8, S7, S6], [a3, a2, a1, a0]);

S5:ratcoeff(LHS,x^5);S4:ratcoeff(LHS,x^4);S3:ratcoeff(LHS,x^3);S2:ratcoeff(LHS,x^2);

solve([S9, S8, S7, S6, S5], [a3, a2, a1, a0]);

solve([S9, S8, S7, S6, S5, S4], [a3, a2, a1, a0]);

solve([S9, S8, S7, S6, S5, S4, S3], [a3, a2, a1, a0]);

solve([S9, S8, S7, S6, S5, S4, S3, S2], [a3, a2, a1, a0]);

S1:ratcoeff(LHS,x^1);S0:ratcoeff(LHS,x, 0);

solve([S9, S8, S7, S6, S5, S4, S3, S2, S1], [a3, a2,a1, a0]);

/* solve([S9, S8, S7, S6, S5, S4, S3, S2, S1, S0], [a3, a2,a1, a0]);*/

solve([S9, S8, S7, S6, S5, S4, S3, S2, S1, S0-1], [a3, a2,a1, a0]);

T1:factor(subst((2*a4^2-25)/5, a3,S0-1 ));

T2:factor(subst((2*a4^3-75*a4)/25, a2,T1 ));

T3:factor(subst((a4^4-75*a4^2+500)/125, a1,T2 ));

T4:factor(subst((a4^5-125*a4^3+2500*a4)/3125, a0,T3 ));

%/* trying d=4 */

lhs(x):=P(x)*P(x-2)*P(x-3)-2*P(x-1)*P(x-1)*P(x-3)+P(x-1)*P(x-2)*P(x-2)
+P(x)*P(x-1)*P(x-3)-2*P(x)*P(x-2)*P(x-2)+P(x-1)*P(x-1)*P(x-2);

P(x):=x^4+a3*x^3+a2*x^2+a1*x+a0;

LHS:expand(lhs(x));

;

S12:ratcoeff(LHS,x^12);S11:ratcoeff(LHS,x^11);S10:ratcoeff(LHS,x^10);

;

S9:ratcoeff(LHS,x^9);S8:ratcoeff(LHS,x^8);S7:ratcoeff(LHS,x^7);S6:ratcoeff(LHS,x^6);

;

S5:ratcoeff(LHS,x^5);S4:ratcoeff(LHS,x^4);S3:ratcoeff(LHS,x^3);S2:ratcoeff(LHS,x^2);

solve([S6, S5], [a2, a1]);

solve([S6, S5, S4], [ a2, a1, a0]);

solve([S6, S5, S4, S3], [a2, a1, a0]);

solve([S6, S5, S4, S3, S2], [ a2, a1, a0]);

S1:ratcoeff(LHS,x^1);S0:ratcoeff(LHS,x, 0);

solve([S6, S5, S4, S3, S2, S1], [a2,a1, a0]);

solve([S6, S5, S4, S3, S2, S1, S0], [ a2,a1, a0]);

/*solve([S6, S5, S4, S3, S2, S1, S0-1], [ a2,a1, a0]);*/

T1:factor(subst((3*a3^2-20)/8, a2, S0-1));

T2:factor(subst((a3^3-20*a3)/16, a1,T1));

T3:factor(subst((a3^4-40*a3^2+144)/256, a0,T2));

%/* trying d=3 */

lhs(x):=P(x)*P(x-2)*P(x-3)-2*P(x-1)*P(x-1)*P(x-3)+P(x-1)*P(x-2)*P(x-2)
+P(x)*P(x-1)*P(x-3)-2*P(x)*P(x-2)*P(x-2)+P(x-1)*P(x-1)*P(x-2);

P(x):=x^3+a2*x^2+a1*x+a0;

LHS:expand(lhs(x));

;

S9:ratcoeff(LHS,x^9);S8:ratcoeff(LHS,x^8);S7:ratcoeff(LHS,x^7);S6:ratcoeff(LHS,x^6);

;

S5:ratcoeff(LHS,x^5);S4:ratcoeff(LHS,x^4);S3:ratcoeff(LHS,x^3);S2:ratcoeff(LHS,x^2);

solve([S3, S2], [a1, a0]);

S1:ratcoeff(LHS,x^1);S0:ratcoeff(LHS,x, 0);

solve([ S3, S2, S1], [a1, a0]);

/* solve([S3, S2, S1, S0], [a1, a0]);*/

T1:factor(subst((a2^2-3)/3,a1,  S0-1));

T2:factor(subst((a2^3-9*a2)/27,a0, T1));

%/* trying d=2 */

lhs(x):=P(x)*P(x-2)*P(x-3)-2*P(x-1)*P(x-1)*P(x-3)+P(x-1)*P(x-2)*P(x-2)
+P(x)*P(x-1)*P(x-3)-2*P(x)*P(x-2)*P(x-2)+P(x-1)*P(x-1)*P(x-2);

P(x):=x^2+a1*x+a0;

LHS:expand(lhs(x));

;

S6:ratcoeff(LHS,x^6);

;

S5:ratcoeff(LHS,x^5);S4:ratcoeff(LHS,x^4);S3:ratcoeff(LHS,x^3);S2:ratcoeff(LHS,x^2);

S1:ratcoeff(LHS,x^1);S0:ratcoeff(LHS,x, 0);

solve([S0-1], [a0]);