## Copyright (c) 2022,  BUAA. All rights reserved.
## Author: Hanqian Fang <hqfang_math@163.com>, Lixin Du <dulixin17@mails.ucas.ac.cn>
## A Maple package for testing whether two given multivariate polynomials are shift equivalent
## if so, then give all shifts
## Modified on April 16, 2022.
###############################################################

SET := module()
  description "polynomial shift equivalent test"; 
  export 
         SET_G,
         SET_KS,
         SET_DOS,
         SET_ADC_expansion,
         SET_ADC_derivation,
         Shifts;
        
  local 
         SET_G_F,
         SelectCoeff,
         HomogeneousTerms,
         DegClassify,
         HomogeneousNoMoreOne,
         IsLeaf,
         LevelPartition,
         SelectCoeff_ADC;

  option package;

######################################################################
# Export
# Name: SET_G
# Calling sequence:
#                          SET_G(p, q, x)
#                          SET_G(p, q, x, output = 'type')
# Input: p, q,  two multivariate polynomials
#        x,  a variable name or a list of variable names
#        output,  (optional) should be a type of solutions in the form:
#                         output = 'Z'.
# Output: L,  a list representing the shift equivalence relation of p and q w.r.t. x;
#         [],  if there is no integer solutions.
#              If there is no condition on 'type', then return L as a general solution of
#                                     q(x) = p(x + a).
#              If 'type' = 'Z', then return L as a general integer solution.
# Remark: A general solution may contains some free variables.
######################################################################

SET_G := proc(p, q, x, {output := 'F'})
    local n, i, j, eqs, L, sol, a;
    n := nops(x);
    i := 'i';
    a := [seq(i[j], j = 1 .. n)];
    L := SET_G_F(p, q, x, a);
    if L = [] then return []; end if;
    if output = 'Z' then
        eqs := {seq(a[j] - L[j], j = 1 .. n)};
        sol := isolve(eqs);
        if evalb({sol} = {}) then
            return [];
        else       
            L := subs(sol, a);
            return L;
        end if;
    end if;
    return L;
end proc:



######################################################################
# Export
# Name : SET_G_F
# Calling sequence:
#                              SET_G_F (p, q, x, a)
# Input: p, q,  two multivariate polynomials
#            x, a,  a variable name or a list of variable names
# Output: k,  a list representation of a general solution of
#                                 q(x) = p(x + a).
#              [],  if there is no solutions.
# Remark: A general solution may contains some free variables in a.
######################################################################
SET_G_F:= proc(p, q, x, a) 
    local n, ly, sy, f, g, sub, d, i, j, k, m, la, sa, eq, eqs, l_y0, d_sub, sz, lz, sub0, c; 
    n := nops(x); 
    ly := convert(x, list); 
    sy := convert(ly, set); 
    f := expand(p); 
    g := expand(q);
    la := convert(a, list); 
    sa := convert(la, set);
    if f = 0 and g = 0 then return la;  end if; 
    sub := expand(f - g); 
    d := degree(f, sy); 
    if d <= degree(sub, sy) then return []; end if;
    l_y0 := [seq(ly[j] = 0, j = 1 .. n)];

#### calculate the situation when d<=1 ###

    if d <= 1 then 
        eqs := {subs([seq(ly[j] = la[j], j = 1 .. n)], f) - subs(l_y0, g)}; 
        eqs := solve(eqs, sa); 
        if eqs = () then return []; end if; 
        eqs := subs(eqs, la); 
        return eqs; 
    end if; 

#### calculate \cap_{j} S^{(j)} ###

    eqs := {}; 
    for j to n do 
        eq := SET_G_F(diff(p, ly[j]), diff(q, ly[j]), x, a); 
        if eq = [] then return []; end if;
        eqs := eqs union {seq(la[k] - eq[k], k = 1 .. n)}; 
    end do; 

#### substitute the solution (A1,...,An) and get the last equation f(A1,...,An)-g(0,...,0)=0###

    eqs := solve(eqs, sa);
    if eqs = () then return []; end if; 
    eq := subs(eqs, la); 
    sub := expand(subs([seq(ly[j] = eq[j], j = 1 .. n)], f) - subs(l_y0, g)); 
    d_sub := degree(sub, sa); 

    if sub = 0 then return eq; end if; 
    if d_sub = 0 then return []; end if; 

#### f(A1,...,An)-g(0,..., 0)= c '(\sum_{1<=j<=m} c_j*λ_j-c_0)^δ ###

    sz := indets(sub) intersect sa; 
    lz := convert(sz, list); 
    m := numelems(lz); 
    sub0 := subs([seq(lz[j] = 0, j = 1 .. m)], sub); 
    if sub0 = 0 then 
        eq := lz[1]; 
        c := SelectCoeff(subs([seq(lz[j] = 0, j = 2 .. m)], sub), [lz[1]], lz[1]^d_sub); 
        for j from 2 to m do 
            eq := eq + lz[j]*SelectCoeff(subs([seq(lz[k] = 0, k = 2 .. j - 1), seq(lz[k] = 0, k = j + 1 .. m)], sub), [lz[1], lz[j]], lz[j]*lz[1]^(d_sub - 1))/(c*d_sub); 
        end do; 
    else 
        eq := 1; 
        c := subs([seq(lz[j] = 0, j = 1 .. m)], sub); 
        for j to m do 
            eq := eq + lz[j]*SelectCoeff(subs([seq(lz[k] = 0, k = 1 .. j - 1), seq(lz[k] = 0, k = j + 1 .. m)], sub), [lz[j]], lz[j])/(c*d_sub); 
        end do; 
    end if; 
    eqs := solve(eqs union {eq = 0}, sa); 
    if eqs = () then return []; end if; 
    return subs(eqs, la); 
end proc;

######################################################################
# Export
# Name : SET_KS
# Calling sequence:
#                          SET_KS(p, q, x)
#                          SET_KS(p, q, x, output = 'type')
# Input: p, q,  two multivariate polynomials
#        x,  a variable name or a list of variable names
#        output,  (optional) should be a type of solutions in the form:
#                         output = 'Z'.
# Output: k,  a list representing the shift equivalence relation of p and q w.r.t. x;
#         [],  if there is no integer solutions.
#              If there is no condition on 'type', then return k as a general solution of
#                                     q(x) = p(x + a).
#              If 'type' = 'Z', then return k as a general integer solution.
# Remark: A general solution may contains some free variables.
######################################################################

SET_KS := proc(p, q, x, {output := 'F'}) 
    local n, a, f, g, S, y, J, G, Eq, sol, k, i, j, vars, globalvars, parameter;
    n := nops(x);
    y := convert(x, list);
    g := expand(p);
    f := expand(q);
    i := 'i';
    a := [seq(i[j], j = 1 .. n)];
    if f = 0 and g = 0 then
        return a;
    end if;
    if degree(f, convert(y, set)) <= degree(expand(f - g), convert(y, set)) then
        return [];
    end if;
#### extract the coefficients of p(x+i) - q(x) w.r.t x, which are polynomial in i and t ###
    g := expand(f - Shifts(g, x, a));
    S := PolynomialIdeals:-PolynomialIdeal({coeffs(g, x)});
#### compute the Groebner basis of radical S w.r.t. graded lex order of i ###
    J := PolynomialIdeals:-Radical(S);
    G := Groebner:-Basis(PolynomialIdeals:-Generators(J), grlex(op(a)));
#### extract the coefficients of every polynomial of G w.r.t parameter t, if exist ###
    Eq := convert(G, set);
    vars := indets(Eq);
    globalvars := convert(a, set);
    parameter := vars minus globalvars;
    Eq := map(coeffs, Eq, parameter);
#### find solutions of a linear system ###
    if output = 'F' then
        sol := solve(Eq, convert(a, set));
    elif output = 'Z' then
        sol := isolve(Eq);
    end if;
    if evalb({sol} = {}) then
        k := [];
    else
        sol := convert(sol, set);
        k := subs(sol, a);
    end if;
    k;
end proc;

#####################################################################
# Export
# Name : SET_DOS
# Calling sequence:
#                          SET_DOS(p, q, x)
#                          SET_DOS(p, q, x, output = 'type')
# Input: p, q,  two multivariate polynomials
#        x,  a variable name or a list of variable names
#        output,  (optional) should be a type of solutions in the form:
#                         output = 'Z'.
# Output: k,  a list representing the shift equivalence relation of p and q w.r.t. x;
#         [],  if there is no integer solutions.
#              If there is no condition on 'type', then return k as a general solution of
#                                     q(x) = p(x + a).
#              If 'type' = 'Z', then return k as a general integer solution.
# Remark: A general solution may contains some free variables.
#####################################################################
SET_DOS := proc(p, q, x, {output := 'F'}) 
    local n, y, f, g, d, hom_f, deg_f, hom_sub, deg_sub, i, j, a, sol_special, eqs, len, index_f, expr, k, dk, lindex, lterms, lindex1, lterms1, l, term, index, dterm, index1, m, o, s, eq, para;
    n := nops(x);
    y := convert(x, list);
    f := expand(p);
    g := expand(q);
    i := 'i';
    a := [seq(i[j], j = 1 .. n)];
    if f = 0 and g = 0 then
        return a;
    end if;
    d := degree(f, convert(y, set));
####同时比较了degree(f)=degree(g)及H(f(x + a))^d = H(g)^d###
    if d <= degree(expand(f - g), convert(y, set)) then
        return [];
    end if;
    hom_f, deg_f := HomogeneousTerms(f, y);
    hom_sub, deg_sub := HomogeneousTerms(f - g, y);
    sol_special := [seq(1, j = 1 .. n)];
    eqs := {};
    len := numelems(deg_f);
    index_f := table([seq(deg_f[j] = [[]], j = 1 .. len)]); #用于记录每个H(f)^k的求偏导过程中，偏导不为0的下标，如n = 2时[[1, 1], [1, 2], [2, 1], [2, 2]]
####以下只需比较H(f(x + a))^(d - j) = H(g)^(d - j)*j and H(g)^(d - j)*j = 1, 2, ... d ###
    for j to d do
####求H(f)^(d - j) - H(g)^(d - j)
        if (d - j) in deg_sub then
            expr := hom_sub[d - j][1];
        else
            expr := 0;
        end if;
####求\sum_i(1/(j!)a_i1...a_ij*\partial_j(H^(d)(f)))+...+\sum_i(a_i*\partial_1(H^(d-j+1)(f)))###
        k := 1;
        dk := deg_f[k];
        lindex := index_f[dk];#目前偏导不为零的下标列表 
        lterms := hom_f[dk];#求了某阶偏导后H(f)^dk列表
        while d - j < dk do
            lindex1 := [];#即将产生的lindex列表
            lterms1 := [];#即将产生的lterms列表
            for l to numelems(lindex) do
####对各lterms中的偏导再求一阶导###
                term := lterms[l];
                index := lindex[l];
####求\sum_i(1/(dk - d + j)!a_i1... a_i(dk - d + j)*\partial_(dk - d + j)(H(f)^dk)), 注意避免重复计算偏导###
                if dk - d + j = 1 then
                    m := 1;
                else
                    m := index[dk - d + j - 1];
                end if;
                for m from m to n do
                    dterm := diff(term, y[m]);
                    if dterm <> 0 then
                        index1 := [op(index), m];
                        lindex1 := [op(lindex1), index1];
                        lterms1 := [op(lterms1), dterm];
                        if 1 < dk - d + j then
                            s := 1;
                            for o to dk - d + j - 1 do
                                dterm := dterm*sol_special[index1[o]]/s;
                                if index1[o + 1] = index1[o] then
                                    s := s + 1;
                                else
                                    s := 1;
                                end if;
                            end do;
                            dterm := dterm*sol_special[index1[o]]/s;
                        else
                            dterm := dterm*i[index1[1]];
                        end if;
                        expr := expr + dterm;
                    end if;
                end do;
            end do;
            index_f[dk] := lindex1;
            hom_f[dk] := lterms1;
            k := k + 1;
            if len < k then
                break;
            end if;
            dk := deg_f[k];
            lindex := index_f[dk];
            lterms := hom_f[dk];
        end do;
#### find solutions of a linear system ###
        eq := solve({coeffs(expand(expr), y)}, convert(a, set));
        if output = 'Z' then 
            if eq = () or isolve(eq)=() then
                 return [];
            end if;
            eqs := solve(eqs union eq, convert(a, set));
            if eqs = () then
                return [];
            else
                 eq:=isolve(eqs);
            end if;
            if eq = () then
                return [];
            end if;    
            eq := subs(eq, a);
        elif output = 'F' then
            if eq = () then
                return [];
            end if;
            eqs := solve(eqs union eq, convert(a, set));
            if eqs = () then
                return [];
            end if;
            eq := subs(eqs, a);
        end if;
        para := convert(indets(eq), list);
        k := nops(para);
        sol_special := subs(zip((x, y) -> x = y, para, convert(LinearAlgebra:-ZeroVector(k), list)), eq);
    end do;
    return eq;
end proc;
#####################################################################
# Export
# Name : SET_ADC_expansion
# Calling sequence:
#                          SET_ADC_expansion(p, q, x)
#                          SET_ADC_expansion(p, q, x, output = 'type')
# Input: p, q,  two multivariate polynomials
#        x,  a variable name or a list of variable names
#        output,  (optional) should be a type of solutions in the form:
#                         output = 'Z'.
# Output: k,  a list representing the shift equivalence relation of p and q w.r.t. x;
#         [],  if there is no integer solutions.
#              If there is no condition on 'type', then return k as a general solution of
#                                     q(x) = p(x + a).
#              If 'type' = 'Z', then return k as a general integer solution.
# Remark: A general solution may contains some free variables.
######################################################################
SET_ADC_expansion := proc(p, q, x, {output := 'F'})
    local n, y, f, g, dl_sub, ht_sub, d_sub, hl, i, j, a, sol_particular, eqs, expr, exprs, tail, eq, para, k;
    n := nops(x);
    y := convert(x, list);
    g := expand(p);
    f := expand(q);
    i := 'i';
    a := [seq(i[j], j = 1 .. n)];
    if f = 0 and g = 0 then
        return a;
    end if;
    if degree(f, convert(y, set)) <= degree(expand(f - g), convert(y, set)) then
        return [];
    end if;
#### extract and classify the coefficients of f(y) - g(y + i) w . r . t y, which are polynomial in i###
    dl_sub, ht_sub := DegClassify([coeffs(expand(f - Shifts(g, y, a)), y)], a);
    sol_particular := [seq(i[j] = 1, j = 1 .. n)];
    eqs := {};
    for d_sub in dl_sub do
#### for each polynomial f in classification, subsubstitute particular solution into f-H^1(f) - H^0(f) in i###
        hl := ht_sub[d_sub];
        exprs := {};
        for expr in hl do
            tail := HomogeneousNoMoreOne(expr, a);
            expr := subs(sol_particular, expr - tail) + tail;
            exprs := {op(exprs), expr};
        end do;
#### solve the linear system ###
        eq := solve(exprs, convert(a, set));
        if output = 'Z' then 
            if eq = () or isolve(eq)=() then
                 return [];
            end if;
            eqs := solve(eqs union eq, convert(a, set));
            if eqs = () then
                return [];
            else
                 eq:=isolve(eqs);
            end if;
            if eq = () then
                return [];
            end if;    
            eq := subs(eq, a);
        elif output = 'F' then
            if eq = () then
                return [];
            end if;
            eqs := solve(eqs union eq, convert(a, set));
            if eqs = () then
                return [];
            end if;
            eq := subs(eqs, a);
        end if;
#### update the particular solution.###
        para := convert(indets(eq), list);
        k := nops(para);
        sol_particular := subs(zip((x, y) -> x = y, para, convert(LinearAlgebra:-ZeroVector(k), list)), eq);
        sol_particular := [seq(i[j] = sol_particular[j], j = 1 .. n)];
    end do;
    return eq;
end proc;
######################################################################
# Export
# Name : SET_ADC_derivation
# Calling sequence:
#                          SET_ADC_derivation(p, q, x)
#                          SET_ADC_derivation(p, q, x, output = 'type')
# Input: p, q,  two multivariate polynomials
#        x,  a variable name or a list of variable names
#        output,  (optional) should be a type of solutions in the form:
#                         output = 'Z'.
# Output: k,  a list representing the shift equivalence relation of p and q w.r.t. x;
#         [],  if there is no integer solutions.
#              If there is no condition on 'type', then return k as a general solution of
#                                     q(x) = p(x + a).
#              If 'type' = 'Z', then return k as a general integer solution.
# Remark: A general solution may contains some free variables.
######################################################################
SET_ADC_derivation := proc(p, q, x, {output := 'F'})
    local n, y, f, g, sub, d, dl, rootterms, llp, llp1, i, j, a, sol_particular, eq, eqs, lr, bool, level, tindex, lindex, term, dterm, df, df1, c, lt, lc, lt1, k, l, lineart, linearc, nlineart, lindex1, m, index, index1, c_l, c_nl, para, s, o, nlinearc, expr_nlinear, t, expr_linear;
    n := nops(x);
    y := convert(x, list);
    f := expand(p);
    g := expand(q);
    i := 'i';
    a := [seq(i[j], j = 1 .. n)];
    if f = 0 and g = 0 then
        return a;
    end if;
    sub := expand(f - g);
    d := degree(f, convert(y, set));
    if d <= degree(sub, convert(y, set)) then
        return [];
    end if;
#### extract and separate the level of the monomials of f(y) w.r.t y    ###
    dl, rootterms, llp := LevelPartition(f, y);
    sol_particular := [seq(i[j] = 0, j = 1 .. n)];
    eqs := {};
    lr := [seq([llp[j + 1][1], 0], j = 1 .. d)];#record the remain terms and its level as [linear,nlinear].
#### solve the equations with degree 0 ###
    lc := [coeffs(g, y, 'lt')];
    lt := [lt];
    g := 0;
    for j to numelems(lt) do
        term := lt[j];
        bool, level := IsLeaf(0, dl, rootterms, term, y);
        if not bool then
            g := g + lc[j]*term;
        else
            lr[level][1] := lr[level][1] - lc[j]*term;
        end if;
    end do;
    if llp[1][1] - g <> 0 then
        return [];
    end if;
#### solve the equations with degree 1,2,...d ###
    tindex := [seq([[]], j = 1 .. d + 1)];
    for k to d do
        expr_nlinear := 0;
        expr_linear := 0;
        lt := {};
        for t to k do
            lindex := tindex[t];
            lindex1 := [];
            df := llp[t];
            df1 := [];
####caculate \sum_{t_1+...t_n=k-t+1}(1/((k-t+1)!)a^(t_1)_1...a^(t_n)_n*\partial_(k-t+1)(L_(t-1)(p))) and avoid repeating ###  
            for l to numelems(lindex) do
#### take the partial derivative of the partial derivatives in df ###
                term := df[l];
                index := lindex[l];
                if k = t then
                    m := 1;
                else
                    m := index[k - t];
                end if;
                for m from m to n do
                    dterm := diff(term, y[m]);
                    if dterm <> 0 then
                        if t = 1 then
                            coeffs(expand(dterm), y, 'lt1');
                            lt := lt union {lt1};
                        end if;
                        index1 := [op(index), m];
                        lindex1 := [op(lindex1), index1];
                        df1 := [op(df1), dterm];
                        s := 1;
                        for o to k - t do
                            dterm := dterm*i[index1[o]]/s;
                            if index1[o + 1] = index1[o] then
                                s := s + 1;
                            else
                                s := 1;
                            end if;
                        end do;
                        dterm := dterm*i[index1[o]]/s;
                        if k = t then
                            expr_linear := expr_linear + dterm;
                        else
                            expr_nlinear := expr_nlinear + dterm;
                        end if;
                    end if;
                end do;
            end do;
            tindex[t] := lindex1;
            llp[t] := df1;
        end do;
        lt := convert(lt, list);
        expr_linear := expr_linear + lr[k][1];
        expr_nlinear := expr_nlinear + lr[k][2];
        linearc := [coeffs(expand(expr_linear), y, 'lineart')];
        lineart := [lineart];
        nlinearc := [coeffs(expand(expr_nlinear), y, 'nlineart')];
        nlineart := [nlineart];
        eq := {};
        for term in lt do
            bool, level := IsLeaf(k, dl, rootterms, term, y);
            c_l, linearc, lineart := SelectCoeff_ADC(linearc, lineart, term);
            c_nl, nlinearc, nlineart := SelectCoeff_ADC(nlinearc, nlineart, term);
            if not bool then
                eq := eq union {subs(sol_particular, c_nl) + c_l};
            else
                lr[level][1] := lr[level][1] + c_l*term;
                lr[level][2] := lr[level][2] + c_nl*term;
            end if;
        end do;
#### solve the linear system ###
        eq := solve(eq, convert(a, set));
       if output = 'Z' then 
            if eq = () or isolve(eq)=() then
                 return [];
            end if;
            eqs := solve(eqs union eq, convert(a, set));
            if eqs = () then
                return [];
            else
                 eq:=isolve(eqs);
            end if;
            if eq = () then
                return [];
            end if;    
            eq := subs(eq, a);
        elif output = 'F' then
            if eq = () then
                return [];
            end if;
            eqs := solve(eqs union eq, convert(a, set));
            if eqs = () then
                return [];
            end if;
            eq := subs(eqs, a);
        end if;
#### update the particular solution.###
        para := convert(indets(eq), list);
        s := nops(para);
        sol_particular := subs(zip((x, y) -> x = y, para, convert(LinearAlgebra:-ZeroVector(s), list)), eq);
        sol_particular := [seq(i[j] = sol_particular[j], j = 1 .. n)];
    end do;
    return eq;
end proc;
#################################################################
# Export
# Name: Shifts
# Calling sequence:
#                             Shifts(f, x, k)
# Input:  f,  a multivariate polynomial or rational function
#            x,  a variable name or a list of variable names
#            k,  an integer or a list of integers
# Output: g,  a rational function such that
#
#                             g(x) = f(x + a).
#
# Remark: If x = [x1, ..., xn] and a = [a1, ..., an] are two lists,
#                then the notation g(x) = f(x + a) means
#                              g(x1, ..., xn) = f(x1 + a1, ..., xn + an).
#################################################################
Shifts := proc(f, x, a) 
    local y, g; 
    y := zip((x, y) -> x = x + y, x, a); 
    g := subs(y, f); 
    g; 
end proc;

#################################################################
# Local
# Name: SelectCoeff
# Calling sequence:
#                                SelectCoeff(p, x, term)
# Input: p, a polynomial.
#            x, a variable name or a list of variable names.
#            term, a monomial.
# Output: c, an expression which is the coefficient of term.
#################################################################
SelectCoeff := proc(p, x, term) 
    local y, f, i, lc, lt; 
    y := convert(x, list); 
    f := expand(p); 
    lc := [coeffs(f, y, 'lt')]; 
    lt := [lt]; 
    for i to numelems(lc) do 
        if lt[i] = term then return lc[i]; end if; 
    end do; 
    return 0; 
end proc;
#################################################################
# Local
# Name: HomogeneousTerms
# Calling sequence:
#                                HomogeneousTerms(f, x)
# Input: f, a multivariate polynomial.
#            x, a variable name or a list of variable names.
# Output: ht,   a table such that
#                         ht = ([d_k = [H^d_k], ... , d_2 = [H^d_2], d_1 = [H^d_1]]);
#              dl, a list such that
#                          dl = [d_k, ...,d_2, d_1];
#              where H^d_i is homogeneous of degree d_i, satsfying
#                    f = H^(d_1)+H^(d_2)+...+H^(d_k) as well as d_k>...>d_2>d_1.   
#################################################################
HomogeneousTerms := proc(f, x)
    local g, y, lc, lm, tl, d, i, ht, dl;
    g := expand(f);
    y := convert(x, set);
#### when f is a monomial, operator "op" may not calculates correctly ###  
    lc, lm := Groebner:-LeadingTerm(g, plex(op(y)));
    if lc*lm - g = 0 then
        return table([degree(g, y)= [g]]), [degree(g, y)];
    end if;
    tl := [op(g)];
    ht := table();
    dl := [];
    for i in tl do
        d := degree(i, y);
        if d in dl then
            ht[d] := [ht[d][1] + i];
        else
            ht[d] := [i];
            dl := [op(dl), d];
        end if;
    end do;
    dl := sort(dl, `>`);
    ht, dl;
end proc;
#################################################################
# Local
# Name: DegClassify
# Calling sequence:
#                                DegClassify(l, x)
# Input: l,  a list of multivariate polynomials
#            x, a variable name or a list of variable names.
# Output: ht,   a table such that
#                          ht = ([d_1 = H^d_1, d_2 = H^d_2, ... , d_k = H^d_k])
#              dl, a list such that
#                          dl = [d_1, d_2,...,d_k];
#              where d_1<d_2<...<d_k;
#                          H^d_i is a sublist of l where its element f satsfying totaldegree(f) = d_i.
#################################################################
DegClassify := proc(l, x)
    local g, y, lc, lm, tl, d, i, ht, dl, hl;
    y := convert(x, set);
    ht := table();
    dl := [];
    for i in l do
        i := expand(i);
        d := degree(i, y);
        if d in dl then
            ht[d] := [op(ht[d]), i];
        else
            ht[d] := [i];
            dl := [op(dl), d];
        end if;
    end do;
    dl := sort(dl, `<`);
    dl, ht;
end proc;
#################################################################
# Local
# Name: HomogeneousNoMoreOne
# Calling sequence:
#                                HomogeneousNoMoreOne(f, x)
# Input: f,  a multivariate polynomial
#            x, a variable name or a list of variable names.
# Output: expr,   an expression such that
#                          expr = H` + H^0
#              where f=H`^(d)+H^(d-1)+...H^(0);
#              H^(i) is homogeneous of degree i. 
#################################################################
HomogeneousNoMoreOne := proc(f, x)
    local g, y, lc, lm, tl, d, i, expr;
    g := expand(f);
    y := convert(x, set);
#### when f is a monomial, operator "op" may not calculates correctly ###  
    lc, lm := Groebner:-LeadingTerm(g, plex(op(y)));
    if lc*lm - g = 0 then
        tl := [g];
    else
        tl := [op(g)];
    end if;
    expr := 0;
    for i in tl do
        d := degree(i, y);
        if d = 0 or d = 1 then
            expr := expr + i;
        end if;
    end do;
    expr;
end proc;
#################################################################
# Local
# Name: LevelPartition
# Calling sequence:
#                                LevelPartition(p, x)
# Input: p,  a multivariate polynomial
#            x, a variable name or a list of variable names.
# Output: ld,    a list such that          
#                          ld = [d_1, ...d_k]; 
#                       where d_1>d_2>...>d_k are the degree of the rootterms.
#               trt,   a table such that 
#                              trt = table([d_1=[x^(alpha_1),x^(beta_1),...],...,d_k=[x^(alpha_k),x^(beta_k),...]]); 
#                       where x^(alpha_i) is the rootterm of p(x) with degree d_i.
#               llp,    a list such that 
#                                llp = [[expr0], [expr1], ... [exprd]]
#                        where expri is the ith level of the p(x).    
#################################################################
LevelPartition := proc(p, x)
    local y, trt, ld, numd, llp, numlp, lc, lterms, l, tag, k, join_tag, index, d, rt, t, i, j, lrt, d_rt;
    y := convert(x, list);
    trt := table();
    ld := [];
    numd := 0;
    lc := [coeffs(sort(expand(p), y), y, 'lterms')];
    lterms := [lterms];
    d := degree(lterms[1], y);
    llp := [seq([0], i = 1 .. d + 1)];
    for l to numelems(lterms) do
        t := lterms[l];
        d := degree(t, y);
        join_tag := 0;
        index := 0;
        for i to numd do
            d_rt := ld[i];
            tag := 0;
            if d < d_rt then
                lrt := trt[d_rt];
                for j to numelems(lrt) do
                    rt := lrt[j];
                    if divide(rt, t) then
                        tag := 1;
                        join_tag := 1;
                        k := d_rt - d;
                        llp[k + 1][1] := llp[k + 1][1] + lc[l]*t;
                        break;
                    end if;
                end do;
            elif d = d_rt then
                index := i;
            end if;
            if tag = 1 then
                break;
            end if;
        end do;
        if join_tag = 0 then
            if index = 0 then
                ld := [op(ld), d];
                numd := numd + 1;
                trt[d] := [t];
            else
                trt[d] := [op(trt[d]), t];
            end if;
            llp[1][1] := llp[1][1] + lc[l]*t;
        end if;
    end do;
    ld, trt, llp;
end proc;
#################################################################
# Local
# Name: IsLeaf
# Calling sequence:
#                                IsLeaf(k,dl,rt,term,x)  
# Input: k,  a nature number which is the degree of the equations we are solving.
#           dl,  a list of the degree of the rootterms. 
#           rt,  a table of the rootterms getting from LevelPartition. 
#           term, a monomial which will be checked is at (d-k)th level or not,  
#               where d is the degree of the p(x).  
#           x,   a variable name or a list of variable names.  
# Output: bool,   a boolean value which is ture or false.  
#               k',     a nature number, the level of which the term is at.    
#################################################################
IsLeaf := proc(k, dl, rt, term, x)
    local y, d, i, rl, r;
    y := convert(x, list);
    d := degree(term, y);
    for i in dl do
        if d < i - k then
            rl := rt[i];
            for r in rl do
                if divide(r, term) then
                    return true, i - d;
                end if;
            end do;
        else
            break;
        end if;
    end do;
    return false, k;
end proc;
#################################################################
# Local
# Name: SelectCoeff_ADC
# Calling sequence:
#                                SelectCoeff_ADC(lc,lt,term) 
# Input: lc,  a list of coefficients.   
#            lt,  a list of monomials. 
#            term, a monomial.
# Output: c,   an expression which is the coefficient of term. 
#               lc,  the new list of coefficients.
#               lt,  the new list of monomials.  
#################################################################
SelectCoeff_ADC := proc(lc, lt, term)
    local c, i, lc1, lt1;
    c := 0;
    for i to numelems(lc) do
        if lt[i] = term then
            c := lc[i];
            lc1 := subsop(i = NULL, lc);
            lt1 := subsop(i = NULL, lt);
            return c, lc1, lt1;
        end if;
    end do;
    return c, lc, lt;
end proc;


end module;

