(*** RSolve.m ***)

(* alpha test version *)

(* last modified: June 7, 1990 *)

(* Copyright C 1990 by Marko Petkovsek *)

BeginPackage["RSolve`"]

RSolve::usage = "RSolve[{eqn1, eqn2, ..}, {a1[n], a2[n], ..}, n,
opts] solves the recurrence equations eqn1, eqn2, ..
for the sequences a1[n], a2[n], ... If there is a single sequence or
equation it need not be given in a list. Equations can either be
recurrences or initial/boundary conditions. An equation may have
the form eqn /; cond where cond is an inequality specifying the range
of values of n for which eqn is valid; when no condition is given
it is assumed to be n >= 0. Alias: RS."

PowerSum::usage = "PowerSum[expr, {z, n, n0:0}] computes
Sum[expr z^n, {n, n0, Infinity}]. If a[n] is a sequence, the name
Gf[a][z] is used by PowerSum to denote Sum[a[n] z^n, {n, 0, Infinity}].
PowerSum is threaded over lists and equations. Alias: PS."

ExponentialPowerSum::usage = "ExponentialPowerSum[expr, {z, n, n0:0}]
computes Sum[expr z^(n-n0) / (n-n0)!, {n, n0, Infinity}]. If a[n] is a
sequence, the name EGf[a][z] is used by ExponentialPowerSum to denote
Sum[a[n] z^n / n!, {n, 0, Infinity}]. ExponentialPowerSum is threaded
over lists and equations. Alias: EPS."

GeneratingFunction::usage = "GeneratingFunction[{eqn1, eqn2, ..},
{a1[n], a2[n], ..}, n, z, opts] computes the ordinary generating
functions Sum[ai[n] z^n, {n, si, Infinity}] for the sequences a1[n],
a2[n], .. defined by the equations eqn1, eqn2, ... Here si denotes the
least value of n such that ai[n] appears in the equations. An equation
may have the form eqn /; cond where cond is an inequality specifying
the range of values of n for which eqn is valid; when no condition
is given it is assumed to be n >= 0. Alias: GF."

ExponentialGeneratingFunction::usage = "ExponentialGeneratingFunction[
{eqn1, eqn2, ..}, {a1[n], a2[n], ..}, n, z] computes the exponential
generating functions Sum[ai[n + si] z^n / n!, {n, 0, Infinity}] for the
sequences a1[n], a2[n], .. defined by the equations eqn1, eqn2, ... Here
si denotes the least value of n such that ai[n] appears in the equations.
An equation may have the form eqn /; cond where cond is an inequality
specifying the range of values of n for which eqn is valid; when no
condition is given it is assumed to be n >= 0. Alias: EGF."

Gf::usage = "Gf[a][z] is used by PowerSum to denote
Sum[a[n] z^n, {n, 0, Infinity}]."

EGf::usage = "EGf[a][z] is used by ExponentialPowerSum to denote
Sum[a[n] z^n / n!, {n, 0, Infinity}]."

HSolve::usage = "HSolve[eq, f, z] returns a hypergeometric formal
series solution to differential equation eq with unknown function
f[z], or {} when it determines that there is no series solution.
It returns Fail when it is unable to do either."

HypergeometricF::usage = "HypergeometricF[{a1, a2, ..}, {b1, b2, ..}, z]
is the generalized hypergeometric function with upper parameters
a1, a2,.., lower parameters b1, b2,.., and argument z."

K::usage = "K is the head used for summation indices by SeriesTerm."

SeriesTerm::usage = "SeriesTerm[expr, {z, a, n}, opts]
returns the n-th coefficient of the Laurent series expansion of expr
around z = a, for symbolic n. Alias: ST."

SymbolicMod::usage = "SymbolicMod[n, k] is equal to Mod[n, k] when n
is a number, and is left unevaluated otherwise."

RealQ::usage = "'RealQ[x] = True' serves to declare x to be real."

ReP::usage = "ReP[x] gives the real part of x."

ImP::usage = "ImP[x] gives the imaginary part of x."

ConjugateQ::usage = "ConjugateQ[x, y] tests whether x and y are complex
conjugates."

SimplifySum::usage = "SimplifySum is a list of rules for simplification
of sums."

SimplifyTrig::usage = "SimplifyTrig is a list of rules which put
trigonometric expressions into canonical form."

SimplifyComplexE1::usage = "SimplifyComplexE1 is one of lists of rules
used to replace sums of conjugate complex quantities by their double
real parts."

SimplifyComplexE2::usage = "SimplifyComplexE2 is one of lists of rules
used to replace sums of conjugate complex quantities by their double
real parts."

SimplifyComplex2::usage = "SimplifyComplex2 is one of lists of rules
used to replace sums of conjugate complex quantities by their double
real parts."

SimplifyComplex3::usage = "SimplifyComplex3 is one of lists of rules
used to replace sums of conjugate complex quantities by their double
real parts."

SimplifyComplex4::usage = "SimplifyComplex4 is one of lists of rules
used to replace sums of conjugate complex quantities by their double
real parts."

SimplifyComplexN::usage = "SimplifyComplexN is one of lists of rules
used to replace sums of conjugate complex quantities by their double
real parts."

FactorialSimplify::usage = "FactorialSimplify[expr] simplifies
occurrences of Factorial[] inside expr."

PartialFractions::usage = "PartialFractions[f, z, opts] gives a partial
fraction decomposition of f[z] over the complex functions."

ArgPi::usage = "When possible ArgPi[z] expresses the complex number z
in the form Abs[z] E^(r Pi I) where r is a rational number."

Entails::usage = "Entails[a, b] determines whether a implies b where
a and b are Boolean combinations of equalities and/or inequalities among
rational functions of a single integer variable."

FirstPos::usage = "FirstPos[l, pattern] returns the position of
the first argument of l which matches pattern, or Null if there is
no such element."

FreeL::usage = "FreeL[expr, list] returns True if none of the
elements of list appears in expr, False otherwise."

Info::usage = "Info[n] is information about n, expressed as a
Boolean combination of inequalities. It should be associated with n."

ISolve::usage = "ISolve[expr] solves a Boolean combination of equalities
and/or inequalities in a single integer variable. The answer is given
as a disjunction of disjoint inclusive inequalities."

LeadingCoef::usage = "LeadingCoef[poly, x] returns the leading
coefficient of the polynomial poly[x]."

ListComplement::usage = "ListComplement[l1_List, l2_List] returns
the multiset-difference of l1 and l2."

MakeList::usage = "MakeList[expr, head:List] returns the list of
arguments of expr if the head of expr matches the given one, or else
it returns {expr}."

MakeTrinomial::usage = "MakeTrinomial[expr_] tries to write expr
in the form Multinomial[a, b, c]."

PatternList::usage = "PatternList[expr, pattern] returns the list of
all subexpressions of expr which match the given pattern, provided
that expr does not contain Rule[]."

PoleMultiplicity::usage = "PoleMultiplicity[f, {z, a}] computes
the multiplicity of the pole of f at z = a."

Reset::usage = "Reset[recur, conds, unknowns, startValues, s0] returns
the list {recur, conds} with starting values of the index reset so
that it starts with s0 rather than with startValues, for all unknown
sequences."

SafeFirst::usage = "SafeFirst[l] returns the first argument of l,
or Null if the length of l is zero."

SafeSeries::usage = "SafeSeries[f, {z, a, n}] overcomes a bug in
Series[] which fails when n is less than the degree of zero of the
denominator of f at a."

TTQ::usage = "TTQ[cond] returns True if it determines that the
information given about the variables implies that cond is true,
and False otherwise."

UserSymbols::usage = "UserSymbols[expr] returns a list of maximal
subexpressions of expr with the property that their head contains
symbols defined in a context different from System`. "

When::usage = "When[cond, a, b] returns a if it determines that the
information given about the variables implies that cond is true,
b if it determines that the information given about the variables
implies that cond is false, and If[cond, a, b] otherwise."

Methods::usage = "Methods is an optional parameter for RSolve.
It specifies which methods and in what order are to be used by RSolve.
Default: Methods -> {MethodGF, MethodEGF}. If a single method is given
it need not be enclosed in a list."

MethodGF::usage = "MethodGF is a possible value for Methods. It denotes
the method of ordinary generating functions."

MethodEGF::usage = "MethodEGF is a possible value for Methods. It denotes
the method of exponential generating functions."

PrecisionHS::usage = "PrecisionHS is an optional parameter for RSolve,
GeneratingFunction, ExponentialGeneratingFunction, and HSolve. It
specifies the precision with which the roots of polynomials are to be
computed by HSolve. When successful, the result of HSolve is a
generalized hypergeometric function, and these roots are its parameters.
Possible values: PrecisionHS -> n, PrecisionHS -> Infinity,
PrecisionHS -> Automatic (default). Infinity means exact roots, but
those roots which cannot be found by Solve[] are computed numerically.
Automatic means that only those roots which Solve[] can find without
recourse to general formulas for solving cubics and quartics are found
exactly, while others are computed numerically."

PrecisionST::usage = "PrecisionST is an optional parameter for RSolve
and SeriesTerm. It specifies the precision with which the roots of
denominators of rational functions to be expanded into power series
are to be computed by SeriesTerm. Possible values: PrecisionST -> n,
PrecisionST -> Infinity, PrecisionST -> Automatic (default). Infinity
means exact roots, but those roots which cannot be found by Solve[]
are computed numerically. Automatic means that only those roots which
Solve[] can find without recourse to general formulas for solving cubics
and quartics are found exactly, while others are computed numerically."

UseApart::usage = "UseApart is an optional parameter for RSolve
and SeriesTerm. It specifies whether or not Apart should be used
before attempting expansion into power series, and in what form.
Possible values: UseApart -> None (don't use Apart at all),
UseApart -> Automatic (default; use Apart without second argument),
UseApart -> All (use Apart with second argument). By default, MethodGF
uses UseApart -> Automatic and MethodEGF uses UseApart -> None.
UseApart -> All is useful when the function to be expanded contains
parameters."

UseMod::usage = "UseMod is an optional parameter for RSolve and
SeriesTerm. Possible values: UseMod -> True (default; use EvenQ, OddQ
and/or SymbolicMod to express the result), UseMod -> False (the
opposite)."

MakeReal::usage = "MakeReal is an optional parameter for RSolve and
SeriesTerm. Possible values: MakeReal -> True (default; replace pairs
of conjugate complex quantities by their doble real parts),
MakeReal -> False (the opposite)."

SpecialFunctions::usage = "SpecialFunctions is an optional parameter
for RSolve and SeriesTerm. Possible values: SpecialFunctions -> True
(default; use special functions such as Legendre polynomials to express
the result), SpecialFunctions -> False (the opposite)."

Begin["`Private`"]

(** ALIASES **)

RS ::= RSolve

PS ::= PowerSum

EPS ::= ExponentialPowerSum

GF ::= GeneratingFunction

EGF ::= ExponentialGeneratingFunction

ST ::= SeriesTerm

(** ATTRIBUTES **)

Attributes[RSolve] = {HoldFirst}

Attributes[GeneratingFunction] = {HoldFirst}

Attributes[ExponentialGeneratingFunction] = {HoldFirst}

Attributes[Parse] = {HoldFirst}

(** OPTIONS **)

Options[RSolve] := Join[{Methods ->{MethodGF, MethodEGF}},
Options[HSolve],
Options[SeriesTerm] ]

Options[GeneratingFunction] := Options[HSolve]

Options[ExponentialGeneratingFunction] := Options[HSolve]

Options[HSolve] = {PrecisionHS -> Automatic}

Options[SeriesTerm] = {PrecisionST -> Automatic,
UseApart -> Automatic,
MakeReal -> True,
UseMod -> True,
SpecialFunctions -> True}

Options[PartialFractions] := Options[SeriesTerm]

(** MESSAGES **)

ST::badopt = "`` -> `` is not a valid option."

PF::badopt = "`` -> `` is not a valid option."

Parse::eqn = "`` is not an equation or a system of equations."

(** REDEFINITIONS OF SYSTEM FUNCTIONS **)

Unprotect[Binomial, Power]

Binomial[n_, n_] := When[n >= 0, 1, 0]

0^0 = 1

Protect[Binomial, Power]

(** (RSolve.m) **)

(* RSOLVE *)

RSolve[list1_, list2_, n_, opts___Rule] :=
Block[{recur, conds, unknowns, startValues, methods, result, prs},
prs = Parse[list1, list2, n];
If[prs === Fail,
Return[Fail],
{recur, conds, unknowns, startValues} = prs ];
methods = MakeList[Methods /.{opts} /.Options[RSolve]];
result = Scan[Block[{temp},
temp = # @@ {recur, conds, unknowns, startValues, n, opts};
If[temp =!= Fail, Return[temp]]]&, methods];
Which[result === Null,
Return[Fail],
True,
Return[result] ]]

(* METHOD GF *)

MethodGF[recur_, conds_, unknowns_, startValues_, n_, opts___Rule] :=
Block[{genf, mm, z, rec, con, i, temp, start, optsList = {opts}},
{rec, con} = Reset[recur, conds, unknowns, startValues, 0];
genf = FunctionSolve[rec, con, unknowns, n, PowerSum,
z, opts];
If[genf =!= Fail,
(Release[start /@ unknowns]) = startValues;
If[FreeQ[optsList, UseApart],
AppendTo[optsList, UseApart -> Automatic] ];
temp = (Table[mm /: Info[mm] = mm >= startValues[[i]];
SeriesTerm[#[[i]],
{z, 0, mm - startValues[[i]]}, Sequence @@ optsList],
{i, Length[unknowns]} ]& /@ genf) /. mm -> n;
Return[temp /. aa_?(MemberQ[unknowns, #]&)[kk_] :>
aa[Expand[kk + start[aa]]] ],
Return[Fail] ]]

(* METHOD EGF *)

MethodEGF[recur_, conds_, unknowns_, startValues_, n_, opts___Rule] :=
Block[{genf, mm, z, rec, con, i, temp, start, optsList = {opts}},
{rec, con} = Reset[recur, conds, unknowns, startValues, 0];
genf = FunctionSolve[rec, con, unknowns, n, ExponentialPowerSum,
z, opts];
If[genf =!= Fail,
(Release[start /@ unknowns]) = startValues;
If[FreeQ[optsList, UseApart],
AppendTo[optsList, UseApart -> None] ];
temp = (Table[mm /: Info[mm] = mm >= startValues[[i]];
FactorialSimplify[
(mm - startValues[[i]])! SeriesTerm[#[[i]],
{z, 0, mm - startValues[[i]]}, Sequence @@ optsList] ],
{i, Length[unknowns]} ]& /@ genf) /. mm -> n;
Return[temp /. aa_?(MemberQ[unknowns, #]&)[kk_] :>
aa[Expand[kk + start[aa]]] ],
Return[Fail] ]]

(** (GeneratingFunctions.m) **)

(* ROOTS OF UNITY *)

Omega[n_] := Exp[2 Pi I / n]

(* SIMPLIFICATION OF EXP[LOG[_]] AND I*SQRT[_] *)

SimplifyExpLog = {E^((a_. Log[b_] + c_.)d_.) -> b^(a d) E^(c d)}

SimplifyI = {Complex[0, a_] Power[b_, Rational[n_, 2]] :>
a (-1)^((n-1)/2) Expand[-b]^(n/2)}

(* SOLVING ROUTINE *)

FunctionSolve[recur_, conds_, unknowns_, n_, Transform_, z_,
opts___Rule] :=
Block[{eqns, extra, range, start, aa, kk, nm,
lo, hi, functions, fnHeads, initials, constants, temp,
series, order, xx, ii, bb, initvals, soln, times,
arg, G, init, subst, pinitials},

(* transform recurrence equations into functional equations *)

eqns = Expand[#[[1]] - #[[2]]]& /@ recur;
eqns = Function[xx,
Block[{hom, inhom, yy},
yy = MakeList[xx, Plus];
inhom = Select[yy, FreeL[#, unknowns]&];
hom = Together[Plus @@ Complement[yy, inhom]];
Numerator[hom] + Expand[(Plus @@
inhom) Denominator[hom]] ]] /@ eqns;
eqns = Transform[eqns, z, n];
If[!FreeL[PatternList[eqns, Transform[__]], unknowns],
Return[Fail]];
Which[Transform === PowerSum, G = Gf,
Transform === ExponentialPowerSum, G = EGf];

(* substitute values which are explicitly given by initial
conditions *)

init = Select[conds, Function[xx, MatchQ[xx,
a_?(MemberQ[unknowns, #]&)[_?NumberQ] == _?NumberQ]]];
subst = Solve[init, Union @@ (PatternList[init, #[_]]& /@
unknowns)];
If[subst === {}, Return[Fail], eqns = eqns /. First[subst]];

(* remember the unknown functions *)

functions = G[#][z]& /@ unknowns;
fnHeads = G /@ unknowns;

(* solve the equations using Solve, DSolve, or HypergeometricSolve *)

If[FreeQ[eqns, Derivative],

temp = Solve[Thread[eqns == 0], functions],

(* before using differential equation solvers, integrate
equations of the form lhs' == rhs' *)

times = Integrable[#, fnHeads, z]& /@ eqns;
eqns = PreIntegrate[#[[1]], fnHeads, {z, #[[2]]}]& /@
Thread[List[eqns, times]];
eqns = Table[Sum[C[ii, jj] z^(jj - 1), {jj, times[[ii]]}] +
eqns[[ii]], {ii, Length[eqns]} ];
eqns = eqns /. Derivative[kk_][ff_] :> D[ff, {z, kk}];
eqns = Thread[eqns == 0];

(* use HSolve then DSolve *)

temp = HSolve1[eqns, functions, z, opts];
temp = temp /. C -> First[unknowns];
If[!FreeL[temp, {HSolve, Fail}],
temp = DSolve[eqns, functions, z] ]];
If[temp === {}, Return[{}]];
If[!FreeL[temp, {Solve, DSolve, Fail}], Return[Fail]];
If[!FreeL[functions /. temp, fnHeads], Return[Fail]];
temp = temp //. SimplifyExpLog;

(* find all initial terms occurring in conditions, equations
or solutions *)

initials = Union @@ (PatternList[{conds, eqns,
functions /. temp}, #[_]]& /@ unknowns);

(* sort initials by falling indices *)

initials = Join @@ (Append[Cases[initials, _[#]]& /@
Reverse[Union[PatternList[initials, _?NumberQ]]], {}]);
pinitials = Select[initials, (#[[1]] >= 0)&];

(* How many series terms will be needed to determine
the initials? *)

(order[#] = Max[Append[First /@ PatternList[initials,
#[_]], 0 ]])& /@ unknowns;

(* determine values of initials and eliminate constants
of integration *)

constants = PatternList[functions /. temp, C[__]] // Union;
series = temp /. (G[aa_?(MemberQ[unknowns, #]&)][z] -> bb_
) :> (aa -> Taylor[bb, z, 0, order[aa]]);
If[!FreeL[series, {Series, Fail}], Return[Fail]];
If[G === EGf,
series = series /. (aa_?(MemberQ[unknowns, #]&) -> bb_
) :> (aa -> Table[nm!, {nm, 0, order[aa]}] bb) ];
initvals = Function[xx, Solve[Union[conds,
(# == (Head[#] /. xx)[[First[#] + 1]])& /@ pinitials,
Thread[Join @@ Rest /@ CList[#[[2,1]]& /@ xx /.
z -> z^(-1), z ] == 0 ] ],
Join[constants, initials] ]] /@ series;

(* plug in initial values *)

soln = Select[Union @@ Table[(Function[xx,
G[xx][z] /. temp[[ii]] ] /@
unknowns) /. initvals[[ii]],
{ii, Length[temp]} ], (Head[#] === List)& ];
Return[soln] ]

(* TAYLOR and CLIST *)

Taylor[f_, z_, a_, n_] :=
Block[{tmp, kk},
tmp = CoefficientList[SafeSeries[f, {z, a, n}], z];
If[Length[tmp] > n + 1, tmp = Take[tmp, n + 1]];
Return[Join[tmp, Table[0,
{kk, n + 1 - Length[tmp]} ]]]]

CList[l_List, x_] := CList[#, x]& /@ l

CList[0, x_] := {0}

CList[a_, x_] := CoefficientList[a, x] /; a =!= 0

(* INTEGRABLE and PREINTEGRATE *)

Integrable[a_Plus, b___] := Min[Integrable[#, b]& /@ (List @@ a)]

Integrable[a_, functions_, z_] := Infinity /; FreeL[a, functions]

Integrable[a_ b_, functions_, z_] :=
Integrable[b, functions, z] /; FreeQ[a, z]

Integrable[a_, functions_, z_] :=
FirstArg[a] /; FirstHead[a] === Derivative

Integrable[a___] := 0

PreIntegrate[a_, functions_, {z_, 0}] := a

PreIntegrate[a_Plus, b___] := PreIntegrate[#, b]& /@ a

PreIntegrate[a_, functions_, {z_, n_}] :=
Nest[Integrate[#, z]&, a, n] /; FreeL[a, functions]

PreIntegrate[a_ b_, functions_, {z_, n_}] :=
a PreIntegrate[b, functions, {z, n}] /; FreeQ[a, z]

PreIntegrate[a_, functions_, {z_, n_}] :=
MapAt[(# - n)&, a,
Append[First[Position[a, Derivative[_]]], 1]] /;
FirstHead[a] === Derivative

(* FIRSTHEAD, FIRSTARG and HEADCOUNT *)

FirstHead[a_] := Nest[Head, a, HeadCount[a]]

FirstArg[a_] := Which[Length[a] == 0,
Null,
True,
First[Nest[Head, a, HeadCount[a] - 1]] ]

HeadCount[a_] := Which[Length[a] == 0,
0,
True,
HeadCount[Head[a]] + 1 ]

(* POWER SUM *)

(* lists and equations *)

PowerSum[expr_List, rest__] := PowerSum[#, rest]& /@ expr

PowerSum[expr_Equal, rest__] := PowerSum[#, rest]& /@ expr

(* revert to inner syntax *)

PowerSum[expr_, {z_Symbol, n_, n0_:0}] :=
Block[{e = expr //. {Sum[a_, {k_, m_}] :> Sum[a, {k, 1, m}],
EvenQ[a_] :> SymbolicMod[a, 2] == 0,
OddQ[a_] :> SymbolicMod[a, 2] == 1,
Mod -> SymbolicMod} },
PowerSum[e, z, n, n0] /.
Derivative[kk_][ff_] :> D[ff, {z, kk}] ];

(* nonzero starting point *)

PowerSum[expr_, z_, n_, n0_] :=
z^n0 PowerSum[Expand[expr /. n -> n + n0] /.
a_Symbol?(Context[#] =!= "System`" && # =!= K &)[m_] :>
a[Expand[m]], z, n] /; NumberQ[N[n0]]

(* the geometric series *)

PowerSum[1, z_, n_] := 1/(1 - z)

(* linearity *)

PowerSum[c_ expr_., z_, n_] :=
c PowerSum[expr, z, n] /; FreeL[c, {n, Blank}]

PowerSum[c_.(expr1_ + expr2_), rest__] :=
PowerSum[c expr1, rest] + PowerSum[c expr2, rest]

PowerSum[Sum[a_, {k_, lo_, hi_}] b_., z_, n_] :=
Sum[PowerSum[Expand[a b], z, n], {k, lo, hi}] /;
FreeQ[{lo, hi}, n] && FreeQ[b, k]

(* powers *)

PowerSum[c_^((a_ + b_)d_) e_., rest__] :=
PowerSum[c^(a d + b d) e, rest]

PowerSum[c_^(a_ + b_) e_., z_, n_] :=
c^a PowerSum[c^b e, z, n] /;
FreeL[c, {n, Blank}] && FreeL[a, {n, Blank}]

PowerSum[(a_ + b_)^c_Integer?Positive d_., rest__] :=
PowerSum[Expand[(a + b)^c d], rest]

(* exponentials *)

PowerSum[e_. c_^((a_. n_ + b_.)d_.), z_, n_] :=
(c^(b d) PowerSum[e, z, n] /. z -> c^(a d) z) /;
(FreeL[#, {n, Blank}]& /@ And[a, b, c, d]) && FreeL[e, {z, Blank}]

(* binomials *)

PowerSum[Binomial[a_, n_ + k_.], z_, n_] :=
z^(-k)((1 + z)^a - Sum[Binomial[a, m] z^m, {m, 0, k - 1}]) /;
FreeQ[a, n] && IntegerQ[k]

PowerSum[Binomial[a_, b_], z_, n_] :=
Block[{e = Expand[n - b], r = Expand[a - 2 b], m},
z^e (((1 - Sqrt[1 - 4 z])/(2 z))^r / Sqrt[1 - 4 z] -
Sum[Binomial[2 m + r, m] z^m, {m, 0, -(e + 1)}] ) /;
IntegerQ[e] && IntegerQ[r]]

(* trinomials *)

PowerSum[Sum[a_. b_Binomial c_Binomial, {k_, 0, n_}],
z_, n_] :=
Block[{aa = a^(1/k)},
1/Sqrt[1 - 2 (2 aa + 1) z + z^2] /;
FreeL[aa, {k, n}] &&
MakeTrinomial[b c] == Sort[Multinomial[k, k, n - k]] ]

PowerSum[Sum[a_. b_Binomial c_Binomial, {k_, 0, n_}],
z_, n_] :=
Block[{aa = a^(1/(n - k))},
1/Sqrt[1 - 2 (2 aa + 1) z + z^2] /;
FreeL[aa, {k, n}] &&
MakeTrinomial[b c] == Sort[Multinomial[k, n - k, n - k]] ]

PowerSum[Sum[a_. b_Binomial c_Binomial, {k_, 0, n_}],
z_, n_] :=
Block[{aa = a^(1/k)},
1/Sqrt[1 - 2 z + (1 - 4 aa) z^2] /;
FreeL[aa, {k, n}] &&
MakeTrinomial[b c] == Sort[Multinomial[k, k, n - 2 k]] ]

PowerSum[Sum[a_. b_Binomial c_Binomial, {k_, 0, n_}],
z_, n_] :=
Block[{aa = a^(1/(n - k))},
1/Sqrt[1 - 2 z + (1 - 4 aa) z^2] /;
FreeL[aa, {k, n}] &&
MakeTrinomial[b c] == Sort[Multinomial[n - k, n - k, 2 k - n]] ]

(* PowerSum is the inverse of SeriesTerm *)

PowerSum[Literal[SeriesTerm[f_, z_, n_ + m_.]],
z_, n_] :=
Block[{k},
z^(-m) (f - Sum[SeriesTerm[f, {z, 0, k}] z^k,
{k, 0, m - 1}]) /;
PoleMultiplicity[f, {z, 0}] == 0 && IntegerQ[m] ]

(* trigonometric inhomogeneity *)

PowerSum[Cos[n_], rest__] :=
PowerSum[(E^(I n) + E^(-I n))/2, rest]

PowerSum[Sin[n_], rest__] :=
PowerSum[(E^(I n) - E^(-I n))/2/I, rest]

(* factorial inhomogeneity *)

PowerSum[expr_./(n_ + a_.)!, z_, n_] :=
Block[{k},
(ExponentialPowerSum[Expand[expr /. n -> n - a], {z, n}] -
Sum[(expr /. n -> k - a) z^k / k!, {k, 0, a - 1}]) / z^a] /;
IntegerQ[a]

(* generating functions of sequences

(here we assume that a[n] is defined for n >= 0) *)

PowerSum[n_^m_. a_[n_ + d_.], z_, n_] :=
Block[{j, k},
Sum[StirlingS2[m, j] z^j
Derivative[j][(PowerSum[a[k], z, k] -
Sum[a[k] z^k, {k, 0, d - 1}])/z^d ],
{j, 0, m}] // Factor] /;
IntegerQ[d] && TTQ[d >= 0] && IntegerQ[m] && TTQ[m >= 0]

PowerSum[a_[k_. n_ + d_.], z_, n_] :=
Block[{j, m},
(z^(-d/k)/k Sum[Omega[k]^(-d j) (Gf[a][Omega[k]^j z^(1/k)] -
Sum[Omega[k]^(m j) a[m] z^(m/k), {m, 0, d-1}]),
{j, 0, k - 1} ] ) // Factor] /;
IntegerQ[d] && TTQ[d >= 0] && IntegerQ[k] && TTQ[k > 0]

(* rational inhomogeneity *)

PowerSum[n_^m_., z_, n_] :=
Block[{j},
Factor[Sum[(StirlingS2[m, j] j! z^j)/(1 - z)^(j + 1),
{j, 0, m}] ] ] /; IntegerQ[m] && TTQ[m >= 0]

PowerSum[(n_ + a_)^m_Integer?Negative, z_, n_] :=
Block[{j},
z^(-a) PolyLog[-m, z] - Sum[z^(j-a) j^m, {j, 1, a-1}] ] /;
IntegerQ[a] && TTQ[a > 0]

(* convolutions *)

PowerSum[Sum[expr_, {k_, lo_, hi_}], z_, n_] :=
Block[{nfree, nfull, s, aa = Expand[lo], bb = Expand[hi - n]},
{nfree, nfull} = If[Head[expr] === Times,
{Select[expr, FreeQ[#, n]&], Select[expr, !FreeQ[#, n]&]},
If[FreeQ[expr, n], {expr, 1}, {1, expr}] ];
Sum[z^k nfree PowerSum[Expand[nfull /. n -> s + k],
z, s, -k], {k, aa, bb - 1}] +
PowerSum[nfree PowerSum[Expand[nfull /. n -> s + k],
z, s, -bb], z, k, Max[aa, bb] ] /;
FreeQ[aa, n] && FreeQ[bb, n] ]

(* multisection of series *)

PowerSum[If[SymbolicMod[n_ + d_., m_] == k_, a_, b_] c_., z_, n_] :=
Block[{l = Mod[k - d, m], jj,
fps = PowerSum[Expand[(a - b)c], z, n]},
PowerSum[Expand[b c], z, n] +
1/m Sum[Omega[m]^(jj(m - l)) fps /. z -> z Omega[m]^jj,
{jj, 0, m - 1}] // Expand ]

(* other conditionals *)

PowerSum[If[cond_, a_, b_] c_., z_, n_] :=
Block[{t, tt, s, k},
t = IneqSolve[cond && n >= 0, n];
tt = IneqSolve[!cond && n >= 0, n];
Plus @@ (Sum[Expand[a c] z^n, {n, First[#], Last[#]}]& /@ t) +
Plus @@ (Sum[Expand[b c] z^n, {n, First[#], Last[#]}]& /@
tt) /.
{Sum[s_. z^n, {n, k_, Infinity}] :>
PowerSum[s, {z, n, k}],
Sum[0, {__}] -> 0} ]

(* PowerSum under Series *)

PowerSum /:
Series[PowerSum[a_, z_, n_], {z_, 0, m_}] :=
Block[{k},
Sum[(a /. n -> k) z^k, {k, 0, m}] + O[z]^(m+1)]

(* EXPONENTIAL POWER SUM *)

(* lists and equations *)

ExponentialPowerSum[expr_List, rest__] :=
ExponentialPowerSum[#, rest]& /@ expr

ExponentialPowerSum[expr_Equal, rest__] :=
ExponentialPowerSum[#, rest]& /@ expr

(* revert to inner syntax *)

ExponentialPowerSum[expr_, {z_Symbol, n_, n0_:0}] :=
Block[{e = expr //. {Sum[a_, {k_, m_}] :> Sum[a, {k, 1, m}],
EvenQ[a_] :> SymbolicMod[a, 2] == 0,
OddQ[a_] :> SymbolicMod[a, 2] == 1,
Mod -> SymbolicMod} },
ExponentialPowerSum[e, z, n, n0] /.
Derivative[kk_][ff_] :> D[ff, {z, kk}] ];

(* nonzero starting point *)

ExponentialPowerSum[expr_, z_, n_, n0_] :=
ExponentialPowerSum[expr /. n -> n + n0, z, n] /; n0 >= 0

(* the exponential series *)

ExponentialPowerSum[1, z_, n_] := E^z

(* linearity *)

ExponentialPowerSum[c_ expr_., z_, n_] :=
c ExponentialPowerSum[expr, z, n] /; FreeL[c, {n, Blank}]

ExponentialPowerSum[c_.(expr1_ + expr2_), rest__] :=
ExponentialPowerSum[c expr1, rest] +
ExponentialPowerSum[c expr2, rest]

ExponentialPowerSum[Sum[a_, {k_, lo_, hi_}] b_., z_, n_] :=
Sum[ExponentialPowerSum[Expand[a b], z, n], {k, lo, hi}] /;
FreeQ[{lo, hi}, n] && FreeQ[b, k]

(* powers *)

ExponentialPowerSum[c_^((a_ + b_)d_), rest__] :=
ExponentialPowerSum[c^(a d + b d), rest]

ExponentialPowerSum[c_^(a_ + b_), z_, n_] :=
c^a ExponentialPowerSum[c^b, z, n] /;
FreeL[c, {n, Blank}] && FreeL[a, {n, Blank}]

(* exponentials *)

ExponentialPowerSum[e_. c_^((a_. n_ + b_.)d_.), z_, n_] :=
(c^(b d) ExponentialPowerSum[e, z, n] /. z -> c^(a d) z) /;
(FreeL[#, {n, Blank}]& /@ And[a, b, c, d]) && FreeL[e, {z, Blank}]

(* polynomial inhomogeneity *)

ExponentialPowerSum[n_^m_., z_, n_] :=
Block[{j},
E^z Factor[Sum[StirlingS2[m, j] z^j,
{j, 0, m}] ] ] /; IntegerQ[m] && TTQ[m >= 0]

(* factorial inhomogeneity *)

ExponentialPowerSum[(n_ + a_.)!, z_, n_] :=
PowerSum[Pochhammer[n + 1, a], z, n] /; IntegerQ[a]

(* convolutions *)

ExponentialPowerSum[Sum[expr_ Binomial[n_ + cc_, k_], {k_, lo_, hi_}],
z_, n_] :=
Block[{aa = Expand[lo], bb = Expand[hi - n]},
Derivative[cc][ExponentialPowerSum[Sum[(expr /. n -> n - cc)
Binomial[n, k], {k, aa, n + bb - cc} ], z, n ] ] /;
IntegerQ[cc] && TTQ[cc > 0] && FreeQ[aa, n] && FreeQ[bb, n] ]

ExponentialPowerSum[Sum[expr_ Binomial[n_, k_], {k_, lo_, hi_}],
z_, n_] :=
Block[{nfree, nfull, s, aa = Expand[lo], bb = Expand[hi - n]},
{nfree, nfull} = If[Head[expr] === Times,
{Select[expr, FreeQ[#, n]&], Select[expr, !FreeQ[#, n]&]},
If[FreeQ[expr, n], {expr, 1}, {1, expr}] ];
Sum[z^k nfree/k! ExponentialPowerSum[Expand[nfull /. n -> s + k],
z, s, -k], {k, aa, bb - 1}] +
ExponentialPowerSum[nfree ExponentialPowerSum[Expand[nfull /.
n -> s + k], z, s, -bb], z, k, Max[aa, bb] ] /;
FreeQ[aa, n] && FreeQ[bb, n] ]

(* multisection of series *)

ExponentialPowerSum[If[SymbolicMod[n_ + d_., m_] == k_, a_, b_] c_.,
z_, n_] :=
Block[{l = Mod[k - d, m], jj,
fps = ExponentialPowerSum[Expand[(a - b)c], z, n]},
ExponentialPowerSum[Expand[b c], z, n] +
1/m Sum[Omega[m]^(jj(m - l)) fps /. z -> z Omega[m]^jj,
{jj, 0, m - 1}] // Expand ]

(* other conditionals *)

ExponentialPowerSum[If[cond_, a_, b_] c_., z_, n_] :=
Block[{t, tt, s, k},
t = IneqSolve[cond && n >= 0, n];
tt = IneqSolve[!cond && n >= 0, n];
Plus @@ (Sum[Expand[a c] z^n / n!, {n, First[#], Last[#]}]& /@
t) +
Plus @@ (Sum[Expand[b c] z^n / n!, {n, First[#], Last[#]}]& /@
tt) /.
{Sum[s_. z^n / n!, {n, k_, Infinity}] :>
ExponentialPowerSum[s, {z, n, k}],
Sum[0, {__}] -> 0} ]

(* generating functions of sequences

(here we assume that a[n] is defined for n >= 0) *)

ExponentialPowerSum[a_[n_ + d_.], z_, n_] :=
Derivative[d][EGf[a][z]] /; IntegerQ[d] && TTQ[d >= 0]

ExponentialPowerSum[n_^m_. a_[n_ + d_.], z_, n_] :=
Block[{j, k},
Sum[StirlingS2[m, j] z^j
Derivative[d + j][EGf[a][z]], {j, 0, m}] ] /;
IntegerQ[d] && TTQ[d >= 0] && IntegerQ[m] && TTQ[m >= 0]

ExponentialPowerSum[n_^m_. a_[n_ - 1], z_, n_] :=
Block[{j, k},
Sum[StirlingS2[m, j] z^j
Derivative[j - 1][EGf[a][z]], {j, m}] ] /;
IntegerQ[m] && TTQ[m > 0]

(* GENERATING FUNCTION *)

GeneratingFunction[list1_, list2_, n_, z_, opts___Rule] :=
Block[{temp, recur, conds, unknowns, startValues, start,
rec, con, i, prs},
prs = Parse[list1, list2, n];
If[prs === Fail,
Return[Fail],
{recur, conds, unknowns, startValues} = prs ];
{rec, con} = Reset[recur, conds, unknowns, startValues, 0];
(Release[start /@ unknowns]) = startValues;
temp = FunctionSolve[rec, con, unknowns, n, PowerSum,
z, opts];
If[temp === Fail,
Return[Fail],
temp = Table[#[[i]] z^startValues[[i]],
{i, Length[unknowns]} ]& /@ temp;
Return[temp/.aa_?(MemberQ[unknowns,#]&)[kk_] :>
aa[Expand[kk + start[aa]]] ]]]

(* EXPONENTIAL GENERATING FUNCTION *)

ExponentialGeneratingFunction[list1_, list2_, n_, z_, opts___Rule] :=
Block[{temp, recur, conds, unknowns, startValues, start,
rec, con, prs},
prs = Parse[list1, list2, n];
If[prs === Fail,
Return[Fail],
{recur, conds, unknowns, startValues} = prs ];
{rec, con} = Reset[recur, conds, unknowns, startValues, 0];
(Release[start /@ unknowns]) = startValues;
temp = FunctionSolve[rec, con, unknowns, n, ExponentialPowerSum,
z, opts];
If[temp === Fail,
Return[Fail],
Return[temp/.aa_?(MemberQ[unknowns,#]&)[kk_] :>
aa[Expand[kk + start[aa]]] ]]]

(** (Hypergeometric.m) **)

(* HSOLVE *)

HSolve[eq_, f_, x_Symbol, opts___Rule] :=
Block[{hsol = HSolve1[eq, f, x, opts], ll},
ll = PList[hsol, C[_]];
Return[hsol /. Thread[ll -> Table[C[k], {k, Length[ll]}]]] ]

HSolve1[eq_, f_, x_Symbol, opts___Rule] :=

Block[{lhs, rhs, terms, weights, inhom, minw, theta, ll, llc, rlc,
aList, bList, badB, nnHS, n0, subst, eqn, y, k, pr, temp,
pm, nneg, n1, badB1},

pr = PrecisionHS /.{opts} /.Options[HSolve];
If[Head[eq] === List,
If[Length[eq] > 1, Return[Fail]];
eqn = First[eq],
eqn = eq ];
If[Head[f] === List,
If[Length[f] > 1, Return[Fail]];
y = Head[First[f]],
y = Head[f] ];
lhs = Expand[Numerator[Together[First[eqn] - Last[eqn]]]];
terms = Select[MakeList[lhs, Plus], !FreeQ[#, y]&];
weights = Union[Weight[#, y[x], x]& /@ terms];
minw = Min[weights];
If[!MemberQ[{{0},{0,1}}, weights - minw], Return[Fail]];
lhs = Expand[lhs/x^minw] /. x^k_.Derivative[n_][y][x] ->
x^(k - n) Sum[StirlingS1[n, j] theta^j y[x], {j,0,n}]
// Expand;
inhom = Plus @@ Select[MakeList[lhs, Plus], FreeQ[#, y]&]
// Expand;
pm = PoleMultiplicity[inhom, {x, 0}];
If[!(IntegerQ[pm] && pm >= 0), Return[Fail]];
If[!PolynomialQ[x^pm inhom, x], Return[Fail]];
lhs = (lhs - inhom) / y[x] // Expand;
rhs = Plus @@ Select[MakeList[lhs, Plus], FreeQ[#, x]&];
lhs = Expand[(rhs - lhs)/x];
{llc, rlc} = {LeadingCoef[lhs, theta], LeadingCoef[rhs, theta]};
aList = If[FreeQ[lhs, theta], {}, -#[[1,2]]& /@
Which[pr === Automatic,
{ToRules[Roots[lhs == 0, theta,
Cubics -> False, Quartics -> False] ]} /.
Literal[Roots[a_, b_, c___]] :> NRoots[a, b],
pr === Infinity,
Solve[lhs == 0, theta] /.
Roots -> NRoots,
NumberQ[pr],
Solve[lhs == 0, theta] /.
Literal[Roots[a_, b_]] :> NRoots[a, b, pr] ]];
bList = If[FreeQ[rhs, theta], {}, (1 - #[[1,2]])& /@
Which[pr === Automatic,
{ToRules[Roots[rhs == 0, theta,
Cubics -> False, Quartics -> False] ]} /.
Literal[Roots[a_, b_, c___]] :> NRoots[a, b],
pr === Infinity,
Solve[rhs == 0, theta] /.
Roots -> NRoots,
NumberQ[pr],
Solve[rhs == 0, theta] /.
Literal[Roots[a_, b_]] :> NRoots[a, b, pr] ]];
badB = Select[bList, (IntegerQ[#] && # <= 0)&];
nnHS = If[badB != {}, 1 - Min[badB], 0];
n0 = Max[nnHS, Exponent[inhom, x]];
badB1 = Select[bList, (IntegerQ[#] && # > 0)&];
nneg = If[badB1 != {}, Max[badB1] - 1, 0];
n1 = Max[nneg, pm];
Clear[c];
c[-n1 - 1] = 0;
subst = {ToRules[Reduce[Table[If[n + 1 == 0,
Coefficient[inhom, x, 0] /. x -> 1/x /. x -> 0,
Coefficient[inhom, x^(n+1)]] +
rlc Times @@ (bList + n) c[n+1]
== llc Times @@ (aList + n) c[n], {n, -n1 - 1, n0 - 1}],
Table[c[n], {n, n0, -n1, -1}] ]]};
If[subst === {}, Return[{}] ];

If[Times @@ (Pochhammer[# + nnHS, n0 - nnHS]& /@
aList) =!= 0 && llc =!= 0,

(* then *)

{aList, bList} = {aList, bList} + nnHS;
If[!MemberQ[bList, 1], AppendTo[aList, 1];
AppendTo[bList, 1] ];
bList = Drop[bList, {FirstPos[bList, 1]}];
If[MemberQ[aList, #], aList = Drop[aList,
{FirstPos[aList, #]}];
bList = Drop[bList,
{FirstPos[bList, #]}]
]& /@ bList;
aList = Sort[aList];
bList = Sort[bList];
If[NumberQ[pr],
aList = N[aList, pr];
bList = N[bList, pr] ];

temp = (Sum[c[k] x^k, {k, -n1, n0 - 1}] +
c[n0] (llc/rlc)^(nnHS - n0) (Times @@
(Pochhammer[#, n0 - nnHS]&
/@ bList)) / (Times @@ (Pochhammer[#, n0 - nnHS]& /@
aList)) (n0 - nnHS)!
x^nnHS (HypergeometricF[aList, bList, llc/rlc x] -
Sum[(Times @@ (Pochhammer[#, k]& /@ aList))/
(Times @@ (Pochhammer[#, k]& /@ bList)) (llc/
rlc x)^k / k!,
{k, 0, n0 - nnHS - 1}] ) ) /. subst // Simplify,

(* else *)

{aList, bList} = {aList, bList} + n0;
If[!MemberQ[bList, 1], AppendTo[aList, 1];
AppendTo[bList, 1] ];
bList = Drop[bList, {FirstPos[bList, 1]}];
If[MemberQ[aList, #], aList = Drop[aList,
{FirstPos[aList, #]}];
bList = Drop[bList,
{FirstPos[bList, #]}]
]& /@ bList;
aList = Sort[aList];
bList = Sort[bList];
If[NumberQ[pr],
aList = N[aList, pr];
bList = N[bList, pr] ];

temp = (Sum[c[k] x^k, {k, -n1, n0 - 1}] +
c[n0] x^n0 HypergeometricF[aList, bList, llc/rlc x])
/. subst // Expand ];

temp = temp /. c -> C;
Return[Thread[y[x] -> #]& /@ {temp}] ]

(* WEIGHT *)

Weight[a_. x_^k_. Derivative[n_][y_][x_], y_[x_], x_] := k - n /;
FreeL[a, {x, y}]

Weight[a_. Derivative[n_][y_][x_], y_[x_], x_] := -n /; FreeL[a, {x, y}]

Weight[a_. x_^k_. y_[x_], y_[x_], x_] := k /; FreeL[a, {x, y}]

Weight[a_. y_[x_], y_[x_], x_] := 0 /; FreeL[a, {x, y}]

Weight[a_, y_[x_], x_] := Infinity /; FreeQ[a, y]

(* HYPERGEOMETRIC F *)

HypergeometricF[aList_List, bList_List, c_.z_] :=
Block[{max, kk, neg = Select[aList, (IntegerQ[#] && # <= 0)&]},
(max = - Max[neg];
Sum[Times @@ (Pochhammer[#, kk]& /@ aList) /
Times @@ (Pochhammer[#, kk]& /@ bList) / kk! (c z)^kk,
{kk, 0, max}]) /; Length[neg] > 0 ]

HypergeometricF[a_List, b_List, 0] := 1

HypergeometricF[{}, {}, z_] := E^z

HypergeometricF[{a_}, {}, z_] := 1/(1 - z)^a

HypergeometricF[{}, {c_}, z_] := Hypergeometric0F1[c, z]

HypergeometricF[{a_}, {c_}, z_] := Hypergeometric1F1[a, c, z]

HypergeometricF[{a_, b_}, {c_}, z_] := Hypergeometric2F1[a, b, c, z]

HypergeometricF /:
Series[HypergeometricF[aList_List, bList_List, c_.z_],
{z_, 0, n_}] :=
Block[{kk},
Sum[Times @@ (Pochhammer[#, kk]& /@ aList) /
Times @@ (Pochhammer[#, kk]& /@ bList) / kk! (c z)^kk,
{kk, 0, n}] + O[z]^(n + 1) ]

Derivative[0, 0, n_Integer?Positive][HypergeometricF][aList_List,
bList_List, c_.z_] :=
c^n Times @@ (Pochhammer[#, n]& /@ aList) /
Times @@ (Pochhammer[#, n]& /@ bList) HypergeometricF[
aList + n, bList + n, c z]

(** (SeriesExpansions.m) **)

(* SERIES TERM *)

SeriesTerm[expr_List, rest__] := SeriesTerm[#, rest]& /@ expr

SeriesTerm[f_, {z_, a_, n_}, opts___Rule] :=
Block[{ff = f /. z -> z + a, nnST, useApart, n0, temp0, temp,
precision, makeReal, useMod, specialFunctions},
IntegerQ[nnST] ^= True;
Info[nnST] ^= Info[n] /. Solve[n == nnST, PatternList[n,
_Symbol?(Context[#] =!= "System`"&) ]][[1]];
sumDepth = 1;
useApart = UseApart /.{opts} /.Options[SeriesTerm];
precision = PrecisionST /.{opts} /.Options[SeriesTerm];
makeReal = MakeReal /.{opts} /.Options[SeriesTerm];
useMod = UseMod /.{opts} /.Options[SeriesTerm];
specialFunctions = SpecialFunctions /.{opts} /.
Options[SeriesTerm];
Which[useApart === Automatic,
ff = Apart[ff],
useApart === All,
ff = Apart[ff, z],
useApart =!= None,
Message[ST::badopt, UseApart, useApart ] ];
temp = SeriesTerm[ff, z, nnST] /. {nnST -> n};
temp = Chop[temp] /. SimplifyFloat;
temp = temp //. SimplifyTrig;
temp = temp /. SimplifyBinomial;
temp = temp /.
{Power[aa_, b_] :> Power[aa, Expand[b]],
Binomial[aa_, b_] :> Binomial[Expand[aa], Expand[b]],
If[aa_, b_, c_] :> If[ISolve[aa], b, c] /;
FreeL[aa, {OddQ, EvenQ, SymbolicMod}] };
temp = temp /. aa_If :> Release //@ aa;
temp = temp /.
{If[OddQ[m_], If[m_ >= aa_, b_, c_] d_., e_] :>
If[OddQ[m], Release[Expand[b d]], e] /;
Entails[Info[m], m >= aa - 1] && OddQ[aa],
If[EvenQ[m_], If[m_ >= aa_, b_, c_] d_., e_] :>
If[EvenQ[m], Release[Expand[b d]], e] /;
Entails[Info[m], m >= aa - 1] && EvenQ[aa] };
temp = temp //. SimplifySum;
If[NumberQ[precision],
temp = N[temp, precision] /. SimplifyFloat ];
If[MatchQ[Info[n], n >= _],
n0 = Info[n][[2]];
temp0 = temp /. If[n >= aa_, b_, _] :> b /; aa == n0 + 1;
If[(temp0 /. n -> n0) == (temp /. n -> n0), temp = temp0] ];
Return[temp] ]

(* finite n *)

SeriesTerm[f_, z_Symbol, 0] := SeriesTerm[z f, z, 1] /;
FreeQ[f, HypergeometricF]

SeriesTerm[f_, z_Symbol, n_Integer] :=
Coefficient[SafeSeries[f, {z, 0, n}], z^n] /;
n != 0 && FreeQ[f, HypergeometricF]

(* linearity *)

SeriesTerm[f_Plus, rest__] := SeriesTerm[#, rest]& /@ f

SeriesTerm[c_ f_, z_Symbol, rest__] := c SeriesTerm[f, z, rest] /;
FreeQ[c, z]

(* SeriesTerm is inverse of Gf and of PowerSum *)

SeriesTerm[Gf[a_][z_], z_Symbol, n_] := a[n]

SeriesTerm[PowerSum[a_, z_, n_], z_Symbol, m_] :=
a /. n -> m

(* SeriesTerm is close to the inverse of EGf and of
ExponentialPowerSum *)

SeriesTerm[EGf[a_][z_], z_Symbol, n_] := a[n] / n!

SeriesTerm[ExponentialPowerSum[a_, z_, n_], z_Symbol, m_] :=
(a / n!) /. n -> m

(* binomials *)

SeriesTerm[(a_.z_ + b_)^p_, z_Symbol, n_] :=
Block[{k},
(b^p When[p < 0,
(-a/b)^n Binomial[n - p - 1, n],
(a/b)^n Binomial[p, n]]) /;
(FreeQ[#, z]&) /@ (a && b && p)]

SeriesTerm[(a_.z_^m_. + b_)^p_, z_Symbol, n_] :=
Cancel[a^p SeriesTerm[z^(m p)(1 + b/a z^(-m))^p, z, n]] /;
(FreeQ[#, z]&) /@ (a && b && m && p) && TTQ[m < 0]

SeriesTerm[(a_.z_^m_ + b_)^p_, z_Symbol, n_] :=
If[SymbolicMod[n, m] == 0,
Cancel[b^p SeriesTerm[(1 + a/b z)^p, z, n/m]],
0] /; (FreeQ[#, z]&) /@ (a && b && m && p) && TTQ[m >= 0] &&
Head[p] =!= Integer

SeriesTerm[(a_.z_^m_. + b_.z_^k_.)^p_, z_Symbol, n_] :=
Cancel[b^p SeriesTerm[z^(k p)(1 + a/b z^(m-k))^p, z, n]] /;
(FreeQ[#, z]&) /@ (a && b && m && k && p)

(* hypergeometrics *)

SeriesTerm[HypergeometricF[aList_, bList_, c_.z_], z_Symbol, n_] :=
When[n >= 0,
(Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
Pochhammer[#, n]]& /@ aList))/
(Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
Pochhammer[#, n]]& /@ bList))/n! c^n, 0] /; FreeQ[c, z]

SeriesTerm[Hypergeometric2F1[a_, b_, c_, d_.z_], z_Symbol, n_] :=
When[n >= 0,
(Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
Pochhammer[#, n]]& /@ {a, b}))/
(Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
Pochhammer[#, n]]& /@ {c, 1})) d^n, 0] /; FreeQ[d, z]

SeriesTerm[Hypergeometric1F1[a_, b_, c_.z_], z_Symbol, n_] :=
When[n >= 0,
If[IntegerQ[a], (a + n - 1)!/(a - 1)!, Pochhammer[a, n]]/
(Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
Pochhammer[#, n]]& /@ {b, 1})) c^n, 0] /; FreeQ[c, z]

SeriesTerm[Hypergeometric0F1[a_, b_.z_], z_Symbol, n_] :=
When[n >= 0,
1/(Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
Pochhammer[#, n]]& /@ {a, 1})) b^n, 0] /; FreeQ[b, z]

(* exponentials *)

SeriesTerm[E^(k_.z_ + m_.), z_Symbol, n_] := E^m k^n/n! /;
FreeQ[k, z] && FreeQ[m, z]

SeriesTerm[a_^b_, z_Symbol, n_] := SeriesTerm[E^(b Log[a]), z, n] /;
a =!= E && FreeQ[a, z]

(* logarithms *)

SeriesTerm[Log[(a_ + b_.z_)^k_.], z_Symbol, n_] :=
When[n == 0, Release[k Log[a]], Release[-k(-b/a)^n / n]] /;
FreeQ[k, z] && FreeQ[a, z] && FreeQ[b, z]

SeriesTerm[Log[a_Times], rest___] :=
Plus @@ (SeriesTerm[Log[#], rest]& /@ a)

(* powers and constants *)

SeriesTerm[z_^k_.f_, z_Symbol, n_] := SeriesTerm[f, z, n-k] /;
IntegerQ[k] && FreeQ[k, z] && !FreeQ[f, z]
SeriesTerm[z_^k_., z_Symbol, n_] := When[n == k, 1, 0] /;
IntegerQ[k] && FreeQ[k, z]
SeriesTerm[f_, z_Symbol, n_] := f When[n == 0, 1, 0] /; FreeQ[f, z]

(* rational functions *)

SeriesTerm[f_. g_^k_Integer?Negative, z_Symbol, n_] :=
Block[{ff = Apart[f g^k, z]},
If[Head[ff] === Plus,
Return[SeriesRat[#, z, n]& /@ ff],
Return[SeriesRat[ff, z, n]] ]] /;
PolynomialQ[f, z] && PolynomialQ[g, z] && !FreeQ[g, z]

SeriesRat[f_, z_Symbol, n_] := SeriesTerm[f, z, n] /; PolynomialQ[f, z]

SeriesRat[c_. z_^k_Integer?Negative, z_Symbol, n_] :=
SeriesTerm[c z^k, z, n] /; FreeQ[c, z]

SeriesRat[f_. g_^k_Integer?Negative, z_Symbol, n_] :=

Block[{ff = Expand[f], gg = Expand[g], poleList, poleSet,
mp, denom = Expand[g^(-k)], i, deg, mod, modList,
displace, temp, l, pr = precision, inverseRoot, kk},

deg = Exponent[gg, z];
If[useMod,
l = CoefficientList[ff, z];
mod = Exponent[gg, z, GCD];
modList = Union[Mod[#, mod]& /@ (Select[Range[Length[l]],
(l[[#]] =!= 0)&] - 1)];
If[Length[modList] == 1,
displace = First[modList],
{displace, mod} = {0, 1}],
{displace, mod} = {0, 1}];
{ff, gg, denom} = {Expand[ff/z^displace], gg, denom} /.
z -> z^(1/mod);
deg = deg/mod;
Which[precision === Automatic && deg > 2,
pr = Indeterminate,
precision === Automatic && deg <= 2,
pr = Infinity];
If[pr == Infinity && !FreeQ[Roots[gg == 0, z], Roots],
pr = Indeterminate];
Which[
pr === Indeterminate,
poleList = Cancel[#[[1,2]]& /@ {ToRules[
NRoots[gg == 0, z] ]}];
poleSet = Union[poleList];
mp = -k Count[poleList, #]& /@ poleSet;
poleUse = poleSet;
temp = SeriesTerm[Plus @@ Table[SeriesDivide[Horner[
ff, {z, poleSet[[i]], mp[[i]]}], Drop[Horner[denom, {z,
poleSet[[i]], 2 mp[[i]]}], mp[[i]]], mp[[i]] ] .
Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
{i, Length[poleSet]}], z, Expand[(n-displace)/mod]],

pr =!= Infinity,
poleList = Cancel[#[[1,2]]& /@ {ToRules[
NRoots[gg == 0, z, pr] ]}];
poleSet = Union[poleList];
mp = -k Count[poleList, #]& /@ poleSet;
poleUse = poleSet;
temp = SeriesTerm[Plus @@ Table[SeriesDivide[Horner[
ff, {z, poleSet[[i]], mp[[i]]}], Drop[Horner[denom, {z,
poleSet[[i]], 2 mp[[i]]}], mp[[i]]], mp[[i]] ] .
Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
{i, Length[poleSet]}], z, Expand[(n-displace)/mod]],

pr === Infinity && deg <= 2,
poleList = Cancel[#[[1,2]]& /@ {ToRules[
Roots[gg == 0, z] ]}];
poleSet = Union[poleList];
mp = -k Count[poleList, #]& /@ poleSet;
inverseRoot = Expand[(1 - gg/(gg /. z -> 0))/z];
poleUse = ArgPi[Apart[
inverseRoot /. z -> #]]^(-1)& /@ poleSet;
temp = SeriesTerm[Plus @@ Table[SeriesDivide[Horner[
ff, {z, poleSet[[i]], mp[[i]]}], Drop[Horner[denom, {z,
poleSet[[i]], 2 mp[[i]]}], mp[[i]]], mp[[i]] ] .
Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
{i, Length[poleSet]}], z, Expand[(n-displace)/mod]],

pr === Infinity && deg > 2,
poleList = Cancel[#[[1,2]]& /@ {ToRules[
Roots[CoefficientList[gg, z].
Table[z^kk, {kk, deg, 0, -1}] == 0, z] ]}];
poleSet = Union[poleList];
mp = -k Count[poleList, #]& /@ poleSet;
poleUse = Complex[Together[ReP[#]//.SimplifyCubic],
Together[ImP[#]//.SimplifyCubic]]& /@ poleSet;
poleUse = poleUse /. Complex[a_, 0] -> a;
If[FreeQ[poleUse, Complex], makeReal = False];
temp = SeriesTerm[Plus @@ Table[SeriesDivide[Horner[
ff /. z -> z / poleUse[[i]], {z, 1, mp[[i]]}],
Drop[Horner[denom /. z -> z / poleUse[[i]], {z,
1, 2 mp[[i]]}], mp[[i]]], mp[[i]] ] .
Table[(z poleUse[[i]] - 1)^(-j), {j, mp[[i]], 1, -1}],
{i, Length[poleSet]}], z, Expand[(n-displace)/mod]]
];
temp = temp //. a_ c_If + b_ c_ -> (a + b) c;
((temp = temp /. SimplifyComplexE1) /; pr === Infinity;
(temp = temp /. SimplifyComplex2) /; pr === Infinity && deg == 2;
(temp = temp /. SimplifyComplex3) /; pr === Infinity && deg == 3;
(temp = temp //. SimplifyComplex4) /; pr === Infinity && deg == 4;
(temp = temp /. SimplifyComplexE2) /; pr === Infinity;
(temp = temp //. SimplifyComplexN) /; pr =!= Infinity) /;
makeReal;
temp = temp /. Complex[a_, b_] :> a + b I;
Return[If[Release[SymbolicMod[n, mod] == Mod[displace, mod]],
Release[temp], 0]]
] /; PolynomialQ[f, z] && PolynomialQ[g, z] && !FreeQ[g, z]

(* special functions *)

SeriesTerm[1 / Sqrt[a_ + b_.z_ + c_.z_^2], z_Symbol, n_] :=
When[n >= 0,
Sqrt[c/a]^n / Sqrt[a] LegendreP[n, -Sgn[a] b /(2 Sqrt[a c])],
0] /; specialFunctions && FreeQ[{a,b,c}, z]

SeriesTerm[Power[a_ + b_.z_ + c_.z_^2, Rational[p_, 2]], z_Symbol, n_] :=
SeriesTerm[Expand[(a + b z + c z^2)^((p + 1)/2)] /
MySqrt[a + b z + c z^2], z, n] /;
p != -1 && specialFunctions && FreeQ[{a,b,c}, z]

SeriesTerm[1 / MySqrt[a_], z_Symbol, n_] :=
SeriesTerm[1 / Sqrt[a], z, n]

SeriesTerm[(a_ + b_.z_)^m_ (c_ + d_.z_)^m_, z_Symbol, n_] :=
SeriesTerm[Expand[(a + b z) (c + d z)]^m, z, n] /;
specialFunctions && FreeQ[{a,b,c,d}, z] &&
Expand[a c]^m == Expand[a^m c^m]

SeriesTerm[(a_ + b_.z_)^m_ (c_ + d_.z_)^m_, z_Symbol, n_] :=
-SeriesTerm[Expand[(a + b z) (c + d z)]^m, z, n] /;
specialFunctions && FreeQ[{a,b,c,d}, z] &&
Expand[a c]^m == -Expand[a^m c^m]

(* products with a rational function *)

SeriesTerm[a_ b_, z_Symbol, n_] :=
Block[{temp, deg, i, m, bn = SeriesTerm[b, z, m], sd = sumDepth},
Increment[sumDepth];
Off[General::itervar];
If[PolynomialQ[a, z],
deg = Exponent[a, z];
temp = CoefficientList[a, z] .
Table[bn /. m -> i, {i, n, n - deg, -1}] /. If -> When,
temp = Sum[SeriesTerm[Apart[a, z], z, n - K[sd]] bn /.
m -> K[sd], {K[sd],
-PoleMultiplicity[b, {z, 0}], n + PoleMultiplicity[a,
{z, 0}] }] ];
On[General::itervar];
Return[temp] ] /; SeparateProduct[a b, (PolynomialQ[#, z] !!
PolynomialQ[#^(-1), z])& ] == {a, b}

(* general products *)

SeriesTerm[a_ b_, z_Symbol, n_] :=
Block[{temp, sd = sumDepth},
Increment[sumDepth];
Off[General::itervar];
temp = Sum[SeriesTerm[a, z, n - K[sd]] *
SeriesTerm[b, z, K[sd]],
{K[sd], -PoleMultiplicity[b, {z, 0}],
n + PoleMultiplicity[a, {z, 0}]}];
On[General::itervar];
Return[temp] ] /; SeparateProduct[a b, (PolynomialQ[#, z] !!
PolynomialQ[#^(-1), z])& ] == {1, a b}

SeriesTerm[a_^k_Integer?Positive, z_Symbol, n_] :=
Block[{temp, sd = sumDepth},
Increment[sumDepth];
Off[General::itervar];
temp = Sum[SeriesTerm[a, z, n - K[sd]] *
SeriesTerm[a^(k - 1), z, K[sd]],
{K[sd], -PoleMultiplicity[a^(k - 1), {z, 0}],
n + PoleMultiplicity[a, {z, 0}]}];
On[General::itervar];
Return[temp] ] /; !FreeQ[a, z]

(* trinomials *)

SeriesTerm[(a_ + b_.z_ + c_.z_^2)^alpha_, z_Symbol, n_] :=
Block[{temp},
Off[General::itervar];
temp = Sum[Binomial[alpha, K[sumDepth]] Binomial[K[sumDepth],
n - K[sumDepth]] a^(alpha - K[sumDepth]) b^(2 K[sumDepth] -
n) c^(n - K[sumDepth]), {K[sumDepth], 0, n}];
On[General::itervar];
++sumDepth;
Return[temp] ] /; FreeQ[{a,b,c}, z]

(* derivatives *)

SeriesTerm[Derivative[k_][f_][z_], z_Symbol, n_] :=
Pochhammer[n + 1, k] SeriesTerm[f[z], z, n + k]

(* HORNER and SERIES DIVIDE *)

Horner[p_, {z_, alpha_, m_}] :=
Block[{qh = Expand[p], ih, kh, nh, ah},
{nh, ah} = {Exponent[qh, z], CoefficientList[qh, z]};
Do[ah[[ih]] = alpha ah[[ih + 1]] + ah[[ih]], {kh, 1, Min[nh, m]},
{ih, nh, kh, -1}];
Do[AppendTo[ah, 0], {ih, m - nh - 1}] /; m > nh + 1;
Return[Take[ah, m]] ]

SeriesDivide[a_, b_, n_] :=
Block[{is, cs = Range[n], ks},
Do[cs[[is]] = (a[[is]] - Sum[cs[[ks]] b[[is - ks + 1]],
{ks, 1, is - 1}])/b[[1]],
{is, 1, n}];
Return[cs] ]

(* SYMBOLIC MOD *)

SymbolicMod[n_, 1] = 0

SymbolicMod[n_?NumberQ, k_] := Mod[n, k]

SymbolicMod /: (SymbolicMod[n_Symbol + a_., 2] == b_) := EvenQ[n] /;
EvenQ[b - a]

SymbolicMod /: (SymbolicMod[n_Symbol + a_., 2] == b_) := OddQ[n] /;
OddQ[b - a]

(* PARTIAL FRACTIONS *)

PartialFractions[f_, z_, opts___Rule] :=
Block[{ff = f, useApart},
useApart = UseApart /.{opts} /.Options[PartialFractions];
Which[useApart == Automatic,
ff = Apart[ff],
useApart == All,
ff = Apart[ff, z],
useApart =!= None,
Message[PF::badopt, UseApart, useApart ] ];
If[Head[ff] === Plus,
Return[PartFrac[#, z, opts]& /@ ff],
Return[PartFrac[f, z, opts]] ]];

PartFrac[f_, z_, opts___Rule] :=
Block[{num = Numerator[f], den = Denominator[f], poleList, poleSet,
mp, i, j, k, deg, temp, precision},
precision = PrecisionST /.{opts} /. Options[PartialFractions];
(deg = Exponent[den, z];
Which[precision === Automatic && deg > 2,
precision = Indeterminate,
precision === Automatic && deg <= 2,
precision = Infinity];
poleList = Cancel[#[[1,2]]& /@ {ToRules[Which[
precision == Indeterminate, NRoots[den == 0, z],
precision =!= Infinity, NRoots[den == 0, z, precision],
precision === Infinity, Roots[den == 0, z]] ]}];
poleSet = Union[poleList];
mp = Count[poleList, #]& /@ poleSet;
If[precision === Infinity,
poleUse = ArgPi[#]& /@ poleSet,
poleUse = poleSet];
temp = Plus @@ Table[SeriesDivide[Horner[num, {z,
poleSet[[i]], mp[[i]]}], Drop[Horner[den, {z,
poleSet[[i]], 2 mp[[i]]}], mp[[i]]], mp[[i]] ] .
Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
{i, Length[poleSet]}]
) /; PolynomialQ[num, z] && PolynomialQ[den, z] ];

(* REALQ, REP, IMP, ANGLE, SGN, ABSV, and CONJUGATEQ *)

RealQ[x_] := True /; IntegerQ[x]
RealQ[x_Rational] = True
RealQ[x_Real] = True
RealQ[Pi] = True
RealQ[E] = True
RealQ[n_!] := True /; RealQ[n]

ReP[x_] := x /; RealQ[x]
ImP[x_] := 0 /; RealQ[x]

ReP[k_Integer x_] := k ReP[x]
ReP[k_Rational x_] := k ReP[x]
ReP[k_Real x_] := k ReP[x]

ImP[k_Integer x_] := k ImP[x]
ImP[k_Rational x_] := k ImP[x]
ImP[k_Real x_] := k ImP[x]

ReP[Complex[x_, y_]] := x
ImP[Complex[x_, y_]] := y

ReP[x_ + y_] := ReP[x] + ReP[y]
ImP[x_ + y_] := ImP[x] + ImP[y]

ReP[x_ Sqrt[y_?Positive]] := ReP[x] Sqrt[y]
ImP[x_ Sqrt[y_?Positive]] := ImP[x] Sqrt[y]
ReP[x_ Sqrt[y_?Negative]] := -ImP[x] Sqrt[-y]
ImP[x_ Sqrt[y_?Negative]] := ReP[x] Sqrt[-y]

ReP[x_ y_] := ReP[x] ReP[y] - ImP[x] ImP[y]
ImP[x_ y_] := ReP[x] ImP[y] + ImP[x] ReP[y]

ReP[x_?Positive ^n_] := x^n /; ImP[n] == 0
ImP[x_?Positive ^n_] := 0 /; ImP[n] == 0

ReP[x_?Negative ^n_Rational] := 0 /; IntegerQ[2n]
ImP[x_?Negative ^n_Rational] := (-x)^n (-1)^(n-1/2) /;
IntegerQ[2n]

ReP[1/x_] := ReP[x] / (ReP[x]^2 + ImP[x]^2)
ImP[1/x_] := -ImP[x] / (ReP[x]^2 + ImP[x]^2)

ReP[x_^Rational[p_,q_]] :=
ReP[x^Quotient[p,q] x^(Mod[p,q]/q)] /;
AbsV[p] >= AbsV[q]

ReP[x_^q_Rational] := AbsV[x]^q Cos[q Angle[x]] /;
AbsV[q] <= 1
ImP[x_^q_Rational] := AbsV[x]^q Sin[q Angle[x]] /;
AbsV[q] <= 1

ReP[E^x_] := Cos[ImP[x]] Exp[ReP[x]]
ImP[E^x_] := Sin[ImP[x]] Exp[ReP[x]]

ReP[x_^2] := ReP[x]^2 - ImP[x]^2
ImP[x_^2] := 2 ReP[x] ImP[x]

ReP[x_^3] := ReP[x]^3 - 3 ReP[x] ImP[x]^2
ImP[x_^3] := 3 ReP[x]^2 ImP[x] - ImP[x]^3

ReP[x_^n_Integer] := Block[{a, b},
a = Round[n/2];
b = n-a;
ReP[x^a] ReP[x^b] - ImP[x^a] ImP[x^b]
] /; n != -1

ImP[x_^n_Integer] := Block[{a, b},
a = Round[n/2];
b = n-a;
ReP[x^a] ImP[x^b] + ImP[x^a] ReP[x^b]
] /; n != -1

Angle[z_] :=
Block[{x = ReP[z], y = ImP[z]},
Which[x == 0, Pi/2 Sgn[y],
y == 0, Pi/2 (1 - Sgn[x]),
True,
ArcTan[y/x] - Pi/2 (Sgn[x]-1) Sgn[y]] /;
FreeQ[{x,y}, ReP] && FreeQ[{x,y}, ImP]]

Sgn[0] = 0
Sgn[a_?Positive] = 1
Sgn[a_?Negative] = -1
Sgn[a_ b_] := Sgn[a] Sgn[b]
Sgn[a_^b_] := Sgn[a]^b
Sgn[Sqrt[a_]] = 1
Sgn[a_ + b_] := Sgn[Chop[N[a + b]]]

AbsV[0] = 0
AbsV[x_?Positive] := x
AbsV[x_?Negative] := -x
AbsV[x_ y_] := AbsV[x] AbsV[y]
AbsV[x_^n_] := AbsV[x]^n /; Head[n] =!= Complex
AbsV[x_] :=
Block[{a = ReP[x], b = ImP[x]},
Which[b == 0, Expand[a Sgn[a]],
a == 0, Expand[b Sgn[b]],
True, Sqrt[Expand[a^2 + b^2]]]]

ConjugateQ[a_, b_] :=
Block[{l = PatternList[{a, b},
_Symbol?(Context[#] =!= "System`"&) ]},
Which[
l === {},
Chop[N[ReP[a]] - N[ReP[b]]] == 0 &&
Chop[N[ImP[a]] + N[ImP[b]]] == 0,
l =!= {},
Expand[ReP[a] - ReP[b]] == 0 &&
Expand[ImP[a] + ImP[b]] == 0 ] ]

(* SIMPLIFICATION RULES *)

SimplifyTrig = {
Sin[x_. + r_Rational Pi] :> Cos[x] /; EvenQ[r-1/2],
Sin[x_. + Pi] :> -Sin[x],
Sin[x_. + n_Integer?OddQ Pi] :> -Sin[x],
Sin[x_. + r_Rational Pi] :> -Cos[x] /; EvenQ[r-3/2],
Sin[x_. + n_Integer?EvenQ Pi] :> Sin[x],
Cos[x_. + r_Rational Pi] :> -Sin[x] /; EvenQ[r-1/2],
Cos[x_. + Pi] :> -Cos[x],
Cos[x_. + n_Integer?OddQ Pi] :> -Cos[x],
Cos[x_. + r_Rational Pi] :> Sin[x] /; EvenQ[r-3/2],
Cos[x_. + n_Integer?EvenQ Pi] :> Cos[x],
Sin[a_?Negative b_.] :> -Sin[-a b],
Cos[a_?Negative b_.] :> Cos[-a b],
Sin[a_Plus] :> -Sin[Expand[-a]] /; MatchQ[a[[1]], _?Negative _.],
Cos[a_Plus] :> Cos[Expand[-a]] /; MatchQ[a[[1]], _?Negative _.],
ArcTan[a_?Negative b_.] :> -ArcTan[-a b],
ArcTan[a_Plus] :> -ArcTan[Expand[-a]] /;
MatchQ[a[[1]], _?Negative _.],
ArcTan[Sqrt[3]] -> Pi/3,
ArcTan[Sqrt[3]/3] -> Pi/6,
ArcTan[1/Sqrt[3]] -> Pi/6}

SimplifySqrt = {a_^Rational[b_,2] :> a^((b-1)/2) Sqrt[a]}

CanonicRadicals = {
(a_ + b_)^Rational[k_,n_] :> a (a + b)^(Mod[k, n] / n) +
b (a + b)^(Mod[k, n] / n) /; Quotient[k, n]==1,
a_^Rational[k_,n_] :> a^(Mod[k, n]/n) Expand[a^Quotient[k, n]]}

SimplifyCubic = Join[SimplifyTrig, SimplifySqrt]

SimplifyBinomial = {
Binomial[1/2, k_ + 1] :> (-4)^(-k) Binomial[2 k, k]/(2(k + 1)) /;
Entails[Info[k], k != -1],
Binomial[1/2, k_] -> (-4)^(-k) Binomial[2 k, k]/(1 - 2 k),
Binomial[-1/2, k_] -> (-4)^(-k) Binomial[2 k, k],
Binomial[k_ - 1/2, k_] -> 4^(-k) Binomial[2 k, k]}

SimplifySum = {
Sum[0, {__}] -> 0,
Sum[1, {k_, a_, b_}] :> When[a <= b, b - a + 1, 0],
Sum[a_. b_, {k_, c__}] :> b Sum[a, {k, c}] /;
FreeQ[b, k],
Sum[If[k_ == m_, a_, b_] c_., {k_, m_, n_}] :>
((a c) /. k -> m) When[n >= m, 1, 0] +
Sum[b c, {k, m + 1, n}],
Sum[(If[-k_ + n_ >= 0, a_, b_] c_. + d_.) e_., {k_, m_, n_}] :>
Sum[Release[Expand[a c e]], {k, m, n}] + Expand[d e],
Sum[If[k_ >= m_, a_, b_] c_., {k_, m_, n_}] :>
Sum[Release[Expand[a c]], {k, m, n}] }

SimplifyComplexE1 = {
a_.I^c_ + d_.(-I)^c_ :> ArgPi[a] E^Expand[I Pi c/2] +
ArgPi[d] E^Expand[-I Pi c/2] /; ConjugateQ[a, d],
a_.E^b_ + d_.E^e_ :> 2 a Cos[ImP[b]] /;
(ImP[a] == a - d == ReP[b] == b + e == 0),
a_.E^b_ + d_.E^e_ :> -2 ImP[a] Sin[ImP[b]] /;
(ReP[a] == a + d == ReP[b] == b + e == 0),
a_.E^b_ + d_.E^e_ :> ArgPi[a] E^b + ArgPi[d] E^e}

SimplifyComplex2 = {
a_.b_^c_ + d_.e_^c_ :>
2 AbsV[a]AbsV[b]^c SimplifyCos[Cos[Angle[a] + c Angle[b]]] /;
(ConjugateQ[a, d] && ConjugateQ[b, e] && ImP[c] == 0)}

SimplifyComplex3 = {
a_.b_Complex^c_ + d_.e_Complex^c_ :>
2 AbsV[a]AbsV[b]^c Cos[Angle[a] + c Angle[b]] /;
(ConjugateQ[a, d] && ConjugateQ[b, e] && ImP[c] == 0)}

SimplifyComplex4 = {
a_.b_Complex^c_ + d_.e_Complex^c_ :>
2 AbsV[a]AbsV[b]^c SimplifyCos[Cos[Angle[a] + c Angle[b]]] /;
(ConjugateQ[a, d] && ConjugateQ[b, e] && ImP[c] == 0)}

SimplifyComplexE2 = {
a_.E^b_ + d_.E^e_ :> Block[{c = ImP[b]},
2 ReP[a] Cos[c] - 2 ImP[a] Sin[c] /;
(ConjugateQ[a, d] && ConjugateQ[b, e] && ReP[b] == 0)]}

SimplifyComplexN = {
a_.b_Complex^c_ + d_.e_Complex^c_ :>
2 Abs[a]Abs[b]^c Cos[Arg[a] + c Arg[b]] /;
(ConjugateQ[a, d] && ConjugateQ[b, e] && ImP[c] == 0)}

SimplifyCos[x_] :=
Block[{x1 = (x /. Cos[a_] :> Cos[Expand[a]]) //. SimplifyTrig,
x2 = x //. SimplifyTrig },
If[LeafCount[x1] <= LeafCount[x2], x1, x2 ]]

SimplifyFloat = {1. -> 1, -1. -> -1, 0. -> 0}

(* FACTORIAL SIMPLIFY *)

FactorialSimplify[expr_] := PostSimplify[DeFact[expr]]

DeFact[expr_] :=

Block[{kin, j, sub, diff, t, tt, arg = {}, ex = expr},

While[!FreeQ[ex, Factorial],
t = First[Position[ex, Factorial[_]]];
tt = First[Part @@ Prepend[t, ex]];
kin = Scan[If[IntegerQ[Expand[# - tt]],
Return[Fact[#]] ]&, arg];
If[kin === Null,
AppendTo[arg, tt]; sub = Fact[tt],
diff = Expand[First[kin] - tt];
If[diff > 0,
sub = kin / Product[tt + j, {j, diff}],
sub = kin Product[tt - j + 1, {j, -diff}] ] ];
ex = MapAt[sub&, ex, {t}] ];
Factor[ex] /. Fact -> Factorial ]

PostSimplify[expr_] :=

Block[{ex = expr /. n! :> Expand[n]!, i},

ex //.
{n_!^k_. m_!^l_.:> Product[m + i, {i, n - m}]^k /;
Expand[n - m] > 0 && k + l == 0,

n_!^k_. m_^l_. :> Expand[n + 1]!^k Expand[n + 1]^(l - k) /;
Expand[m - n - 1] == 0 && Sign[k] == Sign[l],

n_!^k_. m_^l_. :> (-1)^k Expand[n + 1]!^k m^(l - k) /;
Expand[m + n + 1] == 0 && Sign[k] == Sign[l],

n_!^k_. n_^l_. :> Expand[n - 1]!^k /; Expand[k + l] == 0,

n_!^k_. m_^l_. :> (-1)^k Expand[n - 1]!^k /;
Expand[m + n] == Expand[k + l] == 0 } ]

(** (Utilities.m) **)

(* ARGPI *)

ArgPi[z_] :=
Block[{fraction = Rationalize[N[Arg[z]/Pi]]},
If[Head[fraction] === Rational,
Return[AbsV[z] E^(fraction Pi I)],
Return[z]] ]

(* ENTAILS *)

Entails[a_, a_] = True

Entails[a___ && b_ && c___, b_] = True

Entails[False, a_] = True

Entails[a_, b_] :=
Block[{l = Union[PatternList[{a, b},
_Symbol?(Context[#] =!= "System`"&) ]]},
IneqSolve[Implies[a, b], First[l]] ===
{Interval[-Infinity, Infinity]} /; Length[l] == 1]

(* FIRST POS and FREE L *)

FirstPos[l_, pattern_] :=
SafeFirst[Select[Range[Length[l]], MatchQ[l[[#]], pattern]& ]]

FreeL[expr_, l_] := And @@ (FreeQ[expr, #]& /@ l)

(* ISOLVE, INEQSOLVE, LEQ *)

IneqSolve[x_ + a_ > b_, x_] := IneqSolve[x > b - a, x]

IneqSolve[x_ + a_ >= b_, x_] := IneqSolve[x >= b - a, x]

IneqSolve[x_ > n_Integer, x_] := {Interval[n+1, Infinity]}

IneqSolve[x_ >= n_Integer, x_] := {Interval[n, Infinity]}

IneqSolve[a_, x_] := a /; Length[Union[PatternList[a,
_Symbol?(Context[#] =!= "System`"&) ]]] > 1

IneqSolve[a_ && b_, x_] :=
Block[{aa = IneqSolve[a, x],
bb = IneqSolve[b, x],
merge, accum, begin, nn},
If[FreeQ[aa, IneqSolve] && FreeQ[bb, IneqSolve],
aa = Append[aa, {}];
bb = Append[bb, {}];
merge = Sort[Join[
Thread[{Join @@ (aa /. Interval -> List),
Join @@ (aa /. Interval[_,_] -> {-1,1}) }],
Thread[{Join @@ (bb /. Interval -> List),
Join @@ (bb /. Interval[_,_] -> {-1,1}) }] ],
LEQ];
accum = Thread[{First /@ merge,
Accumulate[Plus, Last /@ merge]}];
begin = Join @@ Append[Position[accum, {_,-2}], {}];
Return[(Interval @@ # & /@ Thread[{begin, begin + 1}]) /.
nn_Integer :> First[accum[[nn]]]],
Return[aa && bb] ]]

IneqSolve[!a_, x_] :=
Block[{aa = IneqSolve[a, x], merge, c, d},
If[FreeQ[aa, IneqSolve],
aa = Append[aa, {}];
merge = Join @@ (aa /. Interval[c_, d_] -> {c-1, d+1});
If[merge === {},
merge = {-Infinity, Infinity},
If[First[merge] === -Infinity,
merge = Rest[merge],
PrependTo[merge, -Infinity] ];
If[Last[merge] === Infinity,
merge = Drop[merge,-1],
AppendTo[merge, Infinity] ] ];
Return[Interval @@ # & /@ Partition[merge, 2]],
Return[!aa]] ]

IneqSolve[a_ !! b_, x_] :=
Block[{aa = IneqSolve[!(!a && !b), x]},
If[FreeQ[aa, IneqSolve],
Return[aa],
Return[IneqSolve[a, x] !! IneqSolve[b, x]]] ]

IneqSolve[Inequality[a_, b_, c_, d_, e__], x_] :=
IneqSolve[And @@ (Inequality @@ # & /@
Partition[{a,b,c,d,e},3,2]), x]

IneqSolve[h_[a_, b_, c__], x_] :=
IneqSolve[And @@ (h @@ # & /@
Partition[{a,b,c},2,1]), x] /;
MemberQ[{Less, LessEqual, Greater, GreaterEqual, Equal,
Unequal}, h]

IneqSolve[ineq: h_[a_, b_], x_] :=
Block[{lhs = Together[a - b],
num, den, zeroList, poleList, points, merge, c, d},

{num, den} = {Numerator[lhs], Denominator[lhs]};
Off[Roots::neq];
zeroList = If[FreeQ[num, x], {},
#[[1,2]]& /@ {ToRules[NRoots[num == 0, x]]} ];
poleList = If[FreeQ[den, x], {},
#[[1,2]]& /@ {ToRules[NRoots[den == 0, x]]} ];
On[Roots::neq];
zeroList = Select[zeroList, (Im[#] == 0)&] // Union;
poleList = Select[poleList, (Im[#] == 0)&] // Union;
points = Union[zeroList, poleList];
merge = Partition[Append[Prepend[points,-Infinity],Infinity],
2, 1];
merge = Select[merge, (ineq /. x :>
Which[# === {-Infinity, Infinity}, 0,
#[[1]] === -Infinity, #[[2]] - 1,
#[[2]] === Infinity, #[[1]] + 1,
True, (Plus @@ #)/2 ] )&];
merge = (Interval @@ # & /@ merge) /.
Interval[c_, d_] :> Interval[Ceiling[c], Floor[d]];
merge = merge /. {Ceiling[-Infinity] -> -Infinity,
Floor[Infinity] -> Infinity};
merge = Select[merge, #[[1]] <= #[[2]]&];
zeroList = Sort[Join[Floor /@ zeroList, Ceiling /@ zeroList]];
zeroList = Select[zeroList, ((den /. x->#) != 0)&];
zeroList = Select[zeroList, (ineq /. x->#)&];
merge = Sort[Join[merge, Thread[Interval[zeroList, zeroList]]],
LEQ];
merge = merge /.
Interval[c_?((NumberQ[#] && ((den /. x:>#) == 0 !!
!ineq/.x->#))&), d_] :> Interval[c+1, d];
merge = merge /.
Interval[c_, d_?((NumberQ[#] && ((den /. x:>#) == 0 !!
!ineq/.x->#))&)] :> Interval[c, d-1];
merge = Select[merge, #[[1]] <= #[[2]]&];
args = List @@ (Join @@ Append[merge, Interval[]]);
If[args === {}, {},
{first, last} = {First[args], Last[args]};
Interval @@ # & /@ Partition[Append[Prepend[Join @@
Append[Select[Partition[Drop[Rest[args], -1], 2],
(#[[2]] - #[[1]] > 1)&], {}], first], last], 2] ] ] /;
MemberQ[
{Less, LessEqual, Greater, GreaterEqual, Equal, Unequal}, h]

ISolve[expr_] :=
Block[{l = UserSymbols[expr], n, temp, p, nn},
If[Length[l] == 1,
n = First[l];
temp = IneqSolve[expr /. n -> nn, nn] /. nn -> n;
temp = temp /. List[p___Interval] :> Or[p];
temp = temp /.
{Interval[-Infinity, Infinity] -> True,
Interval[-Infinity, b_] -> n <= b,
Interval[a_, Infinity] -> n >= a,
Interval[a_, a_] -> n == a,
Interval[a_, b_] -> a <= n <= b};
temp = temp /. IneqSolve[a_, n_] -> a;
Return[temp],
Return[expr]] ]

LEQ[h_[a_,b_], h_[c_,d_]] := (a < c) !! (a === c && b <= d)

(* INFO *)

Info[_?NumberQ] = True
Info[_Symbol] = True
Info[_[a___]] := And @@ (Info /@ {a})

(* LEADING COEF *)

LeadingCoef[poly_, x_] := If[poly === 0, 0,
Coefficient[poly, x, Exponent[poly, x]] ]

(* MAKE LIST *)

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

(* MAKE TRINOMIAL and LIST COMPLEMENT *)

MakeTrinomial[a: Binomial[b_, c_] Binomial[d_, e_]] :=
Block[{l1 = {b, d},
l2 = {c, b - c, e, d - e}, l},
l = Intersection[l1, l2];
If[Length[l] != 1, Return[a]];
l1 = ListComplement[l1, l];
l2 = ListComplement[l2, l];
If[{Plus @@ l2} != l1, Return[a]];
Return[Sort[Multinomial @@ l2]] ]

MakeTrinomial[a_] := a

ListComplement[l1_List, l2_List] :=
Block[{l = l1},
If[MemberQ[l, #],
l = Drop[l, {FirstPos[l, #]}] ]& /@ l2;
Return[l] ]

(* PARSING ROUTINE *)

Parse[list1_, list2_, n_] :=
Block[{l1 = MakeList[Release @@ (Hold[list1] /. Condition -> Pair)],
unknowns = MakeList[list2],
eqns, conds, recur, extra, range, startValues, aa, kk,
lo, hi, check},

unknowns = Head /@ unknowns;

check = If[Head[#] === Pair, #[[1]], #]& /@ l1;
If[!MatchQ[check, {__Equal}],
Message[Parse::eqn, check];
Return[Fail] ];

eqns = Select[l1, !FreeQ[#, n]&];
conds = Select[l1, FreeQ[#, n]&];

(* make n >= 0 the default range of validity *)

eqns = If[Head[#] =!= Pair, Pair[#, n >= 0], #]& /@ eqns;

(* solve inequalities, and separate recurrences from
initial conditions *)

eqns = Pair[#[[1]], IneqSolve[#[[2]], n]]& /@ eqns;
recur = Select[eqns, !FreeQ[#, Infinity]&];
extra = Union[Select[eqns, FreeQ[#, Infinity]&],
Pair[#[[1]], Drop[#[[2]], -1]]& /@ recur];
recur = Pair[#[[1]], Last[#[[2]]][[1]]]& /@ recur;
extra = Select[extra, #[[2]] =!= {} &];
extra = Union @@ (Thread /@ extra);
conds = Union[conds, Union @@ (Table[#[[1]], {n, #[[2,1]],
#[[2,2]]}]& /@ extra) ];

(* shift n so that all recurrences will be valid
for n >= 0 *)

recur = (#[[1]] /. (n -> n + #[[2]]) )& /@ recur;

(* determine starting values of n for all sequences *)

{recur, conds} = {recur, conds} //.
Sum[a_, {k_, m_}] :> Sum[a, {k, 1, m}];
{recur, conds} = {recur, conds} /.
{EvenQ[a_] :> SymbolicMod[a, 2] == 0,
OddQ[a_] :> SymbolicMod[a, 2] == 1,
Mod -> SymbolicMod};
range = {conds, recur /. (Sum[aa_, {kk_, lo_, hi_}] :>
{aa /. kk -> lo, aa /. kk -> hi} )};
startValues = Min[First /@ (PatternList[range, #[_]] /.
n -> 0)]& /@ unknowns;

Return[{recur, conds, unknowns, startValues}] ]

(* PATTERN LIST *)

PatternList[expr_?(FreeQ[#, Rule]&), pattern_] :=
Apply[Part, Prepend[#, expr]]& /@ Position[expr, pattern]

PList[expr_, pattern_] :=
Block[{xpr = expr /. ((a_ -> b_) -> {a,b})},
Apply[Part, Prepend[#, xpr]]& /@ Position[xpr, pattern]]

(* POLE MULTIPLICITY *)

PoleMultiplicity[f_, {z_, a_}] :=
Block[{pom = SafeSeries[f /. z -> z + a, {z, 0, 0}]},
If[!FreeL[pom, {Fail, Series}], Return[Infinity],
pom = Normal[pom];
If[pom == 0, Return[0]];
Return[Max[0, Exponent[Expand[pom /. z -> 1/z], z]]] ]]

(* RESET *)

Reset[recur_, conds_, unknowns_, startValues_, s0_] :=
Block[{start},
(Release[start /@ unknowns]) = startValues;
Return[{recur, conds} /.
aa_?(MemberQ[unknowns,#]&)[kk_] :>
aa[kk + s0 - start[aa]] ]]

(* SAFE FIRST, SAFE SERIES *)

SafeFirst[l_] :=
If[Length[l] == 0, Null, First[l]]

SafeSeries[f_, {z_, a_, n_Integer}] :=
Block[{poles = z /. Solve[Denominator[f] == 0, z,
InverseFunctions -> True], ord, tem},
ord = Max[n, Count[poles, a]];
tem = Check[Series[f, {z, a, ord}], Fail, Series::esss];
If[tem =!= Fail, tem = Series[tem, {z, a, n}]];
Return[tem] ]

(* SEPARATE PRODUCT *)

SeparateProduct[p_Times, predicate_] :=
{Select[p, predicate], Select[p, !predicate[#]&]}

(* USER SYMBOLS *)

UserSymbols[expr_] :=
Block[{h = If[Length[expr] == 0, expr, Head[expr]]},
Which[!FreeQ[h,
_Symbol?(Context[#] =!= "System`"&)],
{expr},
Length[expr] == 0,
{},
True,
Union @@ (UserSymbols /@ (List @@ expr))]]

(* TTQ, WHEN *)

TTQ[True] := True
TTQ[a_] := Entails[Info[a], a]

When[cond_, a_, b_] :=
Which[Entails[Info[cond], cond], a,
Entails[Info[cond], !cond], b,
True, If[cond, a, b]]

End[]

(* Protect[RSolve, PowerSum, ExponentialPowerSum, GeneratingFunction,
ExponentialGeneratingFunction, Gf, EGf, HSolve, HypergeometricF, K,
SeriesTerm, SymbolicMod, RealQ, ReP, ImP, ConjugateQ, SimplifySum,
SimplifyTrig, SimplifyComplexE1, SimplifyComplexE2, SimplifyComplex2,
SimplifyComplex3, SimplifyComplex4, SimplifyComplexN, FactorialSimplify,
PartialFractions, ArgPi, Entails, FirstPos, FreeL, Info, ISolve,
LeadingCoef, ListComplement, MakeList, MakeTrinomial, PatternList,
PoleMultiplicity, Reset, SafeFirst, SafeSeries, TTQ, UserSymbols, When,
Methods, MethodGF, MethodEGF, PrecisionHS, PrecisionST, UseApart, UseMod,
MakeReal, SpecialFunctions] *)

EndPackage[]