######################################################################
# INPUT :
# eq - a homogeneous linear ordinary differential equation with 
#      polynomial coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x)
# fn - the unknown function in the form y(x)
# C - name of arbitrary constants (optional argument) (not _C),
#     _B by default 
# 
# OUTPUT :
#  formal power series solution with hypergeometric coefficients;
#  points of these series are chose automatically;
#  solutions may involve arbitrary constants in form _C[i], C[j]
#  (where i,j are some integer).
#
`Slode/hypercoeffs` := proc( eq :: {algebraic, alpgebraic=algebraic},
			fn::anyfunc(name),
			CC::name)
local C, res, Eq, p_r, sng, ons, x, x0 ;
global _B;

	C := `if`( nargs > 2, CC, _B );
	Eq := `Slode/DEdetermine`(eq, fn);

	if not Eq['linear'] then
		ERROR(`The differential equation must be linear on`, Eq['fun'] );
	fi;
	if not Eq['polycfs'] then
		ERROR(`The linear differential equation must be `, 
                      `with polynomial coefficients on`, Eq['var']);
	fi;
	if Eq['rhs'] <> 0 then
		ERROR( `eq must be a homogeneous equation` );
	fi;

	x := Eq['var'];

	p_r := factor(Eq['coeffs'][Eq['ord']]);
	p_r := `if`( type(p_r,`*`), convert(p_r,set), {p_r} );

	sng := map(RootOf,select(has,p_r,x),x);
	ons := map(`+`,sng, -1) minus sng;

	if ons <> {} then 
		res := {`Slode/hyper_at_point`( eq, fn, C, ons[1])};
	elif sng = {} then
		res := {`Slode/hyper_at_point`( eq, fn, C )};
	fi;

	for x0 in sng do
		res := res union {`Slode/hyper_at_point`(eq,fn,C,x0)};
	od;

	res minus {0};

end:
######################################################################
# INPUT :
# eq - a homogeneous linear ordinary differential equation with 
#      polynomial coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x)
# fn - the unknown function in the form y(x)
# C - name of arbitrary constants (optional argument) (not _C)
# x0 - point (optional argument)
# 
`Slode/hyper_at_point` := proc( eq :: {algebraic, algebraic=algebraic},
			fn::anyfunc(name),
			CC::name,
			x0::algebraic )
local xx0, C, v, n, Ser, sln, Rec;
global _B;

	C := `if`( nargs > 2, CC, _B );
	xx0 := `if`( nargs > 3, x0, 0 );

	Ser := `Slode/series`(eq,fn,v(n),C,xx0);
	Rec := `Slode/REdetermine`( Ser[2], v(n));
	if Rec[ord] = 0 then
		sln := 0;
	else
		sln := LREtools['hypergeomsols'](Ser[2],v(n),{},output=basis);
	fi;

	if sln = NULL then
		RETURN( NULL );
	fi;

	`Slode/somecoeff`(Ser, sln, `Slode/defined`(sln,n)+Rec['ord']);
end:
######################################################################
# INPUT: 
#	expr - any expression in x ( rational, with GAMMA-function,
#                                    product )
# OUTPUT:
#	x0 - nonnegint such that expr(x) determine in x0,x0+1,...
`Slode/defined` := proc( expr, x )
local op_gamma, x0s ;

	x0s := select(type,{solve(denom(expr),n)},nonnegint);

	op_gamma := map(
	  (f,xx)->if op(0,f)=`GAMMA` and 
		     type(op(f),linear(xx))
		  then op(f) 
		  fi,
	  indets(expr,function),x);
	x0s := x0s union 
	       map((ln,xx)-> -coeff(ln,xx,0)/coeff(ln,xx),op_gamma,x)
		   union
	       map((f,xx)->
		  if op(0,f)=`product` then
			max(solve(op(1,rhs(op(2,f)))=
			    op(2,rhs(op(2,f)))))-1
		  fi,indets(expr,function),x);
	max(-1,op(select(type, x0s, integer))) + 1 ;
end:
######################################################################
# INPUT :
# eq - a homogeneous linear ordinary differential equation with 
#      polynomial coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x)
# fn - the unknown function in the form y(x)
# C - name of arbitrary constants (optional argument) (not _C),
#     _B by default 
# 
# OUTPUT :
#  formal power series solution with rational coefficients;
#  points of these series are chose automatically;
#  solutions may involve arbitrary constants in form _C[i], C[j]
#  (where i,j are some integer).
#
`Slode/ratcoeffs` := proc( eq :: {algebraic, alpgebraic=algebraic},
			fn::anyfunc(name),
			CC::name)
local C, res, Eq, p_r, sng, addnl, Rec, x, x0, v, n, M, nEq, p0, k ;
global _B;

	C := `if`( nargs > 2, CC, _B );
	Eq := `Slode/DEdetermine`(eq, fn);

	if not Eq['linear'] then
		ERROR(`The differential equation must be linear on`, Eq['fun'] );
	fi;
	if not Eq['polycfs'] then
		ERROR(`The linear differential equation must be `, 
                      `with polynomial coefficients on`, Eq['var']);
	fi;
	if Eq['rhs'] <> 0 then
		ERROR( `eq must be a homogeneous equation` );
	fi;

	x := Eq['var'];

	p_r := factor(Eq['coeffs'][Eq['ord']]);
	p_r := `if`( type(p_r,`*`), convert(p_r,set), {p_r} );

	sng := map(RootOf,select(has,p_r,x),x);

	Rec := `Slode/recurrfordifft_p`(Eq,v(n),min(sng[],1)-1);
	M := max(select( type, 
		{solve(subs( n=n+Rec['b'],Rec['coeffs'][-Rec['b']]),n)},
		nonnegint )[],-1) +1 ;

	nEq := eval(Eq);
	for k to M do
	   p0 := nEq['coeffs'][k-1] ;
	   if p0 <> 0 then
		nEq := DEdetermine(
		    numer(simplify(diff(nEq['L']/p0,x))), fn) ;
	   fi;
	od;
	if M > 0 then
		p_r := factor(nEq['coeffs'][nEq['ord']]);
		p_r := `if`( type(p_r,`*`), convert(p_r,set), {p_r} );
		addnl := map(RootOf,select(has,p_r,x),x);
		addnl := map(`+`,addnl,-1) union addnl;
	else
		addnl := {};
	fi;

	res := {};
	for x0 in sng union addnl do
	   res := res union {`Slode/rat_at_point`(eq,fn,C,x0)};
	od;

	res minus {0};

end:
######################################################################
# INPUT :
# eq - a homogeneous linear ordinary differential equation with 
#      polynomial coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x)
# fn - the unknown function in the form y(x)
# C - name of arbitrary constants (optional argument) (not _C)
# x0 - point (optional argument)
# 
`Slode/rat_at_point` := proc( eq :: {algebraic, algebraic=algebraic},
			fn::anyfunc(name),
			CC::name,
			x0::algebraic )
local xx0, C, v, n, Ser, Rec, sln, x0s;
global _B;

	C := `if`( nargs > 2, CC, _B );
	xx0 := `if`( nargs > 3, x0, 0 );

	Ser := `Slode/series`(eq,fn,v(n),C,xx0);
	Rec := `Slode/REdetermine`( Ser[2], v(n));
	if Rec['ord'] = 0 then
		sln := 0
	else
		sln := `solrec/ratsolrec`(Ser[2],v(n));
	fi;

	if sln = 0 then
		RETURN( NULL );
	fi;

	x0s := select(type,{solve(denom(normal(sln)),n)},nonnegint);	

	`Slode/somecoeff`(Ser, sln, max(x0s[],-1) + 1 + Rec['ord']);
end:
######################################################################
# INPUT :
# eq - a homogeneous linear ordinary differential equation with 
#      polynomial coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x)
# fn - the unknown function in the form y(x)
# C - name of arbitrary constants (optional argument) (not _C),
#     _B by default 
# 
# OUTPUT :
#  formal power series solution with polynomial coefficients;
#  points of these series are chose automatically;
#  solutions may involve arbitrary constants in form _C[i], C[j]
#  (where i,j are some integer).
#
`Slode/polycoeffs` := proc( eq :: {algebraic, alpgebraic=algebraic},
			fn::anyfunc(name),
			CC::name)
local C, res, Eq, p_r, sng, x, x0 ;
global _B;

	C := `if`( nargs > 2, CC, _B );
	Eq := `Slode/DEdetermine`(eq, fn);

	if not Eq['linear'] then
		ERROR(`The differential equation must be linear on`, Eq['fun'] );
	fi;
	if not Eq['polycfs'] then
		ERROR(`The linear differential equation must be `, 
                      `with polynomial coefficients on`, Eq['var']);
	fi;
	if Eq['rhs'] <> 0 then
		ERROR( `eq must be a homogeneous equation` );
	fi;

	x := Eq['var'];

	p_r := factor(Eq['coeffs'][Eq['ord']]);
	p_r := `if`( type(p_r,`*`), convert(p_r,set), {p_r} );

	sng := map(RootOf,select(has,p_r,x),x);

	res := {};
	for x0 in sng do
	   res := res union {`Slode/poly_at_point`(eq,fn,C,x0-1)};
	od;

	res minus {0};

end:
######################################################################
# INPUT :
# eq - a homogeneous linear ordinary differential equation with 
#      polynomial coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x)
# fn - the unknown function in the form y(x)
# C - name of arbitrary constants (optional argument) (not _C)
# x0 - point (optional argument)
# 
`Slode/poly_at_point` := proc( eq :: {algebraic, algebraic=algebraic},
			fn::anyfunc(name),
			CC::name,
			x0::algebraic )
local xx0, C, v, n, Ser, sln;
global _B;

	C := `if`( nargs > 2, CC, _B );
	xx0 := `if`( nargs > 3, x0, 0 );

	Ser := `Slode/series`(eq,fn,v(n),C,xx0);
	sln := LREtools['polysols'](Ser[2],v(n),{},output=basis);

	if sln = NULL then
		RETURN( NULL );
	fi;

	`Slode/somecoeff`(Ser, sln, 0);
end:
######################################################################
# INPUT :
#   FS - formal power series in form
#	[ v(0)+...+v(M)*(x-x0)^M+Sum(v(n)*(x-x0)^n,n=M+1..infinity),
#         recurrence]
#   fn - anything function of n that satisfies the recurrence for
#		n >= M1
#   
# OUTPUT :
#      v~(0)+...+v~(K)*(x-x0)^K+Sum(fn~*(x-x0)^n,n=K+1..infinity)
#
`Slode/somecoeff` := proc( FS :: list, fn, M1::integer )
local FSn, v, n, d, sys, K, i, j, solsys, f, hdr, cons, ri, res, res1;

	FSn := `Slode/FSdetermine`( FS ) ;
	
	n := FSn['rec']['var'];
	v := FSn['rec']['fun'];
	d := FSn['rec']['ord'];

	hdr := eval(FSn['header']);
	if FSn['M'] < M1-1 then
	   for i from hdr['max_i']+1 to M1-1 do
		ri := add(subs(n=i,FSn['rec']['coeffs'][j])*
	    	      hdr['coeffs'][i+j], j=-d..-1) + 
            	      subs(n=i,FSn['rec']['coeffs'][0])*v(i);
            	hdr['coeffs'][i] := 
		      normal(-coeff(ri,v(i),0)/coeff(ri,v(i),1));
	   od;
	   hdr['max_i'] := max(hdr['max_i'], M1-1);
	   K := M1 - 1 - d;
	else
	   K := FSn['M'] - d;
	fi;

	sys := {seq(simplify(subs(n=i,fn))=
		   hdr['coeffs'][i],i=K+1..K+d)};
	solsys := solve(sys);
	if solsys = NULL then
		RETURN( NULL );
	fi;

	f := factor(subs(solsys, fn));
	hdr := subs(solsys, eval(hdr));

	if normal(f) = 0 then 
		RETURN(add(hdr['coeffs'][i]*
			FSn['P_1']^i,i=0..K));
	fi ;

	if type(f,`*`) then
		cons := simplify(remove(depends,f,n)) ;
		f := select(depends,f,n) ;
	elif not depends(f,n) then
		cons := f ;
		f := 1 ;
	else
		cons := 1;
	fi;

	res1 := cons*Sum(f*FSn['P_n'],n=K+1..infinity);
	res := 0;
	while K>= max(M1-d,0) do
	   sys := hdr['coeffs'][K]=subs(n=K,f*cons);
	   solsys := solve({sys});
	   if testeq(sys) then ## f(K) equiv v(K)
		K := K-1;
		res1 := cons*Sum(f*FSn['P_n'],n=K+1..infinity);

	   elif solsys = NULL then ## f(K) <> v(K)
		RETURN( add(hdr['coeffs'][i]*FSn['P_1']^i,i=0..K)+
                        res + res1);
	   else ## f~~(K) = v~~(K)
     		res := res + res1;
	    	hdr := subs(solsys,eval(hdr));
	    	f := factor(subs(solsys,f));
	    	if normal(f) = 0 then 
		   RETURN( add(hdr['coeffs'][i]*FSn['P_1']^i,i=0..K)+
                        res) ;
		fi;
	    	if type(f,`*`) then
		   cons := simplify( cons*remove(depends,f,n) );
		   f := select(depends,f,n) ;
	    	elif not depends(f,n) then
		   cons := simplify(cons*f) ;
		   f := 1 ;
	    	fi;
	    	K := K-1 ;
		res1 := cons*Sum(f*FSn['P_n'],n=K+1..infinity);
		
	   fi;
	od;

	add(hdr['coeffs'][i]*FSn['P_1']^i,i=0..K) + res1 + res ;
end:
######################################################################
# INPUT :
# eq - a linear ordinary differential equation with polynomial
#     coefficients and polynomial right hand side: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x) = f(x)
#
# fn - the unknown function in the form y(x)
# vn - the new function of the recurrence in the form v(n)
# C  - name of arbitrary constants (optional)
# x0 - point (optional argument)
# M  - the required number of initial elements of a power 
#      series solutions (optional argument)
# OUTPUT :
#   formal power series solution:
#	v(0)+v(1)*(x-x0)+...+v(M)*(x-x0)^M + 
#	   Sum(v(n)*(x-x0)^n,n=M+1..infinity)
#   and recurrence for v(n), n>M.
`Slode/series` := proc( eq::{algebraic,algebraic=algebraic},
			fn::anyfunc(name),
			vn::anyfunc(name),
			CC::name )
local  C, x0, Eq, Rec, S, N, M, alpha_A, n, k, 
       V, lf, v, e_n, i, j, _Ck, Cs, M1 ;
global _C;

	x0 := `if`(nargs>4,args[5],0);
	C := `if`(nargs>3,CC,_C); 
	M1 := `if`(nargs>5 and type(args[6],nonnegint),args[6], -1);

	Eq := `Slode/DEdetermine`(eq,fn);

	if not Eq['linear'] then
		ERROR(`The differential equation must be linear on`, Eq['fun'] );
	fi;
	if not Eq['polycfs'] then
		ERROR(`The linear differential equation must be `, 
                      `with polynomial coefficients on`, Eq['var']);
	fi;
	Rec := `Slode/recurrfordifft_p`(Eq,vn,x0);

	n := Rec['var'];
	v := Rec['fun'];
	alpha_A := normal(subs( n=n-Rec['A'],Rec['coeffs'][Rec['A']]));
	S := select( type, map2(op,1,roots(alpha_A,n)),
		     nonnegint);

	N := {S[],seq(k,k=0..Rec['A']-1)};
	M := max( S[], nops(Rec['rhs'])-1+Rec['A'], Rec['ord']-1 );
	M := max( M, M1 );

	V := table(sparse);
	Cs := [];
	lf := table(sparse, [seq(k=Rec['rhs'][k+1],
				k=0..nops(Rec['rhs'])-1)] );
	for k from 0 to M do
	    e_n := subs(n=k-Rec['A'],Rec['L'])-lf[k-Rec['A']];
	    e_n := subs(seq(v(k-i)=V[k-i],i=1..Rec['A']+Rec['b']), e_n );

	    if not member( k, N ) then  ## lcoeff(e_n,v(k))<>0
		V[k] := simplify(-coeff(e_n,v(k),0)/coeff(e_n,v(k),1));
	    else
		V[k] := C[k];
		Cs := [Cs[],C[k]];
		if k >= Rec['A'] and not testeq(e_n) then
		   _Ck := indets(e_n) intersect {Cs[]} ;
		   if _Ck = {} then
			RETURN(NULL);
		   fi;
		   i := max(seq(op(j),j=_Ck));
		   member(C[i],Cs,'j');
		   e_n := collect(e_n,Cs);
		   _Ck := simplify(
				-coeff(e_n,Cs[j],0)/coeff(e_n,Cs[j],1));
		   V := subs(Cs[j]=_Ck,eval(V));
		   Cs := [Cs[1..j-1][],Cs[j+1..-1][]];
		fi;
	    fi;
	od;

	[add(V[k]*(Eq['var']-x0)^k,k=0..M)+
	 Sum( vn*(Eq['var']-x0)^Rec['var'],Rec['var']=M+1..infinity),
	 subs(Rec['var']=Rec['var']-Rec['A'],Rec['L'])]
end:
######################################################################
# INPUT :
# eq - a linear ordinary differential equation with polynomial
#     coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x) = f(x)
#
# fn - the unknown function in the form y(x)
# vn - the new function of the recurrence in the form v(n)
# x0 - point (option argument)
# OUTPUT :
#	recurrence for differential equation in basis (x-x0)^n
`Slode/recurrfordifft_p` := proc( eq )
local  n, v, k, i, j, res, x ;

  if not type(eq, table) then
    if nargs = 1 then
      ERROR(`The 2nd argument fn (of type anyfunc(name)) is missing`);

    elif nargs = 2 then
      ERROR(`The 3rd argument vn (of type anyfunc(name)) is missing`);

    elif not type(args[3],anyfunc(name)) then
      ERROR(`The 3rd argument vn (of type anyfunc(name)) are expected`, 
	    ` but received`, args[3]);

    elif nargs = 3  then		## recurrfordifft_p(eq,fn,vn)
      res := `Slode/recurrfordifft_p`(
			`Slode/DEdetermine`(eq,args[2]), args[3]);
      if res['rhs'] = [] then 
	RETURN(res['L']);
      else
	RETURN([res['L'],res['rhs']]);
      fi;
		
    else 				## recurrfordifft_p(eq,fn,vn,x0)
      RETURN(`Slode/recurrfordifft_p`(
			DEtools['translate'](eq,op(args[2]),args[4],args[2]),
			args[2], args[3] ));
    fi;

  elif nargs = 1 then
    ERROR(`The 2nd argument vn (of type anyfunc(name)) is missing`);

  elif not type(args[2],anyfunc(name)) then
    ERROR(`The 2nd argument vn (of type anyfunc(name)) is expected`,
          ` but received`, args[2]);

  elif nargs > 2 and args[3] <> 0 then	## recurrfordifft_p(eq,vn,x0)
    RETURN( `Slode/recurrfordifft_p`(
            `Slode/DEdetermine`(
		DEtools['translate'](eq['L']-eq['rhs'],eq['var'],args[3],eq['fun'](eq['var'])),
		eq['fun'](eq['var'])),args[2]) );
  fi;

  if not eq['linear'] then
    ERROR(`The differential equation must be linear on`, eq['fun'] );
  fi;
  if not eq['polycfs'] then
    ERROR(`The linear differential equation must be with polynomial coefficients on`, 
	  eq['var']);
  fi ;

  v := op(0,args[2]);
  n := op(args[2]);
  x := eq['var'];

  res := `Slode/REdetermine`(
	    simplify(
	     add( add(eq[_coeffs][i][j]*v(n+i-j)*mul(n+k-j,k=1..i),
	         		j=0..eq['d_max']),
	          i=0..eq['ord']), 
	     RootOf), v(n) );
  if eq['rhs'] <> 0 then
    res['rhs'] := [seq(coeff(eq['rhs'],x,i),i=0..degree(eq['rhs'],x))];
  fi;

  eval( res );
end:
######################################################################
# INPUT :
# eq - a linear ordinary differential equation with polynomial
#     coefficients and polynomial right hand side: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x) = f(x)
#
# fn  - the unknown function in the form y(x)
# 
`Slode/DEdetermine` := proc( eq :: {algebraic, algebraic = algebraic},
                           fn :: anyfunc(name) )
local  Result, cfs, i, p1, pp, k, x;
option remember ;

  Result := table() ;
  Result['fun'] := op(0,fn);
  x := op(fn) ;
  Result['var'] := x;

  # construct coefficient list for the input linear ODE
  cfs := DEtools['convertAlg'](eq,fn);
  if cfs = FAIL then
    # the input ``eq'' is not a linear ODE
    Result['linear'] := false;
    Result['L'] := eq;
    RETURN(eval(Result));
  elif nops(cfs[1]) = 1 and cfs[1][] = 0 then
    # the input ``eq'' is not a DE
    Result['linear'] := false;
    Result['L'] := 0;
    Result['rhs'] := cfs[2];
    RETURN(eval(Result));
  fi;

  Result['linear'] := true;
  Result['ord'] := nops(cfs[1])-1;

  # Normalization
  cfs := DEtools['DEnormal'](cfs,x,fn);
  if denom(cfs[2]) <> 1 then
    cfs := [map(`*`,cfs[1],denom(cfs[2])),numer(cfs[2])];
  fi;

  # Check if the linear ODE has polynomial coefficients
  Result['polycfs'] := type([cfs[1][],cfs[2]], 
			    list(polynom(anything,x))) ;

  if not Result['polycfs'] then
    Result['rhs'] := cfs[2];
    Result['coeffs'] := array(0..Result['ord'],
		[seq(cfs[1][i],i=1..Result['ord']+1)]);
  else
    Result['rhs'] := collect(cfs[2],x);
    Result['coeffs'] := array(0..Result['ord'],
		[seq(collect(cfs[1][i],x),i=1..Result['ord']+1)]);
    Result['d_max']:= max(seq(degree(Result[coeffs][i],x),
			      i=0..Result['ord']));
    Result['_coeffs'] := array(0..Result['ord']);

    for i from 0 to Result['ord'] do
      if Result['coeffs'][i] = 0 then
        Result['_coeffs'][i] := table(sparse,['ord'=-infinity]);
      else
	p1 := [coeffs(Result['coeffs'][i],x,'pp')];
	pp:= [pp];
	Result['_coeffs'][i] := 
	  table(sparse,['ord' = degree(Result['coeffs'][i],x),
			  seq(degree(pp[k],x)=p1[k],k=1..nops(p1))]);
      fi;
    od;
  fi;

  Result['L'] := Result['coeffs'][0]*fn + 
		add(Result['coeffs'][i]*diff(fn,x$i),
			i=1..Result['ord']);
  eval(Result)

end :
######################################################################
# INPUT :
# rec - a linear recurrence with polynomial coefficients: 
#   a[-A](n)v(n+A)+...+a[0](n)v(n)+...+a[b](n)v(n-b) 
# if rec is not homogeneous then right hand side must be
#  finite sequence f(n)
#   [a[-A](n)v(n+A)+...+a[0](n)v(n)+...+a[b](n)v(n-b),
#	[f(0),f(1),...]
# fn  - the unknown function in the form v(n)
# 
`Slode/REdetermine` := proc( 
	eq :: {algebraic, algebraic = algebraic,
	      [{algebraic,algebraic=algebraic},list(anything)]},
        fn :: anyfunc(name) )
local  Result, sh_o, sh, cfs, t, tt, i, k, n;
option remember ;

	Result := table() ;
	Result['L'] := `if`(type(eq,list),eq[1],eq);
	Result['L'] := `if`(type(Result['L'],`=`), 
		lhs(Result['L'])-rhs(Result['L']), Result['L']);
	Result['fun'] := op(0,fn);
	n := op(fn) ;
	Result['var'] := n;

	sh_o := map(`Slode/isshift`,indets(Result['L']),fn);
	sh := map2(op,1,sh_o);

	if sh_o = {} then ## recurrence do not depent of fn
		Result['ord'] := -1;
		Result['A'] := -1;
		Result['b'] := 0;
	else
		Result['A'] := max( map2(op,2,sh_o)[] );
		Result['b'] := -min( map2(op,2,sh_o)[] );
		Result['ord'] := Result['A'] + Result['b'];
	fi;

	Result['linear'] := type(Result['L'], linear(sh));
	if not Result['linear'] then
		RETURN(eval(Result));
	fi;

	Result['rhs'] := `if`(type(eq,list),eq[2],[]) ;
	if not testeq(subs( {seq( t=0, t=sh )}, Result['L'] )) then
		ERROR(`the right hand side of recurrence must`,
		      ` be finite sequence (in form list)`,
		      ` but not`, 
		      subs( {seq( t=0, t=sh )}, Result['L'] ));
	fi;

	Result['L'] := collect(Result['L'], sh);
	cfs := [coeffs(Result['L'],sh,'tt')] ; 

	Result['polycfs'] := 
		type(cfs, list(polynom(anything,n)));

	Result['coeffs'] := table(sparse);
	for t in sh_o do
	    if member(t[1],[tt],'i') then
		Result['coeffs'][t[2]] := collect(cfs[i],n);
	    fi;
	od;
	Result['L'] := add(Result['coeffs'][i]*subs(n=n+i,fn),
			   i=-Result['b']..Result['A']);
	eval(Result)

end :
######################################################################
# INPUT : z, y(x)
# OUTPUT : true if z is y(x+k)
`Slode/isshift` := ( z, y ) ->
   if (z = y) then
	[z,0]
   elif type(z,op(0,y)(anything)) and type(op(z),`+`) 
	and degree(op(z),op(y)) = 1 and lcoeff(op(z),op(y)) = 1 then
		[z, tcoeff(op(z),op(y))]
   fi :
######################################################################
# INPUT :
#   FS - formal power series in form
#	[ v(0)+...+v(M)*(x-x0)^M+Sum(v(n)*(x-x0)^n,n=M+1..infinity),
#         recurrence]
`Slode/FSdetermine` := proc( FS :: list )
local  Result, v, vn, n, h, cf, tt, xx, i, hedr, kn_t, t, Rec, k;
option remember ;

	Result := table() ;

	Result['sum'] := indets(FS[1],'Sum'(anything,anything))[];
	Result['header'] := table(sparse);
	Result['header']['all'] := subs(Result['sum']=0,FS[1]);

	v := op([1,1,0],Result['sum']);
	kn_t := op([1,2,2],Result['sum']);
      n := indets(kn_t,name)[];
	Result['sparse'] := lcoeff(kn_t,n);
	t := coeff(kn_t,n,0);

      if Result['sparse'] = 1 then
	    Result['rec'] := `Slode/REdetermine`( FS[2], v(n) );
	else
	    vn := indets(FS[2],specfunc(anything,v));
          Rec := add(coeff(FS[2],k)*v((op(k)-t)/Result['sparse']),
			k=vn);
	    Result['rec'] := `Slode/REdetermine`(Rec,v(n));
        fi;

	Result['P_n'] := op([1,2],Result['sum']);
	Result['P_1'] := op(1,Result['P_n']);
	if type(Result['P_1'],name) then
		Result['mvar'] := Result[P_1];
		Result['point'] := 0;
	else
		Result['mvar'] := op(1,Result[P_1]);
		Result['point'] := -op(2,Result[P_1]);
	fi;
	Result['n_l'] := op(1,rhs(op(2,Result['sum']))) ;
	Result['M'] := Result['n_l']*Result['sparse'] + t - 1;;

	if Result['M'] >=0 then
		hedr := collect(subs(Result['mvar']=xx+Result['point'],
			Result['header']['all']), xx);
	      cf := [coeffs(hedr,xx,tt)];
		tt := [tt] ;
		Result['header']['coeffs'] := 
		    table(sparse,
		          [seq(degree(tt[i],xx)=cf[i],i=1..nops(tt))] 
			  );
		Result['header']['max_i'] := Result['M'];
	fi;

	eval(Result);
end :
######################################################################
######################################################################
# INPUT :
# eq - a homogeneous linear ordinary differential equation with 
#      polynomial coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x)
# fn - the unknown function in the form y(x)
# C - name of arbitrary constants (optional argument) (not _C),
#     _Const by default 
# 
# OUTPUT :
#  formal power series solution with m-hypergeometric coefficients;
#  points of these series are chose automatically;
#  solutions may involve arbitrary constants in form _C[i], C[j]
#  (where i,j are some integer).
#
`Slode/mhypercoeffs` := proc( eq :: {algebraic, alpgebraic=algebraic},
            fn::anyfunc(name),
            CC::name)
local mSers, v, n  ;

    mSers := `Slode/msparse`(eq,fn,v(n),args[3..-1]);
    if mSers = {} then
	RETURN({});
    fi;
    map(`Slode/mhyper_`,mSers,v,n);

end:

######################################################################
# INPUT :
# eq - a homogeneous linear ordinary differential equation with 
#      polynomial coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x)
# fn - the unknown function in the form y(x)
# m - order of m-hypergeometric solution
# C - name of arbitrary constants (optional argument) (not _C)
# x0 - algebraic number
# 
`Slode/mhyper_at_point` := proc( eq :: {algebraic, algebraic=algebraic},
            fn::anyfunc(name),
	    m :: posint,
            s0::algnum,
	    CC::name )
local res, v, n, mSers, mS, Rec, S, vn, t, vns, k, h, nRec;
    mSers := `Slode/msparse_at_point`(args[1..2],v(n),args[3..-1]);
    if mSers = {} then
	RETURN({});
    fi;
    res := map(`Slode/mhyper_`,mSers,v,n);
end:
######################################################################
`Slode/mhyper_` := proc(mSer,v,n)
local  fn, mS, sys, K, re, d, hdr, i, solsys, cons;
	mS := `Slode/FSdetermine`(mSer);
      re := eval(mS['rec']);
      fn := LREtools[hypergeomsols](re['L'],v(n),{},output=basis);
      if fn = NULL then
	   RETURN(NULL);
      fi;
	
    	d := re['ord'];
    	hdr := eval(mS['header']);

{seq(simplify(subs(n=mS[n_l]-i,f(n))) =
		    v(mS['M']+1-mS['sparse']*i),i=1..d)};
	sys := {seq(simplify(subs(n=mS[n_l]-i,fn)) =
		    hdr['coeffs'][mS['M']+1-mS['sparse']*i],i=1..d)};
	solsys := solve( sys );
    	if solsys = NULL then
          RETURN( NULL );
	fi;

        fn := factor(subs(solsys, fn));
        hdr := subs(solsys, eval(hdr));

        if normal(fn) = 0 then 
            RETURN(add(hdr['coeffs'][i]*mS['P_1']^i,i=0..mS['M']));
        fi ;

        if type(fn,`*`) then
            cons := simplify(remove(depends,fn,n)) ;
            fn := select(depends,fn,n) ;
    	elif not depends(fn,n) then
            cons := fn ;
            fn := 1 ;
    	else
            cons := 1;
        fi;

        add(hdr['coeffs'][i]*mS['P_1']^i,i=0..mS['M']-mS['sparse']*d) + 
	cons*Sum(fn*mS['P_n'],n=mS['n_l']-d..infinity);
end: 
######################################################################
# INPUT :
# eq - a homogeneous linear ordinary differential equation with 
#      polynomial coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x)
# fn - the unknown function in the form y(x)
# CC - name of arbitrary constants (optional argument) 
`Slode/msparse` := proc( eq :: {algebraic, algebraic=algebraic},
            fn::anyfunc(name),
            vn::anyfunc(name),
            CC::name )
local C, R, Rec, v, n, Eq, m, res, f_h, F_H, Rfh, a;
global _B;

    C := `if`( nargs > 3, CC, _B );
    Eq := `Slode/DEdetermine`(eq,fn);

    if not Eq['linear'] then
        ERROR(`The differential equation must be linear in`, Eq['fun'] );
    fi;
    if not Eq['polycfs'] then
        ERROR(`The linear differential equation must have polynomial coefficients on`, 
             Eq['var']);
    fi;
    if Eq['rhs'] <> 0 then
        ERROR( `eq must be a homogeneous equation` );
    fi;
    Rec := `Slode/recurrfordifft_p`(Eq,vn);
    v := Rec['fun'];
    n := Rec['var'];
    res := {};
    for m from 2 to Rec['ord'] do
       f_h := `Slode/first_half`( Eq['L'],m,Eq['var'],Eq['fun'],
		       -Rec['b'],a);
       if f_h[2] <> {} then
    	   res := res union {seq(`Slode/msparse_at_point`(Eq['L'],fn, 
			v(n),m, F_H, C)[], 
			F_H = map(`Slode/eq2algnum`,f_h[2],a) )};
       elif f_h[1] <> 1 then
    	   res := res union {`Slode/msparse_at_point`(f_h[1],fn,
			v(n),m,0,C)[]};
       fi;
    od;
    res;
end:
`Slode/eq2algnum` := proc(eq,x)
local p_r;
    p_r := factor(eq);
    p_r := `if`( type(p_r,`*`), convert(p_r,set), {p_r} );

    map(RootOf,select(has,p_r,x),x)[];
end:
######################################################################
# INPUT :
# eq - a homogeneous linear ordinary differential equation with 
#      polynomial coefficients: 
#   p[r](x)diff(y(x),x$r)+...+p[1](x)diff(y(x),x)+p[0](x)y(x)
# fx - the unknown function in the form y(x)
# vn - the function for recurrence in the form v(n)
# m - integer
# x0 - algebraic number (optional argument)
# CC - name of arbitrary constants (optional argument) 
`Slode/msparse_at_point` := proc( eq :: {algebraic, algebraic=algebraic},
            fx::anyfunc(name),
            vn::anyfunc(name),
	    m::posint,
            x0::algnum,
            CC::name )
local C, v, n, Rec, Rt, j, t, E, V, T, w, w1, tau, W, u, flag, Inw_,
      Ci, i, slv, cf, eqj, res, k, a, Ser, par;
global _B;

    C := `if`( nargs > 5, CC, _B );
    a := `if`( nargs > 4, x0, 0 );

    Rec := `Slode/REdetermine`(`Slode/recurrfordifft_p`( eq, fx, vn, a ),vn);
    v := Rec['fun'];
    n := Rec['var'];

    Rt := array(sparse,0..m-1);
    for j from 0 to Rec['ord'] do
      Rt[j mod m]:=Rt[j mod m]+subs(n=n+Rec['b'],Rec['coeffs'][j-Rec['b']])*E^j;
    od;
    V := Rt[0];
    T := {};
    for t to m-1 while degree(V,E)>0 do
      w:=`Slode/eGCDpar`(V,Rt[t],0,E,n,par);
      if w = [] then RETURN({}) fi;
      V:=w[1][2];
      T:=w[1][1];
    od;

    if degree(V,E) = 0 then RETURN({}) fi;
    
    tau := degree(V,E);
    cf := subs(n=n - Rec['b'],[seq(coeff(V,E,i),i=0..tau)]);
    tau := tau - Rec['b'];
    w := cf[-1];
    w1 := cf[1];
    W := subs( n=n-Rec['b'], `Slode/convert_to_rec`(V,E,v,n) );

    u := max( `Slode/iota2`(w1,n,-Rec['b']), 
	      `Slode/iota2`(w*Rec['coeffs'][-Rec['b']],n,-Rec['b']),
	      `Slode/iota1`(w,n,Rec['A']),
	      `Slode/iota1`(w*Rec['coeffs'][Rec['A']],n,Rec['A']));
    w := u + Rec['ord'];
    Ser := `Slode/FSdetermine`(`Slode/series`( eq, fx, vn, C, a, w ));
    res := {};
    for t from 0 to m-1 do
      t := t;
      Inw_ := copy(Ser['header']['coeffs']);
      for i from u+1 to w while flag <> NULL do
        i := i;
	if (i-t) mod m <> 0 and Inw_[i] <> 0 then
            Ci := indets(Inw_[i],indexed);
	    if Ci = {} then flag := NULL ;
            else
	       Ci := Ci[1];
	       Inw_[i] := collect(Inw_[i],Ci);
	       slv := simplify( -coeff(Inw_[i],Ci,0)/
			coeff(Inw_[i],Ci,1));
	       Inw_ := subs(Ci =slv,eval(Inw_));
            fi;
        fi;
      od;
      for j from u+Rec['b']+1 to w-tau 
      while flag <> NULL do
        subs(n=j,W);
	eqj := simplify(add(subs(n=j,cf[i+1+Rec['b']])*Inw_[j+i],
		i=max(-Rec['b'],-j)..tau));
        if not testeq(eqj) then
          Ci := indets(eqj,indexed);
	  if Ci = {} then flag := NULL ;
          else
	    Ci := Ci[1];
	    eqj := collect(eqj,Ci);
	    slv := simplify( -coeff(eqj,Ci,0)/
			coeff(eqj,Ci,1));
	    Inw_ := subs(Ci = slv,eval(Inw_));
          fi;
        fi;
      od;
      if convert([seq(testeq(i[]),i=entries(Inw_))],`and`) and Rec['ord']<>0 then 
	flag := NULL;
      fi;

      if flag <> NULL then
	w := u + tau + Rec['b'] ;
	res := res union {[add(Inw_[k]*(op(fx)-a)^k,k=0..w)+
     		Sum( v(m*n+t)*(op(fx)-a)^(m*n+t),
                     n=(w+1+(m+t-(w+1 mod m)mod m)-t)/m..infinity),
		     subs(n=m*n+t,subs(n=n-tau,W))]}
      fi;      
    od;
    res;
end:
`Slode/convert_to_rec` := proc(R, Op, v, n)
local cf, t, i;
    cf := factor([coeffs(R,Op,'t')]);
    t := [t];
    convert([seq(cf[i]*v(n+degree(t[i],Op)),i=1..nops(t))],`+`);
    
end:
#####################################################################
`Slode/iota1` := ( qd, n, d ) ->
    d + max(-1, select( type, 
			map2(op,1,roots(qd,n)), nonnegint)[]):

`Slode/iota2` := ( ql, n, l ) ->
    max(l + max(-1, select( type, 
			map2(op,1,roots(ql,n)), nonnegint)[]),
		 -1):
#####################################################################
#####################################################################
`Slode/dGCDpar`:=proc(A1,B1,p1,X,Var,Par)
local A,B,p,r,w,h,res,S,sr,Lc1,Lc2,Lc12,C1,C2,s1,s2,m,n;

A:=A1;B:=B1;
p:=primpart(p1,Par);
if `dGCDpar/my_degree`(p,Par)=0 then RETURN([]) fi;

if p<>0 then A:=`dGCDpar/aldiv`(A,1,p,Par); B:=`dGCDpar/aldiv`(B,1,p,Par) fi;

A:=collect(A,X);
B:=collect(B,X);

# deg(A)>=deg(B) ?:
m:=degree(A,X);
n:=degree(B,X);
if m<n then RETURN(`Slode/dGCDpar`(B,A,p,X,Var,Par)) fi;

if A=0 then RETURN([`dGCDpar/filter`([p,B],X,Var,Par)]) fi;
if B=0 then RETURN([`dGCDpar/filter`([p,A],X,Var,Par)]) fi;

C1:=content(A,[X,Var]);
C2:=content(B,[X,Var]);
Lc1:=content(lcoeff(A,X),Var);
Lc2:=content(lcoeff(B,X),Var);

if p<>0 then
  C1:=gcd(C1,p);
  C2:=gcd(C2,p);
  Lc1:=gcd(Lc1,p);
  Lc2:=gcd(Lc2,p)
else
  C1:=`dGCDpar/sqrfree2`(C1,Par);
  C2:=`dGCDpar/sqrfree2`(C2,Par);
  Lc1:=`dGCDpar/sqrfree2`(Lc1,Par);
  Lc2:=`dGCDpar/sqrfree2`(Lc2,Par)
fi;

Lc12:=gcd(Lc1,Lc2);
divide(C1,gcd(C1,Lc12),'s1');
divide(C2,gcd(C2,Lc12),'s2');
res:=`dGCDpar/filter`([s1,B],X,Var,Par),`dGCDpar/filter`([s2,A],X,Var,Par),
     op(`Slode/dGCDpar`(A,B,Lc12,X,Var,Par));
if p<>0 then divide(p,gcd(p,Lc12*s1*s2),'p') fi;

S:=`dGCDpar/regsreschain`(A,B,X,Var);

if S=[] then
  RETURN([res,`dGCDpar/filter`([p,B],X,Var,Par)])
fi;

r:=content(lcoeff(S[nops(S)],X),Var);

if p=0 then
        r:=`dGCDpar/sqrfree2`(r,Par);
        divide(r,Lc12*s1*s2,'r')
       else r:=gcd(r,p)
fi;
divide(p,r,'p');
res:=res,`dGCDpar/filter`([p,S[nops(S)]],X,Var,Par);
S:=subsop(nops(S)=NULL,S);

while S<>[] do
   sr:=S[nops(S)]; S:=subsop(nops(S)=NULL,S);
   h:=gcd(r,content(lcoeff(sr,X),Var));
   divide(r,h,'w');
   r:=h;
   res:=`dGCDpar/filter`([w,sr],X,Var,Par),res
od;
res:=`dGCDpar/filter`([r,B],X,Var,Par),res;
[res]
end:
#######################################
`dGCDpar/regsreschain`:=proc(A,B,X,Var)
local m,n,A_,B_,a_,b_,l_,i,e_i,temp,u;

m:=degree(A,X);
n:=degree(B,X);
A_:=[A,B];
B_:=[A,B];
a_:=[1,lcoeff(A_[2],X)];
b_:=[1,lcoeff(B_[2],X)^(m-n)];
l_:=[_,m-n+1];

i:=2;
while A_[i]<>0 do
 i:=i+1;
 e_i:=(-1)^l_[i-1]*b_[i-2]^(l_[i-1]-1)*a_[i-2];

 divide(`dGCDpar/my_skew_prem`(A_[i-2],A_[i-1],X,Var),e_i,'temp');
 A_:=[op(A_),collect(temp,X)];
 l_:=[op(l_),degree(A_[i-1],X)-degree(A_[i],X)+1];
 a_:=[op(a_),lcoeff(A_[i],X)];
 if a_[i]<>0 then
   divide(A_[i]*a_[i]^(l_[i]-2),
                b_[i-1]^(l_[i]-2),'temp');
   B_:=[op(B_),collect(temp,X)];
   b_:=[op(b_),lcoeff(B_[i],X)]
 fi
od;
B_:=subsop(1=NULL,2=NULL,B_);
B_
end:
#######################################
`dGCDpar/sqrfree2`:=proc(p,v)
local t;

t:=sqrfree(p,v);
collect(mul(op([2,i,1],t),i=1..nops(op(2,t))),v)
end:
#######################################
`dGCDpar/filter`:=proc(d,X,Var,Par)
local p,q;

p:=op(1,d);q:=op(2,d);
if `dGCDpar/my_degree`(p,Par)=0 then RETURN(NULL)
elif q=0 then okay
else
  q:=`dGCDpar/primpart3`(`dGCDpar/aldiv`(q,content(lcoeff(q,X),Var),p,Par),X,[X,Var,Par])
fi;
if lcoeff(p,`dGCDpar/vlist`(p,[Par]))<0 then p:= -p fi;
if lcoeff(q,`dGCDpar/vlist`(q,[X,Var,Par]))<0 then q:= -q fi;
[`dGCDpar/canform`(p,[Par]),`dGCDpar/canform`(q,[X,Var,Par])]
end:
#######################################
`dGCDpar/vlist`:=proc(p,O)

[op(O),op(sort(convert(indets(p) minus convert(O,set),list),lexorder))]
end:
#######################################
`dGCDpar/primpart3`:=proc(p,X,l)

`dGCDpar/canform`(primpart(p,X),l)
end:
#######################################
`dGCDpar/canform`:=proc(p,O)
local vl;

vl:=`dGCDpar/vlist`(p,O);
sort(collect(expand(p),vl),vl,plex) end:
#######################################
`dGCDpar/aldiv`:=proc(c,q,p,Par)
local res;

if degree(p,Par)=0 then RETURN(c) fi;
res:=evala(Expand(subs(Par=RootOf(p,Par),c/q)));
if degree(p,Par)=1
then
  RETURN(res)
else
  RETURN(subs(RootOf(p,Par)=Par,res))
fi
end:
#######################################
`dGCDpar/my_skew_prem`:=proc(A,B,X,Var)
local m,n,i,res,T;

m:=degree(A,X); n:=degree(B,X);
if m<n then RETURN(A) fi;
T:=array(1..m-n+1);
T[m-n+1]:=B;
for i from m-n by -1 to 1
do
  T[i]:=`dGCDpar/d_skew_product`(X,T[i+1],X,Var);
od;
res:=A;
for i from 1 to m-n+1
do
  res:=expand(lcoeff(T[i],X)*res-coeff(res,X,m-i+1)*T[i])
od;
res
end:
#######################################
`dGCDpar/my_degree`:=proc(p,x)
local t;

t:=collect(p,x);
if t=0 then RETURN(-1)
else RETURN(degree(t,x)) fi
end:
#######################################
`dGCDpar/d_skew_product`:=proc(L1,L2,X,Var)
local i,p,res;

p:=L2; res:=0;
for i from 0 to degree(L1,X)
do
  res:=res+coeff(L1,X,i)*p;
  p:=`dGCDpar/dop`(p,X,Var)
od;
`dGCDpar/canform`(res,[X,Var])
end:
#######################################
`dGCDpar/dop`:=proc(L,D,V)
local i,res,c;

res:=0;
for i from 0 to degree(L,D)
do
  c:=coeff(L,D,i);
  res:=res+c*D^(i+1)+diff(c,V)*D^i
od;
res
end:
#######################################
#####################################################################
`Slode/eGCDpar`:=proc(A1,B1,p1,X,Var,Par)
local A,B,p,r,w,h,res,S,sr,Lc1,Lc2,Lc12,C1,C2,s1,s2,m,n;

A:=A1;B:=B1;
p:=primpart(p1,Par);
if `eGCDpar/my_degree`(p,Par)=0 then RETURN([]) fi;


if p<>0 then A:=`eGCDpar/aldiv`(A,1,p,Par); B:=`eGCDpar/aldiv`(B,1,p,Par) fi;

A:=collect(A,X);
B:=collect(B,X);

# deg(A)>=deg(B) ?:
m:=degree(A,X);
n:=degree(B,X);
if m<n then RETURN(`Slode/eGCDpar`(B,A,p,X,Var,Par)) fi;

if A=0 then RETURN([`eGCDpar/filter`([p,B],X,Var,Par)]) fi;
if B=0 then RETURN([`eGCDpar/filter`([p,A],X,Var,Par)]) fi;

C1:=content(A,[X,Var]);
C2:=content(B,[X,Var]);
Lc1:=content(lcoeff(A,X),Var);
Lc2:=content(lcoeff(B,X),Var);

if p<>0 then
  C1:=gcd(C1,p);
  C2:=gcd(C2,p);
  Lc1:=gcd(Lc1,p);
  Lc2:=gcd(Lc2,p)
else
  C1:=`eGCDpar/sqrfree2`(C1,Par);
  C2:=`eGCDpar/sqrfree2`(C2,Par);
  Lc1:=`eGCDpar/sqrfree2`(Lc1,Par);
  Lc2:=`eGCDpar/sqrfree2`(Lc2,Par)
fi;

Lc12:=gcd(Lc1,Lc2);
divide(C1,gcd(C1,Lc12),'s1');
divide(C2,gcd(C2,Lc12),'s2');
res:=`eGCDpar/filter`([s1,B],X,Var,Par),`eGCDpar/filter`([s2,A],X,Var,Par),
     op(`Slode/eGCDpar`(A,B,Lc12,X,Var,Par));
if p<>0 then divide(p,gcd(p,Lc12*s1*s2),'p') fi;

S:=`eGCDpar/regsreschain`(A,B,X,Var);

if S=[] then
  RETURN([res,`eGCDpar/filter`([p,B],X,Var,Par)])
fi;

r:=content(lcoeff(S[nops(S)],X),Var);

if p=0 then
        r:=`eGCDpar/sqrfree2`(r,Par);
        divide(r,Lc12*s1*s2,'r')
       else r:=gcd(r,p)
fi;
divide(p,r,'p');
res:=res,`eGCDpar/filter`([p,S[nops(S)]],X,Var,Par);
S:=subsop(nops(S)=NULL,S);

while S<>[] do
   sr:=S[nops(S)]; S:=subsop(nops(S)=NULL,S);
   h:=gcd(r,content(lcoeff(sr,X),Var));
   divide(r,h,'w');
   r:=h;
   res:=`eGCDpar/filter`([w,sr],X,Var,Par),res
od;
res:=`eGCDpar/filter`([r,B],X,Var,Par),res;
[res]
end:
#######################################
`eGCDpar/regsreschain`:=proc(A,B,X,Var)
local m,n,A_,B_,a_,b_,l_,i,e_i,temp,u;

m:=degree(A,X);
n:=degree(B,X);
A_:=[A,B];
B_:=[A,B];
a_:=[1,lcoeff(A_[2],X)];
b_:=[1,`eGCDpar/sigfact`(lcoeff(B_[2],X),m-n,Var)];
l_:=[_,m-n+1];

i:=2;
while A_[i]<>0 do
 i:=i+1;
 e_i:=(-1)^l_[i-1]*`eGCDpar/sigfact`(`eGCDpar/sigma`(b_[i-2],Var),l_[i-1]-1,Var)*a_[i-2];

 divide(`eGCDpar/my_skew_prem`(A_[i-2],A_[i-1],X,Var),e_i,'temp');
 A_:=[op(A_),collect(temp,X)];
 l_:=[op(l_),degree(A_[i-1],X)-degree(A_[i],X)+1];
 a_:=[op(a_),lcoeff(A_[i],X)];
 if a_[i]<>0 then
   divide(A_[i]*`eGCDpar/sigfact`(`eGCDpar/sigma`(a_[i],Var),l_[i]-2,Var),
                `eGCDpar/sigfact`(`eGCDpar/sigma`(b_[i-1],Var),l_[i]-2,Var),'temp');
   B_:=[op(B_),collect(temp,X)];
   b_:=[op(b_),lcoeff(B_[i],X)]
 fi
od;
B_:=subsop(1=NULL,2=NULL,B_);
B_
end:
#######################################
`eGCDpar/sigma`:=proc(b,Var)
subs(Var=Var+1,b)
end:
#######################################
`eGCDpar/sigfact`:=proc(b,n,Var)
local i,res,bs;

if n=0
then 1
else
  res:=1;
  bs:=b;
  for i from 0 to n-1 do
    res:=res*bs;
    bs:=`eGCDpar/sigma`(bs,Var)
  od;
  res
fi
end:
#######################################
`eGCDpar/sqrfree2`:=proc(p,v)
local t;

t:=sqrfree(p,v);
collect(mul(op([2,i,1],t),i=1..nops(op(2,t))),v)
end:
#######################################
`eGCDpar/filter`:=proc(d,X,Var,Par)
local p,q;

p:=op(1,d);q:=op(2,d);
if `eGCDpar/my_degree`(p,Par)=0 then RETURN(NULL)
elif q=0 then okay
else
  q:=`eGCDpar/primpart3`(`eGCDpar/aldiv`(q,content(lcoeff(q,X),Var),p,Par),X,[X,Var,Par])
fi;
if lcoeff(p,`eGCDpar/vlist`(p,[Par]))<0 then p:= -p fi;
if lcoeff(q,`eGCDpar/vlist`(q,[X,Var,Par]))<0 then q:= -q fi;
[`eGCDpar/canform`(p,[Par]),`eGCDpar/canform`(q,[X,Var,Par])]
end:
#######################################
`eGCDpar/vlist`:=proc(p,O)

[op(O),op(sort(convert(indets(p) minus convert(O,set),list),lexorder))]
end:
#######################################
`eGCDpar/primpart3`:=proc(p,X,l)

`eGCDpar/canform`(primpart(p,X),l)
end:
#######################################
`eGCDpar/canform`:=proc(p,O)
local vl;

vl:=`eGCDpar/vlist`(p,O);
sort(collect(expand(p),vl),vl,plex) end:
#######################################
`eGCDpar/aldiv`:=proc(c,q,p,Par)
local res;

if degree(p,Par)=0 then RETURN(c) fi;
res:=evala(Expand(subs(Par=RootOf(p,Par),c/q)));
if degree(p,Par)=1
then
  RETURN(res)
else
  RETURN(subs(RootOf(p,Par)=Par,res))
fi
end:
#######################################
`eGCDpar/my_skew_prem`:=proc(A,B,X,Var)
local m,n,i,res,T;

m:=degree(A,X); n:=degree(B,X);
if m<n then RETURN(A) fi;
T:=array(1..m-n+1);
T[m-n+1]:=B;
for i from m-n by -1 to 1
do
  T[i]:=`eGCDpar/e_skew_product`(X,T[i+1],X,Var);
od;
res:=A;
for i from 1 to m-n+1
do
  res:=expand(lcoeff(T[i],X)*res-coeff(res,X,m-i+1)*T[i])
od;
res
end:
#######################################
`eGCDpar/my_degree`:=proc(p,x)
local t;

t:=collect(p,x);
if t=0 then RETURN(-1)
else RETURN(degree(t,x)) fi
end:
#######################################
`eGCDpar/e_skew_product`:=proc(L1,L2,X,Var)
local i,p,res;

p:=L2; res:=0;
for i from 0 to degree(L1,X)
do
  res:=res+coeff(L1,X,i)*p;
  p:=`eGCDpar/eop`(p,X,Var)
od;
`eGCDpar/canform`(res,[X,Var])
end:
#######################################
`eGCDpar/eop`:=proc(L,X,Var)
local i,res,c;

res:=0;
for i from 0 to degree(L,X)
do
  c:=coeff(L,X,i);
  res:=res+`eGCDpar/sigma`(c,Var)*X^(i+1)
od;
res
end:
#######################################################################
#######################################################################
#
#    Title:     first_half
#    Created:   Summer 1997
#    Author:    Eugene Zima
#
###################################################################
#
#
# Temporal patch should be removed in Rmaple with changing mydegree
# to degree
#
#
`Slode/mydegree`:=proc(y,x) if y=0 then -infinity else degree(y,x) fi end:
#
#
`Slode/myldegree`:=proc(y,x) if y=0 then -infinity else ldegree(y,x) fi end:
#
#
###################################################################
#
`Slode/cor_sub` := 
  proc(expr, x, v, fname)

    if not has(expr,diff) and not has(expr,fname)
      then 
        subs(x=v,expr)
      elif
        type(expr, function) and (op(0, expr)='diff' or op(0, expr)=fname )
      then
        expr
      else
        op(0,expr)(seq(`Slode/cor_sub`(op(j,expr),x,v,fname),j=1..nops(expr)))
    fi
end:
#
#
`Slode/l_und_d` :=
  proc(L, x, y)
    local mt, l, d, w, is_infin, j;
#
  mt:=`Slode/DEdetermine`(L ,y(x) );
  l:=min(seq(j-mydegree(mt[coeffs][j],mt[var]) ,j=0..mt[ord]));
  w:=[seq(j-myldegree(mt[coeffs][j],mt[var]) ,j=0..mt[ord])];
  is_infin:=proc(u) evalb(u=infinity) end:
  d:=max(op(remove(is_infin,w)));
  [l,d]
end:
#
#
`Slode/m_split`:= 
  proc(Lw, m, x, y, l)
    local j,i,cf, mta, La, t;
#
  mta:=`Slode/DEdetermine`(Lw ,y(x) );
  La:=array(0..m-1,[seq(0,j=1..m)]);
  cf:=mta[_coeffs];
  for j from 0 to mta[ord] do
    for i from 0 to cf[j][ord] do
      t:= (j-i-l) mod m;
      La[t]:=La[t]+cf[j][i]*x^i*D^j
    od
  od;
  eval(La)
end:
#
#
`Slode/S_und_T`:=
  proc(M_s,x,a,m)
    local w, T, S, i, g, ww, t;
#
  S:=M_s[0];
  T:={};
#
  for t from 1 to m-1 while `Slode/mydegree`(S,D)<>0 do
    w:=`Slode/dGCDpar`(S,M_s[t],0,D,x,a);
    S:=op(nops(w),w)[2];
    T:={seq(w[i][1],i=1..nops(w)-1),op(T)};
    ww:=[coeffs(lcoeff(S,D),x)];
    g:=gcd(ww[1],0);
    for i from 2 to nops(ww) do
      g:=gcd(g,ww[i])
    od;
    if degree(g,a) <> 0
      then T:={g,op(T)};
    fi;
  od;
#
#
S,T
end:
#
#
`Slode/select_point`:=
  proc(u)
    local p,w,i;

  p:=0; w:=0;
  while w =0 do
    w:= map(abs, subs(a=p, u));
    w:= mul(w[i], i=1..nops(w));
    p:=p+1
  od;
#  print(p-1);
  p-1
end:
#
#
`Slode/first_half`:=
  proc(L,m,x,y,l,a)
    local Lw, M_s,w, ordL, a0, cfs, g, u, i, S, SS;

  ordL:=`Slode/DEdetermine`(L ,y(x) )[ord];
  if ordL = 0
    then
      1,{}
    else
      Lw:=`Slode/cor_sub`(L,x,x+a,y);
      M_s:=`Slode/m_split`(Lw,m,x,y,l);
      w:=`Slode/S_und_T`(M_s,x,a,m);
      if ordL=degree(w[1],D)
        then
          RETURN(L,w[2])
      fi;
      a0:=`Slode/select_point`(w[2]);
      SS:=subs(a=a0,w[1]);
      cfs:=[coeffs(SS,D)];
      g:=gcd(cfs[1],0);
      for i from 2 to nops(cfs) do
        g:=gcd(g,cfs[i])
      od;
      SS:=SS/g;
#
      cfs:=seq(coeff(SS,D,i),i=0..degree(SS,D));
      S:=cfs[1]*y(x);
      for i from 2 to nops([cfs]) do
        S:=S+cfs[i]*diff(y(x),x$(i-1))
      od;      
#
      u:=`Slode/first_half`(S,m,x,y,l,a);
      u[1],{op(subs(a=a-a0,u[2])),op(w[2])}
  fi
end:
#
#################################################################
`LREtools/ratpolysols` := proc()
local reqn, info, func_name, var,polys,i,neweq,newsols,n,r,d,
    dis1,h,A,l0,l,j,c,m,G,num,pbasis, opt, u;
option `Copyright 1995 by Jacques Carette, Waterloo Maple Software`;

# Do some argument checking and processing
    if type(args[-1],identical('output')={identical('gensol'),
	identical('basis'),identical('onesol')}) then
        opt := rhs(args[-1]);
        reqn := LREtools[REcreate](args[1..-2]);
    else
        opt := 'gensol';
        reqn := LREtools[REcreate](args)
    fi;
    info := eval(op(4,reqn));

    if not info[type]='linear' then ERROR(`Only linear equations are handled`) fi;
    if not type(info[vars],name) then 
	ERROR(`Only single equations are handled`) 
    fi;

    if select(type, info[coeffs][2..-1], Non(polynom(anything, info[vars])))<>
	[] then
        ERROR(`Only polynomial coefficients are handled`)
    fi;

    if not type(info[coeffs][1], polynom(anything, info[vars])) then
	ERROR(`Only polynomial inhomogeneities are handled`)
    fi;

    # Get function name and variable name
    func_name := info[functions];
    n := info[vars];
    r := info[coeffs];

    d := info[order];
    # if the solution is t_n/s_n, h-1 is an upper bound on k such that
    # degree(gcd(s_{n+k},s_n),n)>0.
    dis1 := LREtools[dispersion](r[2],r[d+2],n);
    if dis1=FAIL or dis1<d then 
        RETURN(LREtools[polysols](reqn, 'output'=opt)) 
    fi;
    h := dis1-d+1;
    userinfo(5,'rsolve',`value of h is `,h);

    # Find a recurrence relating (u_n,u_{n+h},u_{n+2h},...,u_{n+dh})
    # neglecting the inhomogeneous part.
    A:=array(1..d+1,1..d+1,sparse);
    for i from 0 to d do A[i+1,d+1]:=X^i od;
    for i from 0 while i*h<d do A[i+1,i*h+1]:=1 od;
    l0:=[seq(-r[i]/r[d+2],i=2..d+1)];
    l:=l0;
    for i from d to d*h do
        if irem(i,h,'p')=0 then
            for j to d do A[p+1,j]:=l[j] od;
            if p=h then break fi;
        fi;
        c:=normal(subs(n=n+1,l[d]));
        l:=normal([c*l0[1],seq(normal(subs(n=n+1,l[j-1])+c*l0[j]),j=2..d)]);
    od;
    # this gives the equation:
    userinfo(5,'rsolve',`matrix built `,op(A));
    l:=collect(numer(normal(mygausselim(A,d+1,d+1))),X);
    # its order:
    m:=degree(l,X);
    # G is the potential denominator of the solution
    G:=primpart(normal(subs(n=n-d*h,coeff(l,X,m))));
    for i to m while degree(G,n)>0 do
        G:=primpart(gcd(G,normal(subs(n=n-(i-1)*h,coeff(l,X,i-1))))) od;
    if degree(G,n)=0 then RETURN(LREtools[polysols](reqn, 'output'=opt)) fi;

    # Build the equation satisfied by the numerator
    # don't worry about the initial conditions now
    # with the new initial conditions!
    l := normal(`+`(seq(r[i+2]/subs(n=n+i,G)*u(n+i),i=0..d)));
    l:=LREtools[REcreate](denom(l)*r[1]+numer(l),u(n),{});
    # old code for initials conditions:
    #{seq(op(1,i)=op(2,i)/subs(n=op(op(1,i)),G), i=op(3,reqn))});

    num := LREtools[polysols](l, 'output'=`if`(opt='gensol',basis,opt));
    if num=NULL then 
	userinfo(2,{rsolve,LREtools}, `Warning: no rational solutions found`);
	RETURN(NULL) 
    elif opt='gensol' then
	`LREtools/resolveinits`(reqn, pol);
    else 
	RETURN(num/G) 
    fi
end:
#####################################################################
# INPUT :
# eq = p[n](x) y(x+n) + ... + p[1](x) y(x+1) + p[0](x) y(x) = b(x)
#  -- a linear recurrence with polynomial coefficients 
#   p[k](x) - polynomial over a simple algebraic extension of the 
#             rational number field  
# y -- an unknown function in x
# OUTPUT : 
#    - the rational solution and set of constants in form _C[i] :
#                                [y(x),{...} ]
#    - NULL iff no rational solution
`solrec/ratsolrec` := proc( eq, y ) 
local neq, cf, n, denominator,numerator, z, x, Ct ;
global _C ;
  if nargs = 3 then Ct := args[3] 
	       else Ct := _C 
  fi ;

  x := op(1,y);
  neq := subs( x=x-`solrec/recurr_tdegree`(eq,y), eq ) ;
  
  cf := `solrec/recurr_coeffs`( neq, y );
  n := nops(cf) ;

  denominator := `solrec/DENOM`( cf[n], cf[1], n-1, x ) ;

  if denominator <> 1 then
    neq := `solrec/subs_to_recurrence`( neq, y, z(x)/denominator ) ;
    if type( neq, `=` ) then neq := lhs(neq) - rhs(neq) fi ;
    neq := numer( normal(neq) ) ;
  else
    neq := `solrec/subs_to_recurrence`( neq, y, z(x) );
  fi ;

  numerator := `solrec/polsolrec`( neq, z(x),  Ct ) ;
  if numerator = NULL then NULL
  elif denominator <> 1 then
	simplify(numerator/denominator);
#       [collect(simplify(numerator[1]/denominator),numerator[2]),
#	    numerator[2] ];
  else  numerator fi ;
end :
#####################################################################
##
## DENOM
##
#  divisors( x )
#  INPUT: x - a positive integer
#  OUTPUT: [x, ..., k, ..., 1] where k - divisor of x
`solrec/divisors` := proc( x ) local div, y, k ;
    div := [] ;
    if x mod 2 = 0 then y := x/2 else y := (x-1)/2 fi ;
    for k to y do if x mod k = 0 then
	 div := [k, op(div) ] ;
    fi od ;
    RETURN( [x,op(div) ] );
end :
# param
# INPUT:
# OUTPUT: hs - array of h such that gcd(A,B(x+h)) is nontrivial
#                  h_1 < h_2 < ...
#         ds - array of gcd(A(x),B(x+h_i)) where h_i in hs
`solrec/param` := proc( A, B, x ) local s, r, h, i, hs, ds ;
   A ;
   subs(x=x+h,B);
   s := resultant( A, subs(x=x+h,B), x ) ;
   hs := [`solrec/nonnegintsolve_a`( s, h )] ;
   hs := sort( hs ) ;
   ds := [seq( evala(Gcd(A,subs(x=x+i,B))),i=hs ) ] ;
   [hs,ds] ;
end :
# DENOM( an, a0, n, x )
# INPUT: an - a leading coefficient of the difference operator
#             ( polynomial in x )
#        a0 - a trailing coefficient (polynomial in x)
#        n - a degree of the operator
#        x - a variable
# OUTPUT: u - polynomial which may be a denominator of a rational
#             solution of the recurrence equation
`solrec/DENOM` := proc( an, a0, n, x ) local u, A, B, c, m, k, arr_h, arr_d,i ;
   A := subs( x=x-n, an ) ;
   B := a0 ;
   u := 1 ;
   arr_h := `solrec/param`( A, B, x ) ;
   arr_d := arr_h[2] ;
   arr_h := arr_h[1] ;
#   for k to nops(arr_h) do
    for k from nops(arr_h) by -1 to 1 do
      c := evala(Gcd( A, arr_d[k] )) ;
      c := evala(Gcd( B, subs(x=x-arr_h[k], c))) ;
      B := evala(Quo(B,c,x)) ;
      c := subs(x=x+arr_h[k],c) ;
      A := evala(Quo(A,c,x)) ;
      u := u * product(subs(x=x-i, c),i=0..arr_h[k]) ;
   od ;
   RETURN(u) ;
end :
######################################################################
# INPUT :
# eq = p[n](x) y(x+n) + ... + p[1](x) y(x+1) + p[0](x) y(x) = g(x)
#  -- a linear recurrence with polynomial coefficients
#   p[k](x) - polynomial over a simple algebraic extension of the 
#             rational number field  
# y -- an unknown function in x
# OUTPUT : 
#    - the polynomial solution and set of constants in form _C[i] :
#                                [y(x),{_C[0],...} ]
#    - NULL iff no polynomial solution
`solrec/polsolrec` := proc( eq, y ) 
local neq,cf, g, b, n, x, N, slv, j, k, s, h, t, _n, E, slv1, Ct, 
      res, sbs, i ;
global _C ;

  if nargs = 3 then Ct := args[3] 
	       else Ct := _C 
  fi ;

  x := op(1,y);
  neq := subs( x=x-`solrec/recurr_tdegree`(eq,y), eq ) ;

  cf := `solrec/recurr_coeffs`( neq, y, g ) ;

  cf := `solrec/recurr_to_diffc`( cf ) ;
  cf := `solrec/recurr_for_diffc`( cf, x, n ) ;
  g := `solrec/ls`( g, x ) ;

  b := `solrec/bound_index`(cf) ;
  N := `solrec/dg_bound`( cf[b], nops(g)-1, b, n ) ;

  slv := `solrec/solve_recurr`( eval(cf), g, n, N, Ct ) ;

  if slv = NULL then RETURN( NULL ) fi ;

  slv := `solrec/next_el_is_null`(slv,eval(cf),g,n,N);
  if slv = NULL then RETURN( NULL ) fi;

  t := nops(slv[3]) ;
  for j to t do
       s[j] := collect(simplify(add(
	  convert(slv[2][j][k+1]*binomial(x,k),factorial),
			   k=0..N)), x ) ;
  od ;
  h := collect(simplify(add(
	  convert(slv[1][k+1]*binomial(x,k),factorial),
			   k=0..N)), x ) ;

  res := add(slv[3][j]*s[j],j=1..t)+h ;

  collect(res,x);
end :
######################################################################
`type/recurrence` := proc( L, y ) local nL ;
   if type(L,equation) then nL := lhs(L) - rhs(L) ;
		       else nL := L fi ;
   if type(nL,algebraic) and type(y, anyfunc(name)) then
      type(nL,mylinear(map(`solrec/isshift`,indets(nL),y))) fi :

end :
`type/difference` := proc( L, y ) local nL ;
   if type(L,equation) then nL := lhs(L) - rhs(L) ;
		       else nL := L fi ;
   if type(nL,algebraic) and type(y, anyfunc(name)) then
      type(nL,mylinear(map(`solrec/isdelta`,indets(nL),y))) fi :

end :
`type/mylinear` := (expr, x) -> type(expr,linear(x)) or
				numboccur(expr,x) = 0 :
######################################################################
# INPUT : p 
# OUTPUT : true iff p is x+k 
`solrec/isshiftarg` := ( p, x ) ->
   (type(p,`+`) or type(p,`-`)) and 
        lcoeff(p,x)=1 and type(tcoeff(p,x),integer) :
# INPUT : z
# OUTPUT : true iff z is y(x+k)
`solrec/isshift` := ( z, y ) ->
   if (z = y) or type(z,op(0,y)(anything)) and `solrec/isshiftarg`(op(1,z),op(1,y))
	then z  fi :
# INPUT : y(x) - any function
#         i - non-negative integer
# OUTPUT : y(x+i)
`solrec/eval_shift` := ( y, i ) -> (op(0,y))(op(1,y)+i) :
#
# L - `p[n](x) y(x+n) + ... + p[1](x) y(x+1) + p[0](x) y(x) = b(x)`
# ==> n
`solrec/recurr_degree` := proc( L, y ) local shiftvars, nL, k, x ;
 if not type( L, recurrence(y) ) then
   ERROR( `1st argument must be a linear recurrence`, L ) fi ;

 if type(L,`=`) then nL := lhs(L) - rhs(L) else nL := L fi ;
   shiftvars := map( `solrec/isshift`, indets(nL), y ) ;
 if shiftvars = {} then
     ERROR( `The recurrence `, nL,` mast be in `, y ) fi ;

 x := op(1,y) ;
 max(seq(coeff(op(1,k),x,0), k=shiftvars)) ;

end :
#
# L - `p[n](x) y(x+n) + ... + p[1](x) y(x+1) + p[0](x) y(x) + ... + p[-m] y(x-m) = b(x)`
# ==> n
`solrec/recurr_ldegree` := ( L, y ) -> `solrec/recurr_degree`( L, y ) :

# ==> m
`solrec/recurr_tdegree` := proc( L, y ) local nL, shiftvars, x, k; 
 if not type( L, recurrence(y) ) then
   ERROR( `1st argument must be a linear recurrence`, L ) fi ;

 if type(L,`=`) then nL := lhs(L) - rhs(L) else nL := L fi ;
   shiftvars := map( `solrec/isshift`, indets(nL), y ) ;
 if shiftvars = {} then
     ERROR( `The recurrence `, nL,` mast be in `, y ) fi ;

 x := op(1,y) ;
 min(seq(coeff(op(1,k),x,0), k=shiftvars)) ;
end :

#
# eq = p[n](x) y(x+n) + ... + p[1](x) y(x+1) + p[0](x) y(x) + b(x)
# => [p[n](x), ..., p[1](x), p[0](x)], g := b(x) ;
#
`solrec/recurr_coeffs` := proc( eq, y, g ) local n, m, neq, shiftvars, i, t ;
 n := `solrec/recurr_degree`( eq, y ) ;
 m := `solrec/recurr_tdegree`( eq, y ) ;
 if type(eq,`=`) then neq := lhs(eq) - rhs(eq) else neq := eq fi ;
 shiftvars := map( `solrec/isshift`,indets(neq),y ) ;
 neq := collect(neq,shiftvars) ;
 if nargs = 3 then
    g := -subs( {seq( t=0, t=shiftvars )}, neq ) fi ;

 [seq( coeff(neq,`solrec/eval_shift`(y,i)), i = m..n )] ;
end :
#
# eq = p[n](x) y(x+n) + ... + p[1](x) y(x+1) + p[0](x) y(x) + b(x)
# => p[n](x)
#
`solrec/recurr_lcoeff` := proc( eq, y, g ) local cf ;
 cf := `solrec/recurr_coeffs`( eq, y ) ;
 cf[nops(cf)] 
end :
`solrec/recurr_rhs` := proc( eq, y ) local n, neq, shiftvars, t ;
 n := `solrec/recurr_degree`( eq, y ) ;
 if type(eq,`=`) then neq := lhs(eq) - rhs(eq) else neq := eq fi ;
 shiftvars := map( `solrec/isshift`,indets(neq),y ) ;

 -subs( {seq( t=0, t=shiftvars )}, neq ) ;
end :
#   max(0,-b-1, deg_f-b, seq(rhs(op(k)),k={isolve(a_b, n)}) )
# y = y(x+i) ==> 'y(x+i) = u(x+i)'
#
`solrec/subs_to_shift` := ( y, u ) ->
      y = subs( op(indets(op(1,y))) = op(1,y), u ) :
# L is a linear recurrence of y(x)
# ==> "subs( y(x) = u(x), L )" ;
`solrec/subs_to_recurrence` := proc( L, y, u ) local rL, lL, shiftvars ;
 if not type( L, recurrence(y) ) then
   ERROR( `1st argument must be a linear recurrence`, L ) fi ;

 if type(L,`=`) then
     lL := lhs(L) ;
     shiftvars :=map(`solrec/isshift`,indets(lL),y) ;
     lL := simplify(subs( map( `solrec/subs_to_shift`, shiftvars, u ), lL )) ;

     rL := rhs(L) ;
     shiftvars :=map(`solrec/isshift`,indets(rL),y) ;
     rL := simplify(subs( map( `solrec/subs_to_shift`, shiftvars, u ), rL )) ;

     lL = rL ;
 else
     shiftvars :=map(`solrec/isshift`,indets(L),y) ;
     simplify(subs( map( `solrec/subs_to_shift`, shiftvars, u ), L )) ;
 fi ;

end :
######################################################################
# L - a linear recurrence of y(x)
# => equal difference operator
`solrec/recurrence_to_difference` := proc( L, y ) local nL, b, n ;
 nL := `solrec/recurr_coeffs`( L, y, b ) ;
 nL := collect(`solrec/recurr_to_diffc`( nL ),op(1,y)) ;
 nL := subs(Delta(op(0,y),op(1,y),1)=Delta(op(0,y),op(1,y)),
      add( Delta(op(0,y),op(1,y),n-1)*nL[n], n=2..nops(nL) ) + y*nL[1]) ;
 if type(L,`=`) then nL = b else nL - b fi ;
end :
# L - coefficients of a linear recurrence
# => coefficients of a equal difference operator
`solrec/recurr_to_diffc` := proc( L ) local n, k, L1, L2 ;
   n := nops(L)-1 ;
   if n = 0 then L
   else
      L1 := [seq( L[n+1]*binomial(n,k), k=0..n )] ;
      L2 := `solrec/recurr_to_diffc`([seq(L[k],k=1..n)]) ;
      [seq(L1[k]+L2[k], k=1..n),L1[n+1]]
   fi ;
end :
#
######################################################################
#
# z is Delta(y,x,k), k >= 0
#
`solrec/isdelta` := ( z, y ) ->
   if (z = y) or
      (type(z,Delta(name,name)) or type(z,Delta(name,name,integer)))
      and op(1,z)=op(0,y) and op(2,z)=op(1,y) then z  fi :
#
#  y(x) --> Delta(y,x,i)
#
`solrec/eval_delta` := ( y, i ) -> if i>1 then Delta(op(0,y),op(1,y),i)
			  elif i=1 then Delta(op(0,y),op(1,y))
			  else y fi :
#
# L = `p[n](x) Delta(y,x,n) + ... + p[1](x) Delta(y,x) + p[0](x) y(x) = b(x)`
# ==> n
`solrec/diffc_degree` := proc( L, y ) local deltavars, nL, k, x ;
 if not type( L, difference(y) ) then
   ERROR( `1st argument must be a linear difference equation`, L ) fi ;

 if type(L,`=`) then nL := lhs(L) - rhs(L) else nL := L fi ;
   deltavars := map( `solrec/isdelta`, indets(nL), y ) ;
 if deltavars = {} then
     ERROR( `The operator `, nL,` mast be in `, y ) fi ;

 x := op(1,y) ;
 max(seq(`solrec/delta_degree`(k), k=deltavars)) ;

end :
#
# z = Delta(y,x,n) or Delta(y,x) or y(x) == n or 1 or 0
#
`solrec/delta_degree` := (z) -> if type(z,Delta(name,name)) then 1
		       elif type(z,Delta(name,name,integer)) then op(3,z)
		       else 0 fi :
#
# eq = p[n](x) Delta(y,x,n) + ... + p[1](x) Delta(y,x) + p[0](x) y(x) + b(x)
# => [p[n](x), ..., p[1](x), p[0](x)], g := b(x) ;
#
`solrec/diffc_coeffs` := proc( eq, y, g ) local n, neq, deltavars, i, t ;
 n := `solrec/diffc_degree`( eq, y ) ;
 if type(eq,`=`) then neq := lhs(eq) - rhs(eq) else neq := eq fi ;
 deltavars := map( `solrec/isdelta`,indets(neq),y ) ;
 neq := collect(neq,deltavars) ;
 if nargs = 3 then
    g := -subs( {seq( t=0, t=deltavars )}, neq ) fi ;

 [seq( coeff(neq,`solrec/eval_delta`(y,i)), i = 0..n )] ;
end :
`solrec/diffc_rhs` := proc( eq, y ) local n, neq, deltavars, t ;
 n := `solrec/delta_degree`( eq, y ) ;
 if type(eq,`=`) then neq := lhs(eq) - rhs(eq) else neq := eq fi ;
 deltavars := map( `solrec/isdelta`,indets(neq),y ) ;

 -subs( {seq( t=0, t=deltavars )}, neq ) ;
end :
#
#
#
`solrec/subs_to_difference` := proc( L, y, u ) local cf, res,  delta_u, i, g ;
   cf := `solrec/diffc_coeffs`( L, y, g ) ;
   delta_u := `solrec/deltas`( u, op(1,y), nops(cf)-1 ) ;

   res := 0 ;
   for i to nops(cf) do
	 res := simplify(cf[i]*delta_u[i] + res) ;
   od ;

   if type( L, `=` ) then collect(res,x)=collect(g,x)
		     else collect(res-g,x)  fi ;
end :
######################################################################
#
# L = p[n](x) Delta(y,x,n) + ... + p[1](x) Delta(y,x) + p[0](x) y(x)
#        
# y(x) = ly(0)P[0] + ly(1)P[1] + ... + ly(n)P[n] +...
# where P[n] = binomial(x,n)
#
# OUTPUT : recurrence for {y[n]}
#     alpha[-A](n+A)ly(n+A)+...+alpha[0](n)ly(n)+...+alpha[b](n-b)ly(n-b)
`solrec/recurrence_for_difference` := proc( L, y, ly ) 
local cf, n, Alph, A,b, g, i ;

  cf := `solrec/diffc_coeffs`( L, y, g ) ;
  if type(L,`=`) or g <> 0 then  
	ERROR( `1st argument must be a linear difference operator`, L ) 
  fi ;
  if not type(ly,anyfunc(name)) then
	ERROR( `3st argument is not right`, ly ) ;
  fi ;
  n := op(ly) ;
  Alph :=  `solrec/recurr_for_diffc`( cf, op(1,y), n ) ;
	# Alph=array[-A..b, [alpha[-A](n),...,alpha[b](n)]]

  A := -`solrec/near_index`(Alph) ;
  b := `solrec/bound_index`(Alph) ;
  add(subs(n=n-i,Alph[i])*op(0,ly)(op(1,ly)-i),i=-A..b) ;
end :
# cf - coefficients of a difference operator at x
# ==> array of coefficients of recurrence at n
# [-A .. b, [alpha(n),...,alpha(n)]]
#
`solrec/recurr_for_diffc` := proc( cf, x, n) 
local alp, r, d, ncf, p, i, k, j, tmp, tmp1, Arr;
   r := nops(cf) - 1 ;
   d := max( seq(degree(p,x), p = cf )) ;
   ncf := array( 0..r, 0..d, map(`solrec/ls`,cf,x,d) ) ;

   Arr := table([(-r,0)=1]);
   alp := array( -r..d, sparse);
   alp[-r] := ncf[r,0];
   for j to d do
	Arr[-r,j] := expand(Arr[-r,j-1]*(n-r-j+1)/j);
	alp[-r]:=alp[-r]+Arr[-r,j]*ncf[r,j];
   od;
   for i from -r+1 to d do
     Arr[i-1,-1] := 0;
     for j from 0 to d do
	Arr[i,j] := Arr[i-1,j-1] + Arr[i-1,j];
	tmp1 := binomial(j,max(0,i));
	tmp := 0;
	for k from max(0,-i) to min(r,j-i) do
	    tmp := tmp + tmp1*ncf[k,j];
	    tmp1 := tmp1*(j-(i+k))/(i+k+1);
	od;
        alp[i]:=alp[i]+Arr[i,j]*tmp;
     od;
   od;
  
   while alp[d] = 0 do d:=d-1 od ;
   array(-r..d,[seq(alp[i], i=-r..d )] ) ;
end :
#
# p -- polynomial degree d in x ==> [ p,Delta(p,x),...,Delta(p,x,d) ]
#
`solrec/deltas` := proc( p,x ) local d, i, cf ;
   if nargs = 2 then d := degree(p,x)
		else d := args[3] ;
   fi ;
   cf := array(0..d,[p]) ;
   for i to d do
      cf[i] := collect(subs( x= x+1, cf[i-1] ) - cf[i-1], x )
   od ;
   convert(cf,list) ;
end :
#
# ๐ ญญ polynomial in x
#
`solrec/ls` := ( p, x ) -> subs( x=0, `solrec/deltas`(args) ) :
######################################################################
# a[b](n)*y(n-b)+...+a[-i](n)*y(n+i)+..+a[-A](n)*y(n+A) = f(n)
#     a[-A](n+A)ly(n+A)+...+a[0](n)ly(n)+...+a[b](n-b)ly(n+b)
# ==> solutions of recurrence and set of parameters C[i_1],...,C[i_t]
# [[y(0),...,y(M)],[C[i_1],...,C[i_t]]
`solrec/solve_recurrence` := proc( rec, f, y, M, C ) 
local b, A, cf, g, i, n, slv, k ;
	A := `solrec/recurr_ldegree`( rec, y ) ;
	b := -`solrec/recurr_tdegree`( rec, y ) ;
	cf := `solrec/recurr_coeffs`( rec, y, g ) ;
	if g <> 0 then ERROR() fi ;
	# cf = [a[b](n-b),...,a[-A](n+A)]

	n := op(y) ;
	slv := `solrec/solve_recurr`(array( -A..b, 
		[seq(collect(subs(n=n+i,cf[-i+b+1]),n),i=-A..b)]),
		f, op(y), M, C ) ;
	[[seq(slv[1][k]+
	         add(slv[2][i][k]*slv[3][i],i=1..nops(slv[3])),
              k=1..M+1)],
         {op(slv[3])} ] ;
end :
# a[b](n-b)*y(n-b)+...+a[-i](n+i)*y(n+i)+..+a[-A](n+A)*y(n+A) = f(n)
# ==> solutions of recurrence
# INPUT : cf = array( -A..b, [a[-A],...,a[b]] )
#         f = [ f(0), f(1), ... ]
# OUTPUT :
#  g =   [ g[0], ..., g[M] ]
#  V = [ [v1[0], ..., v1[M],
#    ... [vt[0], ..., vt[M] ]
#  L = [ _C[i1], ..., _C[it] ]
`solrec/solve_recurr` := proc( cf, f, n, M, _C )
local A, b, S, k, j, _n, t, alpha, V, L, E, e_n, g, lf, v,
      i, m, _Ck ;

   b := `solrec/bound_index`( cf ) ;
   A := -`solrec/near_index`( cf ) ;

   if type(op(2,cf),range) then alpha := table(sparse,op(3,cf));
			   else alpha := table(sparse,op(cf));
   fi;

   g := table(sparse) ;

   lf := table(sparse, [seq(k=f[k+1],k=0..nops(f)-1)] ) ;

   E:= {} ;

   t := min(A-1,M)+1 ;
   V := [seq(table(sparse,[k=1]),k=0..t-1)];
   L := [seq(_C[k],k=0..t-1)];

   for _n from min(A,M+1) to M do
      _n ;
      S := subs(n=_n,alpha[-A]);
      if  not testeq(S) then
	for v in V do
	  v[_n]:= simplify(-( add(subs(n=_n-i,alpha[i-A])*v[_n-i],i=1..A+b))/
			    S);
	od ;
	g[_n] :=simplify((lf[_n-A]-add(subs(n=_n-i,alpha[i-A])*g[_n-i],i=1..A+b)) /
			S) ;
      else
	t := t+1 ;
	V := [op(V),table(sparse,[_n = 1])] ;
	L := [op(L), _C[_n]] ;
	e_n := simplify(add(subs(n=_n+i-A,alpha[-i])*(g[_n+i-A]+
			 add(L[j]*V[j][_n+i-A],j=1..t-1)),i=-b..A)
		  -lf[_n-A]) ;
	if e_n <> 0 then
	      if indets(e_n) intersect {op(L)} = {} 
		then RETURN( NULL ) fi ;
	      i := max(seq(op(k),k=indets(e_n) intersect {op(L)})) ;
	      member( _C[i],L,'k') ;
	      _Ck := collect(simplify(-coeff(e_n,L[k],0)/coeff(e_n,L[k],1)),L);
	      for j from i to _n-1 do
		  for m to k-1 do
		     V[m][j] := V[m][j] + V[k][j]*coeff(_Ck,L[m]) ;
		  od ;
		  g[j] := g[j] + V[k][j]*subs(seq(v=0,v=L),_Ck) ;
	      od ;
	      V := [op(V[1..k-1]),op(V[k+1..t])] ;
	      L := [op(L[1..k-1]),op(L[k+1..t])] ;
	      t := t-1 ;
	fi ;
      fi ;
    od ; k := 'k';
    [ [seq(g[k],k=0..M)],
      [seq([seq(V[i][k],k=0..M)],i=1..t)],
      L ]
end :
#############################################################################
`solrec/next_el_is_null` := proc(S,cf,f,n,N)
local A, b, k, j, _n, t, alpha, V, L, E, e_n, g, lf, v,
      i, m, _Ck ;

   b := `solrec/bound_index`( cf ) ;
   A := -`solrec/near_index`( cf ) ;

   if type(op(2,cf),range) then alpha := table(sparse,op(3,cf));
			   else alpha := table(sparse,op(cf));
   fi;

   g := table(sparse,[seq(i=S[1][i+1],i=0..N)]) ;
   lf := table(sparse, [seq(k=f[k+1],k=0..nops(f)-1)] ) ;

   E:= {} ;

   t := nops(S[3]) ;
   V := [seq(table(sparse,[seq(i=S[2][k+1][i+1],i=0..N)]),k=0..t-1)];
   L := S[3];

   for _n from max(A,N+1) to N+A+b do
      _n ;
	e_n := simplify(add(subs(n=_n+i-A,alpha[-i])*(g[_n+i-A]+
			 add(L[j]*V[j][_n+i-A],j=1..t)),i=-b..A)
		  -lf[_n-A]) ;
	if e_n <> 0 then
	      if indets(e_n) intersect {L[]} = {} 
		then RETURN( NULL ) fi ;
	      i := max(seq(op(k),k=indets(e_n) intersect {op(L)})) ;
	      member( _C[i],L,'k') ;
	      k;
	      _Ck := collect(simplify(-coeff(e_n,L[k],0)/coeff(e_n,L[k],1)),L);
	      for j from i to N do
		  for m to k-1 do
		     V[m][j] := V[m][j] + V[k][j]*coeff(_Ck,L[m]) ;
		  od ;
		  g[j] := g[j] + V[k][j]*subs(seq(v=0,v=L),_Ck) ;
	      od ;
	      V := [op(V[1..k-1]),op(V[k+1..t])] ;
	      L := [op(L[1..k-1]),op(L[k+1..t])] ;
	      t := t-1 ;
	fi ;
    od ; k := 'k';
    [ [seq(g[k],k=0..N)],
      [seq([seq(V[i][k],k=0..N)],i=1..t)],
      L ]
end :
#############################################################################
# arr -- array( n..m, [...] ) or [ ... ]
# ==> m  - bound index
#     n or 1 - near index
`solrec/bound_index` := ( arr ) ->
	       if type(arr,vector) then nops(arr)
		   elif type(arr,array) then op(2,op(2,eval(arr)))
		       fi :
`solrec/near_index` := ( arr ) ->
	       if type(arr,vector) then 1
		   elif type(arr,array) then op(1,op(2,eval(arr)))
		       fi :

#####################################################################
# INPUT :
# eq = p[n](x) y(x+n) + ... + p[1](x) y(x+1) + p[0](x) y(x) = g(x)
#  -- a linear recurrence with polynomial coefficients
#   p[k](x) - polynomial over a simple algebraic extension of the 
#             rational number field  
# y -- an unknown function in x
# OUTPUT : 
#    - a degree bound for polynomial solutions of eq
`solrec/degree_bound` := proc(eq,yx)
local cf, g, b, k, j, alpha_b, l_cf, r, d, n;
    cf := `solrec/recurr_coeffs`(eq,yx,g);
    cf := `solrec/recurr_to_diffc`(cf);
    r := nops(cf)-1;
    d := max(seq(degree(cf[k+1],op(yx)),k=0..r));
    b :=max(seq(degree(cf[k+1],op(yx))-k,k=0..r));
    l_cf := array( 0..r, 0..d, map(ls,cf,x,d));
    alpha_b := add(add(simplify(convert(binomial(n+b,j),factorial))
                       *binomial(j,b+k)*l_cf[k,j],
                       j=0..d),k=0..r);
    `solrec/dg_bound`( expand(alpha_b), degree(g,op(yx)), b, n);
end:
#
#
`solrec/dg_bound` := proc( a_b, deg_f, b, n ) local k;
   max(0,-b-1, deg_f-b, `solrec/nonnegintsolve_a`(a_b, n) )
end :
######################################################################
# P = a[n] x^n + ... + a[1] x + a[0],
#   a[k] = b[m-1] r^(m-1) + ... + b[1] r + b[0], k = 0..n
#      where r is an algebraic number (irreducible polynomial)
# OUTPUT : nonnegative integer roots of P  r1, ..., rk 
`solrec/nonnegintsolve_a` := ( P, x ) ->
  select( type, map2(op,1,roots(P,x)),
		     nonnegint)[]:
######################################################################