LRS := module()
export HypergeometricSolution,
       Resolving,
       SelectIndicator,
       ByLinearSolve,
       ReductionStep,
       ByRationalSolution,
       ByMatrixInverse,
       CyclicVector;
local  select_not_similar;

# Input: A -- square matrix of rational functions in n
#        n -- independent variable
# Output: hypergeometric solutions of the system
#        Y(n+1) = A Y(n)
#  where Y is a vector of unknown functions

HypergeometricSolution := proc(A::Matrix(square),
                               n::name,
                               {select_indicator:=SelectIndicator},
                               {the_case_of_cyclic_vector:=ByLinearSolve},
                               {hypersols:=false},
                               {by_cyclic_vector:=false},
                               $)
local m, st, st_total, L, B, Hs, ell, t, Ctfs, Res, RatSols, AA;
global _c;

    st_total := time();
    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;

    if select_indicator::integer and select_indicator > m then
        error cat("the option argument 'select_indicator' is expected <= ", m)
    end if;

    if by_cyclic_vector = true then
        L, B := CyclicVector(A, n, y);
    else
        L, B := Resolving(A, n, y, 'indicator' = select_indicator);
    end if;

    if hypersols = false then
        st := time():
        Hs := LREtools['hypergeomsols'](L[1], y(n), {}, output = basis);
        if Hs = 0 then
            Hs := [];
        end if;
        userinfo(3, LRS, nops(Hs), "dimension hypergeometric solution space is found:", time()-st);
    else
        Hs := hypersols;
    end if;

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

    if nops(L) = 1 then
        Res := the_case_of_cyclic_vector(A, n, [], Hs, B, m)[1];
        userinfo(3, LRS, "all time:", time()-st_total);

        return Res;
    end if;

    Res, ell := ByRationalSolution(A, n, [], Hs);

    for t from 2 to nops(L) while nops(Res) < m and hypersols = false do

        st := time():
        Hs := LREtools['hypergeomsols'](L[t], y(n), {}, 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(A, n, 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(square), x, y, {indicator:=1}, $)
local A, m, mm, B, NS, t, d, i, k, L, l, st;

    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 indicator::integer and t = 1 then
            l := indicator;
        elif indicator::integer then
            l := 1;
        else
            l := 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, "Ã 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, Ã, 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);
    if convert(sol, 'set') = {0} then
        return [Res];
    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);
    Ã := 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'));
        Ã[i] := Vector['row']([seq(coeff(tmp, funs[j]), j = 1 .. mm)]);
    end do;

    return Ã, 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(A, n, Ell, Hs)
local st, ell, Ctfs, Res, i, j, j1, Ai, RatSols, indts;
global _c;

    st := time():
    ell := select_not_similar(Ell, Hs, n);
    userinfo(4, LRS, "fill ell: ", time()-st);

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

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

        Ai := map(radnormal, map(`*`, A, 1/Ctfs[i][1]));

        st := time():
        RatSols :=
                 LinearFunctionalSystems:-RationalSolution(Ai, n, 'difference');
        indts := indets(RatSols, _c['integer']);
        if nops(indts) = 0 then
            userinfo(3, LRS, "there are no rational solutions: ", time()-st);
        else
            userinfo(3, LRS, nops(indts), "dimension rational solution space are found: ", time()-st);
        end if;

        st := time():
        RatSols := map(`*`, RatSols, radnormal(ell[i]/Ctfs[i][2]));
        Res := Res, seq(Vector(radnormal(eval(eval(RatSols, {j = 1}),
                                        {seq(j1 = 0,j1 = indts)}))), j = indts);
        userinfo(4, LRS, "hypergeometric solutions: ", time()-st);
    end do;

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

end proc:
ByMatrixInverse := proc(A, n, Ell, Hs, Bb, m)
local st, B, Res, Ctf, Rhs, i, j, j1;

    st := time():
    B := Bb^(-1):
    userinfo(3, LRS, "B^(-1) is computed:", time()-st);

    Res := NULL;
    for j to nops(Hs) do
        RationalNormalForms:-IsHypergeometricTerm(Hs[j], n, 'Ctf');
        Rhs := Vector([seq(mul(eval(Ctf, n = n+j1), j1 = 0..i-1),
                                                                i = 0 .. m-1)]);
        st := time();
        Res := Res,
               map(radnormal, map(simplify,
                               Hs[j]*(B.Rhs), {GAMMA, RootOf, power, symbolic}));
        userinfo(3, LRS, "compute solutions: ", time()-st);
    end do;

    return [Res], [];

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, st_total;
global _c;

    m := LinearAlgebra:-RowDimension(AA);
    A := copy(AA);
    d := 0;
    mm := m;
    for t while d < m do

        st := time():
        B := Matrix(1, mm);
        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[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, "tilde M is constructed:", time()-st);
        end if;
    end do;

    return [seq(L[i], i = 1 .. t-1)], B[1..-2];
end proc:

end module: