(* :Title: Symbolics *)

(* :Context: Symbolics` *)

(* :Author: Eric Rowland *)
(*          http://math.tulane.edu/~erowland *)

(* :Date: {2009, 9, 20} *)

(* :Package Version: 1.01 *)

(* :Mathematica Version: 7.0 *)

(* :Discussion:
	This package is a collection of functions for computing with polynomials and formal power series.
*)


BeginPackage["Symbolics`"]

{
	CoefficientList,
	FragileTogether,
	FromCoefficientList,
	GeneratingFunctionRelations,
	ImplicitSeries,
	PolynomialSolve,
	SeriesCoefficientList
}


Unprotect["Symbolics`*"]

Begin["`Private`"]

var[variables___String] :=
	StringJoin[Riffle[("\*StyleBox[\"" <> # <> "\", \"TI\"]" &) /@ {variables}, ", "]]

sub[variable_String, indices : {___}] :=
StringJoin[Riffle[
	Replace[indices, {
		"..." -> "\*StyleBox[\"\[Ellipsis]\", \"TI\"]",
		s_String :> "\!\(\*SubscriptBox[StyleBox[\"" <> variable <> "\", \"TI\"], StyleBox[\"" <> s <> "\", \"TI\"]]\)",
		s_ :> "\!\(\*SubscriptBox[StyleBox[\"" <> variable <> "\", \"TI\"], StyleBox[\"" <> ToString[s] <> "\", \"TR\"]]\)"
	},
	{1}
], ", "]]
sub["...", _] := var["..."]
sub[variable_String, index_] := sub[variable, {index}]
sub[variables_List, index : Except[_List]] := sub @@ ({#, index} &) /@ variables
sub[variables : {_String, _} ...] := StringJoin[Riffle[sub @@@ {variables}, ", "]]


CoefficientList::usage =
"CoefficientList[" <> var["poly", "var"] <> "] gives a list of coefficients of powers of " <> var["var"] <> " in " <> var["poly"] <> ", starting with power 0.
CoefficientList[" <> var["poly"] <> ", {" <> sub["var", {1, 2, "..."}] <> "}] gives an array of coefficients of the " <> sub["var", "i"] <> ".
CoefficientList[" <> var["poly", "var", "polynomial"] <> "] gives the list of coefficients of " <> var["poly"] <> " in the basis specified by " <> var["polynomial"] <> ", which is a one\[Hyphen]parameter family of polynomials in " <> var["var"] <> " or one of \"BernoulliB\", \"Binomial\", \"ChebyshevT\", \"ChebyshevU\", \"EulerE\", \"Fibonacci\", \"HermiteH\", \"LaguerreL\", \"LegendreP\", \"NorlundB\", \"Pochhammer\", \"Power\"."

FragileTogether::usage =
"FragileTogether[expr] puts terms in a sum over a common denominator, without expanding the numerator."

FromCoefficientList::usage =
"FromCoefficientList[" <> var["list", "var"] <> "] constructs a polynomial in " <> var["var"] <> " from its coefficient list " <> var["list"] <> ".
FromCoefficientList[" <> var["array"] <> ", {" <> sub["var", {1, 2, "..."}] <> "}] constructs a multivariate polynomial from its coefficient array " <> var["array"] <> ".
FromCoefficientList[" <> var["list", "var", "polynomial"] <> "] constructs a polynomial in " <> var["var"] <> " from its coefficient list " <> var["list"] <> " in the basis specified by " <> var["polynomial"] <> ", which is a one\[Hyphen]parameter family of polynomials in " <> var["var"] <> " or one of \"BernoulliB\", \"Binomial\", \"ChebyshevT\", \"ChebyshevU\", \"EulerE\", \"Fibonacci\", \"HermiteH\", \"LaguerreL\", \"LegendreP\", \"NorlundB\", \"Pochhammer\", \"Power\"."

GeneratingFunctionRelations::usage =
"GeneratingFunctionRelations[{" <> sub["eqn", {1, 2, "..."}] <> "}, " <> var["a"] <> "[" <> sub["n", {1, 2, "...", "k"}] <> "], " <> var["f"] <> "[" <> sub["x", {1, 2, "...", "k"}] <> "]] gives equations satisfied by the multivariate generating function " <> var["f"] <> "[" <> sub["x", {1, 2, "...", "k"}] <> "] of the sequence " <> var["a"] <> "[" <> sub["n", {1, 2, "...", "k"}] <> "] satisfying a linear system of partial recurrence equations."

ImplicitSeries::usage =
"ImplicitSeries[" <> var["eqn", "f[x]"] <> ", {" <> var["x"] <> ", " <> sub["x", 0] <> ", " <> var["n"] <> "}] generates a power series expansion for the implicit function " <> var["f[x]"] <> " in " <> var["eqn"] <> " about the point " <> var["x"] <> " = " <> sub["x", 0] <> " to order \!\(\*SuperscriptBox[RowBox[{\"(\", RowBox[{StyleBox[\"x\", \"TI\"], \"-\", SubscriptBox[StyleBox[\"x\", \"TI\"], StyleBox[\"0\", \"TR\"]]}], \")\"}], StyleBox[\"n\", \"TI\"]]\)."

PolynomialSolve::usage =
"PolynomialSolve[" <> var["eqn", "vars"] <> "] solves the system of equations obtained by equating two polynomials in " <> var["vars"] <> ".
PolynomialSolve[" <> var["eqn", "vars", "solvevars"] <> "] solves for the variables " <> var["solvevars"] <> "."

SeriesCoefficientList::usage =
"SeriesCoefficientList[" <> var["f"] <> ", {" <> var["x"] <> ", " <> sub["x", 0] <> ", " <> var["n"] <> "}] gives a list of coefficients of the power series expansion for " <> var["f"] <> " about the point " <> var["x"] <> " = " <> sub["x", 0] <> " to order \!\(\*SuperscriptBox[RowBox[{\"(\", RowBox[{StyleBox[\"x\", \"TI\"], \"-\", SubscriptBox[StyleBox[\"x\", \"TI\"], StyleBox[\"0\", \"TR\"]]}], \")\"}], StyleBox[\"n\", \"TI\"]]\).
SeriesCoefficientList[" <> var["f"] <> ", {" <> var["x"] <> ", " <> sub["x", 0] <> ", " <> var["m"] <> "}, {" <> var["y"] <> ", " <> sub["y", 0] <> ", " <> var["n"] <> "}, \[Ellipsis]] gives an array of coefficients of the series expansion with respect to " <> var["x"] <> ", then " <> var["y"] <> ", etc."


Classify[l : head_[___], f_ : Identity, g_ : Identity, h_ : Identity] :=
Module[{list = {f[#], #} & /@ List @@ l, vals = {}, classes = {}},
	Function[{image, element},
		Occurrence[
			vals,
			image -> (AppendTo[classes[[#]], g[element]] &),
			AppendTo[vals, image]; AppendTo[classes, {g[element]}]
		]
	] @@@ list;
	Transpose[{h /@ head @@@ classes, vals}]
]

SetAttributes[Occurrence, HoldRest]
Occurrence[list : _[___], Rule[p_, f_], alt_ : Null] :=
If[Length[#] >= 1,
	f[#[[-1,1]]],
	alt
] & @ Position[list, p, {1}, 1, Heads -> False]


changebasis[poly_, var_, matrix_] :=
With[{e = Exponent[poly, var] + 1},
	LinearSolve[Transpose[Take[matrix, e, e]], CoefficientList[poly, var]]
]
Unprotect[CoefficientList]
CoefficientList[0, _, _] = {}
CoefficientList[
	poly_,
	var_,
	polynomial : "BernoulliB" | "ChebyshevT" | "ChebyshevU" | "EulerE" | "GegenbauerC" | "HermiteH" | "LaguerreL" | "LegendreP" | "NorlundB"
] :=
	CoefficientList[poly, var, ToExpression[polynomial][#, var] &]
CoefficientList[poly_, var_, "Binomial"] :=
	changebasis[poly, var,
		Table[
			StirlingS1[i, Range[0, Exponent[poly, var]]] / i!,
			{i, 0, Exponent[poly, var]}
		]
	]
CoefficientList[poly_, var_, "Fibonacci"] := CoefficientList[poly, var, Fibonacci[# + 1, var] &]
CoefficientList[poly_, var_, "Pochhammer"] := CoefficientList[poly, var, Pochhammer[var, #] &]
CoefficientList[poly_, var_, "Power"] := CoefficientList[poly, var]
CoefficientList[poly_, var_, basispoly_] /; PolynomialQ[poly, var] :=
	changebasis[poly, var,
		PadRight[#, Exponent[poly, var] + 1] & /@
			CoefficientList[basispoly /@ Range[0, Exponent[poly, var]], var]
	]
SyntaxInformation[CoefficientList] = {"ArgumentsPattern" -> {_, _, _., OptionsPattern[]}}
Protect[CoefficientList]

FragileTogether[expr_Plus] :=
With[{numerators = Numerator[List @@ expr], denominators = Denominator[List @@ expr]},
	Total[MapIndexed[#1 Times @@ Delete[denominators, First[#2]] &, numerators]] / Times @@ denominators
]
FragileTogether[expr_] := expr
SyntaxInformation[FragileTogether] = {"ArgumentsPattern" -> {_}}

FromCoefficientList[coeffs_List, vars_List] /; ArrayDepth[coeffs] >= Length[vars] :=
	Total[coeffs (Outer[Times, ##] &) @@ (vars^(Range[Dimensions[coeffs]] - 1)), Infinity]
FromCoefficientList[coeffs_List, 0] := First[coeffs]
FromCoefficientList[coeffs_List, var_] := coeffs.var^Range[0, Length[coeffs] - 1]
FromCoefficientList[
	coeffs_List,
	var_,
	polynomial : "BernoulliB" | "ChebyshevT" | "ChebyshevU" | "EulerE" | "GegenbauerC" | "HermiteH" | "LaguerreL" | "LegendreP" | "NorlundB"
] :=
	coeffs . ToExpression[polynomial][Range[0, Length[coeffs] - 1], var]
FromCoefficientList[coeffs_List, var_, polynomial : "Binomial" | "Pochhammer"] :=
	coeffs . ToExpression[polynomial][var, Range[0, Length[coeffs] - 1]]
FromCoefficientList[coeffs_List, var_, "Fibonacci"] :=
	coeffs . Fibonacci[Range[Length[coeffs]], var]
FromCoefficientList[coeffs_List, var_, "Power"] := FromCoefficientList[coeffs, var]
SyntaxInformation[FromCoefficientList] = {"ArgumentsPattern" -> {_, _, _.}}

GeneratingFunctionRelations[eqn_Equal, expr_, genfun_] :=
	GeneratingFunctionRelations[{eqn}, expr, genfun]
GeneratingFunctionRelations[eqns_And, expr_, genfun_] :=
	GeneratingFunctionRelations[List @@ eqns, expr, genfun]
GeneratingFunctionRelations[
	eqns : {HoldPattern[Equal][_, _] ..},
	(a : Except[List])[v__],
	(f : Except[List])[x__]
] /; Length[{v}] == Length[{x}] :=
Module[
	{
		vars = {v}, formalvars = {x}, MySum,
		monomiallists = MonomialList /@ Subtract @@@ eqns,
		linearsystemq, lowindices, highindices, expr, shiftrules, shiftedlowindices
	},
	linearsystemq = And @@ Flatten[Map[
		MatchQ[#, c_ | c_. a @@ (_Integer | # + Optional[_Integer] &) /@ vars /; FreeQ[c, a | Alternatives @@ vars]] &,
		monomiallists,
		{2}
	]];
	Function[monomiallist,
		lowindices = (Max[Cases[monomiallist, a[___, # + shift_., ___] :> -shift, {0, Infinity}]] &) /@ vars;
		highindices = Replace[lowindices, {-Infinity -> 0, _ -> Infinity}, {1}];
		lowindices =  lowindices /. -Infinity -> 0;
		Collect[Numerator[FragileTogether[Total[Apply[
			Function[{coeff, label},
				expr = Replace[label, {{l_} :> l, _ -> 1}];
				shiftrules = Select[
					Thread[vars -> vars - Replace[List @@ expr - vars, Except[_Integer] -> Null, {1}]],
					FreeQ[#, Null] &
				];
				shiftedlowindices = lowindices + vars - (vars /. shiftrules);
				coeff Times @@ (formalvars^-shiftedlowindices) (
					MySum[
						Times @@ (formalvars^vars) (expr /. shiftrules),
						Sequence @@ Transpose[{
							vars,
							shiftedlowindices,
							highindices
						}]
					] //. {
						MySum[arg_, i___, {s_, l_?Positive, Infinity}, j___] :>
							MySum[arg, i, {s, 0, Infinity}, j] - MySum[arg, i, {s, 0, l - 1}, j],
						MySum[arg_, i___, {s_, l_, u_Integer}, j___] :>
							Sum[MySum[arg, i, j], {s, l, u}]
					} //. {
						MySum[arg_] :> arg,
						MySum[c_ arg_., i__] /; FreeQ[c, Alternatives @@ First /@ {i}] :> c MySum[arg, i]
					} /. {
						MySum[arg_, i__] :>
							With[{sumvars = Replace[vars, Except[Alternatives @@ First /@ {i}] -> 0, {1}]},
								If[MatchQ[arg, Times @@ (formalvars^sumvars) a @@ sumvars],
									f @@ Replace[Transpose[{formalvars, sumvars}], {{_, 0} -> 0, {y_, _} :> y}, {1}],
									MySum[arg, i]
								]
							]
					} /.
						MySum -> Sum
				)
			],
			Classify[monomiallist, Cases[#, a[__], {0, Infinity}] &, # /. a[__] -> 1 &, Total],
			{1}
		]]]], _f] == 0
	] /@ monomiallists /; linearsystemq
]
SyntaxInformation[GeneratingFunctionRelations] = {"ArgumentsPattern" -> {_, _, _}}

ImplicitSeries[lhs_ == rhs_, f_[x_], {x_, x0_, n_Integer?NonNegative}] /; !FreeQ[lhs == rhs, f[_?(!FreeQ[#, x] &)]] :=
Module[
	{
		coefficients = {{}},
		expr = (Denominator[#2] Numerator[#1] - Denominator[#1] Numerator[#2] &)
			[Together[lhs], Together[rhs]]
	},
	While[Min[Length /@ coefficients] <= n,
		coefficients =
			If[Length[#] > n,
				#,
				With[{sols = Solve[(expr /. x -> x0 /. #) == 0, Derivative[Length[#]][f][x0]]},
					Switch[sols,
						{{}}, #,
						_, Sequence @@ Thread[Append[#, First /@ sols]]
					]
				]
			] & /@ coefficients;
		expr = D[expr, x]
	];
	(f[x] -> FromCoefficientList[Last /@ # / Range[0, n]!, x] + O[x - x0]^(n + 1) &) /@
		coefficients
]
SyntaxInformation[ImplicitSeries] = {"ArgumentsPattern" -> {_, _, _}}

PolynomialSolve[lhs_ == rhs_, vars_, solvevars_ : {}] :=
With[{coeffs = CoefficientList[lhs - rhs, vars]},
	Solve[coeffs == ConstantArray[0, Dimensions[coeffs]], solvevars]
]
SyntaxInformation[PolynomialSolve] = {"ArgumentsPattern" -> {_, _, _.}}

SeriesCoefficientList[expr_, varspecs : ({_, _, _Integer?NonNegative} ..)] :=
With[{poly = Normal[Series[expr, varspecs]]},
	PadRight[CoefficientList[poly, First /@ {varspecs}], Last /@ {varspecs} + 1] /;
		PolynomialQ[poly, First /@ {varspecs}]
]
SyntaxInformation[SeriesCoefficientList] = {"ArgumentsPattern" -> {_, __}}


End[]

Protect["Symbolics`*"]

EndPackage[]