RS := module()
export EG, ResolvingSequence, Sdetermine, ResolvingEquation, ResolvingMatrix, Indicator,
       ResolvingDependence, ResolvingSystems, CyclicVector;
option package;
local check_arguments, EG_doit, Row_Length, Shift_Row, ResolvingStep, doit,
      indicator, Bs, Ds, Ss, ReduceNormalForm;

# ResolvingSequence
# Input list/set form:
#    Sys  -- homogeneous system
#    vars -- list/set of unknown functions, [y1(x),y2(x)] for example
#    Ring -- Ore polynomial ring, given by OreTools:-SetOreRing
#   'select_indicator' = si -- (option) si is a name from vars,
#                                   or a procedure to select an unknown with arguments
#                                       (system-in-set-form, list of vars, Ring)
# Output: a resolving sequence of equations for the system Sys
# Input (matrix-form):
#    Sys  -- linear homogeneous equation with matrix coefficients
#       A_n . ξ^n(y(x),x) + ... A_1 . ξ(y(x),x) + A_0 . y(x) = 0
#    y(x) -- unknown function
#    Ring -- Ore polynomial ring, given by OreTools:-SetOreRing
#    'select_indicator' = si -- (option) si is an integer,
#                                   or a procedure to select integer with arguments
#                                       (system-in-matrix-form, var, Ring)
# Output: a resolving sequence of equations for the system Sys
# Input (normal-form):
#    A  -- square matrix
#    y(x) -- unknown function (to construct y[1](x), y[2](x),...)
#    Ring -- Ore polynomial ring, given by OreTools:-SetOreRing
#    'select_indicator' = si -- (option) si is an integer,
#                                   or a procedure to select integer with arguments
#                                       (matrix, var, Ring)
# Output: a resolving sequence of equations for the system ξy=Ay

ResolvingSequence := proc()
local pp, si, s_f, st, i, result;

    st := time():
    pp, si, s_f := check_arguments(args);

    if s_f <> false then
        error "Invalid optional arguments; first one is '%1'", 'start_from' = s_f;
    elif si::'name' then
        i := si;
    elif args[1]::'Matrix' then
        i := si(args[1], args[3]);
    else
        i := si(args[1], args[2], args[3]);
    end if;

    Bs := NULL;
    Ds := NULL;
    Ss := pp['source_system'];
    result := doit(pp, i, si);

    userinfo(3, ResolvingSequence, "all time:", time()-st);

    indicator := map2(op, 2, result);
    if op(0,args[1]) = 'OrePoly' then
        map2(op, 1, result);
    else
        map(res -> OreTools:-Converters:-FromOrePolyToLinearEquation(
                         res[1], op(0,pp['functions'][res[2]]), pp['Ring'])=0,
            result);
    end if;

end proc:
doit := proc(pp, i, si)
local eq, B, m, n, d, st, vars, k, j, sys, slv, indts, S, A, tmp, u, sbst;

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

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

    d := LinearAlgebra:-RowDimension(B);
    if m > 1 then
        if pp['source_system']::{'set','list'} then
            Bs := Bs, ('`.`')(copy(B[1..-2]),
               Vector([seq(map(el->
                          OreTools:-Apply(OrePoly(0$k,1), el, pp['Ring']),
                     pp['functions'])[], k = 0 .. n-1)]));
        else
            Bs := Bs, copy(B[1..-2]);
        end if;
    end if;

    if d = n*m+1 then
        return [[eq,i]];
    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 := subsop(i = NULL, pp['functions']);

    sbst := [seq(map(el-> op(0,el)[k] =
                          OreTools:-Apply(OrePoly(0$k,1), el, pp['Ring']),
                     vars)[], k = 0 .. n-1)];
    vars := [seq(map(el-> op(0,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))];

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

    if pp['source_system']::'Matrix' then
        S := <S[1..i-1], Matrix(1,LinearAlgebra:-ColumnDimension(S)),S[i..-1]>;
        sys := <sys[1 .. i-1], Vector[row](m-1), sys[i .. -1]>;
        sys := <sys[1 .. -1, 1 .. i-1]| Vector(m)| sys[1 .. -1, i .. -1]>;
        Ds := Ds, sys;
        tmp := ReduceNormalForm(pp, S, eval(indts, sbst));
    else
        if pp['source_system']::{'list','set'} then
            Ds := Ds, Vector(pp['functions']) =
                               <S[1..i-1],0,S[i..-1]>.eval(Vector(indts), sbst);
        else
            Ds := Ds, sys;
        end if;
        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, u(pp['var']), pp['Ring']);
    end if;
    tmp[1]['functions'] := map(el-> el(pp['var']), indts);
    userinfo(3, ResolvingSequence, "reduction is evaluated:", time()-st);
    Ss := Ss, tmp[1]['source_system'];

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

    for k from 2 to nops(eq) do
        d := op([0,1],tmp[1]['functions'][eq[k][2]]);
        member(op([0,0],tmp[1]['functions'][eq[k][2]])(pp['var']),
                   pp['functions'],
                   'j');
        eq[k] := [OrePoly(0$d, op(eq[k][1])), j];
    end do;

    return eq;
end proc:
ReduceNormalForm := proc(pp, S, vars)
local A, u, new_pp;

    A := pp['source_system'];

    A := LinearAlgebra:-SubMatrix(A,
             map(proc(el) local k;
               `if`(member(el, pp['functions'], 'k'), k, NULL) end proc, vars),
             [1..-1]);
    A := map(radnormal, A.S);

    new_pp := check_arguments(A, u(pp['var']), pp['Ring']);

    eval(new_pp);
end proc:
ResolvingStep := proc(pp, i, m, n)
local AA, k, B, slv, st, C, indts;

    AA := pp['normal_form'];

    st := time():

    if i = "random" then
        B := Matrix(1,m*n,{seq((1,k)=RandomTools:-Generate(integer(range=-10..10)),
                                                                  k = 1 .. m)});
        userinfo(3, ResolvingSequence, "resolving for ", convert(B[1,1..m], 'list'));
    elif i::'list' then
        B := convert(i, Matrix);
        userinfo(3, ResolvingSequence, "resolving for ", convert(B[1,1..m], 'list'));
    else
        userinfo(3, ResolvingSequence, "resolving for ", [0$(i-1), 1, 0$(m-i)]);
        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)])>>;
    end if;

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

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

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

    end do;

    userinfo(3, ResolvingSequence, 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);

    return OrePoly(convert(slv, 'list')[]), B;

end proc:
EG_doit := proc(ExplicitMatrix, n, m, rdim, Ring)
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, Ring);
        return ExplicitMatrix;
    end if;

    for k to rdim  do
        ExplicitMatrix[k] := Shift_Row(ExplicitMatrix[k], m, Ring);
    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, Ring);
        return ExplicitMatrix;
    else
        ExplicitMatrix[c] := R;
        return procname(ExplicitMatrix, n, m, rdim, Ring);
    end if;
end proc:

Shift_Row := proc(C, m, Ring)
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(OreTools:-Properties:-GetSigma(Ring), c,
                       OreTools:-Properties:-GetVariable(Ring))[m+1..-1]|
               Vector[row](m)>>[1] +
             map(OreTools:-Properties:-Getdelta(Ring), c,
                                OreTools:-Properties:-GetVariable(Ring));
        c := map(radnormal, c);
        l := Row_Length(c);
    end do;
    den := lcm(map(denom,convert(c, 'set'))[]);
    c := map(radnormal, den*c);

    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, q, s_i, n, inverseA_l,
      dummy, dot_indts, indts, tmp, s_f;
option remember;

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

    if op(0,args[1]) = 'OrePoly' then

        if nargs = 1  or
           not args[2]::'UnivariateOreRing' then
               error
                   "The second argument (of type UnivariateOreRing) is expected";
        end if;
        x := OreTools:-Properties:-GetVariable(args[2]);
        q := attributes(args[2]):-qpar;

        pp := table({
             'source_system' = args[1],
             'var' = x,
             'Ring' = args[2]});

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

        S := [op(args[1])];
        n := nops(S) - 1;
        if n = 0 then
            error "A nonzero-order system is expected"
        end if;

        M := select(type, S, 'Matrix');
        if M = {} then
            m1 := 1;
            m := 1;
        else
            m := {map(el->[LinearAlgebra:-Dimension(el)], M)[]};
            if nops(m) <> 1 then
                error "The system's coefficients are expected with the equal size";
            elif m[1][1] < m[1][2] then
                error "A full rank system is expected";
            end if;
            m1, m := m[1][];
        end if;

        if m = 0 or m1 = 0 then
            error "Non empty matrices are expected"
        end if;

        S := map(`.`, S, LinearAlgebra:-IdentityMatrix(m1,m));

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

        S := map(radnormal, S);
        pp['explicit_form'] := Array(0..n, [seq(S[k+1], k = 0..n)]);
        M := S[n+1];
        for k from n to 1 by -1 do
            M := <<M | S[k]>>;
        end do;
        cfs := convert(M, 'set');

    elif 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 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 nargs = 2 or
           not args[3]::'UnivariateOreRing' then
               error "The third argument (of type UnivariateOreRing) is expected";
        end if;

        if x <> OreTools:-Properties:-GetVariable(args[3]) then
            error "The second argument (function(%1)) is expected",
                                     OreTools:-Properties:-GetVariable(args[3]);
        end if;
        q := attributes(args[3]):-qpar;
        if has({q, x}, y) then
            error "%1 cannot be used for an unknown", y(x);
        end if;


        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({
             'source_system' = 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,
             'Ring' = args[3]});

        pp['embracing_system'] := map(radnormal@expand,
             Matrix(m, m,
                {seq((j, j) = denom(radnormal(add(M[j, i]*y^i, i = 1 .. m))),
                     j = 1 .. m)}).
             pp['embracing_system']);

        cfs := convert(args[1], 'set');

        opts := [args[4..-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/set of unknowns) is expected"
        elif not args[2]::{'list','set'}('function'('name')) then
            error "The second argument (a list/set of function(name)) is expected"
        end if;
        vars := [args[2][]];
        y := map2(op, 0, vars);
        m := nops(vars);
        m1 := nops(S);
        if nops({vars[]}) < m then
            error
           "The second argument (a list/set of different unknowns) is expected";
        end if;

        if nargs = 2 or
           not args[3]::'UnivariateOreRing' then
               error
                   "The third argument (of type UnivariateOreRing) is expected";
        end if;

        x := OreTools:-Properties:-GetVariable(args[3]);
        q := attributes(args[3]):-qpar;
        if {map(op, vars)[]} <> {x} then
            error
                 "The second argument (a {list,set}(function(%1))) is expected",
                                                                              x;
        end if;

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

        pp := table({
             'source_system' = args[1],
             'var' = x,
             'functions' = vars,
             'Ring' = args[3]});

        opts := [args[4..-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 system of %1 type is expected", op(2,pp['Ring']);
        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]);
        if nargs = 2 or
           not args[3]::'UnivariateOreRing' then
               error
                   "The third argument (of type UnivariateOreRing) is expected";
        end if;
        x := OreTools:-Properties:-GetVariable(args[3]);
        q := attributes(args[3]):-qpar;
        if op(args[2]) <> x then
            error "The second argument (a function(%1)) is expected", x;
        end if;
        if has({q, x}, y) then
            error "%1 cannot be used for an unknown", y(x);
        end if;

        pp := table({
             'source_system' = args[1],
             'var' = x,
             'Ring' = args[3]});

        opts := [args[4..-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 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;

        M := eval(S, map(`=`, dot_indts, 0));
        if M <> 0 then
            try
                M := OreTools:-Converters:-FromLinearEquationToOrePoly(
                                                            M, y, pp['Ring']);
            catch:
                error "Wrong a system's representation";
            end try;
            tmp := map(`.`, M, Matrix(m1, m, 'shape' = 'identity'));
            M := op(1, tmp);
            for i from 2 to nops(tmp) do
                M := <op(i, tmp)| M>;
            end do;
            for i from nops(tmp) to n do
                M := <Matrix(m1, m)| M>;
            end do;
        else
            M := Matrix(m1,m*(n+1));
        end if;
        try
            dot_indts := map(el -> `if`(op(2, el) = y(x),
               [op(1, el), OrePoly(1), coeff(S,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;
        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;

    q := attributes(pp['Ring']):-qpar;

    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_doit(M, n, m, m1, pp['Ring']);

        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];
        inverseA_l := (M[1 .. m,1 .. m])^(-1);
        pp['normal_form'] := Array(0..n-1,
            [seq(copy(map(radnormal,
                   -inverseA_l.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]::'set', vars[1], 1);
    hasoption(opts, 'select_indicator', 's_i', 'opts');
    if s_i::'posint' then
        if args[2]::'set' then
            error
        "the option argument 'select_indicator' is expected a function from %1",
                                                                           vars;
        elif s_i > m then
            error "The option argument 'select_indicator' is expected <= %1", m;
        end if;
    elif s_i::'function'('name') then
        if not member(s_i, pp['functions'], 's_i') then
            error
        "the option argument 'select_indicator' is expected a function from %1",
                                                                           vars;
        end if;
    end if;

    s_f := false;
    if hasoption(opts, 'start_from'=list, 's_f', 'opts') then
        if nops(s_f) <> m then
            error "the option argument 'start_from' is expected a list of %1 elements",
                                                                           m;
        end if;
    end if;


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

    return pp, s_i, s_f;

end proc:
Sdetermine := proc();
    eval(check_arguments(args)[1]);
end proc:
EG := proc()
local pp, i, res, m;
    pp := check_arguments(args)[1];

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

    if op(0,args[1]) = 'OrePoly' then
        OrePoly(seq( pp['embracing_system'][1..-1,
                                      m*(pp['order']-i)+1..m*(pp['order']-i+1)],
                     i = 0 .. m));
    elif args[2]::'function'('name') then
        add(`if`(convert(pp['embracing_system'][1..-1, i*m+1..m*(i+1)], 'set') =
                     {0},
                 0,
                 pp['embracing_system'][1..-1, i*m+1..m*(i+1)].
                 OreTools:-Converters:-FromOrePolyToLinearEquation(
                       OrePoly(0$pp['order']-i, 1), op(0,args[2]), pp['Ring'])),
            i = 0 .. pp['order']) = 0;
    else
        res := add(pp['embracing_system'][1..-1, i*m+1..m*(i+1)].
                   Vector(map2(OreTools:-Converters:-FromOrePolyToLinearEquation,
                            OrePoly(0$(pp['order']-i),1),
                               map2(op, 0, pp['functions']), pp['Ring'])),
                   i = 0 .. pp['order']);
        if args[1]::'set' then
            map(`=`, convert(res, 'set'), 0);
        else # args[1]::'list' then
            map(`=`, convert(res, 'list'), 0);
        end if;
    end if;
end proc:
CyclicVector := proc()
local pp, res, u, s_f;

    pp := check_arguments(args);
    pp, s_f := pp[1], pp[3];

    if s_f = false then
        res := ResolvingStep(pp, "random", nops(pp['functions']), pp['order']);
    else
        res := ResolvingStep(pp, s_f, nops(pp['functions']), pp['order']);
    end if;
    Bs := copy(res[2][1..-2]);
    indicator := NULL;

    if op(0,args[1]) = 'OrePoly' then
        res[1];
    else
        OreTools:-Converters:-FromOrePolyToLinearEquation(
               res[1], u, pp['Ring'])=0;
    end if;
end proc:
ResolvingEquation := proc()
local pp, si, s_f, m, i, k, res;

    pp, si, s_f := check_arguments(args);

    if s_f <> false then
        si := s_f;
    elif si::'name' then
        i := si;
    elif args[1]::'Matrix' then
        i := si(args[1], args[3]);
    else
        i := si(args[1], args[2], args[3]);
    end if;

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

    res := ResolvingStep(pp, i, m, pp['order']);
    if m > 1 then
        if pp['source_system']::{'set','list'} then
            Bs := ('`.`')(copy(res[2][1..-2][1..-2]),
               Vector([seq(map(el->
                          OreTools:-Apply(OrePoly(0$k,1), el, pp['Ring']),
                     pp['functions'])[], k = 0 .. pp['order']-1)]));
        else
            Bs := copy(res[2][1..-2]);
        end if;
    end if;

    Bs := copy(res[2][1..-2]);
    indicator := [i];

    if op(0,args[1]) = 'OrePoly' then
        res[1];
    else
        OreTools:-Converters:-FromOrePolyToLinearEquation(
               res[1],
               op(0,pp['functions'][i]),
               pp['Ring'])=0;
    end if;
end proc:
ResolvingMatrix := proc() map(copy, [Bs]) end proc:
ResolvingDependence := proc() [Ds] end proc:
ResolvingSystems := proc() [Ss] end proc:
Indicator := proc() indicator end proc:

end module:

with(RS):