############################################################################
##PROCEDURE(help,test=false,label="LiouvillianSolution") 
##      LiouvillianSolution
##HALFLINE find Liouvillian solutions of Linear (q-)Recurrence equations
##AUTHOR Khmelnov Denis
##DATE March, 2008
##AUTHOR for q-case Ryabenko Anna
##DATE for q-case December, 2012
##
##CALLINGSEQUENCE
##- LiouvillianSolution('reqn', 'fcn', `'indices'`=['n','k'], `'low'` = 'lb')
##
##PARAMETERS
##- 'reqn' : linear (q-)recurrence.
##- 'fcn' : function name, 'v(n)' for instance.
##- `'indices'`=['n','k'] : optional; where 'n', 'k'  are names.
##- `'low'`='lb' : optional; where 'lb' is a name
##
##RETURNS the Liouvillian solution of the given linear (q-)recurrence 'reqn'
##        or 0 if there is no such solution. 
##DESCRIPTION
##- The procedure finds the Liouvillian solution of the given linear 
##  (q-)recurrence 'reqn' in 'fcn' with the rational function coefficients
##  using the algorithm by Hendriks & Singer.
##- The Liouvilian solution is a generalization of the (q-)hypergeometric solution.
##  Let _H_ is the set of all (q-)hypergeometric sequences and _L_ is the smallest
##  subring of the ring _S_ of all sequences which contains _H_ and is closed under
##  (q-)shifts, summation and interlacing. The elements of _L_ are called 
##  Liouvillian sequences and 'reqn' has a Liouvillian solution if it 
##  has a nonzero solution in _L_.
##- The interlacing of _m_ sequences _a[1] = a[1,0], a[1,1],`...`_, 
##  _a[2] = a[2,0], a[2,1], `...`_,  _a[m] = a[m,0], a[m,1], `...`_ is the sequence
##  _a[1,0], a[2,0], `...`, a[m,0], a[1,1], a[2,1], `...`, a[m,1], `...`_. 
##- The solution is constructed by factoring the given equation, finding solutions for
##  the factors and combining the bases 'B[i]' of the factor's solutions in the solution of 
##  the given system. 
##- The procedure returns the solution that that is a linear combination of the 
##  independent solutions. The independent solutions will have indexed coefficients of the form 
##  _\_C[1], \_C[2],`...`,\_C[n]_, where _\_C_ is the global name.
##
##OPTIONS
##- `'indices'`=[n,k] Specifies base names for dummy variables. The default values are the global 
##  names _\_n_ and _\_k_, respectively. The _n1_, _n2_, etc. used as summation indices in the Liouvilian solution. 
##  The name _k_ is used as the product index.
##- If the option `'low'`='lb' is given, then 'lb' is used to assign a value
##  of the low bound of the index starting form which the found solution is valid.
##  
##ALGORITHM(nohelp)
##- The algorithm by Hendriks & Singer implemented in LiouvillianSolution reduces the task of
##  finding the Liovillian solution of the general form to one of finding the Liouvillian solutions
##  which are the interlacings of the (q-)hypergeometric solutions for the factors of
##  the given systems, which is in turn is reduced to one of finding the (q-)hypergeometric solutions
##  of a specially constructed equation for a factor.
##- Combination of the bases 'B[i]' is using Casaratian determinant. The implementation computes
##  it using the recurrence for the determinant, rather than computing it by the definition.
##  It allows computing more compact expressions in solutions. In addition the determinant's roots
##  are used for defining the lower summation bounds.
##
##TEST
## kernelopts(opaquemodules=false):
##
##EXAMPLES
##> read "/home/anna/Work/Maple_2013/LiouvillianSolution.mm":
##> rec := y(k)-(1+k-3*k^3)*y(k+2):
##> Res := LiouvillianSolution(rec, y(k));
##> qrec := -(a^2+b*x)*q*y(x)+(-a+a*q)*y(q*x)+y(q^2*x):
##> LiouvillianSolution(qrec, y(x));
##> qrec := (x+4)*(y(x*q^4)+q*y(x*q^3))-(x-1)*(y(x*q)+q*y(x)):
##> LiouvillianSolution(qrec, y(x));
##> rec := (n+1)*n*(n-1)*y(n+4)-n^2*(2*n+3)*(n-1)*y(n+3)+(n-1)*(n^4-8*n^2-9*n-1)*y(n+2)+(n+1)*(2*n^4+4*n^3-n^2-5*n-1)*y(n+1)-n^3*(n+1)*(n^2-2)*y(n):
##> LiouvillianSolution(rec,y(n));
##
##REFERENCES
##- P.A.Hendriks, M.F.Singer. Solving Difference Equations in Finite Terms.
##  J. Symb. Computation, 1999, 27 (3), pp. 239 - 259   
##- S.A.Abramov, M.A.Barkatou, D.E.Khmelnov. On m-Interlacing Solutions of Linear Difference Equations.
##    In: Computer Algebra in Scientific Computing, 11th International Workshop, CASC 2009, Kobe, Japan, 
##    September 2009, Proceedings, pp. 1-17
##
##SEEALSO
##- "LREtools[hypergeomsols]"
##- "LREtools"
##- "QDifferenceEquations[QHypergeometricSolution]"
##- "QDifferenceEquations"
##- "piecewise"
##############################################################################
LiouvillianSolution := module()
local ModuleApply, LS_doit, checkargs, 
      HypergeomSols, 
      EqCreate, Hyper_term, 
      isolveX, EvalHterm, 
      Ring, case, x, y, q, sum_index, prod_index, FirstNonzero,
      Shift_i,  make_hyper;

ModuleApply := proc(rc::algebraic, fn::anyfunc(name))
global _C;
local indts, N, rec, fns, reqn, info, sol, ndx, sbs, i, dummy, lb, n_0, n0_name;

    indts := indets(rc,'name');
    rec, fns := eval([rc, fn],{seq(indts[i]=N[i],i=1..nops(indts))})[];

    reqn := LREtools['REcreate'](rec, fns, {});
    info := op(-1,reqn);
    if not assigned(info['coeffs']) then
        reqn := QDifferenceEquations:-QECreate(rc, fn);
        case := 'qdifference';
    else
        reqn := LREtools['REcreate'](rc, fn, {});
        case := 'difference';
    end if;
    
    try n0_name := checkargs(reqn, args[3..-1]);
    catch:
        error;
    end try:
    unassign('_C');

    sol, lb := LS_doit(reqn, _C, 1, n_0);  

    if sol=FAIL then
        return(0);
    end if;

    # Renumber indetermined constants
    ndx := select(has,select(type,indets(sol),name),_C);
    sbs := seq(ndx[i]=_C[i],i=1..nops(ndx));
    sol := subs(sbs,sol);
    dummy := select(type,indets(sol),'Sum'(0,anything));
    if nops(dummy)>0 then
      sol := subs(seq(i=0,i=dummy),sol);
    end if;
    # Substitute the lower summation bound if any
    sol := eval(sol,n_0=lb);

    sol := subsindets(
            sol, 'specfunc'('anything', 'Hterm'),
            i -> Hyper_term(op(i)));
    if case = 'difference' then
        sol := eval(sol,sum_index=x);
    end if;
    if n0_name<>undefined then
       assign(n0_name,lb);
    end if;
    sol
end proc:
##############################################################################
LS_doit := proc(reqn, CC, ndx, n_0)
local tau_m, L_m, info, d, L, tau_n, m, H, Lm, L1, sol_1, L_1, N, lb_1, sol_2,
      lb_2, R, lb_C, lb, ndts1, ord1, basis1, ndts2, ord2, basis2, DC, a, l,
      s, i, M, k, rh, u, basis2in1, D, sol, B1, B2, quot, eq, r, sys, sbs; 

  tau_m := proc(m, max_shift, tau_n)
    option remember;
    local res;
           
    if m >= 0 and m<=max_shift-1 then
       res := Vector[column](max_shift);
       res[m+1] := 1;
    elif m = max_shift then            
      res := copy(tau_n);
    else
      res := LinearAlgebra:-Map(Shift_i, copy(tau_m(m-1,max_shift, tau_n)), 1);
      res := LinearAlgebra:-Map(normal,
                 Vector[column]([Vector[column]([0]),LinearAlgebra:-SubVector(res,[1..-2])])+
                 tau_n*res[-1]);
    end if;
    return res;
  end proc:

  L_m := proc(m, max_shift, tau_n)
    local sol, k, P, ns, i, A, lm, LL, N;    
    if m = 1 then
        sol, LL, N := HypergeomSols(reqn);
        LL := OreTools:-LCM['left'](OrePoly(1),seq(OrePoly(-k,1),k=LL),Ring);
        return([sol,LL,N]);
    else
       A := tau_m(0, max_shift, tau_n);
       k := 0; 
       ns := {};
       while nops(ns)=0 do
         k := k+1;
         A := <<A|tau_m(k*m, max_shift, tau_n)>>;
         ns := LinearAlgebra:-NullSpace(A);
       end do;
       P := OrePoly(seq(ns[1][i+1],i=0..k));
       P := OreTools:-Primitive(P);
       if case = 'difference' then
           lm := eval(P,x=m*x);
       else
           lm := eval(P,x=x^m);
       end if;
       lm := OreTools:-Converters:-FromOrePolyToLinearEquation(lm,y,Ring);
       sol, LL, N := HypergeomSols(EqCreate(lm,y(x)));
       if sol=0 then 
           return([0, OrePoly(1),0]); 
       end if;
       if case = 'difference' then
           LL := radnormal(eval(LL, x=x/m));
       else
           LL := radnormal(eval(LL, x=x^(1/m)));
       end if;
       sol := map(make_hyper,LL);
       N := max(seq(op(5,k),k=sol));
       LL := OreTools:-LCM['left'](seq(OrePoly(-k,0$(m-1),1),k=LL),Ring);
       [sol,LL,N]      
    end if;        
   end; 

   info := op(-1, reqn);
   d := info['order'];
   L := OrePoly(info['coeffs'][2..-1][]);
   tau_n :=Vector['column']([seq(-info['coeffs'][i+2]/info['coeffs'][-1], i=0..d-1)]);

   m := 0;
   H := OrePoly(1);
   while (m<d) and (H=OrePoly(1)) do
     m := m+1;
     Lm := L_m(m,d,tau_n);
     H := OreTools:-GCD['right'](L,Lm[2],Ring);
   end do;
   forget(tau_m);
   if H=OrePoly(1) then
     # GCRD (L_m, L) is trivial for all 1<=m<=d - no Liouvillian solution
     return(FAIL, 0);
   end if;
   L1 := OreTools:-Remainder['right'](L,H,Ring,'quot');
   L1 := [quot,L1];
   if m>1 then
      sol_1 := piecewise(seq(op([irem(sum_index,m)=k,
                   add(Hyper_term(op(EvalHterm(Lm[1][l],m,k)))*CC[ndx,k,l],
                                     l=1..nops(Lm[1]))]),
                              k=0..m-1));
   else
       sol_1:=add(Hyper_term(op(Lm[1][l]))*CC[ndx,1,l],l=1..nops(Lm[1]));
   end if;
   L_1 := OreTools:-Remainder['right'](Lm[2],H,Ring,'quot');
   L_1 := [quot, L_1];
   N := max(Lm[3], map(FirstNonzero, map(denom, convert(L_1[1], 'set'))),
                   map(FirstNonzero, map(denom, convert(Lm[2], 'set'))),
                   map(FirstNonzero, map(denom, convert(H,'set')))); 
   d := OreTools:-Utility:-Degree(H);
   r := OreTools:-Utility:-Degree(L_1[1]);
   if r > 0 then 
       if case = 'difference' then 
           eq := OreTools:-Apply(H, sol_1, Ring);
       else
           eq := add(eval(OreTools:-Utility:-Coefficient(H, i), x = q^sum_index)*
                     eval(sol_1, sum_index = sum_index+i), i = 0 .. d);
       end if;
       sys := {seq(simplify(eval(eval(eq,sum_index=max(N,0)+k),Product=mul)),k=0..r-1)};
       sbs := solve(sys, select(el-> op(0, el) = CC, indets(sys, indexed)));
       sol_1 := eval(sol_1, sbs);
   end if;


   lb_1 := N;
   sol_2, lb_2 := LS_doit(EqCreate(
          OreTools:-Converters:-FromOrePolyToLinearEquation(L1[1],y,Ring),y(x)), 
                            CC, ndx+1, n_0);
   R := (-1)^d*OreTools:-Utility:-TrailingCoefficient(H)/
                OreTools:-Utility:-LeadingCoefficient(H);
   lb_C := FirstNonzero(numer(R)*denom(R));
   lb := max(lb_2, lb_1, lb_C);
   ndts1 := select(type,select(has,indets(sol_1),CC),name); 
   ord1 := nops(ndts1);
   basis1 := [seq(B1[ndx][l](sum_index)=eval(sol_1,{seq(ndts1[k]=`if`(k=l,1,0),
                                                          k=1..ord1)}),l=1..ord1)];

   if sol_2=FAIL then 
      return(sol_1, lb_1) 
   end if;
   ndts2 := select(type,select(has,indets(sol_2),CC),name); 
   ord2 := nops(ndts2);
   basis2 := [seq(B2[ndx][l](sum_index)=eval(sol_2,{seq(ndts2[k]=`if`(k=l,1,0),k=1..ord2)}),
                    l=1..ord2)];

   R := (-1)^d*OreTools:-Utility:-TrailingCoefficient(H)/
                OreTools:-Utility:-LeadingCoefficient(H);
   DC := Hyper_term(op(make_hyper(R)));
   DC := eval(DC,sum_index=sum_index+1);
   s := {};
   for i to ord1 do
     M := Matrix([seq([seq(B1[ndx][k](sum_index+l),k=1..i-1),
                       `if`(l=ord1,b2(sum_index),0),
                       seq(B1[ndx][k](sum_index+l),k=i+1..ord1)],
                  l=1..ord1)]);
     s := s union {a[i] = LinearAlgebra:-Determinant(M)/D[ndx](sum_index)};
   end do;
     
    for k from 1 to ord1 do
       rh := normal(rhs(op(select(has,s,a[k]))));
       u[k]:= `if`(rh<>0,Sum(subs(sum_index=cat(sum_index,ndx),rh),
                            cat(sum_index,ndx)=n_0..sum_index-1),0);
    end do; 
    d := add(u[k]*B1[ndx][k](sum_index),k=1..ord1);
    basis2in1 := [seq(eval(d,{b2(sum_index)=B2[ndx][k](sum_index),
               b2(cat(sum_index,ndx))=B2[ndx][k](cat(sum_index,ndx))}),k=1..ord2)];

    #resolve
    D[ndx] := unapply(DC, sum_index);
    sol := add(CC[ndx,k,2]*basis2in1[k],k=1..ord2)+
           add(CC[ndx,k+ord2,2]*B1[ndx][k](sum_index),k=1..ord1);
    for k to ord1 do 
       B1[ndx][k] := unapply(rhs(basis1[k]),sum_index); 
    end do;
    for k to ord2 do 
       B2[ndx][k] := unapply(rhs(basis2[k]),sum_index); 
    end do;
    sol:=eval(sol);
    sol, lb;
end proc:
##############################################################################
checkargs := proc(reqn)
local a, nk, info, n0_name;

    info := op(-1,reqn);
    if not assigned(info['coeffs']) then
      error "only single linear equations (on a single variables) are handled";
    elif info['linear']=false then
      error "only linear equations are handled";
    elif nops([info['vars']])<>1 or not type(info['vars'],'name') then
      error "only equations in one variable are handled";
    elif remove(type,info['coeffs'][2..-1],
                                    'polynom'('anything',info['vars']))<>[] then
      error "only polynomial coefficients are handled";
    elif not Testzero(info['coeffs'][1]) then
      error "only homogeneous equations are handled";
    end if;

    a := [args[2..-1]];
    if not a::'list'('name'='anything') then
        error "wrong options"
    end if;

    n0_name := undefined;
    hasoption(a, 'low'='name', 'n0_name', 'a');
    if not hasoption(a,'indices'=['name','name'],'nk','a') then
        if case = 'qdifference' then
            sum_index := `tools/genglobal`('_n',indets([args],'name'));
        else
            sum_index := info['vars'];
        end if;
        prod_index := `tools/genglobal`('_k',indets([args],'name'));
    elif indets(nk) intersect indets(reqn) <> {} then
        error  "invalid right-hand side for option indices=%1", nk;
    else
        sum_index, prod_index := nk[];
    end if;

    if a <> [] then
        error "unknown options '%1'", a[];
    end if;   

   x := info['vars'];
   y := info['functions'];

   if case = 'difference' then 
       Ring := OreTools:-SetOreRing(x,'shift');
       Shift_i := proc(expr,i) eval(expr,x=x+i) end proc;
       FirstNonzero := proc(expr) `LREtools/firstnonzero`(expr,x) end proc;
   else
       q := info['q_par'];
       Ring := OreTools:-SetOreRing([x,q],'qshift');
       Shift_i := proc(expr,i) eval(expr,x=x*q^i) end proc;
       FirstNonzero := proc(expr) max(isolveX(expr,x,q),-1)+1 end proc;
   end if;

   return n0_name;
end proc:
##############################################################################
isolveX := proc(P, X, q)
local s, slv;

    s := sum_index;
    slv := {solve(eval(P, X = q^s), s)};
    slv := select(is, slv, 'integer');
    return slv[];
end proc:
##############################################################################
EqCreate := proc(rec, fns);

    if case = 'difference' then
        LREtools['REcreate'](rec, fns, {});
    else 
        QDifferenceEquations:-QECreate(rec, fns);
    end if;
end proc:
##############################################################################
HypergeomSols := proc(reqn)
local sol, LL, N, k, indts, ctf;
    if type(reqn,'specfunc'('anything','RESol'))  then 
        sol := LREtools['hypergeomsols'](reqn,'output'='basis');
        if sol=0 then 
            return(0, [], 0) 
        end if;
        LL := NULL;
        for k in sol do
            ctf := simplify(eval(k,x=x+1)/k);
            indts := indets(ctf,'specfunc'('anything','product'));
            indts := select((el)->op([2,2,2],el)=x,indts);
            ctf := eval(ctf, 
                   map(el->el=eval(op(1,el),op([2,1],el)=x)*
                        product(op(1,el),op([2,1],el)=op([2,2,1],el)..x-1),indts));
            LL := LL, ctf;
        end do;
        LL := [LL];
    else 
        LL := QDifferenceEquations:-QHypergeometricSolution(
                      op([1,1],reqn), op([2,1],reqn), 'output'='Certificate');
        if op(LL) = NULL then 
            return(0, OrePoly(1), 0) 
        else
            indts := indets(LL,'specindex'('integer',_C));
            if indts <> {} then 
                LL := [{seq(eval(eval(LL,k=1),map(el->el=0,indts))[],k=indts)}[]];
            else 
                LL := [LL[]];
            end if;
         end if;
      end if;
      sol := map(make_hyper, LL);
      N := max(seq(op(5,k),k=sol));
      return sol, LL, N;
end proc:
EvalHterm := proc(HT,m,N)
local n0, R;
      R := op(1,HT);
      n0 := FirstNonzero(eval(numer(R)*denom(R),sum_index=m*op(2,HT)+N));
     'Hterm'(R,sum_index,m,N,n0);
end proc:
make_hyper := proc(R)
local n0;
    n0 := FirstNonzero(numer(R)*denom(R));
    if case = 'qdifference' then
        'Hterm'(eval(R,x=q^sum_index),sum_index,1,0,n0);
    else
        'Hterm'(R,sum_index,1,0,n0);
    end if;
end proc:
Hyper_term := proc(ctf,n,m,N,K)
    local k, h;
        
    k := prod_index;
    h := simplify(Product(subs(n=m*k+N,ctf),k=K..(n-N)/m-1), 'Product');
    if case ='difference' then
        simplify(eval(h, Product=product),GAMMA);
    else
        subs(Product(1,k=K..n-1)=1, h);
    end if
end proc:
##############################################################################
end module: