LRS := module()
option package;
export HypergeometricSolution,
       SelectIndicator,
       Resolving,
       UniversalDenominator,
       RationalSolution, CyclicVector;
local check_arguments, EG, Row_Length, Shift_Row,  ResolvingSequence,
       ResolvingStep, ByLinearSolve,
       ReductionStep,
       ByRationalSolution,
       select_not_similar;

# Input list/set form:
#    Sys  -- homogeneous 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 homogeneous 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, by_cyclic_vector;
global _c;

    st_total := time();
    try
        pp, i, si, by_cyclic_vector := check_arguments(args);
    catch:
        error;
    end try;

    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);
        return [];
    end if;

    m := nops(pp['functions']);
    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, ell := ByRationalSolution(pp, [], Hs);

    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, ell := ByRationalSolution(pp, ell, Hs);
            Res := [Res[], RatSols[]];
        end if;

    end do;

    userinfo(3, 'LRS', "all time:", time()-st_total);
    return Res;

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(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
    sol := LinearAlgebra:-LinearSolve(
                  <<B[2 .. -1, 1 .. StartRow-1]| B[2 .. -1, StartRow+1 .. -1]>>,
                                                         Vector(k-1), 'free' = z);
    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)
local st, ell, Ctfs, Res, i, j, j1, Ai, RatSols;
global _c;

    st := time():
    ell := select_not_similar(Ell, Hs, pp['var']);
    userinfo(4, 'LRS', "fill ell: ", time()-st);

    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(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;
            Ai := convert(Ai, 'list');
            RatSols := RationalSolution(Ai, pp['functions']);
        end if;

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

        st := time():
        Res := Res, map(`*`, RatSols, radnormal(ell[i]/Ctfs[i][2]))[];
        userinfo(4, 'LRS', "hypergeometric solutions: ", time()-st);
    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(A, pp['var']);
    userinfo(3, 'LRS', "reduction is evaluated:", time()-st);
    tmp[1]['functions'] := map(el-> el(pp['var']), indts);

    if si::'integer' or si::'name' then
        eq := eq = 0, procname(tmp[1], 1, 1)[];
    elif pp['input_type'] = 'Matrix' then
        si(A, pp['ring']);
        error "no tests";
    else
        error "no tests";
    end if;

    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 pp, m, m1, y, opts, S, M, vars, cfs, i, k, j, x, s_i, n,
      dummy, dot_indts, indts, cycl_vctr;
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({
             '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),
                 convert(args[1], 'set')) minus {0};
        S := [S[]];
        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);
        x := op(args[2][1]);
        m := nops(vars);
        m1 := nops(S);
        if nops({vars[]}) < m then
            error
           "The second argument (a list of different unknowns) is expected";
        end if;

        if {map(op, vars)[]} <> {x} then
            error
                 "The second argument (a list(function(%1))) is expected",
                                                                              x;
        end if;

        indts := select(el->has(x, el), y);
        if indts <> []  then
             error "%1 cannot be used for an unknown", indts[1](x);
        end if;

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

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

        if {eval(S, map(`=`, y, 0))[]} <> {0} then
            error "A homogeneous system is expected";
        end if;
        M := [seq(eval(S, map(`=`,subsop(i = NULL, y), 0)), i = 1 .. m)];
        try
            M := [seq(map(el ->
                     OreTools:-Converters:-FromLinearEquationToOrePoly(
                                           el+y[i](x)*dummy, y[i], pp['Ring']),
                          M[i]), i = 1 .. m)];
        catch:
            M := FAIL;
        end try;
        if has(M, FAIL) then
            error "A linear differential system is expected";
        end if;
        M := eval(M, dummy = 0);

        n := max(map(nops,map(op, M))) - 1;
        if n = 0 then
            error "A nonzero-order system is expected"
        end if;
        pp['order'] := n;

        if m > nops(S) then
            error "A full rank system is expected";
        end if;

        M := Matrix(nops(M[1]),m*(n+1),
                    {seq(seq(seq((k,(n+1-j+1)*m-(m-i))=op(j,M[i][k]),
                                 j = 1..nops(M[i][k])),
                             k = 1 .. nops(M[i])),
                         i = 1..m)});
        M := map(radnormal, M);
        if radnormal(expand({seq(S[k]-add(add(M[k][m*j-m+i]*
                  OreTools:-Apply(OrePoly(0$(n-j+1),1), vars[i], pp['Ring']),
                             j = 1..n+1),
                         i = 1..m),
                k = 1 .. nops(S))})) <> {0} then
            error "A linear system of %1 type is expected", op(2,pp['Ring']);
        end if;
        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 "Wrong a system's representation";
        end try;

        if S::'Matrix(square)' and LinearAlgebra:-RowDimension(S) = 0 then
            error "Non empty system is expected"
        end if;

        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(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
            M := convert(eval(S, y = 0), 'set');
        catch:
            error "Wrong a system's representation";
        end try;
        if M = {} then
            error "Non empty matrices are expected"
        elif M <> {0} then
            error "A homogeneous linear system with matrix coefficients is expected";
        end if;

        dot_indts := indets(S, `.`);
        dot_indts, indts := selectremove(el -> nops(el) = 2 and
                                               op(1,el)::'Matrix', dot_indts);
        if indts <> {} then
            error "Wrong a system's representation";
        end if;

        M := eval(S, {seq(dot_indts[i] = cfs[i]*op(2,dot_indts[i]),
                                                    i = 1 .. nops(dot_indts))});
        try
            M := OreTools:-Converters:-FromLinearEquationToOrePoly(
                                              M + dummy*y(x), y, pp['Ring']);
        catch:
            M := FAIL;
        end try;
        if has(M, FAIL) then
            error "A linear system of %1 type is expected", op(2,pp['Ring']);
        end if;
        M := eval(M, dummy = 0);
        n := nops(M) - 1;
        if n = 0 then
            error "A nonzero-order system is expected"
        end if;

        m := map(el-> [LinearAlgebra:-Dimension(op(1, el))], dot_indts);
        if nops(m) <> 1 then
            error "Wrong a system's representation";
        elif m[1][1] < m[1][2] then
            error "A full rank system is expected";
        end if;
        m1, m := m[1][];
        if m = 0 then
            error "Non empty matrices are expected"
        elif m = 1 then
            vars := [y(x)];
        else
            vars := [seq(y[i](x), i = 1..m)];
        end if;
        pp['functions'] := vars;
        pp['order'] := n;

        try
            dot_indts := map(el ->
               [op(1, el),
                OreTools:-Converters:-FromLinearEquationToOrePoly(
                                                     op(2, el), y, pp['Ring']),
                coeff(S,el)],
                             dot_indts);
        catch:
            error "Wrong a system's representation";
        end try;

        if select(el -> nops(el[2]) = 0 or
                        op(-1, el[2]) <> 1 or
                        not {op(1..-2, el[2])} in {{},{0}}, dot_indts) <> {} then
            error "Wrong a system's representation";
        end if;
        M := Matrix(m1,m*(n+1));
        for i in dot_indts do
            M[1..-1,m*(n+1)-m*nops(i[2])+1..m*(n+1)-m*nops(i[2])+m] :=
                 M[1..-1,m*(n+1)-m*nops(i[2])+1..m*(n+1)-m*nops(i[2])+m] +
                 i[1]*i[3];
        end do;
        M := map(radnormal, M);
        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');

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

    if has(cfs, y) then
        if y::'name' then
            error "%1 cannot be used for an unknown", y(x);
        else
           y := select(el->has(cfs, el), y);
           if nops(y) = 1 then
                error "The name %1 cannot be used for unknowns", y[];
           else
                error "The names %1 cannot be used for unknowns", y;
           end if;
        end if;
    end if;

    # check system's coefficients
    if not cfs::'set'('ratpoly'('anything',x)) then
        error
         "The system's coefficients of type ratpoly(anything, %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} then
            error "A full rank system is expected";
        end if;
        pp['embracing_system'] := M[1..m];
        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;

    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);
    hasoption(opts, 'select_indicator', 's_i', 'opts');
    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 pp['input_type'] = 'Matrix' then
        hasoption(opts, 'by_cyclic_vector' = truefalse, 'cycl_vctr', 'opts');
    end if;

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

    return pp, i, s_i, cycl_vctr;

end proc:
# RationalSolution
# Input list/set form:
#    Sys  -- homogeneous 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 homogeneous equation with matrix coefficients
#       A_n . y(x+n) + ... A_1 . y(x+1) + A_0 . y(x) = 0
#    y(x) -- unknown function
#    A  -- homogeneous system
#    vars -- list/set of unknown functions, [y[1](n),y[2](n)] for example
# 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;
    try
        pp := check_arguments(args)[1];
    catch:
        error;
    end try;
    if pp['input_type'] = 'Matrix' then
        RatSols :=
                 LinearFunctionalSystems:-RationalSolution(args[1], args[2], 'difference');
    else
        u := UniversalDenominator(args[1], args[2]);
        Syst := add(pp['explicit_form'][j].
            eval(u*Vector(pp['functions']), pp['var'] = pp['var'] + j),
                          j = 0 .. pp['order']);
        Syst := convert(Syst, 'list');
        PolySols := LinearFunctionalSystems:-PolynomialSolution(Syst, pp['functions']);
        RatSols := radnormal(map(`*`, PolySols, u));
    end if;

    indts := indets(RatSols, _c['integer']);
    return [seq(Vector(radnormal(eval(eval(RatSols, {j = 1}),
                                       {seq(j1 = 0,j1 = indts)}))), j = indts)];

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

    try
        pp := check_arguments(args)[1];
    catch:
        error;
    end try;

    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];
    V := lcm(convert(map(denom,pp['leading_invers']), 'set')[]);
    W := lcm(convert(map(denom, M[1..m,-m..-1]^(-1)), 'set')[]);
    return 1/`LREtools/DENOM`(V, W, pp['order'], pp['var']);
end proc:
end module: