(* Multivariate linear recurrences with constant coefficients for Mma 3.0 *)
(* by Marko Petkovsek *)
(* August 22, 1999 *)

BeginPackage["Multivariate`",{"DiscreteMath`RSolve`"}];

Off[General::spell1];


DirectionVector::usage = "DirectionVector[H] finds the computational direction
vector for a set of shifts H, provided that the convexity condition is 
satisfied.";

Apex::usage = "Apex[H] finds the apex for a set of shifts H.";

MakeList::usage = "If the head of expr is h then MakeList[expr,h] replaces it by 
List, otherwise it returns {expr}.";

DiscreteInitialValueProblem::usage = "DiscreteInitialValueProblem[a[i1,12,...,ik]
== rhs, initconds, s] finds the solution generating function of the given linear
recurrence with constant coeficients (valid from the starting point s on) and 
initial conditions provided that the convexity condition is satisfied and the 
apex has at most one positive coordinate.";

FirstTerm::usage = "FirstTerm[f, {x, x0}] finds the lowest-order term(s) in the 
expansion of f around x = x0.";

Subsets::usage = "Subsets[l] generates the set of all subsets of l.";

x::usage = "x[1], x[2], ..., x[k] are the variables of the generating function
returned by DiscreteInitialValueProblem.";


Begin["Private`"];

Share[];

Unprotect[Power];
0^0 = 1;
Protect[Power];

DirectionVector[H_]:=Module[{n=Length[First[H]],out},
		
		(* Checks the convex-
        hull condition and finds a suitable direction vector pointing into 
          the first orthant along which recursive computation can proceed *) 

		Off[ConstrainedMin::nsat];
		out=LinearProgramming[Table[0,{n}],Join[-H,IdentityMatrix[n]],
        Table[1,{Length[H]+n}]];On[ConstrainedMin::nsat];
    Return[If[Head[out]===LinearProgramming,False,out]]];

Apex[H_]:=Max/@Thread[Append[H,Table[0,{Length[First[H]]}]]];

MakeList[expr_,head_]:=If[Head[expr]===head,List@@expr,List@expr];

MakeRule[a_[index__]==y_,f_]:=
  Evaluate[f@@(If[Head[#]===Symbol,Pattern[#,Blank[]],#]& /@ {index})]:=y;

DiscreteInitialValueProblem[a_[index__]==rhs_,initconds_,start_]:=
  Module[{ind={index},n,coeffs,shifts,v,p,Q,K,U,vars,phi,t,pt,roots},
		n=Length[ind];
		terms=MakeList[rhs,Plus];
		vars=x/@Range[n];
		coeffs=terms/._a->1;
		shifts=terms/.c_. a[z__]-> {z};
		shifts=shifts/._?(MemberQ[ind,#]&)->0;
		v=DirectionVector[shifts];
		If[v===False,Print["Convexity condition not satisfied"];Return[$Failed],
          Print["Computational direction vector: ",v]];
		p=Apex[shifts];
		Print["Apex: ", p];
		Q=Factor[Times@@(vars^p)-coeffs.(Times@@(vars^(p-#))&)/@shifts];
		Print["Characteristic polynomial: ",Q];
		ClearAll[phi];
		MakeRule[#,phi]& /@ initconds;
		K=Simplify[
        Times@@(vars^(-start)) (
            coeffs.(Times@@(vars^(p-#))Ph[start,#,phi]&)/@shifts)]/.phi->a;
		Print["Known initial function: ",K];
	    If[Count[p,0]<n-1,Print["APEX HAS MORE THAN ONE POSITIVE COORDINATE"]; 
		  Return[$Failed]];
		
		(* the rational case *)
		
		If[p===Table[0,{n}], Print["RATIONAL CASE"];
          Return[Simplify[Times@@(vars^start)K/Q +Ph[start,-start,phi]/.phi->a]]];
		
		(* the algebraic case *)
		
        Print["ALGEBRAIC CASE"];
        t=First[First[Position[p,_?Positive]]];
		pt=p[[t]];
		others=Complement[vars,{x[t]}];
		roots=Simplify[x[t]/.Solve[Q==0,x[t]]];
			Off[Power::infy,Infinity::indet];
		roots=Select[roots,
        Function[
          u,((Fold[FirstTerm,u,Thread[{vars,Table[0,{n}]}]]
		  /.Thread[vars->0])) == 0]];
		On[Power::infy,Infinity::indet];	  
		values=Simplify[K/.x[t]->#& /@ roots];
		U=Simplify[InterpolatingPolynomial[Thread[{roots,values}],x[t]]];
		Return[Simplify[(Times@@(vars^start)(K-U)/Q +Ph[start,-start,phi]/.
		phi->a) ]]
	];

FirstTerm[f_,{x_,x0_}]:=
  Module[{t=Normal[Series[f,{x,x0,0}]]},
    If[t=!=0,t,Expand[FirstTerm[f/x,{x,x0}]x]]];

Subsets[x_List]:=Distribute[{{},{#}}& /@ x,List,List,List,Union];

Ph[s_,h_,phi_]:=Plus@@(MakeSum[s,h,phi,#]&/@Rest[Subsets[Range[Length[s]]]]);

MakeSum[s_,h_,phi_,subset_]:=
  Module[{n=Length[s],dvars,vars,i,monomial,lower,upper,iterators},
		vars=x/@Range[n];
		dvars=i/@Range[n];
		monomial=Times@@(vars^dvars);
		lower=s+h;
		upper=MapIndexed[If[MemberQ[subset,First[#2]],#1,Infinity]&, s-1];
        iterators=Thread[{dvars,lower,upper}];
		(-1)^(Length[subset]-1)Fold[MySum,phi@@dvars monomial,iterators]//.
		PowerSum[
          expr_,{x_[j_],n_,n0_}]:>
        Module[{z},PowerSum[expr/.x[j]->z,{z,n,n0}]/.z->x[j]]];

MySum[f_,{i_,a_Integer,b_Integer}]:=Sum[f,{i,a,b}];
MySum[f_,{i_[k_],a_Integer,Infinity}]:=Module[{ik},PowerSum[f/x[k]^ik/.i[k]->ik,{x[k],ik,a}]];

End[];

EndPackage[];


On[General::spell1];
