# The third version, May 2018

LRS := module()
option package;
export HypergeometricSolution,
       SelectIndicator,
       Resolving,
       UniversalDenominator,
       RationalSolution, CyclicVector;
local  check_arguments, check_order, RHSdetermine,
       EG, Row_Length, Shift_Row,
       ResolvingSequence,
       ResolvingStep, ByLinearSolve,
       ReductionStep,
       ByRationalSolution,
       select_not_similar,
       HS, RS, UD, R, CV,
       max_isolve;

# Input list/set form:
#    Sys  -- difference system
#    vars -- list of unknown functions, [y1(x),y2(x)] for example
#   'select_indicator' = si -- (option) si is
#                              a name from vars or
#                              a procedure to select an unknown from vars,
#                              with arguments (system-in-set-form,
#                              list of vars)
# Output: a basis of hypergeometric solutions of the system Sys
# Input (matrix-form):
#    Sys  -- linear difference equation with matrix coefficients
#       A_n . y(x+n) + ... A_1 . y(x+1) + A_0 . y(x) = 0
#    y(x) -- unknown function
#    'select_indicator' = si -- (option) si is
#                             an integer or
#                             a procedure to select integer,
#                             with arguments (system-in-matrix-form, var)
# Input (normal form):
#    A  -- square matrix
#    x  -- independent variable
#    'select_indicator' = si -- (option) si is an integer,
#                                   or a procedure to select integer by
#                                      the given matrix and x
#    'by_cyclic_vector' = true/false -- (option) apply a cyclic vector method,
#                                       false by default
# Output: hypergeometric solutions of the system y(x+1) = A y(x)
#  where y is a vector of unknown functions
HypergeometricSolution := proc()
local y, pp, si, m, st, st_total, L, B, Hs, ell, t, Res, RatSols, i, j, by_cyclic_vector, part_res, prt_sol;
global _c;

    st_total := time();
    try
        pp, i, si, by_cyclic_vector, prt_sol := check_arguments[HS](args);
    catch:
        error;
    end try;
    m := nops(pp['functions']);
    part_res := [[], Vector(m)];
    ell := [];
    if not pp['homogeneous'] then
        try RHSdetermine(pp);
        catch:
            error;
        end try;
        for j to nops(pp['RHS_internel']) do
            Res := ByRationalSolution(pp, ell, [pp['RHS_internel'][j][1]], j);
            if [Res] = [] then
                return NULL;
            else
                Res, ell := Res;
            end if;
            part_res := [[part_res[1][], Res[1][]],
                         part_res[2] + Res[2]];
        end do;
    end if;
    if prt_sol then
        return part_res[2];
    end if;

    if by_cyclic_vector = true and pp['input_type'] = 'Matrix' then
        L, B := CyclicVector(pp['normal_form'][0], pp['var'], y);
    elif pp['input_type'] = 'Matrix' then
        L, B := Resolving(pp['normal_form'][0], pp['var'], y, 'select_indicator' = si);
    else
        st := time():
        L := ResolvingSequence(pp, i, si);
        userinfo(3, 'LRS', "ResolvingSequence is constructed:", time()-st);
    end if;

    st := time():
    Hs := LREtools['hypergeomsols'](
                 L[1], op(map2(op, 0, indets(L[1], 'function')) intersect
                                {y, map2(op, 0, pp['functions'])[]})(pp['var']),
                                                       {}, 'output' = 'basis');
    if Hs = 0 then
        Hs := [];
    end if;
    userinfo(3, 'LRS', nops(Hs),
                "dimension hypergeometric solution space is found:", time()-st);

    if nops(L) = 1 and Hs = [] then
        userinfo(3, 'LRS', "all time:", time()-st_total);
        if pp['homogeneous'] then
            return [];
        else
            return part_res;
        end if;
    end if;

    if nops(L) = 1 and pp['input_type'] = 'Matrix' and 'by_linear_solve' = true then
        Res := ByLinearSolve(pp['normal_form'][0], pp['var'], [], Hs, B, m)[1];
        userinfo(3, 'LRS', "all time:", time()-st_total);

        return Res;
    end if;

    Res := ByRationalSolution(pp, ell, Hs);
    if [Res] <> [] then
        Res, ell := Res;
    else
        Res := [];
    end if;

    for t from 2 to nops(L) while nops(Res) < m*pp['order'] do

        st := time():
        Hs := LREtools['hypergeomsols'](
                 L[t], op(map2(op, 0, indets(L[t], 'function')) intersect
                                {y, map2(op, 0, pp['functions'])[]})(pp['var']),
                                                       {}, 'output' = 'basis');
        if Hs = 0 then
            userinfo(3, 'LRS',
              "0 dimension hypergeometric solution space is found:", time()-st);
        else
            userinfo(3, 'LRS', nops(Hs),
                "dimension hypergeometric solution space is found:", time()-st);

            RatSols := ByRationalSolution(pp, ell, Hs);
            if [RatSols] <> NULL then
                RatSols, ell := RatSols;
                Res := [Res[], RatSols[]];
            end if;
        end if;

    end do;

    userinfo(3, 'LRS', "all time:", time()-st_total);
    if not pp['homogeneous'] then
       return [[Res[], part_res[1][]], part_res[2]];
    else
       return Res;
    end if;

end proc:
select_not_similar := proc(ell, HS, n)
local Hs, i;

    Hs := HS;
    for i to nops(ell) do
        Hs := remove(el -> SumTools:-Hypergeometric:-AreSimilar(el, ell[i], n),
                                                                            Hs);
    end do;
    for i while i < nops(Hs) do
        Hs := [Hs[1 .. i][],
               remove(el -> SumTools:-Hypergeometric:-AreSimilar(el, Hs[i], n),
                                                              Hs[i+1 .. -1])[]];
    end do;

    return Hs;

end proc:
Resolving := proc(AA::Matrix, x::name, y::name,
                 {select_indicator:=SelectIndicator}, $)
local A, m, mm, B, NS, t, d, i, k, L, l, st;

    try
        check_arguments[R](AA, x);
    catch:
        error;
    end try;

    m := LinearAlgebra:-RowDimension(AA);
    A := copy(AA);
    B := Matrix(m+1, m);
    d := 0;
    mm := m;
    for t while d < m do
        st := time();
        if select_indicator::integer and t = 1 then
            l := select_indicator;
        elif select_indicator::integer then
            l := 1;
        else
            l := select_indicator(A, x);
        end if;
        B := Matrix(1, mm, {(1,l) = 1});
        userinfo(3, 'LRS', "resolving for ", convert(B[1], 'list'),
                                                               ":", time()-st);

        st := time():
        NS := {};
        for i while NS = {} do
            B := <<B, map(radnormal,eval(B[-1], x = x+1).A)>>;
            NS := LinearAlgebra:-NullSpace(LinearAlgebra:-Transpose(B));
        end do;
        L[t] := collect(numer(radnormal(
                        Vector['row'](i, [seq(y(x+k), k = 0 .. i-1)]).(NS[1]))),
                                                   {seq(y(x+k), k = 0 .. i-1)});
        userinfo(3, 'LRS', i-1, "order resolving equation is constructed:",
                                                                     time()-st);

        d := d + i - 1;
        if d < m then
            st := time():
            A, mm := ReductionStep(B[1..-2], A, x, i-1, l);
            userinfo(3, 'LRS', "AA is constructed:", time()-st);
        end if;
    end do;

    return [seq(L[i], i = 1 .. t-1)], B[1..-2];
end proc:
ReductionStep := proc(B, A, n, k, StartRow)
local z, indts, funs, sol, sys, mm, AA, i, j, tmp;

    # set y[StartRow](n)=0 to  BY = 0 and solve this algebraic system
    if LinearAlgebra:-RowDimension(B) = 1 then
        sol := Vector(LinearAlgebra:-ColumnDimension(B)-1, 'symbol' = z);
    else
        sol := LinearAlgebra:-LinearSolve(
                  <<B[2 .. -1, 1 .. StartRow-1]| B[2 .. -1, StartRow+1 .. -1]>>,
                                                         Vector(k-1), 'free' = z);
    end if;
    indts := [indets(sol, z['integer'])[]];
    funs := map(el -> el(n), indts);
    sol := eval(sol, {seq(indts[i] = funs[i], i = 1 .. nops(indts))});

    sys := radnormal(convert(
         eval(sol, n = n+1) -
         << <<A[1..StartRow-1, 1..StartRow-1] |
                                     A[1..StartRow-1, StartRow+1..-1]>>,
            <<A[StartRow+1..-1, 1..StartRow-1]|
                                     A[StartRow+1..-1, StartRow+1..-1]>> >>.sol,
                                                              'set')) minus {0};
    sys := solve(sys, convert(eval(funs, n = n+1), 'set'));
    mm := nops(indts);
    AA := Matrix(mm, mm);
    for i to mm do
        tmp := select(has, sys, eval(funs[i], n = n+1))[1];
        sys := sys minus {tmp};
        tmp := collect(rhs(tmp), convert(funs, 'set'));
        AA[i] := Vector['row']([seq(coeff(tmp, funs[j]), j = 1 .. mm)]);
    end do;

    return AA, mm;
end proc:
ByLinearSolve := proc(A, n, Ell, Hs, B, m)
local st, Ctf, Res, j, j1, i, slv, Rhs;

    Res := NULL;
    Rhs := Matrix(m, nops(Hs));
    for j to nops(Hs) do
        RationalNormalForms:-IsHypergeometricTerm(Hs[j], n, 'Ctf');
        Rhs[1..-1,j] := Vector([seq(mul(eval(Ctf, n = n+j1), j1 = 0..i-1),
                                                                i = 0 .. m-1)]);
    end do;

    st := time();
    slv := LinearAlgebra:-LinearSolve(B, Rhs);
    userinfo(3, 'LRS', "LinearSolve: ", time()-st);
    for j to nops(Hs) do
        Res := Res,
               map(radnormal,
                    map(simplify, slv[1..-1,j]*Hs[j],
                                  {GAMMA, RootOf, 'power', 'symbolic'}));
    end do;

    return [Res], [];

end proc:
ByRationalSolution := proc(pp, Ell, Hs, N := 0)
local st, ell, Ctfs, Res, i, j, j1, Ai, RatSols;
global _c;
    if N = 0 then
        st := time():
        ell := select_not_similar(Ell, Hs, pp['var']);
        userinfo(4, 'LRS', "fill ell: ", time()-st);
    else
        ell := Hs;
    end if;

    st := time():
    Ctfs := map(
    proc(el)
       local ctf, z, a, b, c, d;
       if el::'ratpoly'('anything', pp['var']) then
           return [1, el];
       end if;
       RationalNormalForms:-IsHypergeometricTerm(el, pp['var'], 'ctf');
       z,a,b,c, d := RationalNormalForms:-RationalCanonicalForm(expand(ctf), pp['var']);
          return [z*a/b, c/d];
    end proc, ell);
    userinfo(4, 'LRS', "compute certificates: ", time()-st);

    Res := NULL;
    for i to nops(Ctfs) do

        if pp['input_type'] = 'Matrix' then
            Ai := map(radnormal, map(`*`, pp['normal_form'][0], 1/Ctfs[i][1]));
            st := time():
            RatSols := RationalSolution(Ai, pp['var']);
        else
            Ai := pp['explicit_form'][0].Vector(pp['functions']);
            j1 := 1;
            for j to pp['order'] do
                j1 := j1*eval(Ctfs[i][1], pp['var']=pp['var']+j-1);
                Ai := Ai + pp['explicit_form'][j]*j1.
                      Vector(eval(pp['functions'], pp['var']=pp['var']+j));
            end do;
            if N <> 0 then
                Ai := Ai-pp['RHS_internel'][N][4];
            end if;
            Ai := convert(Ai, 'list');
            try
                RatSols := RationalSolution(Ai, pp['functions']);
            catch:
                RatSols := NULL;
            end try;
        end if;

        if RatSols = NULL then
            userinfo(3, 'LRS', "there are no rational solutions: ", time()-st);
            return NULL;
        elif N = 0 and nops(RatSols) = 0 then
            userinfo(3, 'LRS', "there are no rational solutions: ", time()-st);
        elif N = 0 then
            userinfo(3, 'LRS', nops(RatSols), "dimension rational solution space is found: ", time()-st);
        end if;

        if N = 0 then
            st := time():
            Res := Res, map(`*`, RatSols, radnormal(ell[i]/Ctfs[i][2]))[];
            userinfo(4, 'LRS', "hypergeometric solutions: ", time()-st);
        else
            Res := map(`*`, RatSols[1], Hs[1]), Hs[1]*RatSols[2];
        end if;
    end do;

    return [Res], [Ell[], ell[]];

end proc:
SelectIndicator := proc(A::Matrix(square), n::name)
local m, num_zero, i, max_num_zero, Res, Res1, min_dgr, dgr;

    if not convert(A, 'set')::'set'(ratpoly(anything,n)) then
        error cat("rational function in ", n,
                                       " of data stored in Matrix is expected");
    end if;

    m := LinearAlgebra:-RowDimension(A);

    if m = 0 then
        error "non empty matrix is expected"
    end if;

    max_num_zero := 0;
    Res := NULL;
    for i to m do
        num_zero := nops(select(`=`, convert(A[i], 'list'), 0));
        if num_zero > max_num_zero then
            max_num_zero := num_zero;
            Res := i;
        elif num_zero = max_num_zero then
            Res := Res, i;
        end if;
    end do;

    if nops([Res]) = 1 then
        return Res;
    end if;

    Res1 := NULL;
    min_dgr := infinity;
    for i in [Res] do
        dgr := convert(map(el -> degree(numer(el), n)+degree(denom(el), n),
                                   remove(`=`, convert(A[i], 'list'), 0)), `+`);
        if min_dgr > dgr then
            min_dgr := dgr;
            Res1 := i;
        elif min_dgr = dgr then
            Res1 := Res1, i
        end if;
    end do;

    if nops([Res1]) = 1 then
        return Res1;
    else
        return Res1[rand(1 .. nops([Res1]))()]
    end if;

end proc:
CyclicVector := proc(AA, x, y)
local A, B, d, m, mm, t, i, NS, L, st, k;
global _c;

    m := LinearAlgebra:-RowDimension(AA);
    A := copy(AA);
    d := 0;

        st := time():
        B := Matrix(1, m);
        B[1] := Vector[row]([seq(RandomTools:-Generate(integer(range=-10..10)),
                                                                  i = 1 .. m)]);
        userinfo(3, 'LRS', "check if ", convert(B[1], 'list'), "is a cyclic vector");
        NS := {};
        for i while NS = {} do
            B := <<B, map(radnormal,eval(B[-1], x = x+1).A)>>;
            NS := LinearAlgebra:-NullSpace(LinearAlgebra:-Transpose(B));
        end do;
        L := collect(numer(radnormal(
                        Vector['row'](i, [seq(y(x+k), k = 0 .. i-1)]).(NS[1]))),
                                                   {seq(y(x+k), k = 0 .. i-1)});
        userinfo(3, 'LRS', i-1, "order resolving equation is constructed:",
                                                                     time()-st);
    return [L], B[1..-2];
end proc:
ResolvingSequence := proc(pp, i, si)
local eq, B, m, n, d, st, vars, k, j, sys, slv, indts, S, A, tmp, u;

    m := nops(pp['functions']);
    n := pp['order'];

    eq, B := ResolvingStep(pp, i, m, n);

    d := LinearAlgebra:-RowDimension(B);
    if d = n*m+1 then
        return [eq=0];
    end if;

    # reduction step
    # solve B(x) Y(x) = 0 with yi(x)=0, ξyi(x)=0, ξ^{n-1}yi(x) = 0
    st := time():

    B := LinearAlgebra:-DeleteColumn(B[n+1..-2], [seq(i+m*k,k=0..n-1)]);
    vars := map2(op, 0, subsop(i = NULL, pp['functions']));
    vars := [seq(map(el-> el[k], vars)[], k = 0 .. n-1)];

    sys := convert(B.Vector(vars), 'set');
    slv := solve(sys, {vars[]});
    indts := indets(map(rhs, slv));
    indts := select(has, vars, indts);

    slv := remove(evalb, slv);
    slv := [seq(eval(vars[k], slv), k = 1 .. nops(vars))];
    S := Matrix([seq([seq(coeff(slv[k], indts[j]), j = 1 .. nops(indts))],
        k = 1 .. nops(vars))]);

    A := Matrix(convert(pp['normal_form'], 'list'));
    A := LinearAlgebra:-DeleteColumn(
           LinearAlgebra:-DeleteRow(A, i), [seq(i+m*k,k=0..n-1)]);
    A := LinearAlgebra:-SubMatrix(A,
             map(proc(el) local k;
               `if`(member(el, vars[-m+1..-1], 'k'), k, NULL) end proc, indts),
             [1..-1]).S;
    A := map(radnormal, A);

    A := <<LinearAlgebra:-SubMatrix(S,
       map(proc(el) local k;
          `if`(member(el, vars[1..-m], 'k'), k+m-1, NULL) end proc, indts),
           [1..-1]), A>>;

    tmp := check_arguments[R](A, pp['var']);
    userinfo(3, 'LRS', "reduction is evaluated:", time()-st);
    tmp['functions'] := map(el-> el(pp['var']), indts);
    eq := eq = 0, procname(tmp, 1, 1)[];

    indts := indets([eq[2..-1]], 'function');
    indts := select(el-> op(0, el) in vars, indts);
    eq := eval([eq], map(el->op(0,el)=op([0,0],el), indts));
    if select(has, indets([eq[2..-1]]), vars[m..-1]) <> {} then
        error "no tests";
    end if;

    return eq;
end proc:
ResolvingStep := proc(pp, i, m, n)
local AA, k, B, slv, eq, st, C, indts;
eval(pp);
    userinfo(3, 'LRS', "resolving for ", [0$(i-1), 1, 0$(m-i)]);

    AA := pp['normal_form'];

    st := time():

    B := Matrix(n,m*n,{seq((k+1,m*k+i)=1,k=0..n-1)});

    B := <<B, Matrix(1, m*n, [seq(AA[k][i], k=0..n-1)])>>;

    slv := LinearAlgebra:-LinearSolve(B^%T, Vector(m*n), 'free' = C);

    while convert(slv, 'set') = {0} do

       B := <<B,
              map(expand@radnormal, Matrix(1, m*n,
                 [map(eval, B[-1,-m .. -1], pp['var'] = pp['var']+1).AA[0],
        seq(
          map(eval, B[-1,(k-1)*m + 1 .. k*m], pp['var']=pp['var']+1)+
          map(eval,  B[-1,-m .. -1],  pp['var']=pp['var']+1).AA[k],
        k = 1..n-1)])) >>;
        slv := LinearAlgebra:-LinearSolve(B^%T, Vector(m*n), 'free' = C);

    end do;

    userinfo(3, 'LRS', LinearAlgebra:-RowDimension(slv)-1,
                              "order resolving equation is constructed:",
                                                                     time()-st);

    indts := [indets(slv, C['integer'])[], C[0]];
    slv := eval(slv, {indts[1] = 1, seq(k = 0, k = indts[2..-1])});
    slv := map(expand@radnormal,
                    lcm(convert(map(denom, map(radnormal, slv)), 'set')[])*slv);

    eq :=  OreTools:-Converters:-FromOrePolyToLinearEquation(
               OrePoly(convert(slv, 'list')[]),
               op(0,pp['functions'][i]),
               pp['Ring']);
    return eq, B;

end proc:
EG := proc(ExplicitMatrix, n, m, rdim, x)
local k, v, c, l_c, l_c1, R, den;

    if convert(ExplicitMatrix[-1], 'set') = {0} then
        ExplicitMatrix[1..-2] :=
                            procname(ExplicitMatrix[1..-2], n, m, rdim-1, x);
        return ExplicitMatrix;
    end if;

    for k to rdim  do
        ExplicitMatrix[k] := Shift_Row(ExplicitMatrix[k], m, x);
    end do;

    v := LinearAlgebra:-NullSpace(ExplicitMatrix[1..-1, 1..m]^%T);
    if v = {} then
        return ExplicitMatrix;
    end if;

    l_c := 0;
    for k to rdim do
        l_c1 := Row_Length(ExplicitMatrix[k]);
        if v[1][k] <> 0 and l_c < l_c1 then
            c := k;
            l_c := l_c1;
        end if;
    end do;

    R := map(radnormal,(v[1]^%T).ExplicitMatrix);
    den := lcm(map(denom,convert(R, 'set'))[]);
    R := map(radnormal, den*R);
    if convert(R, 'set') = {0} then

        if c < rdim then
            ExplicitMatrix[c] := ExplicitMatrix[-1];
        end if;

        ExplicitMatrix[-1] := R;
        ExplicitMatrix[1..-2] :=
                            procname(ExplicitMatrix[1..-2], n, m, rdim-1, x);
        return ExplicitMatrix;
    else
        ExplicitMatrix[c] := R;
        return procname(ExplicitMatrix, n, m, rdim, x);
    end if;
end proc:

Shift_Row := proc(C, m, x)
local l, c, den;

    c := C;
    l := Row_Length(c);
    while l <> 0 and convert(c[1..m], 'set') = {0} do
        den := c[l];
        c := map(radnormal, c/den);
        c := <<map(eval, c, x = x+1)[m+1..-1] |
               Vector[row](m)>>[1];
        c := map(radnormal, c);
        den := lcm(map(denom,convert(c, 'set'))[]);
        c := map(radnormal, den*c);
        l := Row_Length(c);
    end do;

    return c;

end proc:
Row_Length := proc(c)
local k;
    for k from LinearAlgebra:-ColumnDimension(c) by -1 to 1 do
        if c[k] <> 0 then
            return k;
        end if;
    end do;
    return 0;
end proc:

check_arguments := proc()
 local m, x, vars, M, pp, cfs, opts, S, y, m1, n, b, indts, term, dot_indts, k,
       i, j, s_i, cycl_vctr, prt_sol;
option remember;

    if nargs = 0 then
        error "The first argument (a system) is expected";
    end if;

    if args[1]::'Matrix' then
        if not args[1]::'Matrix(square)' then
            error "A square matrix is expected"
        end if;

        m := LinearAlgebra:-RowDimension(args[1]);
        if m = 0 then
            error "Non empty matrix is expected"
        end if;

        if nargs = 1 then
            error "The second argument (an independent variable) is expected"
        elif not args[2]::'name' then
            error "The second argument (a name) is expected"
        end if;
        x := args[2];

        if m = 1 then
            vars := [y(x)];
        else
            vars := [seq(y[i](x), i = 1..m)];
        end if;

        M := map(radnormal, args[1]);
        pp := table({
             'homogeneous' = true,
             'RHS' = Vector(m),
             'embracing_system' = <<Matrix(m, 'shape' = 'identity') | -M>>,
             'normal_form' = Array(0..0, [M]),
             'explicit_form' = Array(0..1, [-M, Matrix(m, 'shape' = 'identity')]),
             'order' = 1,
             'var' = x,
             'functions' = vars,
             'input_type' = 'Matrix',
             'leading_invers' = LinearAlgebra:-IdentityMatrix(m),
             'Ring' = OreTools:-SetOreRing(x, 'shift')});
        cfs := convert(args[1], 'set');

        opts := [args[3..-1]];

    elif args[1]::{'set', 'list'}({'algebraic', '`=`'}) then

        S := map(el->`if`(el::`=`, (lhs-rhs)(el), el), args[1]);
        if args[1]::'set' then
            S := convert(S minus {0}, 'list');
        end if;
        if nops(S) = 0 then
            error "Non empty system is expected"
        end if;

        if nargs = 1 then
            error "The second argument (a list of unknowns) is expected"
        elif not args[2]::'list'('function'('name')) then
            error "The second argument (a list of function(name)) is expected"
        end if;
        vars := args[2];
        y := map2(op, 0, vars);
        m := nops(vars);
        if m = 0 then
            error "Non empty list of unknowns is expected";
        end if;
        m1 := nops(S);
        if nops({vars[]}) < m then
            error
           "The second argument (a list of different unknowns) is expected";
        end if;

        x := {map(op, vars)[]};
        if nops(x) <> 1 then
            error
             "The second argument (a list of function(%1)) is expected", x[1];
        end if;
        x := x[1];
        if has(y, x) then
            error "%1 cannot be used for an unknown",
                               (select(has, y, x)[1])(x);
        elif has(x, y) then
            error "%1 cannot be used for an unknown",
                               (select(el -> has(x, el), y)[1])(x);
        end if;

        pp := table({
             'var' = x,
             'functions' = vars,
             'input_type' = 'list',
             'Ring' = OreTools:-SetOreRing(x, 'shift')});

        opts := [args[3..-1]];

        try
            n := check_order(S, y, x);
        catch:
            error;
        end try;
        pp['order'] := n;

        b := eval(S, map(`=`, y, 0));
        pp['homogeneous'] := evalb(b = [0$m1]);
        pp['RHS'] := <-b>;
        indts := indets(S, 'specfunc'({y[]}));
        S := collect(expand(S), indts, radnormal);
        if select(el -> not el::'linear'(indts), S) <> [] then
            error "A linear system is expected";
        end if;

        M := Matrix([seq([seq(seq(
                                        coeff(S[i], y[k](x+n-j), 1),
                                                 k = 1 .. m), j = 0 .. n)],
                                                                i = 1 .. m1)]);
        pp['explicit_form'] := Array(0..n,
                 [seq(M[1..-1,(n+1)*m-(m*(k+1)-1)..(n+1)*m-m*k],
                                                                   k = 0..n)]);
        cfs := convert(M, 'set');

    elif args[1]::{'algebraic',`=`} then
        try
            S := `if`(args[1]::`=`, (lhs-rhs)(args[1]), args[1]);
        catch:
            error
         "The first argument (a system in matrix or list/set form) is expected";
        end try;

        if nargs = 1 then
            error "The second argument (an unknown) is expected"
        elif not args[2]::'function'('name') then
            error "The second argument (a function(name)) is expected"
        end if;
        y := op(0, args[2]);
        x := op(args[2]);
        if has(y, x) or has(x, y) then
            error "%1 cannot be used for an unknown", y(x);
        end if;

        pp := table({
             'var' = x,
             'input_type' = 'list',
             'Ring' = OreTools:-SetOreRing(x, 'shift')});

        opts := [args[3..-1]];

        try
            n := check_order(S, y, x);
        catch:
            error;
        end try;
        pp['order'] := n;

        indts := indets(S, 'Matrix');
        if indts = {} then
            m := {[1, 1]};
        else
            m := map(el -> [LinearAlgebra:-Dimension(el)], indts);
        end if;
        if nops(m) <> 1 or m[1][1] = 0 or m[1][2] = 0 then
            error "System's coefficients have to be non-empty matrices with equal dimensions";
        end if;
        m := m[1];

        S := convert(S + dummy, list);
        pp[explicit_form] := Array(0 .. n,
                                   [seq(Matrix(m[1], m[2]), k = 0 .. n)]);
        pp['RHS'] := Vector(m[1]);
        pp['homogeneous'] := true;
        for term in S do
            if term = dummy then
                next;
            elif not term::dependent(y) then
                if term::'Vector' and LinearAlgebra:-Dimension(term) = m[1] then
                    pp['RHS'] := pp['RHS'] - term;
                    pp['homogeneous'] := pp['homogeneous'] and
                                         evalb(convert(term, 'set') = {0});
                elif m[1] = 1 then
                    pp['RHS'] := pp['RHS'] - <term>;
                    pp['homogeneous'] := false;
                else
                    error "The system's term %1 is wrong, a Vector with %2 rows is expected", term, m[1];
                end if;
                next;
            end if;
            dot_indts := indets(term, `.`);
            dot_indts := select(has, dot_indts, y);
            if dot_indts = {} then
                indts := indets(term, 'specfunc'(y));
                if nops(indts) <> 1 or
                   not(term = indts[1] or
                       term::`*` and select(`=`, term, indts[1]) = indts[1]) or
                   m[1] <> m[2] then
                    error "The system's term %1 is wrong", term;
                end if;
                indts := indts[1];
                k := op(1, indts) - x;
                pp['explicit_form'][k] :=
                      pp['explicit_form'][k] +
                      coeff(term, indts, 1)*LinearAlgebra:-IdentityMatrix(m[1]);
                next;
            elif nops(dot_indts) > 1 then
                error "The system's term %1 is wrong", term;
            end if;
            dot_indts := dot_indts[1];
            if nops(dot_indts) <> 2 then
                error "The `.`-terms of two arguments is expected, but %1 has %2 arguments",
                      dot_indts, nops(dot_indts);
            elif not op(2, dot_indts)::'specfunc'(y) then
                error "The `.`-terms with the second argument of type 'specfunc'(%1) is expected, but %2 is given",
                      y, op(2, dot_indts);
            end if;
            k := op([2,1], dot_indts)-x;
            if op(1, dot_indts)::Matrix then
                pp['explicit_form'][k] :=
                                     pp['explicit_form'][k] +
                                     coeff(term, dot_indts, 1)*op(1, dot_indts);
            elif m[1] = m[2] then
                pp['explicit_form'][k] :=
                           pp['explicit_form'][k] +
                           coeff(term, dot_indts, 1)*op(1, dot_indts)*
                                            LinearAlgebra:-IdentityMatrix(m[1]);
            else
                error "The system's term %1 is wrong", term;
            end if;
        end do;
        M := Matrix(m[1],m[2]*(n+1));
        for i from 0 to n do
            M[1..-1, i*m[2]+1 .. (i+1)*m[2]] := pp['explicit_form'][n-i];
        end do;
        if m[2] = 1 then
            pp['functions'] := [y(x)];
        else
            pp['functions'] := [seq(y[i](x), i = 1..m[2])];
        end if;
        cfs := convert(M, 'set');
        m1, m := m[1], m[2];

    else
        error
         "The first argument (a system in matrix or list/set form) is expected";
    end if;

    indts := select(el -> has(el, y) or has(y, indets(el, 'name')), cfs);
    if indts <> {} then
        error "The system has a wrong coefficient: %1", indts[1];
    end if;
    # check system's coefficients
    if not cfs::'set'('ratpoly'('anything',x)) then
        error
             "The system's coefficients of type ratpoly in %1 are expected", x;
    elif op(procname) = HS and not cfs::'set'('ratpoly'('rational',x)) then
        error
             "The system's coefficients of type ratpoly(rational(%1)) are expected", x;
    elif not cfs::'set'('ratpoly'('algnum',x)) then
        error
             "The system's coefficients of type ratpoly(algnum(%1)) are expected", x;
    end if;

    if not assigned(pp['embracing_system']) then
        M := EG(M, n, m, m1, pp['var']);
        if m < m1 and convert(M[m+1..-1], 'set') <> {0} then
            error "It's impossible"
        end if;
        if convert(M[m], 'set') = {0} and pp['homogeneous'] then
            error "A full rank system is expected";
        else
            pp['embracing_system'] := M[1..m];
            if convert(M[m], 'set') <> {0} then
                pp['leading_invers'] := (M[1 .. m,1 .. m])^(-1);
                pp['normal_form'] := Array(0..n-1,
                    [seq(copy(map(radnormal,
                             -pp['leading_invers'].
                             M[1 .. m, (n+1)*m-(m*(k+1)-1)..(n+1)*m-m*k])),
                         k=0..n-1)]);
            end if;
        end if;
    end if;

    if not opts::'list'(`=`) then
        error
    "optional arguments are expected to be of type equation, but received '%1'",
                                                      remove(type, opts, `=`)[];
    end if;

    s_i := `if`(args[2]::'list', vars[1], 1);
    if op(procname) = HS then
        hasoption(opts, 'select_indicator', 's_i', 'opts');
    end if;
    if not s_i::'posint' and not s_i::'function'('name') then
        i := s_i(args[1], args[2]);
    else
        i := s_i;
    end if;
    if i::'posint' then
        if args[2]::'set' then
            error "'select_indicator' is expected a function from %1", vars;
        elif i > m then
            error "The option argument 'select_indicator' is expected <= %1", m;
        end if;
    elif i::'function'('name') then
        if not member(i, vars, 'i') then
            error "'select_indicator' is expected a function from %1", vars;
        end if;
    end if;

    cycl_vctr := false;
    if op(procname) = HS and pp['input_type'] = 'Matrix' then
        hasoption(opts, 'by_cyclic_vector' = truefalse, 'cycl_vctr', 'opts');
    end if;

    prt_sol := false;
    if op(procname) = HS and not pp['homogeneous'] and
                              hasoption(opts, 'output', 'prt_sol', 'opts') then
        if prt_sol <> 'partsol' then
            error "optional argument 'output' is expected to be 'partsol'";
        end if;
        prt_sol := true;
    end if;

    if opts <> [] then
        error "Invalid optional arguments; first one is '%1'", opts[1];
    end if;

    if op(procname) = HS then
        return pp, i, s_i, cycl_vctr, prt_sol;
    else
        return pp;
    end if;

end proc:
check_order := proc(S, y, x)
local shift_indts, indt, n;

    shift_indts := indets(S, 'specfunc'(y));
    if shift_indts = {} then
        error "A linear system for %1 is expected", map(el->el(x), y);
    end if;
    indt := select(el -> nops(el) <> 1 or
                         not (op(el)-x)::'nonnegint', shift_indts);
    if nops(indt) = 1 then
        error "The system has a wrong element %1", indt[1];
    elif nops(indt) > 1 then
        error "The system has wrong elements %1", indt;
    end if;
    n := max(map(el -> op(el) - x, shift_indts));
    return n;

end proc:
RHSdetermine := proc(pp)
local m1, m, x, dummy, term, t1, R, Rhs, RHs, k0, i, j;
    if assigned(pp['RHS_k0']) then return end if;

    x := pp['var'];
    m1 := LinearAlgebra:-RowDimension(pp['RHS']);
    Rhs := Array(1..m1, 'fill'=[]);
    k0 := 0;
    for i to m1 do
        if pp['RHS'][i] = 0 then next end if;
        for term in pp['RHS'][i] + dummy do
            if term = dummy then next; end if;
            try
                t1 := SumTools:-Hypergeometric:-IsHypergeometricTerm(
                                                                 term, x, 'R');
            catch:
                error "The right-hand side of sum of hypergeometric terms is expected but %1 is given", term;
            end try;
            if not t1 then
                t1 := expand(term);
                if t1 = term then
                    error "The right-hand side of sum of hypergeometric terms is expected but %1 is given", term;
                end if;
                pp['RHS'][i] := eval(pp['RHS'][i], term = t1);
                try procname(pp);
                    return;
                catch:
                    error "The right-hand side of sum of hypergeometric terms is expected but %1 is given", term;
                end try;
            end if;
            k0 := max(k0, max_isolve(numer(R),x)+1, max_isolve(denom(R),x)+1);
            R := SumTools:-Hypergeometric:-RationalCanonicalForm[1](
                                                                 expand(R), x);
            R := R[1], R[2]/R[3], R[4]/R[5];
            Rhs[i] := [Rhs[i][], Array([term, R[], 0])];
        end do;
    end do;
    RHs := [];
    for i to m1 do
        for term in Rhs[i] do
            try
                term[5] := value(eval(term[1]/term[4], x = k0));
            catch:
                error "The right-hand side has a term %1 which  raises `%2`, for large enough %3", term[1], lastexception[2], x;
            end try;
            term[5] := radnormal(term[5]);
            if term[5] = 0 then next end if;
            for j to nops(RHs) while term[5] <> 0 do
                if RHs[j][2] <> term[2] then
                    next;
                elif RHs[j][3] = term[3] then
                    RHs[j][4][i] := RHs[j][4][i] + term[4]*term[5];
                    break;
                end if;
                R := SumTools:-Hypergeometric:-RationalCanonicalForm[1](
                                                 expand(term[3]/RHs[j][3]), x);
                R := R[1], radnormal(R[2]/R[3]), R[4]/R[5];
                if R[2] = 1 then
                    RHs[j][4][i] := RHs[j][4][i] +
                                    term[4]*term[5]*R[3]/eval(R[3], x = k0);
                    break;
                end if;
            end do;
            if j > nops(RHs) then
                RHs := [RHs[], [radnormal(term[1]/term[5]/term[4]),
                                    term[2], term[3],
                                             Vector(m1, {i=term[5]*term[4]})]];
            end if;
        end do;
    end do;
    pp['RHS'] := add(j[1].j[4], j = RHs);
    pp['RHS_k0'] := k0;
    pp['RHS_internel'] := RHs;
    return;
end proc:
# RationalSolution
# Input list/set form:
#    Sys  -- difference system
#    vars -- list of unknown functions, [y1(x),y2(x)] for example
# Output: a basis of rational solutions of the system Sys
# Input (matrix-form):
#    Sys  -- linear equation with matrix coefficients
#       A_n . y(x+n) + ... A_1 . y(x+1) + A_0 . y(x) = 0
#    y(x) -- unknown function
# Output: a basis of rational solutions, a universal denominator for the system A.
# Input (normal form):
#    A  -- square matrix
#    x  -- independent variable
# Output: a basis of rational solutions, a universal denominator for the system
#        y(x+1) = A y(x)
#  where y(x) is a vector of unknown functions
RationalSolution := proc()
local pp, u, Syst, j, j1, RatSols, PolySols, indts, cfs;
    try
        pp := check_arguments[RS](args);
    catch:
        error;
    end try;
    if pp['input_type'] = 'Matrix' then
        RatSols :=
                 LinearFunctionalSystems:-RationalSolution(args[1], args[2], 'difference');
    else
        if not pp['homogeneous'] then
            cfs := convert(pp['RHS'], 'set');
            if not cfs::'set'('ratpoly'('anything', pp['var'])) then
                return NULL;
            elif not cfs::'set'('ratpoly'('algnum', pp['var'])) then
                error
                   "The system's right-hand sides of type ratpoly(algnum(%1)) are expected", x;
            end if;
        end if;
        try
            u := UniversalDenominator(args[1], args[2]);
        catch:
            return NULL;
        end try;
        if u = FAIL then
            return NULL;
        end if;
        Syst := add(pp['explicit_form'][j].
            eval(u*Vector(pp['functions']), pp['var'] = pp['var'] + j),
                          j = 0 .. pp['order'])-pp['RHS'];
        Syst := convert(Syst, 'list');
        try
            PolySols := LinearFunctionalSystems:-PolynomialSolution(
                                                        Syst, pp['functions']);
        catch:
            return NULL;
        end try;
        if PolySols = FAIL then
            return NULL;
        end if;
        RatSols := radnormal(map(`*`, PolySols, u));
    end if;

    if RatSols = [] then
        return NULL;
    end if;
    indts := indets(RatSols, _c['integer']);

    RatSols := [[seq(Vector(radnormal(eval(eval(RatSols, {j = 1}),
                                       {seq(j1 = 0,j1 = indts)}))), j = indts)],
             <radnormal(eval(RatSols,{seq(j1 = 0, j1 = indts)}))>];
    if convert(RatSols[2], 'set') = {0} then
        return RatSols[1];
    else
        RatSols := [map(`+`, RatSols[1], -RatSols[2]), RatSols[2]];
        return RatSols;
    end if;

end proc:
UniversalDenominator := proc()
local pp, m, i, j, M, V, W, Syst, cfs, newargs;

    try
        pp := check_arguments[UD](args);
    catch:
        error;
    end try;
    if not pp['homogeneous'] then
        cfs := convert(pp['RHS'], 'set');
        if not cfs::'set'('ratpoly'('anything', pp['var'])) then
            error "The right-hand side of type ratpoly in %1 is expected",
                                                                    pp['var'];
        elif not cfs::'set'('ratpoly'('algnum', pp['var'])) then
            error "The system's right-hand sides of type ratpoly(algnum(%1)) are expected", x;
        end if;
        Syst := add(pp['explicit_form'][j].
            eval(Vector(pp['functions']), pp['var'] = pp['var'] + j),
                          j = 0 .. pp['order'])-pp['RHS'];
        Syst := convert(Syst, 'list');
        try
            newargs := LinearFunctionalSystems:-HomogeneousSystem(
                                                'homo', Syst, pp['functions']);
            return procname(newargs);
        catch:
            return FAIL;
        end try;
    end if;

    m := nops(pp['functions']);
    M := LinearFunctionalSystems:-MatrixTriangularization(
            Matrix([seq(pp['explicit_form'][pp['order']-i], i=0..pp['order'])]),
                                          pp['order']+1, pp['var'], 'trail')[1];
    try
        V := lcm(convert(map(denom,pp['leading_invers']), 'set')[]);
        W := lcm(convert(map(denom, M[1..m,-m..-1]^(-1)), 'set')[]);
    catch:
        lasterror;
        return FAIL;
    end try;
    return 1/`LREtools/DENOM`(V, W, pp['order'], pp['var']);
end proc:

###############################################################
# Name: max_isolve                                            #
# Input:                                                      #
#   p -- a polynomial in Q(a)[x] where a is an                #
#        algebraic number represented by RootOf               #
#   x -- a variable name                                      #
# Output: max nonnegative integer root of p                   #
###############################################################

max_isolve := proc(P,x)
local p, indts, g, sols, k;

    p := convert(P,'RootOf');
    indts := indets(p, {'RootOf','name'}) minus {x};
    if nops(indts)<=1 and p::'polynom'('radalgnum', x) then
        g := content(p,indts);

        sols := isolve(g);
        if sols = {} then
            return -1;
        else
            sols := seq(op([1,2],k), k=[sols]);
        end if;

        return max(-1, sols);
    elif type(r, 'polynom'('radalgnum', x)) then
        sols := roots(p, x);
    else
        sols := frontend(roots, [p,x], ['{`+`,`*`,`=`,`..`,radalgnum}',{}]);
    end if;
    return max(-1, op(select(type,map2(op,1,sols),'integer')));
end proc:

end module: