home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / DISCRETE.PAK / RSOLVE.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  85.4 KB  |  2,466 lines
  1.  
  2. (* :Title: RSolve *)
  3.  
  4. (* :Author:  Marko Petkovsek *)
  5.  
  6. (* :Summary: Recurrence (difference) equation solver *)
  7.  
  8. (* :Context: DiscreteMath`RSolve` *)
  9.  
  10. (* :Package Version: 1.0 *)
  11.  
  12. (* :Copyright: Copyright 1990 by Marko Petkovsek *)
  13.  
  14. (* :History:
  15.   Alpha test version by Marko Petkovsek, June 7, 1990
  16.   Modified by ECM, Wolfram Research Inc., Jan., 1991
  17. *)
  18.  
  19. (* :Keywords:
  20.   recurrence, difference, discrete, solve, equation 
  21. *)
  22.  
  23. (* :Warning: Binomial and Power (specifically, 0^0) are redefined. *)
  24.  
  25. (* :Mathematica Version: 2.0 *)
  26.  
  27. (* :Limitation: none known. *)
  28.  
  29. (* :Discussion: *)
  30.  
  31.  
  32. BeginPackage["DiscreteMath`RSolve`"]
  33.  
  34.  
  35. RSolve::usage =
  36. "RSolve[{eqn1, eqn2, ..}, {a1[n], a2[n], ..}, n, opts] solves the recurrence
  37. equations eqn1, eqn2, .. for the sequences a1[n], a2[n], ... If there is a
  38. single sequence or equation, it need not be given in a list. Equations can
  39. either be recurrences or initial/boundary conditions. An equation may have 
  40. the form eqn /; cond where cond is an inequality specifying the range 
  41. of values of n for which eqn is valid; when no condition is given 
  42. it is assumed to be n >= 0."
  43.  
  44. PowerSum::usage =
  45. "PowerSum[expr, {z, n, n0:0}] computes Sum[expr z^n, {n, n0, Infinity}]. 
  46. If a[n] is a sequence, the name Gf[a][z] is used by PowerSum to denote
  47. Sum[a[n] z^n, {n, 0, Infinity}].   PowerSum is threaded over lists and
  48. equations."
  49.  
  50. ExponentialPowerSum::usage =
  51. "ExponentialPowerSum[expr, {z, n, n0:0}] computes Sum[expr z^(n-n0) / (n-n0)!,
  52. {n, n0, Infinity}]. If a[n] is a  sequence, the name EGf[a][z] is used by
  53. ExponentialPowerSum to denote Sum[a[n] z^n / n!, {n, 0, Infinity}].
  54. ExponentialPowerSum is threaded over lists and equations."
  55.  
  56. GeneratingFunction::usage =
  57. "GeneratingFunction[{eqn1, eqn2, ..}, {a1[n], a2[n], ..}, n, z, opts] computes
  58. the ordinary generating functions Sum[ai[n] z^n, {n, si, Infinity}] for the
  59. sequences a1[n], a2[n], .. defined by the equations eqn1, eqn2, ... Here si
  60. denotes the least value of n such that ai[n] appears in the equations. An 
  61. equation may have the form eqn /; cond where cond is an inequality specifying 
  62. the range of values of n for which eqn is valid; when no condition 
  63. is given it is assumed to be n >= 0."
  64.  
  65. ExponentialGeneratingFunction::usage =
  66. "ExponentialGeneratingFunction[{eqn1, eqn2, ..}, {a1[n], a2[n], ..}, n, z]
  67. computes the exponential generating functions Sum[ai[n + si] z^n / n!,
  68. {n, 0, Infinity}] for the sequences a1[n], a2[n], .. defined by the
  69. equations eqn1, eqn2, ... Here si denotes the least value of n such that
  70. ai[n] appears in the equations.   An equation may have the form eqn /; cond
  71. where cond is an inequality specifying the range of values of n for which eqn
  72. is valid; when no condition is given it is assumed to be n >= 0."
  73.  
  74. Gf::usage =
  75. "Gf[a][z] is used by PowerSum to denote Sum[a[n] z^n, {n, 0, Infinity}]." 
  76.  
  77. EGf::usage = "EGf[a][z] is used by ExponentialPowerSum to denote 
  78. Sum[a[n] z^n / n!, {n, 0, Infinity}]."
  79.  
  80. HSolve::usage = "HSolve[eq, f[z], z] returns a hypergeometric formal 
  81. series solution to differential equation eq with unknown function 
  82. f[z], or {} when it determines that there is no series solution.
  83. It returns $Failed when it is unable to do either."
  84.  
  85. HypergeometricF::usage = "HypergeometricF[{a1, a2, ..}, {b1, b2, ..}, z]
  86. is the generalized hypergeometric function with upper parameters
  87. a1, a2,.., lower parameters b1, b2,.., and argument z."
  88.  
  89. K::usage = "K is the head used for summation indices by SeriesTerm."
  90.  
  91. SeriesTerm::usage = "SeriesTerm[expr, {z, a, n}, opts] 
  92. returns the n-th coefficient of the Laurent series expansion of expr 
  93. around z = a, for symbolic n. Alias: ST."
  94.  
  95. SymbolicMod::usage = "SymbolicMod[n, k] is equal to Mod[n, k] when n
  96. is a number, and is left unevaluated otherwise."
  97.  
  98. RealQ::usage = "'RealQ[x] = True' serves to declare x to be real."
  99.  
  100. ReP::usage = "ReP[x] gives the real part of x."
  101.  
  102. ImP::usage = "ImP[x] gives the imaginary part of x."
  103.  
  104. ConjugateQ::usage = "ConjugateQ[x, y] tests whether x and y are complex
  105. conjugates."
  106.  
  107. SimplifySum::usage = "SimplifySum is a list of rules for simplification 
  108. of sums."
  109.  
  110. SimplifyTrig::usage = "SimplifyTrig is a list of rules that put 
  111. trigonometric expressions into canonical form."
  112.  
  113. SimplifyComplexE1::usage = "SimplifyComplexE1 is one of the lists of rules 
  114. used to replace sums of conjugate complex quantities by their double
  115. real parts."
  116.  
  117. SimplifyComplexE2::usage = "SimplifyComplexE2 is one of the lists of rules 
  118. used to replace sums of conjugate complex quantities by their double
  119. real parts."
  120.  
  121. SimplifyComplex2::usage = "SimplifyComplex2 is one of the lists of rules 
  122. used to replace sums of conjugate complex quantities by their double
  123. real parts."
  124.  
  125. SimplifyComplex3::usage = "SimplifyComplex3 is one of the lists of rules 
  126. used to replace sums of conjugate complex quantities by their double
  127. real parts."
  128.  
  129. SimplifyComplex4::usage = "SimplifyComplex4 is one of the lists of rules 
  130. used to replace sums of conjugate complex quantities by their double
  131. real parts."
  132.  
  133. SimplifyComplexN::usage = "SimplifyComplexN is one of the lists of rules 
  134. used to replace sums of conjugate complex quantities by their double
  135. real parts."
  136.  
  137. FactorialSimplify::usage = "FactorialSimplify[expr] simplifies 
  138. occurrences of Factorial[] inside expr."
  139.  
  140. PartialFractions::usage = "PartialFractions[f, z, opts] gives a partial
  141. fraction decomposition of f[z] over the complex functions."
  142.  
  143. ArgPi::usage = "When possible, ArgPi[z] expresses the complex number z 
  144. in the form Abs[z] E^(r Pi I) where r is a rational number." 
  145.  
  146. Entails::usage = "Entails[a, b] determines whether a implies b, where
  147. a and b are Boolean combinations of equalities and/or inequalities among
  148. rational functions of a single integer variable."
  149.  
  150. FirstPos::usage = "FirstPos[l, pattern] returns the position of
  151. the first argument of l that matches pattern, or Null if there is
  152. no such element."
  153.  
  154. FreeL::usage = "FreeL[expr, list] returns True if none of the 
  155. elements of list appears in expr, False otherwise."
  156.  
  157. Info::usage = "Info[n] is information about n, expressed as a 
  158. Boolean combination of inequalities. It should be associated with n."
  159.  
  160. ISolve::usage = "ISolve[expr] solves a Boolean combination of equalities 
  161. and/or inequalities in a single integer variable. The answer is given 
  162. as a disjunction of disjoint inclusive inequalities."
  163.  
  164. LeadingCoef::usage = "LeadingCoef[poly, x] returns the leading
  165. coefficient of the polynomial poly[x]."
  166.  
  167. ListComplement::usage = "ListComplement[list1, list2] returns 
  168. the multiset-difference of l1 and l2."
  169.  
  170. MakeList::usage = "MakeList[expr, head:List] returns the list of
  171. arguments of expr if the head of expr matches the given one, otherwise,
  172. it returns {expr}."
  173.  
  174. MakeTrinomial::usage = "MakeTrinomial[expr_] tries to write expr
  175. in the form Multinomial[a, b, c]."
  176.  
  177. PatternList::usage = "PatternList[expr, pattern] returns the list of
  178. all subexpressions of expr that match the given pattern, provided
  179. that expr does not contain Rule[]."
  180.  
  181. PoleMultiplicity::usage = "PoleMultiplicity[f, {z, a}] computes
  182. the multiplicity of the pole of f at z = a."
  183.  
  184. Reset::usage = "Reset[recur, conds, unknowns, startValues, s0] returns
  185. the list {recur, conds} with starting values of the index reset so 
  186. that it starts with s0 rather than with startValues, for all unknown 
  187. sequences."
  188.  
  189. SafeFirst::usage = "SafeFirst[l] returns the first argument of l,
  190. or Null if the length of l is zero."
  191.  
  192. SafeSeries::usage = "SafeSeries[f, {z, a, n}] overcomes a bug in 
  193. Series[] that fails when n is less than the degree of zero of the
  194. denominator of f at a."
  195.  
  196. TTQ::usage = "TTQ[cond] returns True if it determines that the
  197. information given about the variables implies that cond is true, 
  198. and False otherwise."
  199.  
  200. UserSymbols::usage = "UserSymbols[expr] returns a list of maximal 
  201. subexpressions of expr with the property that their head contains 
  202. symbols defined in a context different from System`. "
  203.  
  204. When::usage = "When[cond, a, b] returns a if it determines that the
  205. information given about the variables implies that cond is true, 
  206. b if it determines that the information given about the variables 
  207. implies that cond is false, and If[cond, a, b] otherwise."
  208.  
  209. Methods::usage = "Methods is an optional parameter for RSolve.
  210. It specifies which methods and in what order are to be used by RSolve.
  211. Default: Methods -> {MethodGF, MethodEGF}. If a single method is given
  212. it need not be enclosed in a list."
  213.  
  214. MethodGF::usage = "MethodGF is a possible value for the Methods
  215. parameter for RSolve. It denotes the method of ordinary generating functions."
  216.  
  217. MethodEGF::usage = "MethodEGF is a possible value for the Methods parameter 
  218. for RSolve. It denotes the method of exponential generating functions."
  219.  
  220. PrecisionHS::usage = "PrecisionHS is an optional parameter for RSolve,
  221. GeneratingFunction, ExponentialGeneratingFunction, and HSolve. It 
  222. specifies the precision with which the roots of polynomials are to be 
  223. computed by HSolve. When successful, the result of HSolve is a 
  224. generalized hypergeometric function, and these roots are its parameters. 
  225. Possible values: PrecisionHS -> n, PrecisionHS -> Infinity, 
  226. PrecisionHS -> Automatic (default). Infinity means exact roots, but 
  227. those roots that cannot be found by Solve[] are computed numerically. 
  228. Automatic means that only those roots that Solve[] can find without 
  229. recourse to general formulas for solving cubics and quartics are found 
  230. exactly, while others are computed numerically."
  231.  
  232. PrecisionST::usage = "PrecisionST is an optional parameter for RSolve
  233. and SeriesTerm. It specifies the precision with which the roots of 
  234. denominators of rational functions to be expanded into power series
  235. are to be computed by SeriesTerm. Possible values: PrecisionST -> n,
  236. PrecisionST -> Infinity, PrecisionST -> Automatic (default). Infinity 
  237. means exact roots, but those roots that cannot be found by Solve[] 
  238. are computed numerically. Automatic means that only those roots which 
  239. Solve[] can find without recourse to general formulas for solving cubics 
  240. and quartics are found exactly, while others are computed numerically."
  241.  
  242. UseApart::usage = "UseApart is an optional parameter for RSolve
  243. and SeriesTerm. It specifies whether Apart should be used
  244. before attempting expansion into power series, and in what form. 
  245. Possible values: UseApart -> None (don't use Apart at all), 
  246. UseApart -> Automatic (default; use Apart without second argument),
  247. UseApart -> All (use Apart with second argument). By default, MethodGF
  248. uses UseApart -> Automatic and MethodEGF uses UseApart -> None.
  249. UseApart -> All is useful when the function to be expanded contains
  250. parameters."
  251.  
  252. UseMod::usage = "UseMod is an optional parameter for RSolve and
  253. SeriesTerm. Possible values: UseMod -> True (default; use Even, Odd 
  254. and/or SymbolicMod to express the result), UseMod -> False (the 
  255. opposite)."
  256.  
  257. MakeReal::usage = "MakeReal is an optional parameter for RSolve and
  258. SeriesTerm. Possible values: MakeReal -> True (default; replace pairs 
  259. of conjugate complex quantities by their doble real parts), 
  260. MakeReal -> False (the opposite)."
  261.  
  262. SpecialFunctions::usage = "SpecialFunctions is an optional parameter 
  263. for RSolve and SeriesTerm. Possible values: SpecialFunctions -> True
  264. (default; use special functions such as Legendre polynomials to express
  265. the result), SpecialFunctions -> False (the opposite)."
  266.  
  267. Even::usage =
  268. "Even[n] is equivalent to EvenQ[n] except that it remains unevaluated when
  269. the head of n is Symbol."
  270.  
  271. Odd::usage =
  272. "Odd[n] is equivalent to OddQ[n] except that it remains unevaluated when the
  273. head of n is Symbol."
  274.  
  275.  
  276.  
  277. Begin["`Private`"]
  278.  
  279.  
  280.  
  281. (** ALIASES **)
  282.  
  283. RS = RSolve
  284. PS = PowerSum
  285. EPS = ExponentialPowerSum
  286. GF = GeneratingFunction
  287. EGF = ExponentialGeneratingFunction
  288. ST = SeriesTerm
  289.  
  290.  
  291.  
  292. (** ATTRIBUTES **)
  293.  
  294. Attributes[RSolve] = {HoldFirst}
  295. Attributes[iRSolve] = {HoldFirst}
  296. Attributes[GeneratingFunction] = {HoldFirst}
  297. Attributes[ExponentialGeneratingFunction] = {HoldFirst}
  298. Attributes[Parse] = {HoldFirst}
  299.  
  300.  
  301.  
  302. (** OPTIONS **)
  303.  
  304. Options[RSolve] := Join[{Methods ->{MethodGF, MethodEGF}},
  305.                          Options[HSolve],
  306.                          Options[SeriesTerm] ]
  307. Options[GeneratingFunction] := Options[HSolve]
  308. Options[ExponentialGeneratingFunction] := Options[HSolve]
  309. Options[HSolve] = {PrecisionHS -> Automatic}
  310. Options[SeriesTerm] = {PrecisionST -> Automatic, 
  311.                        UseApart -> Automatic,
  312.                        MakeReal -> True,
  313.                        UseMod -> True,
  314.                        SpecialFunctions -> True}
  315. Options[PartialFractions] := Options[SeriesTerm] 
  316.  
  317.  
  318.  
  319.  
  320. (** MESSAGES **)
  321.  
  322.  
  323. ST::badopt = "`` -> `` is not a valid option."
  324. PF::badopt = "`` -> `` is not a valid option."
  325. Parse::eqn = "`` is not an equation or a system of equations."
  326.  
  327.  
  328.  
  329.  
  330. (** REDEFINITIONS OF SYSTEM FUNCTIONS **)
  331.  
  332. Unprotect[Binomial, Power]
  333.     Binomial[n_, n_] := When[n >= 0, 1, 0]
  334.     0^0 = 1
  335. Protect[Binomial, Power]
  336.  
  337.  
  338.  
  339.  
  340. (* RSOLVE *)
  341.  
  342. RSolve[list1_, list2_, n_, opts___Rule] :=
  343.     Block[{result = iRSolve[list1, list2, n, opts]},
  344.         result /; result =!= $Failed
  345.     ]
  346.  
  347. iRSolve[list1_, list2_, n_, opts___Rule] :=
  348.     Block[{recur, conds, unknowns, startValues, methods, result, prs},
  349.         prs = Parse[list1, list2, n];
  350.         If[prs === $Failed,
  351.             Return[$Failed],
  352.           {recur, conds, unknowns, startValues} = prs ];
  353.         methods = MakeList[Methods /.{opts} /.Options[RSolve]];
  354.         result = Scan[
  355.          Block[{temp},
  356.                  temp = # @@ {recur, conds, unknowns, startValues, n, opts};
  357.                  If[temp =!= $Failed,
  358.             Return[ Map[Thread[list2->#]&,temp] ]
  359.         ]
  360.          ]&,
  361.          methods];
  362.         Which[result === Null,  
  363.                   Return[$Failed],
  364.               True, 
  365.         Return[result] ]
  366.     ]
  367.  
  368.  
  369. (* METHOD GF *)
  370.  
  371. MethodGF[recur_, conds_, unknowns_, startValues_, n_, opts___Rule] :=
  372.     Block[{genf, mm, z, rec, con, i, temp, start, optsList = {opts}},
  373.        {rec, con} = Reset[recur, conds, unknowns, startValues, 0];
  374.         genf = FunctionSolve[rec, con, unknowns, n, PowerSum,
  375.                              z, opts];
  376.         If[genf =!= $Failed,
  377.            (Evaluate[start /@ unknowns]) = startValues;
  378.             If[FreeQ[optsList, UseApart],
  379.                  AppendTo[optsList, UseApart -> Automatic]
  380.         ];
  381.             temp = (Table[mm /: Info[mm] = mm >= startValues[[i]];
  382.                           SeriesTerm[#[[i]], 
  383.                                      {z, 0, mm - startValues[[i]]},
  384.                      Sequence @@ optsList
  385.                     ], 
  386.                           {i, Length[unknowns]}
  387.             ]& /@ genf) /. mm -> n; 
  388.             Return[temp /. aa_?(MemberQ[unknowns, #]&)[kk_] :> 
  389.                            aa[Expand[kk + start[aa]]] ],
  390.             Return[$Failed]
  391.     ]
  392.     ]
  393.             
  394.             
  395. (* METHOD EGF *)
  396.  
  397. MethodEGF[recur_, conds_, unknowns_, startValues_, n_, opts___Rule] :=
  398.     Block[{genf, mm, z, rec, con, i, temp, start, optsList = {opts}},
  399.        {rec, con} = Reset[recur, conds, unknowns, startValues, 0];
  400.         genf = FunctionSolve[rec, con, unknowns, n, ExponentialPowerSum,
  401.                              z, opts];
  402.         If[genf =!= $Failed,
  403.            (Evaluate[start /@ unknowns]) = startValues;
  404.             If[FreeQ[optsList, UseApart],
  405.                  AppendTo[optsList, UseApart -> None] ];
  406.             temp = (Table[mm /: Info[mm] = mm >= startValues[[i]];
  407.                 FactorialSimplify[
  408.                (mm - startValues[[i]])! SeriesTerm[#[[i]], 
  409.                {z, 0, mm - startValues[[i]]}, Sequence @@ optsList] ], 
  410.                {i, Length[unknowns]} ]& /@ genf) /. mm -> n;
  411.             Return[temp /. aa_?(MemberQ[unknowns, #]&)[kk_] :> 
  412.                            aa[Expand[kk + start[aa]]] ],
  413.           Return[$Failed] ]]
  414.             
  415.                                            
  416.                                
  417.              
  418. (** (GeneratingFunctions.m) **)
  419.  
  420.            
  421. (* ROOTS OF UNITY *)
  422.  
  423. Omega[n_] := Exp[2 Pi I / n]
  424.  
  425.  
  426.  
  427. (* SIMPLIFICATION OF EXP[LOG[_]] AND I*SQRT[_] *)
  428.  
  429. SimplifyExpLog = {E^((a_. Log[b_] + c_.)d_.) -> b^(a d) E^(c d)}
  430.  
  431. SimplifyI = {Complex[0, a_] Power[b_, Rational[n_, 2]] :> 
  432.                               a (-1)^((n-1)/2) Expand[-b]^(n/2)}
  433.  
  434.  
  435. (* SOLVING ROUTINE *)
  436.  
  437. FunctionSolve[recur_, conds_, unknowns_, n_, Transform_, z_, 
  438.               opts___Rule] :=
  439.     Block[{eqns, extra, range, start, aa, kk, nm,
  440.            lo, hi, functions, fnHeads, initials, constants, temp, 
  441.            series, order, xx, ii, bb, initvals, soln, times,
  442.            arg, G, init, subst, pinitials},
  443.                        
  444. (* transform recurrence equations into functional equations *)      
  445.  
  446.         eqns = Expand[#[[1]] - #[[2]]]& /@ recur;
  447.         eqns = Function[xx,
  448.                    Block[{hom, inhom, yy},
  449.                        yy = MakeList[xx, Plus];
  450.                        inhom = Select[yy, FreeL[#, unknowns]&];
  451.                        hom = Together[Plus @@ Complement[yy, inhom]];
  452.                        Numerator[hom] + Expand[(Plus @@ 
  453.                            inhom) Denominator[hom]] ]] /@ eqns;
  454.         eqns = Transform[eqns, z, n];
  455.         If[!FreeL[PatternList[eqns, Transform[__]], unknowns], 
  456.             Return[$Failed]];
  457.         Which[Transform === PowerSum,            G = Gf,
  458.               Transform === ExponentialPowerSum, G = EGf];
  459.               
  460. (* substitute values explicitly given by initial conditions *)
  461.  
  462.         init = Select[conds, Function[xx, MatchQ[xx,
  463.             a_?(MemberQ[unknowns, #]&)[_?NumberQ] == _?NumberQ]]];
  464.         subst = Solve[init, Union @@ (PatternList[init, #[_]]& /@ 
  465.             unknowns)];
  466.         If[subst === {}, Return[$Failed], eqns = eqns /. First[subst]];
  467.             
  468. (* remember the unknown functions *)
  469.                
  470.         functions = G[#][z]& /@ unknowns;
  471.         fnHeads = G /@ unknowns;
  472.   
  473. (* solve the equations using Solve, DSolve, or HypergeometricSolve *)
  474.     
  475.         If[FreeQ[eqns, Derivative],
  476.         
  477.             temp = Solve[Thread[eqns == 0], functions],
  478.             
  479.     (* before using differential equation solvers, integrate
  480.        equations of the form  lhs' == rhs'  *)
  481.        
  482.             times = Integrable[#, fnHeads, z]& /@ eqns;
  483.             eqns = PreIntegrate[#[[1]], fnHeads, {z, #[[2]]}]& /@
  484.                 Thread[List[eqns, times]];
  485.             eqns = Table[Sum[C[ii, jj] z^(jj - 1), {jj, times[[ii]]}] +
  486.                 eqns[[ii]], {ii, Length[eqns]} ];
  487.             eqns = eqns /. Derivative[kk_][ff_] :> D[ff, {z, kk}];
  488.             eqns = Thread[eqns == 0];
  489.             
  490.     (* use HSolve1 then DSolve *)
  491.             
  492.         temp = HSolve1[eqns, functions, z, opts];
  493.             temp = temp /. C -> First[unknowns];
  494.             If[!FreeL[temp, {HSolve, $Failed}],
  495.                 temp = DSolve[eqns, functions, z] 
  496.         ]
  497.     ];
  498.         If[temp === {}, Return[{}]];
  499.         If[!FreeL[temp, {Solve, DSolve, $Failed}], Return[$Failed]];
  500.         If[!FreeL[functions /. temp, fnHeads], Return[$Failed]];
  501.         temp = temp //. SimplifyExpLog;
  502.  
  503. (* find all initial terms occurring in conditions, equations,
  504.    or solutions *)
  505.                
  506.         initials = Union @@ (PatternList[{conds, eqns, 
  507.             functions /. temp}, #[_]]& /@ unknowns);
  508.             
  509.             (* sort initials by falling indices *)
  510.             
  511.         initials = Join @@ (Append[Cases[initials, _[#]]& /@
  512.             Reverse[Union[PatternList[initials, _?NumberQ]]], {}]);
  513.         pinitials = Select[initials, (#[[1]] >= 0)&];
  514.             
  515.             (* How many series terms will be needed to determine 
  516.                the initials? *)
  517.                
  518.        (order[#] = Max[Append[First /@ PatternList[initials, 
  519.             #[_]], 0 ]])& /@ unknowns;
  520.             
  521. (* determine values of initials and eliminate constants
  522.    of integration *)
  523.                
  524.         constants = PatternList[functions /. temp, C[__]] // Union;
  525.  
  526.         series = temp /. (G[aa_?(MemberQ[unknowns, #]&)][z] -> bb_
  527.             ) :> (aa -> Taylor[bb, z, 0, order[aa]]);
  528.         If[!FreeL[series, {Series, $Failed}], Return[$Failed]];
  529.         If[G === EGf,
  530.             series = series /. (aa_?(MemberQ[unknowns, #]&) -> bb_
  531.                 ) :> (aa -> Table[nm!, {nm, 0, order[aa]}] bb) ];
  532.         initvals = Function[xx, Solve[Union[conds,
  533.            (# == (Head[#] /. xx)[[First[#] + 1]])& /@ pinitials,
  534.             Thread[Join @@ Rest /@ CList[#[[2,1]]& /@ xx /. 
  535.                 z -> z^(-1), z ] == 0 ] ],
  536.             Join[constants, initials] ]] /@ series;
  537.             
  538. (* plug in initial values *)
  539.  
  540.         soln = Select[Union @@ Table[(Function[xx, 
  541.             G[xx][z] /. temp[[ii]] ] /@ 
  542.             unknowns) /. initvals[[ii]],
  543.            {ii, Length[temp]} ], (Head[#] === List)& ];
  544.         Return[soln] ]
  545.         
  546.        
  547. (* TAYLOR and CLIST *)
  548.  
  549. Taylor[f_, z_, a_, n_] :=
  550.     Block[{tmp, kk, series = SafeSeries[f, {z, a, n}]},
  551.       Off[General::poly];
  552.       tmp = CoefficientList[series,z];
  553.       If[SameQ[Head[tmp],CoefficientList],
  554.              tmp = CoefficientList[Normal[series],z];
  555.              If[SameQ[Head[tmp],CoefficientList],
  556.                      tmp = {Normal[series]}
  557.       ]];
  558.       On[General::poly];
  559.        If[Length[tmp] > n + 1, tmp = Take[tmp, n + 1]];
  560.        Return[Join[tmp, Table[0, 
  561.            {kk, n + 1 - Length[tmp]} ]]]]
  562.            
  563. CList[l_List, x_] := CList[#, x]& /@ l
  564.  
  565. CList[0, x_] := {0}
  566.  
  567. CList[a_, x_] :=
  568.   Block[{tmp},
  569.     Off[General::poly];
  570.     tmp = CoefficientList[a, x];
  571.     If[SameQ[Head[tmp],CoefficientList],
  572.       tmp = CoefficientList[Normal[a],x];
  573.       If[SameQ[Head[tmp],CoefficientList],
  574.     tmp = {Normal[a]}
  575.     ]];
  576.     On[General::poly];
  577.     tmp
  578.  ] /; a =!= 0
  579.         
  580.         
  581. (* INTEGRABLE and PREINTEGRATE *)        
  582.                 
  583. Integrable[a_Plus, b___] := Min[Integrable[#, b]& /@ (List @@ a)]
  584.  
  585. Integrable[a_, functions_, z_] := Infinity /; FreeL[a, functions]
  586.  
  587. Integrable[a_ b_, functions_, z_] := 
  588.     Integrable[b, functions, z] /; FreeQ[a, z]
  589.     
  590. Integrable[a_, functions_, z_] := 
  591.     FirstArg[a] /; FirstHead[a] === Derivative
  592.                            
  593. Integrable[a___] := 0
  594.  
  595.  
  596. PreIntegrate[a_, functions_, {z_, 0}] := a
  597.  
  598. PreIntegrate[a_Plus, b___] := PreIntegrate[#, b]& /@ a
  599.  
  600. PreIntegrate[a_, functions_, {z_, n_}] := 
  601.     Nest[Integrate[#, z]&, a, n] /; FreeL[a, functions]
  602.  
  603. PreIntegrate[a_ b_, functions_, {z_, n_}] := 
  604.     a PreIntegrate[b, functions, {z, n}] /; FreeQ[a, z]
  605.     
  606. PreIntegrate[a_, functions_, {z_, n_}] := 
  607.      MapAt[(# - n)&, a, 
  608.          Append[First[Position[a, Derivative[_]]], 1]] /;
  609.                              FirstHead[a] === Derivative
  610.                            
  611.  
  612.  
  613. (* FIRSTHEAD, FIRSTARG and HEADCOUNT *)
  614.  
  615. FirstHead[a_] := Nest[Head, a, HeadCount[a]]
  616.  
  617.  
  618.  
  619. FirstArg[a_] := Which[Length[a] == 0,
  620.                           Null,
  621.                       True,
  622.                           First[Nest[Head, a, HeadCount[a] - 1]] ]
  623.                           
  624.  
  625. HeadCount[a_] := Which[Length[a] == 0,
  626.                            0,
  627.                        True,
  628.                            HeadCount[Head[a]] + 1 ]
  629.                            
  630.                            
  631.                            
  632. (* POWER SUM *)
  633.  
  634.  
  635.     (* lists and equations *)
  636.     
  637. PowerSum[expr_List, rest__] :=
  638.     PowerSum[#, rest]& /@ expr
  639.  
  640. PowerSum[expr_Equal, rest__] :=
  641.     PowerSum[#, rest]& /@ expr
  642.  
  643.  
  644.     (* revert to inner syntax *)
  645.      
  646. PowerSum[expr_, {z_Symbol, n_, n0_:0}] :=
  647.     Block[{e = expr //. {Sum[a_, {k_, m_}] :> Sum[a, {k, 1, m}],
  648.                          Literal[Even[a_]] :> SymbolicMod[a, 2] == 0,
  649.                          Literal[Odd[a_]]  :> SymbolicMod[a, 2] == 1,
  650.                          Mod -> SymbolicMod} },
  651.     PowerSum[e, z, n, n0] /. 
  652.             Derivative[kk_][ff_] :> D[ff, {z, kk}]
  653.     ]
  654.  
  655.  
  656.     (* non-zero starting point *)
  657.    
  658. PowerSum[expr_, z_, n_, n0_] := 
  659.     z^n0 PowerSum[Expand[expr /. n -> n + n0] /. 
  660.             a_Symbol?(Context[#] =!= "System`" && # =!= K &)[m_] :> 
  661.                 a[Expand[m]], z, n
  662.      ]  /; NumberQ[N[n0]]
  663.                             
  664.                            
  665.  
  666.     (* the geometric series *)
  667.    
  668. PowerSum[1, z_, n_] :=
  669.     1/(1 - z)
  670.  
  671.     (* linearity *)
  672.    
  673. PowerSum[c_ expr_., z_, n_] := 
  674.     c PowerSum[expr, z, n]     /; FreeL[c, {n, Blank}]
  675.                                                           
  676. PowerSum[c_.(expr1_ + expr2_), rest__] := 
  677.         PowerSum[c expr1, rest] + PowerSum[c expr2, rest] 
  678.     
  679.                     
  680. PowerSum[Sum[a_, {k_, lo_, hi_}] b_., z_, n_] :=
  681.     Sum[PowerSum[Expand[a b], z, n],
  682.     {k, lo, hi}]     /;      FreeQ[{lo, hi}, n] && FreeQ[b, k]
  683.     
  684.                        
  685.     (* powers *)
  686.     
  687. PowerSum[c_^((a_ + b_)d_) e_., rest__] := 
  688.     PowerSum[c^(a d + b d) e, rest]
  689.  
  690. PowerSum[c_^(a_ + b_) e_., z_, n_] := 
  691.     c^a PowerSum[c^b e, z, n]    /; FreeL[c, {n, Blank}] && FreeL[a, {n, Blank}]
  692.                         
  693. PowerSum[(a_ + b_)^c_Integer?Positive d_., rest__] := 
  694.     PowerSum[Expand[(a + b)^c d], rest]
  695.     
  696.                                           
  697.                                           
  698.     (* exponentials *)
  699.     
  700. PowerSum[e_. c_^((a_. n_ + b_.)d_.), z_, n_] :=
  701.     (c^(b d) PowerSum[e, z, n] /. z -> c^(a d) z)    /;
  702.             ( (FreeL[#, {n, Blank}]& /@ And[a, b, c, d]) &&
  703.               FreeL[e, {z, Blank}]
  704.             )
  705.                     
  706.                     
  707.     (* binomials *)
  708.     
  709. PowerSum[Binomial[a_, n_ + k_.], z_, n_] :=
  710.     z^(-k) *
  711.     ((1 + z)^a - Sum[Binomial[a, m] z^m, {m, 0, k - 1}])    /;
  712.             FreeQ[a, n] && IntegerQ[k]
  713.  
  714. PowerSum[Binomial[a_, b_], z_, n_] :=
  715.     Block[{e = Expand[n - b], r = Expand[a - 2 b], m},
  716.         z^e (   ((1 - Sqrt[1 - 4 z])/(2 z))^r / Sqrt[1 - 4 z] -
  717.                 Sum[Binomial[2 m + r, m] z^m, {m, 0, -(e + 1)}]
  718.         )
  719.         /;    IntegerQ[e] && IntegerQ[r]
  720.     ]
  721.  
  722.  
  723.     (* trinomials *)
  724.     
  725. PowerSum[Sum[a_. b_Binomial c_Binomial, {k_, 0, n_}],
  726.     z_, n_] :=
  727.     Block[{aa = a^(1/k)}, 
  728.         1/Sqrt[1 - 2 (2 aa + 1) z + z^2]
  729.             /;     (FreeL[aa, {k, n}] &&
  730.                 MakeTrinomial[b c] == Sort[Multinomial[k, k, n - k]])
  731.     ] 
  732.  
  733. PowerSum[Sum[a_. b_Binomial c_Binomial, {k_, 0, n_}],
  734.     z_, n_] :=
  735.     Block[{aa = a^(1/(n - k))}, 
  736.         1/Sqrt[1 - 2 (2 aa + 1) z + z^2]
  737.         /;
  738.                FreeL[aa, {k, n}] && 
  739.                MakeTrinomial[b c] == Sort[Multinomial[k, n - k, n - k]]
  740.     ] 
  741.  
  742. PowerSum[Sum[a_. b_Binomial c_Binomial, {k_, 0, n_}],
  743.     z_, n_] :=
  744.     Block[{aa = a^(1/k)}, 
  745.         1/Sqrt[1 - 2 z + (1 - 4 aa) z^2]
  746.         /;
  747.         FreeL[aa, {k, n}] && 
  748.         MakeTrinomial[b c] == Sort[Multinomial[k, k, n - 2 k]]
  749.     ]
  750.  
  751. PowerSum[Sum[a_. b_Binomial c_Binomial, {k_, 0, n_}],
  752.     z_, n_] :=
  753.     Block[{aa = a^(1/(n - k))}, 
  754.         1/Sqrt[1 - 2 z + (1 - 4 aa) z^2]
  755.         /;
  756.         FreeL[aa, {k, n}] && 
  757.         MakeTrinomial[b c] == Sort[Multinomial[n - k, n - k, 2 k - n]] 
  758.     ] 
  759.  
  760.  
  761.  
  762.     (* PowerSum is the inverse of SeriesTerm *)
  763.     
  764. PowerSum[Literal[SeriesTerm[f_, z_, n_ + m_.]], 
  765.     z_, n_] := 
  766.     Block[{k},
  767.         z^(-m) (f - Sum[SeriesTerm[f, {z, 0, k}] z^k, 
  768.                     {k, 0, m - 1}])
  769.         /;
  770.                 PoleMultiplicity[f, {z, 0}] == 0 && IntegerQ[m] 
  771.     ]
  772.  
  773.  
  774.     (* trigonometric inhomogeneity *)
  775.     
  776. PowerSum[Cos[n_], rest__] := 
  777.     PowerSum[(E^(I n) + E^(-I n))/2, rest]
  778.         
  779.                            
  780. PowerSum[Sin[n_], rest__] := 
  781.     PowerSum[(E^(I n) - E^(-I n))/2/I, rest]
  782.  
  783.     
  784.     (* factorial inhomogeneity *)
  785.     
  786. PowerSum[expr_./(n_ + a_.)!, z_, n_] :=
  787.     Block[{k},
  788.        (ExponentialPowerSum[Expand[expr /. n -> n - a], {z, n}] - 
  789.                Sum[(expr /. n -> k - a) z^k / k!, {k, 0, a - 1}]) / z^a
  790.     ] /;  IntegerQ[a]
  791.                                                   
  792.  
  793.     (* generating functions of sequences
  794.     
  795.       (here we assume that a[n] is defined for n >= 0) *)
  796.  
  797.  
  798.  
  799. PowerSum[n_^m_. a_[n_ + d_.], z_, n_] :=
  800.     Block[{j, k},
  801.          Sum[StirlingS2[m, j] z^j * 
  802.             Derivative[j][(PowerSum[a[k], z, k] - 
  803.                            Sum[a[k] z^k, {k, 0, d - 1}])/z^d ],
  804.            {j, 0, m}] // Factor
  805.     ] /; 
  806.      IntegerQ[d] && TTQ[d >= 0] && IntegerQ[m] && TTQ[m >= 0]    
  807.             
  808. PowerSum[a_[k_. n_ + d_.], z_, n_] := 
  809.     Block[{j, m},
  810.         (z^(-d/k)/k Sum[Omega[k]^(-d j) (Gf[a][Omega[k]^j z^(1/k)] -
  811.             Sum[Omega[k]^(m j) a[m] z^(m/k), {m, 0, d-1}]),
  812.             {j, 0, k - 1} ] ) // Factor
  813.     ] /;
  814.          IntegerQ[d] && TTQ[d >= 0] && IntegerQ[k] && TTQ[k > 0]
  815.  
  816.  
  817.  
  818.     (* rational inhomogeneity *)
  819.     
  820. PowerSum[n_^m_., z_, n_] := 
  821.     Block[{j},
  822.         Factor[Sum[(StirlingS2[m, j] j! z^j)/(1 - z)^(j + 1),
  823.                    {j, 0, m}] ]
  824.     ] /; IntegerQ[m] && TTQ[m >= 0]
  825.                    
  826.  
  827. PowerSum[(n_ + a_)^m_Integer?Negative, z_, n_] :=
  828.     Block[{j},
  829.         z^(-a) PolyLog[-m, z] - Sum[z^(j-a) j^m, {j, 1, a-1}]
  830.     ] /;    IntegerQ[a] && TTQ[a > 0]
  831.         
  832.                            
  833.     (* convolutions *)
  834.     
  835. PowerSum[Sum[expr_, {k_, lo_, hi_}], z_, n_] :=
  836.     Block[{nfree, nfull, s, aa = Expand[lo], bb = Expand[hi - n]},
  837.          ({nfree, nfull} = If[Head[expr] === Times, 
  838.                        {Select[expr,  FreeQ[#, n]&],
  839.                     Select[expr, !FreeQ[#, n]&]},
  840.                     If[FreeQ[expr, n],
  841.                     {expr, 1},
  842.                     {1, expr}
  843.                 ]
  844.             ];
  845.         Sum[z^k nfree PowerSum[Expand[nfull /. n -> s + k], z, s, -k],
  846.          {k, aa, bb - 1}] +
  847.         PowerSum[nfree PowerSum[Expand[nfull /. n -> s + k], z, s, -bb],
  848.                   z, k, Max[aa, bb] ]
  849.         )     /;     
  850.                         FreeQ[aa, n] && FreeQ[bb, n] 
  851.     ]
  852.  
  853.  
  854.     (* multisection of series *)
  855.                                     
  856. PowerSum[If[SymbolicMod[n_ + d_., m_] == k_, a_, b_] c_., z_, n_] :=
  857.     Block[{l = Mod[k - d, m], jj,
  858.            fps = PowerSum[Expand[(a - b)c], z, n]}, 
  859.         PowerSum[Expand[b c], z, n] + 
  860.         1/m Sum[Omega[m]^(jj(m - l)) fps /. z -> z Omega[m]^jj,
  861.         {jj, 0, m - 1}] // Expand
  862.     ]
  863.                                     
  864.  
  865.     (* other conditionals *)
  866.                                     
  867. PowerSum[If[cond_, a_, b_] c_., z_, n_] := 
  868.     Block[{t, tt, s, k},
  869.         t = IneqSolve[cond && n >= 0, n];
  870.         tt = IneqSolve[!cond && n >= 0, n];
  871.         Plus @@ (Sum[Expand[a c] z^n, {n, First[#], Last[#]}]& /@ t) +
  872.             Plus @@ (Sum[Expand[b c] z^n, {n, First[#], Last[#]}]& /@ 
  873.             tt) /.
  874.                {Sum[s_. z^n, {n, k_, Infinity}] :> 
  875.                     PowerSum[s, {z, n, k}],
  876.                 Sum[0, {__}] -> 0}
  877.     ]
  878.  
  879.  
  880.     (* PowerSum under Series *)
  881.     
  882. PowerSum /: 
  883.     Series[PowerSum[a_, z_, n_], {z_, 0, m_}] :=
  884.         Block[{k}, 
  885.         Sum[(a /. n -> k) z^k, {k, 0, m}] + O[z]^(m+1)]
  886.     
  887.     
  888. (* EXPONENTIAL POWER SUM *)
  889.  
  890.  
  891.     (* lists and equations *)
  892.     
  893. ExponentialPowerSum[expr_List, rest__] := 
  894.     ExponentialPowerSum[#, rest]& /@ expr
  895.  
  896. ExponentialPowerSum[expr_Equal, rest__] := 
  897.     ExponentialPowerSum[#, rest]& /@ expr
  898.  
  899.  
  900.     (* revert to inner syntax *)
  901.      
  902. ExponentialPowerSum[expr_, {z_Symbol, n_, n0_:0}] :=
  903.     Block[{e = expr //. {Sum[a_, {k_, m_}] :> Sum[a, {k, 1, m}],
  904.                          Literal[Even[a_]] :> SymbolicMod[a, 2] == 0,
  905.                          Literal[Odd[a_]]  :> SymbolicMod[a, 2] == 1,
  906.                          Mod -> SymbolicMod} },
  907.         ExponentialPowerSum[e, z, n, n0] /. 
  908.             Derivative[kk_][ff_] :> D[ff, {z, kk}] ];
  909.  
  910.  
  911.  
  912.     (* non-zero starting point *)
  913.    
  914. ExponentialPowerSum[expr_, z_, n_, n0_] := 
  915.     ExponentialPowerSum[expr /. n -> n + n0, z, n] /; n0 >= 0
  916.                           
  917.  
  918.     (* the exponential series *)
  919.    
  920. ExponentialPowerSum[1, z_, n_] := E^z
  921.  
  922.  
  923.  
  924.     (* linearity *)
  925.    
  926. ExponentialPowerSum[c_ expr_., z_, n_] := 
  927.     c ExponentialPowerSum[expr, z, n] /; FreeL[c, {n, Blank}]
  928.                                                           
  929. ExponentialPowerSum[c_.(expr1_ + expr2_), rest__] := 
  930.     ExponentialPowerSum[c expr1, rest] + 
  931.     ExponentialPowerSum[c expr2, rest]
  932.  
  933. ExponentialPowerSum[Sum[a_, {k_, lo_, hi_}] b_., z_, n_] :=
  934.     Sum[ExponentialPowerSum[Expand[a b], z, n], {k, lo, hi}] /; 
  935.                                    FreeQ[{lo, hi}, n] && FreeQ[b, k]
  936.                     
  937.                     
  938.     (* powers *)
  939.     
  940. ExponentialPowerSum[c_^((a_ + b_)d_), rest__] := 
  941.     ExponentialPowerSum[c^(a d + b d), rest]
  942.  
  943. ExponentialPowerSum[c_^(a_ + b_), z_, n_] := 
  944.     c^a ExponentialPowerSum[c^b, z, n] /; 
  945.                         FreeL[c, {n, Blank}] && FreeL[a, {n, Blank}]
  946.                                           
  947.                                           
  948.     (* exponentials *)
  949.     
  950. ExponentialPowerSum[e_. c_^((a_. n_ + b_.)d_.), z_, n_] :=
  951.     (c^(b d) ExponentialPowerSum[e, z, n] /. z -> c^(a d) z) /;
  952.   (FreeL[#, {n, Blank}]& /@ And[a, b, c, d]) && FreeL[e, {z, Blank}]
  953.                     
  954.                     
  955.     (* polynomial inhomogeneity *)
  956.     
  957. ExponentialPowerSum[n_^m_., z_, n_] := 
  958.     Block[{j},
  959.         E^z Factor[Sum[StirlingS2[m, j] z^j,
  960.                    {j, 0, m}] ] ]  /; IntegerQ[m] && TTQ[m >= 0]
  961.                    
  962.                    
  963.     (* factorial inhomogeneity *)
  964.     
  965. ExponentialPowerSum[(n_ + a_.)!, z_, n_] := 
  966.     PowerSum[Pochhammer[n + 1, a], z, n] /; IntegerQ[a]
  967.                    
  968.                    
  969.     (* convolutions *)
  970.     
  971. ExponentialPowerSum[Sum[expr_ Binomial[n_ + cc_, k_], {k_, lo_, hi_}],
  972.     z_, n_] :=
  973.     Block[{aa = Expand[lo], bb = Expand[hi - n]},
  974.         Derivative[cc][ExponentialPowerSum[Sum[(expr /. n -> n - cc) *
  975.                                       Binomial[n, k],
  976.                         {k, aa, n + bb - cc} ], z, n ] ] /;
  977.       IntegerQ[cc] && TTQ[cc > 0] && FreeQ[aa, n] && FreeQ[bb, n] ]
  978.  
  979.  
  980.  
  981. ExponentialPowerSum[Sum[expr_ Binomial[n_, k_], {k_, lo_, hi_}],
  982.     z_, n_] :=
  983.     Block[{nfree, nfull, s, aa = Expand[lo], bb = Expand[hi - n]},
  984.        {nfree, nfull} = If[Head[expr] === Times, 
  985.            {Select[expr,  FreeQ[#, n]&], Select[expr, !FreeQ[#, n]&]},
  986.             If[FreeQ[expr, n], {expr, 1}, {1, expr}] ];
  987.         Sum[z^k nfree/k! ExponentialPowerSum[Expand[nfull /. n -> s + k],
  988.                  z, s, -k], {k, aa, bb - 1}] +
  989.             ExponentialPowerSum[nfree ExponentialPowerSum[Expand[nfull /.
  990.                 n -> s + k], z, s, -bb], z, k, Max[aa, bb] ] /; 
  991.                                           FreeQ[aa, n] && FreeQ[bb, n] ]
  992.  
  993.  
  994.     (* multisection of series *)
  995.                                     
  996. ExponentialPowerSum[If[SymbolicMod[n_ + d_., m_] == k_, a_, b_] c_., 
  997. z_, n_] :=
  998.     Block[{l = Mod[k - d, m], jj,
  999.            fps = ExponentialPowerSum[Expand[(a - b)c], z, n]}, 
  1000.         ExponentialPowerSum[Expand[b c], z, n] + 
  1001.         1/m Sum[Omega[m]^(jj(m - l)) fps /. z -> z Omega[m]^jj,
  1002.         {jj, 0, m - 1}] // Expand ]
  1003.                                     
  1004.  
  1005.     (* other conditionals *)
  1006.  
  1007. ExponentialPowerSum[If[cond_, a_, b_] c_., z_, n_] := 
  1008.     Block[{t, tt, s, k},
  1009.         t = IneqSolve[cond && n >= 0, n];
  1010.         tt = IneqSolve[!cond && n >= 0, n];
  1011.         Plus @@ (Sum[Expand[a c] z^n / n!, {n, First[#], Last[#]}]& /@ 
  1012.             t) +
  1013.         Plus @@ (Sum[Expand[b c] z^n / n!, {n, First[#], Last[#]}]& /@ 
  1014.             tt) /.
  1015.                {Sum[s_. z^n / n!, {n, k_, Infinity}] :> 
  1016.                     ExponentialPowerSum[s, {z, n, k}], 
  1017.                 Sum[0, {__}] -> 0} ]
  1018.  
  1019.        
  1020.        
  1021.     (* generating functions of sequences 
  1022.     
  1023.       (here we assume that a[n] is defined for n >= 0) *)
  1024.         
  1025. ExponentialPowerSum[a_[n_ + d_.], z_, n_] := 
  1026.     Derivative[d][EGf[a][z]]  /; IntegerQ[d] && TTQ[d >= 0]
  1027.                                                           
  1028. ExponentialPowerSum[n_^m_. a_[n_ + d_.], z_, n_] :=
  1029.     Block[{j, k},
  1030.         Sum[StirlingS2[m, j] z^j * 
  1031.             Derivative[d + j][EGf[a][z]], {j, 0, m}] ] /; 
  1032.      IntegerQ[d] && TTQ[d >= 0] && IntegerQ[m] && TTQ[m >= 0]            
  1033.             
  1034. ExponentialPowerSum[n_^m_. a_[n_ - 1], z_, n_] :=
  1035.     Block[{j, k},
  1036.         Sum[StirlingS2[m, j] z^j * 
  1037.             Derivative[j - 1][EGf[a][z]], {j, m}] ] /; 
  1038.                                     IntegerQ[m] && TTQ[m > 0]            
  1039.             
  1040.  
  1041.  
  1042.  
  1043.  
  1044.  
  1045. (* GENERATING FUNCTION *)
  1046.  
  1047. GeneratingFunction[list1_, list2_, n_, z_, opts___Rule] := 
  1048.     Block[{temp, recur, conds, unknowns, startValues, start, 
  1049.            rec, con, i, prs},
  1050.         prs = Parse[list1, list2, n];
  1051.         If[prs === $Failed,
  1052.             Return[$Failed],
  1053.           {recur, conds, unknowns, startValues} = prs ];
  1054.        {rec, con} = Reset[recur, conds, unknowns, startValues, 0];
  1055.        (Evaluate[start /@ unknowns]) = startValues;
  1056.         temp = FunctionSolve[rec, con, unknowns, n, PowerSum, 
  1057.                              z, opts];
  1058.         If[temp === $Failed,
  1059.             Return[$Failed],
  1060.           temp = Table[#[[i]] z^startValues[[i]], 
  1061.              {i, Length[unknowns]} ]& /@ temp;
  1062.           Return[temp/.aa_?(MemberQ[unknowns,#]&)[kk_] :> 
  1063.                            aa[Expand[kk + start[aa]]] ]]]
  1064.  
  1065.  
  1066. (* EXPONENTIAL GENERATING FUNCTION *)
  1067.  
  1068. ExponentialGeneratingFunction[list1_, list2_, n_, z_, opts___Rule] := 
  1069.     Block[{temp, recur, conds, unknowns, startValues, start,
  1070.            rec, con, prs},
  1071.         prs = Parse[list1, list2, n];
  1072.         If[prs === $Failed,
  1073.             Return[$Failed],
  1074.           {recur, conds, unknowns, startValues} = prs ];
  1075.        {rec, con} = Reset[recur, conds, unknowns, startValues, 0];
  1076.        (Evaluate[start /@ unknowns]) = startValues;
  1077.         temp = FunctionSolve[rec, con, unknowns, n, ExponentialPowerSum, 
  1078.                              z, opts];
  1079.         If[temp === $Failed,
  1080.             Return[$Failed],
  1081.             Return[temp/.aa_?(MemberQ[unknowns,#]&)[kk_] :> 
  1082.                          aa[Expand[kk + start[aa]]] ]]]
  1083.  
  1084.  
  1085. (** (Hypergeometric.m) **)
  1086.  
  1087.  
  1088. (* HSOLVE *)
  1089.  
  1090. HSolve[eq_, f_, x_Symbol, opts___Rule] :=
  1091.     Block[{hsol = HSolve1[eq, f, x, opts], ll,k},
  1092.         (ll = PList[hsol, C[_]];
  1093.         Return[hsol /. Thread[ll -> Table[C[k], {k, Length[ll]}]]])
  1094.     /; hsol =!= $Failed
  1095.     ]
  1096.  
  1097. HSolve1[eq_, f_, x_Symbol, opts___Rule] :=
  1098.  
  1099.     Block[{lhs, rhs, terms, weights, inhom, minw, theta, ll, llc, rlc,
  1100.            aList, bList, badB, nnHS, n0, subst, eqn, y, k, pr, temp,
  1101.            pm, nneg, n1, badB1},
  1102.             
  1103.         pr = PrecisionHS /.{opts} /.Options[HSolve];
  1104.         If[Head[eq] === List,
  1105.             If[Length[eq] > 1, Return[$Failed]];
  1106.             eqn = First[eq],
  1107.             eqn = eq ];
  1108.         If[Head[f] === List,
  1109.             If[Length[f] > 1, Return[$Failed]];
  1110.             y = Head[First[f]],
  1111.             y = Head[f] ];
  1112.         lhs = Expand[Numerator[Together[First[eqn] - Last[eqn]]]];
  1113.         terms = Select[MakeList[lhs, Plus], !FreeQ[#, y]&];
  1114.         weights = Union[Weight[#, y[x], x]& /@ terms];
  1115.         minw = Min[weights];
  1116.         If[!MemberQ[{{0},{0,1}}, weights - minw],
  1117.         Return[$Failed]];
  1118.     Off[StirlingS1::intnm];
  1119.         lhs = Expand[lhs/x^minw] /. x^k_.Derivative[n_][y][x] -> 
  1120.             x^(k - n) Sum[StirlingS1[n, j] theta^j y[x], {j,0,n}]
  1121.             // Expand;
  1122.     On[StirlingS1::intnm];
  1123.         inhom = Plus @@ Select[MakeList[lhs, Plus], FreeQ[#, y]&] 
  1124.             // Expand;
  1125.         pm = PoleMultiplicity[inhom, {x, 0}];
  1126.         If[!(IntegerQ[pm] && pm >= 0), Return[$Failed]];
  1127.         If[!PolynomialQ[Expand[x^pm inhom], x], Return[$Failed]];
  1128.         lhs = (lhs - inhom) / y[x] // Expand;
  1129.         rhs = Plus @@ Select[MakeList[lhs, Plus], FreeQ[#, x]&];
  1130.         lhs = Expand[(rhs - lhs)/x];
  1131.         {llc, rlc} = {LeadingCoef[lhs, theta], LeadingCoef[rhs, theta]};
  1132.         aList = If[FreeQ[lhs, theta], {}, -#[[1,2]]& /@ 
  1133.             Which[pr === Automatic,
  1134.                       {ToRules[Roots[lhs == 0, theta,
  1135.                        Cubics -> False, Quartics -> False] ]} /.
  1136.                        Literal[Roots[a_, b_, c___]] :> NRoots[a, b],
  1137.                   pr === Infinity,
  1138.                        Solve[lhs == 0, theta] /.
  1139.                        Roots -> NRoots,
  1140.                   NumberQ[pr],
  1141.                        Solve[lhs == 0, theta] /.
  1142.                        Literal[Roots[a_, b_]] :> NRoots[a, b, pr] ]];
  1143.         bList = If[FreeQ[rhs, theta], {}, (1 - #[[1,2]])& /@ 
  1144.             Which[pr === Automatic,
  1145.                       {ToRules[Roots[rhs == 0, theta,
  1146.                        Cubics -> False, Quartics -> False] ]} /.
  1147.                        Literal[Roots[a_, b_, c___]] :> NRoots[a, b],
  1148.                   pr === Infinity,
  1149.                        Solve[rhs == 0, theta] /.
  1150.                        Roots -> NRoots,
  1151.                   NumberQ[pr],
  1152.                        Solve[rhs == 0, theta] /.
  1153.                        Literal[Roots[a_, b_]] :> NRoots[a, b, pr] ]];
  1154.         badB = Select[bList, (IntegerQ[#] && # <= 0)&];
  1155.         nnHS = If[badB != {}, 1 - Min[badB], 0];        
  1156.         n0 = Max[nnHS, Exponent[inhom, x]];        
  1157.         badB1 = Select[bList, (IntegerQ[#] && # > 0)&];        
  1158.         nneg = If[badB1 != {}, Max[badB1] - 1, 0]; 
  1159.         n1 = Max[nneg, pm];   
  1160.         Clear[c];    
  1161.         c[-n1 - 1] = 0;
  1162.         subst = {ToRules[Reduce[Table[If[n + 1 == 0, 
  1163.             Coefficient[inhom, x, 0] /. x -> 1/x /. x -> 0,
  1164.             Coefficient[inhom, x^(n+1)]] + 
  1165.             rlc Times @@ (bList + n) c[n+1]
  1166.              == llc Times @@ (aList + n) c[n], {n, -n1 - 1, n0 - 1}],
  1167.             Table[c[n], {n, n0, -n1, -1}] ]]};
  1168.         If[subst === {}, Return[{}] ];
  1169.         
  1170.         If[Times @@ (Pochhammer[# + nnHS, n0 - nnHS]& /@ 
  1171.            aList) =!= 0 && llc =!= 0,
  1172.             
  1173.         (* then *)
  1174.         
  1175.         {aList, bList} = {aList, bList} + nnHS;
  1176.         If[!MemberQ[bList, 1], AppendTo[aList, 1]; 
  1177.                                AppendTo[bList, 1] ];
  1178.         bList = Drop[bList, {FirstPos[bList, 1]}];
  1179.         If[MemberQ[aList, #], aList = Drop[aList, 
  1180.                                  {FirstPos[aList, #]}];
  1181.                               bList = Drop[bList, 
  1182.                                  {FirstPos[bList, #]}]
  1183.           ]& /@ bList;
  1184.         aList = Sort[aList];
  1185.         bList = Sort[bList];
  1186.         If[NumberQ[pr],
  1187.             aList = N[aList, pr];
  1188.             bList = N[bList, pr] ];
  1189.         
  1190.         temp = (Sum[c[k] x^k, {k, -n1, n0 - 1}] +
  1191.             c[n0] (llc/rlc)^(nnHS - n0) (Times @@ 
  1192.             (Pochhammer[#, n0 - nnHS]&
  1193.             /@ bList)) / (Times @@ (Pochhammer[#, n0 - nnHS]& /@ 
  1194.             aList)) (n0 - nnHS)!
  1195.             x^nnHS (HypergeometricF[aList, bList, llc/rlc x] -
  1196.             Sum[(Times @@ (Pochhammer[#, k]& /@ aList))/
  1197.                 (Times @@ (Pochhammer[#, k]& /@ bList)) (llc/
  1198.                 rlc x)^k / k!,
  1199.                 {k, 0, n0 - nnHS - 1}] ) ) /. subst // Simplify,
  1200.         
  1201.         (* else *)
  1202.         
  1203.         {aList, bList} = {aList, bList} + n0;
  1204.         If[!MemberQ[bList, 1], AppendTo[aList, 1]; 
  1205.                                AppendTo[bList, 1] ];
  1206.         bList = Drop[bList, {FirstPos[bList, 1]}];
  1207.         If[MemberQ[aList, #], aList = Drop[aList, 
  1208.                                  {FirstPos[aList, #]}];
  1209.                               bList = Drop[bList, 
  1210.                                  {FirstPos[bList, #]}]
  1211.           ]& /@ bList;
  1212.         aList = Sort[aList];
  1213.         bList = Sort[bList];
  1214.         If[NumberQ[pr],
  1215.             aList = N[aList, pr];
  1216.             bList = N[bList, pr] ];
  1217.         
  1218.         temp = (Sum[c[k] x^k, {k, -n1, n0 - 1}] +
  1219.             c[n0] x^n0 HypergeometricF[aList, bList, llc/rlc x])
  1220.             /. subst // Expand ];
  1221.             
  1222.         temp = temp /. c -> C;
  1223.         Return[Thread[y[x] -> #]& /@ {temp}] ]
  1224.  
  1225.  
  1226.  
  1227. (* WEIGHT *)
  1228.  
  1229.  
  1230. Weight[a_. x_^k_. Derivative[n_][y_][x_], y_[x_], x_] := k - n /;
  1231.                                                  FreeL[a, {x, y}]
  1232.                                                         
  1233. Weight[a_. Derivative[n_][y_][x_], y_[x_], x_] := -n /; FreeL[a, {x, y}]
  1234.  
  1235. Weight[a_. x_^k_. y_[x_], y_[x_], x_] := k  /; FreeL[a, {x, y}]
  1236.  
  1237. Weight[a_. y_[x_], y_[x_], x_] := 0  /; FreeL[a, {x, y}]
  1238.  
  1239. Weight[a_, y_[x_], x_] := Infinity  /; FreeQ[a, y]   
  1240.  
  1241.  
  1242.  
  1243. (* HYPERGEOMETRIC F *)                                  
  1244.             
  1245.                                                  
  1246. HypergeometricF[aList_List, bList_List, c_.z_] :=
  1247.     Block[{max, kk, neg = Select[aList, (IntegerQ[#] && # <= 0)&]},
  1248.        (max = - Max[neg];
  1249.         Sum[Times @@ (Pochhammer[#, kk]& /@ aList) /
  1250.             Times @@ (Pochhammer[#, kk]& /@ bList) / kk! (c z)^kk,
  1251.                 {kk, 0, max}]) /; Length[neg] > 0 ]
  1252.  
  1253. HypergeometricF[a_List, b_List, 0] := 1
  1254.  
  1255. HypergeometricF[{}, {}, z_] := E^z
  1256.  
  1257. HypergeometricF[{a_}, {}, z_] := 1/(1 - z)^a
  1258.  
  1259. HypergeometricF[{}, {c_}, z_] := Hypergeometric0F1[c, z]
  1260.  
  1261. HypergeometricF[{a_}, {c_}, z_] := Hypergeometric1F1[a, c, z]
  1262.  
  1263. HypergeometricF[{a_, b_}, {c_}, z_] := Hypergeometric2F1[a, b, c, z]
  1264.  
  1265. HypergeometricF /: 
  1266.     Series[HypergeometricF[aList_List, bList_List, c_.z_], 
  1267.     {z_, 0, n_}] :=
  1268.         Block[{kk},
  1269.             Sum[Times @@ (Pochhammer[#, kk]& /@ aList) /
  1270.                 Times @@ (Pochhammer[#, kk]& /@ bList) / kk! (c z)^kk,
  1271.                 {kk, 0, n}] + O[z]^(n + 1) ]
  1272.                 
  1273. Unprotect[ Derivative]
  1274.  
  1275. Derivative[0, 0, n_Integer?Positive][HypergeometricF][aList_List, 
  1276.     bList_List, c_.z_] :=
  1277.         c^n Times @@ (Pochhammer[#, n]& /@ aList) /
  1278.             Times @@ (Pochhammer[#, n]& /@ bList) HypergeometricF[
  1279.             aList + n, bList + n, c z]
  1280.             
  1281. Unprotect[ Derivative]
  1282.  
  1283.  
  1284.  
  1285. (** (SeriesExpansions.m) **)
  1286.  
  1287.  
  1288. (* SERIES TERM *)
  1289.  
  1290. SeriesTerm[expr_List, rest__] :=
  1291.     SeriesTerm[#, rest]& /@ expr
  1292.  
  1293. SeriesTerm[f_, {z_, a_, n_}, opts___Rule] :=
  1294.     Block[{ff = f /. z -> z + a, nnST, useApart, n0, temp0, temp,
  1295.            precision, makeReal, useMod, specialFunctions},
  1296.         IntegerQ[nnST] ^= True;
  1297.         Info[nnST] ^= Info[n] /. Solve[n == nnST, PatternList[n, 
  1298.                 _Symbol?(Context[#] =!= "System`"&) ]][[1]];
  1299.         sumDepth = 1;
  1300.     {useApart,precision,makeReal,useMod,specialFunctions} =
  1301.         {UseApart,PrecisionST,MakeReal,UseMod,SpecialFunctions} /.
  1302.             {opts} /.Options[SeriesTerm];
  1303.         Which[useApart === Automatic,
  1304.                   ff = Apart[Factor[ff]];
  1305.           ff = If[Head[ff]===Plus,
  1306.             Map[If[PolynomialQ[Denominator[#],z],
  1307.                 ExpandDenominator[#], #]&, ff],
  1308.             If[PolynomialQ[Denominator[ff],z],
  1309.                 ExpandDenominator[ff], ff]
  1310.             ],
  1311.               useApart === All,
  1312.                   ff = Apart[Factor[ff], z],
  1313.               useApart =!= None,
  1314.                   Message[ST::badopt, UseApart, useApart ] ];
  1315.         temp = SeriesTerm[ff, z, nnST] /. {nnST -> n};
  1316.         temp = Chop[temp] /. SimplifyFloat;
  1317.         temp = temp //. SimplifyTrig;
  1318.         temp = temp /. SimplifyBinomial;
  1319.         temp = temp /. 
  1320.                 {Power[aa_, b_] :> Power[aa, Expand[b]],
  1321.                  Binomial[aa_, b_] :> Binomial[Expand[aa], Expand[b]],
  1322.                  If[aa_, b_, c_] :> If[ISolve[aa], b, c] /;
  1323.                      FreeL[aa, {Odd, Even, SymbolicMod}] };
  1324.         temp = temp /. aa_If :> Release //@ aa;
  1325.         temp = temp /. 
  1326.                 {If[Literal[Odd[m_]], If[m_ >= aa_, b_, c_] d_., e_] :>
  1327.                      If[Odd[m], ReleaseHold[Expand[b d]], e] /;
  1328.                      Entails[Info[m], m >= aa - 1] && OddQ[aa],
  1329.                  If[Literal[Even[m_]], If[m_ >= aa_, b_, c_] d_., e_] :>
  1330.                      If[Even[m], ReleaseHold[Expand[b d]], e] /;
  1331.                      Entails[Info[m], m >= aa - 1] && EvenQ[aa] };
  1332.         temp = temp //. SimplifySum;
  1333.         If[NumberQ[precision], 
  1334.             temp = N[temp, precision] /. SimplifyFloat ];
  1335.         If[MatchQ[Info[n], n >= _],
  1336.             n0 = Info[n][[2]];
  1337.             temp0 = temp /. If[n >= aa_, b_, _] :> b /; aa == n0 + 1;
  1338.             If[(temp0 /. n -> n0) == (temp /. n -> n0), temp = temp0] ];
  1339.         Return[temp] ]
  1340.         
  1341.  
  1342.     (* finite n *)
  1343.         
  1344. SeriesTerm[f_, z_Symbol, 0] := 
  1345.     SeriesTerm[z f, z, 1] /; FreeQ[f, HypergeometricF]
  1346.  
  1347. SeriesTerm[f_, z_Symbol, n_Integer] := 
  1348.     Coefficient[SafeSeries[f, {z, 0, n}], z^n]  /; 
  1349.                                 n != 0 && FreeQ[f, HypergeometricF]
  1350.  
  1351.     (* linearity *)
  1352.     
  1353. SeriesTerm[f_Plus, rest__] :=
  1354.     SeriesTerm[#, rest]& /@ f
  1355.  
  1356. SeriesTerm[c_ f_, z_Symbol, rest__] :=
  1357.     c SeriesTerm[f, z, rest] /;      FreeQ[c, z] && !FreeQ[f, z]
  1358.  
  1359.     (* SeriesTerm is inverse of Gf and of PowerSum *)
  1360.     
  1361. SeriesTerm[Gf[a_][z_], z_Symbol, n_] := a[n]
  1362.  
  1363. SeriesTerm[PowerSum[a_, z_, n_], z_Symbol, m_] :=  a /. n -> m
  1364.     
  1365.  
  1366.     (* SeriesTerm is close to the inverse of EGf and of 
  1367.        ExponentialPowerSum *)
  1368.     
  1369. SeriesTerm[EGf[a_][z_], z_Symbol, n_] := a[n] / n!
  1370.  
  1371. SeriesTerm[ExponentialPowerSum[a_, z_, n_], z_Symbol, m_] :=  (a / n!) /. n->m
  1372.  
  1373.  
  1374.  
  1375.     (* binomials *)
  1376.     
  1377. SeriesTerm[(a_.z_ + b_)^p_, z_Symbol, n_] :=
  1378.     Block[{k},
  1379.        b^p When[p < 0,
  1380.                 PowerExpand[(-a/b)^n] Binomial[n - p - 1, n],
  1381.                 PowerExpand[ (a/b)^n] Binomial[p, n]
  1382.        ]
  1383.     
  1384.         /;    (FreeQ[#, z]&) /@ (a && b && p)
  1385.     ]
  1386.         
  1387. SeriesTerm[(a_.z_^m_. + b_)^p_, z_Symbol, n_] :=
  1388.          Cancel[a^p SeriesTerm[z^(m p)(1 + b/a z^(-m))^p, z, n]]  /;
  1389.                (FreeQ[#, z]&) /@ (a && b && m && p) && TTQ[m < 0]  
  1390.  
  1391. SeriesTerm[(a_.z_^m_ + b_)^p_, z_Symbol, n_] :=
  1392.     If[SymbolicMod[n, m] == 0,
  1393.         Cancel[b^p SeriesTerm[(1 + a/b z)^p, z, n/m]],
  1394.       0]  /; (FreeQ[#, z]&) /@ (a && b && m && p) && TTQ[m >= 0] &&
  1395.                                                Head[p] =!= Integer
  1396.     
  1397. SeriesTerm[(a_.z_^m_. + b_.z_^k_.)^p_, z_Symbol, n_] :=
  1398.     Cancel[b^p SeriesTerm[z^(k p)(1 + a/b z^(m-k))^p, z, n]]  /;
  1399.                          (FreeQ[#, z]&) /@ (a && b && m && k && p)   
  1400.  
  1401.  
  1402.     (* hypergeometrics *)
  1403.     
  1404. SeriesTerm[HypergeometricF[aList_, bList_, c_.z_], z_Symbol, n_] :=
  1405.      When[n >= 0,
  1406.        (Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
  1407.           Pochhammer[#, n]]& /@ aList))/                      
  1408.        (Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
  1409.           Pochhammer[#, n]]& /@ bList))/n! c^n,
  1410.        0]  /; FreeQ[c, z]
  1411.         
  1412. SeriesTerm[Hypergeometric2F1[a_, b_, c_, d_.z_], z_Symbol, n_] :=
  1413.      When[n >= 0,
  1414.        (Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
  1415.            Pochhammer[#, n]]& /@ {a, b}))/                      
  1416.        (Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
  1417.            Pochhammer[#, n]]& /@ {c, 1})) d^n,
  1418.        0]  /; FreeQ[d, z]
  1419.  
  1420. SeriesTerm[Hypergeometric1F1[a_, b_, c_.z_], z_Symbol, n_] :=
  1421.      When[n >= 0, 
  1422.          If[IntegerQ[a],
  1423.        (a + n - 1)!/(a - 1)!,
  1424.        Pochhammer[a, n]
  1425.     ] / (Times @@ (If[IntegerQ[#],
  1426.               (# + n - 1)!/(# - 1)!,
  1427.                   Pochhammer[#, n]
  1428.                ]& /@ {b, 1})) c^n,
  1429.     0]  /; FreeQ[c, z]
  1430.  
  1431.  
  1432. SeriesTerm[Hypergeometric0F1[a_, b_.z_], z_Symbol, n_] :=
  1433.      When[n >= 0,
  1434.         1/(Times @@ (If[IntegerQ[#], (# + n - 1)!/(# - 1)!,
  1435.            Pochhammer[#, n]]& /@ {a, 1})) b^n,
  1436.     0]  /; FreeQ[b, z]
  1437.  
  1438.         
  1439.      (* exponentials *)
  1440.      
  1441. SeriesTerm[E^(k_.z_ + m_.), z_Symbol, n_] :=
  1442.     E^m k^n/n!  /;    FreeQ[k, z] && FreeQ[m, z]
  1443.                                 
  1444. SeriesTerm[a_^b_, z_Symbol, n_] :=
  1445.     SeriesTerm[E^(b Log[a]), z, n] /;    a =!= E && FreeQ[a, z]
  1446.                                     
  1447.     (* logarithms *)
  1448.     
  1449. SeriesTerm[Log[(a_ + b_.z_)^k_.], z_Symbol, n_] :=
  1450.     When[n == 0, Release[k Log[a]], Release[-k(-b/a)^n / n]]  /;
  1451.                    FreeQ[k, z] && FreeQ[a, z] && FreeQ[b, z]
  1452.                    
  1453. SeriesTerm[Log[a_Times], rest___] :=
  1454.     Plus @@ (SeriesTerm[Log[#], rest]& /@ a) 
  1455.     
  1456.                                     
  1457.     (* powers and constants *)
  1458.     
  1459. SeriesTerm[z_^k_.f_, z_Symbol, n_] :=
  1460.     SeriesTerm[f, z, n-k]  /; IntegerQ[k] && FreeQ[k, z] && !FreeQ[f, z]                                                         
  1461. SeriesTerm[z_^k_., z_Symbol, n_] :=
  1462.     When[n == k, 1, 0]  /; IntegerQ[k] && FreeQ[k, z]                                                        
  1463. SeriesTerm[f_, z_Symbol, n_] :=
  1464.     f When[n == 0, 1, 0]  /; FreeQ[f, z]
  1465.  
  1466.  
  1467.     (* rational functions *)
  1468.     
  1469. SeriesTerm[f_. g_^k_Integer?Negative, z_Symbol, n_] := 
  1470.     Block[{ff = Apart[Factor[f g^k], z]},
  1471.         If[Head[ff] === Plus,
  1472.             Return[SeriesRat[#, z, n]& /@ ff] ,
  1473.             Return[SeriesRat[ff, z, n]]
  1474.     ]
  1475.     ] /;  PolynomialQ[Expand[f], z] && PolynomialQ[Expand[g], z] &&
  1476.         !FreeQ[g, z]
  1477.  
  1478.     
  1479. SeriesRat[f_, z_Symbol, n_] :=
  1480.     SeriesTerm[f, z, n]  /; PolynomialQ[Expand[f], z]
  1481.          
  1482. SeriesRat[c_. z_^k_Integer?Negative, z_Symbol, n_] := 
  1483.     SeriesTerm[c z^k, z, n]  /; FreeQ[c, z]
  1484.          
  1485.  
  1486. SeriesRat[f_. g_^k_Integer?Negative, z_Symbol, n_] :=
  1487.   
  1488.     Block[{ff = Expand[f], gg = Expand[g], poleList, poleSet,
  1489.            mp, denom = Expand[g^(-k)], i, deg, mod, modList,
  1490.            displace, temp, l, pr = precision, inverseRoot, kk},
  1491.         
  1492.          deg = Exponent[gg, z];
  1493.          If[useMod, 
  1494.              l = CoefficientList[ff, z];
  1495.              mod = Exponent[gg, z, GCD];
  1496.              modList = Union[Mod[#, mod]& /@ (Select[Range[Length[l]],
  1497.                  (l[[#]] =!= 0)&] - 1)];
  1498.              If[Length[modList] == 1,
  1499.                  displace = First[modList],
  1500.                 {displace, mod} = {0, 1}
  1501.          ],
  1502.              {displace, mod} = {0, 1}
  1503.      ];     
  1504.          {ff, gg, denom} = {Expand[ff/z^displace], gg, denom} /.
  1505.              z -> z^(1/mod);
  1506.          deg = deg/mod;
  1507.          Which[precision === Automatic && deg > 2, 
  1508.                    pr = Indeterminate,
  1509.                precision === Automatic && deg <= 2, 
  1510.                    pr = Infinity
  1511.      ];
  1512.          If[pr == Infinity && !FreeQ[Roots[gg == 0, z], Roots],
  1513.              pr = Indeterminate];
  1514.          Which[
  1515.              pr === Indeterminate, 
  1516.              poleList = Cancel[#[[1,2]]& /@ {ToRules[
  1517.                             NRoots[gg == 0, z] ]}];
  1518.              poleSet = Union[poleList];
  1519.              mp = -k Count[poleList, #]& /@ poleSet;
  1520.              poleUse = poleSet;
  1521.              temp = SeriesTerm[Plus @@
  1522.             Table[
  1523.               SeriesDivide[
  1524.                 Horner[ff, {z, poleSet[[i]], mp[[i]]}],
  1525.                      Drop[Horner[denom, {z, poleSet[[i]], 2 mp[[i]]}],
  1526.                        mp[[i]]
  1527.                 ],
  1528.                 mp[[i]]
  1529.               ] .
  1530.                           Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
  1531.                           {i, Length[poleSet]}
  1532.             ],
  1533.             z, Expand[(n-displace)/mod]
  1534.             ],
  1535.              
  1536.              pr =!= Infinity, 
  1537.                  poleList = Cancel[#[[1,2]]& /@ {ToRules[
  1538.                             NRoots[gg == 0, z, pr] ]}];
  1539.                  poleSet = Union[poleList];
  1540.                  mp = -k Count[poleList, #]& /@ poleSet;
  1541.                  poleUse = poleSet;
  1542.                  temp = SeriesTerm[Plus @@
  1543.               Table[
  1544.                 SeriesDivide[
  1545.                   Horner[ff, {z, poleSet[[i]], mp[[i]]}],
  1546.                        Drop[Horner[denom, {z, poleSet[[i]], 2 mp[[i]]}],
  1547.                      mp[[i]]
  1548.                   ],
  1549.                   mp[[i]]
  1550.                 ] .
  1551.                          Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
  1552.                          {i, Length[poleSet]}
  1553.               ],
  1554.               z, Expand[(n-displace)/mod]
  1555.             ],
  1556.              
  1557.              pr === Infinity && deg <= 2, 
  1558.                  poleList = Cancel[#[[1,2]]& /@ {ToRules[
  1559.                             Roots[gg == 0, z] ]}];
  1560.                  poleSet = Union[poleList];
  1561.                  mp = -k Count[poleList, #]& /@ poleSet;
  1562.                  inverseRoot = Expand[(1 - gg/(gg /. z -> 0))/z];
  1563.                  poleUse = ArgPi[Apart[
  1564.                      inverseRoot /. z -> #]]^(-1)& /@ poleSet;
  1565.                  temp = SeriesTerm[Plus @@
  1566.               Table[
  1567.                 SeriesDivide[
  1568.                   Horner[ff, {z, poleSet[[i]], mp[[i]]}],
  1569.                        Drop[Horner[denom, {z, poleSet[[i]], 2 mp[[i]]}],
  1570.                          mp[[i]]
  1571.                   ],
  1572.                   mp[[i]]
  1573.                 ] .
  1574.                          Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
  1575.                          {i, Length[poleSet]}
  1576.               ],
  1577.               z, Expand[(n-displace)/mod]
  1578.             ],
  1579.              
  1580.              pr === Infinity && deg > 2, 
  1581.                  poleList = Cancel[#[[1,2]]& /@ {ToRules[
  1582.                             Roots[CoefficientList[gg, z].
  1583.                             Table[z^kk, {kk, deg, 0, -1}] == 0, z] ]}];
  1584.                  poleSet = Union[poleList];
  1585.                  mp = -k Count[poleList, #]& /@ poleSet;
  1586.                  poleUse = Complex[Together[ReP[#]//.SimplifyCubic],
  1587.                            Together[ImP[#]//.SimplifyCubic]]& /@ poleSet;
  1588.                  poleUse = poleUse /. Complex[a_, 0] -> a;
  1589.                  If[FreeQ[poleUse, Complex],
  1590.             makeReal = False
  1591.          ];
  1592.                  temp = SeriesTerm[Plus @@
  1593.               Table[
  1594.                 SeriesDivide[
  1595.                   Horner[ff /. z -> z / poleUse[[i]],
  1596.                 {z, 1, mp[[i]]}],
  1597.                               Drop[Horner[denom /. z -> z / poleUse[[i]],
  1598.                 {z, 1, 2 mp[[i]]}],
  1599.                     mp[[i]]
  1600.                   ],
  1601.                   mp[[i]]
  1602.                 ] .
  1603.                          Table[(z poleUse[[i]] - 1)^(-j),
  1604.                 {j, mp[[i]], 1, -1}],
  1605.                          {i, Length[poleSet]}
  1606.               ],
  1607.               z, Expand[(n-displace)/mod]
  1608.             ]
  1609.               ];
  1610.  
  1611.          temp = temp //. a_ c_If + b_ c_ -> (a + b) c;
  1612.     If[makeReal,
  1613.       If[pr === Infinity,
  1614.         temp = temp /. SimplifyComplexE1];
  1615.       If[pr === Infinity && deg == 2,
  1616.         temp = temp /. SimplifyComplex2];
  1617.       If[pr === Infinity && deg == 3,
  1618.         temp = temp /. SimplifyComplex3];
  1619.       If[pr === Infinity && deg == 4,
  1620.         temp = temp //. SimplifyComplex4];
  1621.       If[pr === Infinity,
  1622.         temp = temp /. SimplifyComplexE2];
  1623.       If[pr =!= Infinity,
  1624.         temp = temp //. SimplifyComplexN]
  1625.     ];
  1626.         temp = temp /. Complex[a_, b_] :> a + b I;
  1627.         Return[If[Evaluate[SymbolicMod[n, mod] == Mod[displace, mod]], 
  1628.                 Evaluate[temp], 0]]
  1629.         ] /; PolynomialQ[Expand[f], z] && PolynomialQ[Expand[g], z] &&
  1630.         !FreeQ[g, z]
  1631.          
  1632.  
  1633.     (* special functions *)
  1634.  
  1635. SeriesTerm[1 / Sqrt[a_ + b_.z_ + c_.z_^2], z_Symbol, n_] := 
  1636.     When[n >= 0, 
  1637.         Sqrt[c/a]^n / Sqrt[a] LegendreP[n, -Sgn[a] b /(2 Sqrt[a c])], 
  1638.             0]  /; specialFunctions && FreeQ[{a,b,c}, z]
  1639.  
  1640. SeriesTerm[Power[a_ + b_.z_ + c_.z_^2, Rational[p_, 2]], z_Symbol, n_] := 
  1641.     SeriesTerm[Expand[(a + b z + c z^2)^((p + 1)/2)] / 
  1642.         MySqrt[a + b z + c z^2], z, n]  /; 
  1643.                        p != -1 && specialFunctions && FreeQ[{a,b,c}, z]
  1644.                
  1645. SeriesTerm[1 / MySqrt[a_], z_Symbol, n_] := 
  1646.     SeriesTerm[1 / Sqrt[a], z, n] 
  1647.     
  1648. SeriesTerm[(a_ + b_.z_)^m_ (c_ + d_.z_)^m_, z_Symbol, n_] :=
  1649.     SeriesTerm[Expand[(a + b z) (c + d z)]^m, z, n]  /;
  1650.         specialFunctions && FreeQ[{a,b,c,d}, z] &&
  1651.         Expand[a c]^m == Expand[a^m c^m]
  1652.  
  1653. SeriesTerm[(a_ + b_.z_)^m_ (c_ + d_.z_)^m_, z_Symbol, n_] :=
  1654.    -SeriesTerm[Expand[(a + b z) (c + d z)]^m, z, n]  /;
  1655.         specialFunctions && FreeQ[{a,b,c,d}, z] &&
  1656.         Expand[a c]^m == -Expand[a^m c^m]
  1657.  
  1658.     
  1659.  
  1660.     (* products with a rational function *)
  1661.     
  1662. SeriesTerm[a_ b_, z_Symbol, n_] :=
  1663.     Block[{temp, deg, i, m, bn = SeriesTerm[b, z, m], sd = sumDepth},
  1664.          Increment[sumDepth];
  1665.          Off[General::itervar];
  1666.          If[PolynomialQ[Expand[a], z],
  1667.              deg = Exponent[a, z];
  1668.              temp = CoefficientList[a, z] . 
  1669.                  Table[bn /. m -> i, {i, n, n - deg, -1}] /. If -> When,
  1670.              temp = Sum[SeriesTerm[Apart[a, z], z, n - K[sd]] bn /. m -> K[sd], 
  1671.                        Evaluate[{K[sd], -PoleMultiplicity[b, {z, 0}],
  1672.                       n + PoleMultiplicity[a, {z, 0}]}] ]
  1673.          ];
  1674.          On[General::itervar];
  1675.          Return[temp] ] /; SeparateProduct[a b, (PolynomialQ[Expand[#], z] || 
  1676.                            PolynomialQ[Expand[#^(-1)], z])& ] == {a, b}
  1677.         
  1678.     (* general products *)
  1679.     
  1680. SeriesTerm[a_ b_, z_Symbol, n_] := 
  1681.     Block[{temp, sd = sumDepth},
  1682.          Increment[sumDepth];
  1683.          Off[General::itervar];
  1684.          temp = Sum[SeriesTerm[a, z, n - K[sd]] * 
  1685.              SeriesTerm[b, z, K[sd]],
  1686.             Evaluate[{K[sd], -PoleMultiplicity[b, {z, 0}],
  1687.                    n + PoleMultiplicity[a, {z, 0}]}] ];
  1688.          On[General::itervar];
  1689.          Return[temp] ] /; SeparateProduct[a b, (PolynomialQ[Expand[#], z] || 
  1690.                            PolynomialQ[Expand[#^(-1)], z])& ] == {1, a b}
  1691.         
  1692. SeriesTerm[a_^k_Integer?Positive, z_Symbol, n_] := 
  1693.     Block[{temp, sd = sumDepth},
  1694.          Increment[sumDepth];
  1695.          Off[General::itervar];
  1696.          temp = Sum[SeriesTerm[a, z, n - K[sd]] * 
  1697.              SeriesTerm[a^(k - 1), z, K[sd]],
  1698.             Evaluate[{K[sd], -PoleMultiplicity[a^(k - 1), {z, 0}],
  1699.                    n + PoleMultiplicity[a, {z, 0}]}] ];
  1700.          On[General::itervar];
  1701.          Return[temp] ] /; !FreeQ[a, z]
  1702.  
  1703.  
  1704.    (* trinomials *)
  1705.    
  1706. SeriesTerm[(a_ + b_.z_ + c_.z_^2)^alpha_, z_Symbol, n_] := 
  1707.     Block[{temp},
  1708.          Off[General::itervar];
  1709.          temp = Sum[Binomial[alpha, K[sumDepth]] Binomial[K[sumDepth],
  1710.              n - K[sumDepth]] a^(alpha - K[sumDepth]) b^(2 K[sumDepth] -
  1711.              n) c^(n - K[sumDepth]), Evaluate[{K[sumDepth], 0, n}]];
  1712.          On[General::itervar];
  1713.          ++sumDepth;
  1714.          Return[temp] ] /; FreeQ[{a,b,c}, z]
  1715.          
  1716.          
  1717.    (* derivatives *)
  1718.    
  1719. SeriesTerm[Derivative[k_][f_][z_], z_Symbol, n_] :=
  1720.     Pochhammer[n + 1, k] SeriesTerm[f[z], z, n + k]
  1721.       
  1722.  
  1723.  
  1724. (* HORNER and SERIES DIVIDE *)
  1725.  
  1726. Horner[p_, {z_, alpha_, m_}] :=
  1727.     Block[{qh = Expand[p], ih, kh, nh, ah},
  1728.         {nh, ah} = {Exponent[qh, z], CoefficientList[qh, z]};
  1729.         Do[ah[[ih]] = alpha ah[[ih + 1]] + ah[[ih]], {kh, 1, Min[nh, m]},
  1730.             {ih, nh, kh, -1}];
  1731.     If[m > nh + 1,
  1732.             Do[AppendTo[ah, 0], {ih, m - nh - 1}]
  1733.     ];
  1734.         Return[Take[ah, m]]
  1735.     ]
  1736.         
  1737.         
  1738. SeriesDivide[a_, b_, n_] :=
  1739.     Block[{is, cs = Range[n], ks},
  1740.         Do[cs[[is]] = (a[[is]] - Sum[cs[[ks]] b[[is - ks + 1]], 
  1741.             {ks, 1, is - 1}])/b[[1]], 
  1742.             {is, 1, n}];
  1743.         Return[cs] ]
  1744.         
  1745.         
  1746.  
  1747. (* SYMBOLIC MOD *)
  1748.  
  1749. SymbolicMod[n_, 1] = 0
  1750.  
  1751. SymbolicMod[n_?NumberQ, k_] := Mod[n, k]
  1752.  
  1753. SymbolicMod /: (SymbolicMod[n_ + a_., 2] == b_) := Even[n] /;
  1754.                EvenQ[b - a] && (SameQ[Head[n],Symbol] || MatchQ[n,K[_Integer]])
  1755.  
  1756. SymbolicMod /: (SymbolicMod[n_ + a_., 2] == b_) := Odd[n] /;
  1757.                OddQ[b - a] && (SameQ[Head[n],Symbol] || MatchQ[n,K[_Integer]])
  1758.  
  1759.      
  1760. (* PARTIAL FRACTIONS *)
  1761.  
  1762. PartialFractions[f_, z_, opts___Rule] :=
  1763.     Block[{ff = f, useApart},
  1764.         useApart = UseApart /.{opts} /.Options[PartialFractions];
  1765.         Which[useApart == Automatic,
  1766.                   ff = Apart[ff],
  1767.               useApart == All,
  1768.                   ff = Apart[ff, z],
  1769.               useApart =!= None,
  1770.                   Message[PF::badopt, UseApart, useApart ] ];
  1771.         If[Head[ff] === Plus,
  1772.             Return[PartFrac[#, z, opts]& /@ ff],
  1773.             Return[PartFrac[f, z, opts]] ]];
  1774.  
  1775.     
  1776. PartFrac[f_, z_, opts___Rule] :=
  1777.     Block[{num = Numerator[f], den = Denominator[f], poleList, poleSet,
  1778.            mp, i, j, k, deg, temp, precision},         
  1779.          precision = PrecisionST /.{opts} /. Options[PartialFractions];
  1780.         (deg = Exponent[den, z];
  1781.          Which[precision === Automatic && deg > 2, 
  1782.                    precision = Indeterminate,
  1783.                precision === Automatic && deg <= 2, 
  1784.                    precision = Infinity];
  1785.          poleList = Cancel[#[[1,2]]& /@ {ToRules[Which[
  1786.              precision == Indeterminate, NRoots[den == 0, z],
  1787.              precision =!= Infinity, NRoots[den == 0, z, precision],
  1788.              precision === Infinity, Roots[den == 0, z]] ]}];
  1789.          poleSet = Union[poleList];
  1790.          mp = Count[poleList, #]& /@ poleSet;
  1791.          If[precision === Infinity, 
  1792.              poleUse = ArgPi[#]& /@ poleSet,
  1793.              poleUse = poleSet];
  1794.          temp = Plus @@ Table[SeriesDivide[Horner[num, {z, 
  1795.              poleSet[[i]], mp[[i]]}], Drop[Horner[den, {z, 
  1796.              poleSet[[i]], 2 mp[[i]]}], mp[[i]]], mp[[i]] ] .
  1797.              Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
  1798.             {i, Length[poleSet]}]
  1799.          ) /; PolynomialQ[Expand[num], z] && PolynomialQ[Expand[den], z] ];
  1800.          
  1801.  
  1802. (* REALQ, REP, IMP, ANGLE, SGN, ABSV, and CONJUGATEQ *)
  1803.  
  1804. RealQ[x_] := True /; IntegerQ[x]
  1805. RealQ[x_Rational] = True
  1806. RealQ[x_Real] = True
  1807. RealQ[Pi] = True
  1808. RealQ[E] = True
  1809. RealQ[n_!] := True /; RealQ[n]
  1810.  
  1811. ReP[x_] := x /; RealQ[x]
  1812. ImP[x_] := 0 /; RealQ[x]
  1813.  
  1814. ReP[k_Integer x_] := k ReP[x]
  1815. ReP[k_Rational x_] := k ReP[x]
  1816. ReP[k_Real x_] := k ReP[x]
  1817.  
  1818. ImP[k_Integer x_] := k ImP[x]
  1819. ImP[k_Rational x_] := k ImP[x]
  1820. ImP[k_Real x_] := k ImP[x]
  1821.  
  1822. ReP[Complex[x_, y_]] := x
  1823. ImP[Complex[x_, y_]] := y
  1824.  
  1825. ReP[x_ + y_] := ReP[x] + ReP[y]
  1826. ImP[x_ + y_] := ImP[x] + ImP[y]
  1827.  
  1828. ReP[x_ Sqrt[y_?Positive]] :=  ReP[x] Sqrt[y]
  1829. ImP[x_ Sqrt[y_?Positive]] :=  ImP[x] Sqrt[y]
  1830. ReP[x_ Sqrt[y_?Negative]] := -ImP[x] Sqrt[-y]
  1831. ImP[x_ Sqrt[y_?Negative]] :=  ReP[x] Sqrt[-y]
  1832.  
  1833. ReP[x_ y_] := ReP[x] ReP[y] - ImP[x] ImP[y]
  1834. ImP[x_ y_] := ReP[x] ImP[y] + ImP[x] ReP[y]
  1835.  
  1836. ReP[x_?Positive ^n_] := x^n /; ImP[n] == 0
  1837. ImP[x_?Positive ^n_] := 0   /; ImP[n] == 0
  1838.  
  1839. ReP[x_?Negative ^n_Rational] := 0 /; IntegerQ[2n]
  1840. ImP[x_?Negative ^n_Rational] := (-x)^n (-1)^(n-1/2) /; 
  1841.                                           IntegerQ[2n]
  1842.  
  1843. ReP[1/x_] := ReP[x] / (ReP[x]^2 + ImP[x]^2)
  1844. ImP[1/x_] := -ImP[x] / (ReP[x]^2 + ImP[x]^2)
  1845.  
  1846. ReP[x_^Rational[p_,q_]] := 
  1847.   Module[{u = x^Quotient[p,q], v = x^(Mod[p,q]/q)},
  1848.     ReP[u] ReP[v] - ImP[u] ImP[v]
  1849.   ]    /;     AbsV[p] >= AbsV[q]
  1850.  
  1851. ImP[x_^Rational[p_,q_]] :=
  1852.   Module[{u = x^Quotient[p,q], v = x^(Mod[p,q]/q)},
  1853.     ReP[u] ImP[v] + ImP[u] ReP[v]
  1854.   ]    /;     AbsV[p] >= AbsV[q]
  1855.     
  1856.  
  1857.  
  1858. ReP[x_^q_Rational] := AbsV[x]^q Cos[q Angle[x]] /; 
  1859.                                       AbsV[q] <= 1
  1860. ImP[x_^q_Rational] := AbsV[x]^q Sin[q Angle[x]] /; 
  1861.                                       AbsV[q] <= 1
  1862.  
  1863. ReP[E^x_] := Cos[ImP[x]] Exp[ReP[x]]
  1864. ImP[E^x_] := Sin[ImP[x]] Exp[ReP[x]]
  1865.  
  1866. ReP[x_^2] := ReP[x]^2 - ImP[x]^2
  1867. ImP[x_^2] := 2 ReP[x] ImP[x]
  1868.  
  1869. ReP[x_^3] := ReP[x]^3 - 3 ReP[x] ImP[x]^2
  1870. ImP[x_^3] := 3 ReP[x]^2 ImP[x] - ImP[x]^3
  1871.  
  1872. ReP[x_^n_Integer] := Block[{a, b},
  1873.     a = Round[n/2];
  1874.     b = n-a;
  1875.     ReP[x^a] ReP[x^b] - ImP[x^a] ImP[x^b]
  1876.     ] /; n != -1
  1877.  
  1878. ImP[x_^n_Integer] := Block[{a, b},
  1879.     a = Round[n/2];
  1880.     b = n-a;
  1881.     ReP[x^a] ImP[x^b] + ImP[x^a] ReP[x^b]
  1882.     ] /; n != -1
  1883.  
  1884.  
  1885.  
  1886. Angle[z_] := 
  1887.     Block[{angle =
  1888.          Block[{x = ReP[z], y = ImP[z]},
  1889.                  Which[
  1890.            TrueQ[Expand[x] == 0], Pi/2 Sgn[y],
  1891.                    TrueQ[Expand[y] == 0], Pi/2 (1 - Sgn[x]),
  1892.                    True, ArcTan[y/x] - Pi/2 (Sgn[x]-1) Sgn[y]
  1893.              ]
  1894.          ]},
  1895.       angle /; FreeQ[angle,Sgn]
  1896.     ] /;  FreeQ[{ReP[z],ImP[z]},ReP] && FreeQ[{ReP[z],ImP[z]},ImP]
  1897.                   
  1898. Sgn[0] = 0
  1899. Sgn[a_?Positive] = 1
  1900. Sgn[a_?Negative] = -1
  1901. Sgn[a_ b_] := Sgn[a] Sgn[b]
  1902. Sgn[a_^b_] := Sgn[a]^b
  1903. Sgn[Sqrt[a_]] = 1
  1904. Sgn[a_ + b_] := Sgn[Chop[N[a + b]]] /; !SameQ[a + b, Chop[N[a + b]]]
  1905.  
  1906.     
  1907. AbsV[0] = 0
  1908. AbsV[x_?Positive] := x
  1909. AbsV[x_?Negative] := -x
  1910. AbsV[x_ y_] :=
  1911.    Block[{a = AbsV[x], b = AbsV[y]},
  1912.     a b /; FreeQ[{a,b},AbsV]
  1913.    ]
  1914. AbsV[x_^n_] :=
  1915.    Block[{a = AbsV[x]},
  1916.     a^n /; FreeQ[a,AbsV]
  1917.    ] /; Head[n] =!= Complex
  1918. AbsV[x_] :=
  1919.    Block[{absv = Block[{a = ReP[x], b = ImP[x]},
  1920.             Which[
  1921.               TrueQ[Expand[b] == 0], Expand[a Sgn[a]],
  1922.               TrueQ[Expand[a] == 0], Expand[b Sgn[b]],
  1923.               True, Sqrt[Expand[a^2 + b^2]]
  1924.                         ]
  1925.                  ]},
  1926.          absv  /; FreeQ[absv, Sgn]
  1927.   ] /; FreeQ[{ReP[x],ImP[x]}, ReP] && FreeQ[{ReP[x],ImP[x]}, ImP]
  1928.  
  1929. ConjugateQ[a_, b_] :=
  1930.     Block[{l = PatternList[{a, b}, 
  1931.                _Symbol?(Context[#] =!= "System`"&) ]},
  1932.         Which[
  1933.             l === {},
  1934.                 Chop[ N[ReP[a] - ReP[b]] ] == 0 &&
  1935.                 Chop[ N[PolynomialMod[ImP[a] + ImP[b],2 Pi]] ] == 0,
  1936.             l =!= {},
  1937.                 Expand[ReP[a] - ReP[b]] == 0 &&
  1938.                 Expand[PolynomialMod[ImP[a] + ImP[b],2 Pi]] == 0
  1939.     ]
  1940.     ]
  1941.             
  1942.         
  1943. (* SIMPLIFICATION RULES *)
  1944.  
  1945. SimplifyTrig = {
  1946.     Sin[x_. + r_Rational Pi] :> Cos[x] /; EvenQ[r-1/2],
  1947.     Sin[x_. + n_Integer?OddQ Pi] :> -Sin[x],
  1948.     Sin[x_. + r_Rational Pi] :> -Cos[x] /; EvenQ[r-3/2],
  1949.     Sin[x_. + n_Integer?EvenQ Pi] :> Sin[x],
  1950.     Cos[x_. + r_Rational Pi] :> -Sin[x] /; EvenQ[r-1/2],
  1951.     Cos[x_. + n_Integer?OddQ Pi] :> -Cos[x],
  1952.     Cos[x_. + r_Rational Pi] :> Sin[x] /; EvenQ[r-3/2],
  1953.     Cos[x_. + n_Integer?EvenQ Pi] :> Cos[x],
  1954.     Sin[a_?Negative b_.] :> -Sin[-a b],
  1955.     Cos[a_?Negative b_.] :>  Cos[-a b],
  1956.     Sin[a_Plus] :> -Sin[Expand[-a]] /; MatchQ[a[[1]], _?Negative _.],
  1957.     Cos[a_Plus] :>  Cos[Expand[-a]] /; MatchQ[a[[1]], _?Negative _.],
  1958.     ArcTan[a_?Negative b_.] :> -ArcTan[-a b],
  1959.     ArcTan[a_Plus] :> -ArcTan[Expand[-a]] /; 
  1960.                                        MatchQ[a[[1]], _?Negative _.],
  1961.     ArcTan[Sqrt[3]] -> Pi/3,
  1962.     ArcTan[Sqrt[3]/3] -> Pi/6,
  1963.     ArcTan[1/Sqrt[3]] -> Pi/6}
  1964.     
  1965.     
  1966. SimplifySqrt = {a_^Rational[b_,2] :> a^((b-1)/2) Sqrt[a]}
  1967.  
  1968. CanonicRadicals = {
  1969.    (a_ + b_)^Rational[k_,n_] :> a (a + b)^(Mod[k, n] / n) + 
  1970.         b (a + b)^(Mod[k, n] / n)   /; Quotient[k, n]==1,
  1971.     a_^Rational[k_,n_] :> a^(Mod[k, n]/n) Expand[a^Quotient[k, n]]}
  1972.  
  1973.  
  1974. SimplifyCubic = Join[SimplifyTrig, SimplifySqrt]
  1975.     
  1976.  
  1977. SimplifyBinomial = {
  1978.     Binomial[1/2, k_ + 1] :> (-4)^(-k) Binomial[2 k, k]/(2(k + 1)) /;
  1979.         Entails[Info[k], k != -1],
  1980.     Binomial[1/2, k_] -> (-4)^(-k) Binomial[2 k, k]/(1 - 2 k),
  1981.     Binomial[-1/2, k_] -> (-4)^(-k) Binomial[2 k, k],
  1982.     Binomial[k_ - 1/2, k_] -> 4^(-k) Binomial[2 k, k]}
  1983.  
  1984. SimplifySum = {
  1985.     Sum[0, {__}] -> 0,
  1986.     Sum[1, {k_, a_, b_}] :> When[a <= b, b - a + 1, 0],
  1987.     Sum[a_. b_, {k_, c__}] :> b Sum[a, {k, c}] /; 
  1988.         FreeQ[b, k],
  1989.     Sum[If[k_ == m_, a_, b_] c_., {k_, m_, n_}] :>
  1990.        ((a c) /. k -> m) When[n >= m, 1, 0] + 
  1991.         Sum[b c, {k, m + 1, n}],
  1992.     Sum[(If[-k_ + n_ >= 0, a_, b_] c_. + d_.) e_., {k_, m_, n_}] :>
  1993.          Sum[Evaluate[Expand[a c e + d e]], {k, m, n}],
  1994.     Sum[(If[k_ >= m_, a_, b_] c_. + d_.) e_., {k_, m_, n_}] :>
  1995.          Sum[Evaluate[Expand[a c e + d e]], {k, m, n}] }
  1996.                   
  1997. SimplifyComplexE1 = {
  1998.     a_.I^c_ + d_.(-I)^c_ :> ArgPi[a] E^Expand[I Pi c/2] +
  1999.         ArgPi[d] E^Expand[-I Pi c/2] /;
  2000.       (ConjugateQ[a, d] && FreeQ[{ArgPi[a],ArgPi[d]},ArgPi]),
  2001.     a_.E^b_ + d_.E^e_ :> 2 a Cos[ImP[b]] /; 
  2002.        (ImP[a] == a - d == ReP[b] == PolynomialMod[b + e, 2 Pi I] == 0),     
  2003.     a_.E^b_ + d_.E^e_ :> -2 ImP[a] Sin[ImP[b]] /;  
  2004.        (ReP[a] == a + d == ReP[b] == PolynomialMod[b + e, 2 Pi I] == 0),     
  2005.     a_.E^b_ + d_.E^e_ :> ArgPi[a] E^b + ArgPi[d] E^e /;
  2006.        FreeQ[{ArgPi[a],ArgPi[d]},ArgPi]
  2007.             }
  2008.  
  2009.  
  2010. SimplifyComplex2 = {
  2011.     a_.b_^c_ + d_.e_^c_ :> 
  2012.         2 AbsV[a]AbsV[b]^c SimplifyCos[Cos[Angle[a] + c Angle[b]]] /; 
  2013.         (ConjugateQ[a, d] && ConjugateQ[b, e] && ImP[c] == 0 &&
  2014.      FreeQ[{AbsV[a],AbsV[b]},AbsV] && FreeQ[{Angle[a],Angle[b]},Angle])
  2015.                }
  2016.         
  2017. SimplifyComplex3 = {
  2018.     a_.b_Complex^c_ + d_.e_Complex^c_ :> 
  2019.         2 AbsV[a]AbsV[b]^c Cos[Angle[a] + c Angle[b]] /;
  2020.         (ConjugateQ[a, d] && ConjugateQ[b, e] && ImP[c] == 0 &&
  2021.      FreeQ[{AbsV[a],AbsV[b]},AbsV] && FreeQ[{Angle[a],Angle[b]},Angle])}
  2022.         
  2023. SimplifyComplex4 = {
  2024.     a_.b_Complex^c_ + d_.e_Complex^c_ :> 
  2025.         2 AbsV[a]AbsV[b]^c SimplifyCos[Cos[Angle[a] + c Angle[b]]] /; 
  2026.         (ConjugateQ[a, d] && ConjugateQ[b, e] && ImP[c] == 0 &&
  2027.      FreeQ[{AbsV[a],AbsV[b]},AbsV] && FreeQ[{Angle[a],Angle[b]},Angle])}
  2028.  
  2029. SimplifyComplexE2 = {
  2030.     a_.E^b_ + d_.E^e_ :>
  2031.         Block[{c = ImP[b]},
  2032.             2 ReP[a] Cos[c] - 2 ImP[a] Sin[c]
  2033.         ] /; (ConjugateQ[a, d] && ConjugateQ[b, e] && ReP[b] == 0 &&
  2034.           FreeQ[{ReP[a],ImP[a]},ReP] && FreeQ[{ReP[a],ImP[a]},ImP])
  2035.             }
  2036.  
  2037.         
  2038. SimplifyComplexN = {
  2039.     a_.b_Complex^c_ + d_.e_Complex^c_ :> 
  2040.         2 Abs[a]Abs[b]^c Cos[Arg[a] + c Arg[b]] /; 
  2041.         (ConjugateQ[a, d] && ConjugateQ[b, e] && ImP[c] == 0)}
  2042.         
  2043. SimplifyCos[x_] :=
  2044.     Block[{x1 = (x /. Cos[a_] :> Cos[Expand[a]]) //. SimplifyTrig,
  2045.            x2 =  x //. SimplifyTrig },
  2046.         If[LeafCount[x1] <= LeafCount[x2], x1, x2 ]]
  2047.         
  2048. SimplifyFloat = {1. -> 1, -1. -> -1, 0. -> 0}
  2049.  
  2050.  
  2051. (* FACTORIAL SIMPLIFY *)
  2052.  
  2053. FactorialSimplify[expr_] := PostSimplify[DeFact[expr]]
  2054.  
  2055.  
  2056. DeFact[expr_] :=
  2057.  
  2058.     Block[{kin, j, sub, diff, t, tt, arg = {}, ex = expr},
  2059.     
  2060.         While[!FreeQ[ex, Factorial],
  2061.             t = First[Position[ex, Factorial[_]]];
  2062.             tt = First[Part @@ Prepend[t, ex]];
  2063.             kin = Scan[If[IntegerQ[Expand[# - tt]],
  2064.                 Return[Fact[#]] ]&, arg];
  2065.             If[kin === Null,
  2066.                 AppendTo[arg, tt]; sub = Fact[tt],
  2067.                 diff = Expand[First[kin] - tt];
  2068.                 If[diff > 0,
  2069.                     sub = kin / Product[tt + j, {j, diff}],
  2070.                     sub = kin Product[tt - j + 1, {j, -diff}] ] ];
  2071.             ex = MapAt[sub&, ex, {t}] ];
  2072.         Factor[ex] /. Fact -> Factorial ]
  2073.         
  2074.         
  2075. PostSimplify[expr_] :=
  2076.     (* DEBUG; added "_" to n! below... blank was missing I think! *)
  2077.     Block[{ex = expr /. n_! :> Expand[n]!, i},
  2078.  
  2079.         ex //.
  2080.            {n_!^k_. m_!^l_.:> Product[m + i, {i, n - m}]^k /;
  2081.                                      Expand[n - m] > 0  &&  k + l == 0,
  2082.                                                  
  2083.             n_!^k_. m_^l_. :> Expand[n + 1]!^k Expand[n + 1]^(l - k) /;
  2084.                         Expand[m - n - 1] == 0  &&  Sign[k] == Sign[l],
  2085.                                     
  2086.             n_!^k_. m_^l_. :> (-1)^k Expand[n + 1]!^k m^(l - k) /; 
  2087.                          Expand[m + n + 1] == 0  && Sign[k] == Sign[l],
  2088.                          
  2089.             n_!^k_. n_^l_. :> Expand[n - 1]!^k  /;  Expand[k + l] == 0,
  2090.             
  2091.             n_!^k_. m_^l_. :> (-1)^k Expand[n - 1]!^k /;
  2092.                                Expand[m + n] == Expand[k + l] == 0 } ]
  2093.                                
  2094.  
  2095. (** (Utilities.m) **)
  2096.  
  2097.  
  2098. (* ARGPI *)
  2099.  
  2100. ArgPi[z_] :=
  2101.     Block[{fraction = Rationalize[N[Arg[z]/Pi]]},
  2102.         If[Head[fraction] === Rational, 
  2103.             Return[AbsV[z] E^(fraction Pi I)],
  2104.             Return[z]
  2105.     ]
  2106.     ]
  2107.                              
  2108.  
  2109. (* ENTAILS *)
  2110.  
  2111. Entails[a_, a_] = True
  2112.  
  2113. Entails[a___ && b_ && c___, b_] = True
  2114.  
  2115. Entails[False, a_] = True
  2116.             
  2117. Entails[a_, b_] :=
  2118.     Block[{l = Union[PatternList[{a, b}, 
  2119.                      _Symbol?(Context[#] =!= "System`"&) ]]},
  2120.     IneqSolve[Implies[a, b], First[l]] === 
  2121.                  {Interval[-Infinity, Infinity]}
  2122.       /; Length[l] == 1
  2123.     ]
  2124.              
  2125.  
  2126. (* FIRST POS and FREE L *)
  2127.  
  2128. FirstPos[l_, pattern_] := 
  2129.     SafeFirst[Select[Range[Length[l]], MatchQ[l[[#]], pattern]& ]]
  2130.  
  2131. FreeL[expr_, l_] := And @@ (FreeQ[expr, #]& /@ l)
  2132.  
  2133.  
  2134. (* ISOLVE, INEQSOLVE, LEQ *)
  2135.  
  2136. IneqSolve[x_ + a_ > b_, x_] := IneqSolve[x > b - a, x]
  2137.  
  2138. IneqSolve[x_ + a_ >= b_, x_] := IneqSolve[x >= b - a, x]
  2139.  
  2140. IneqSolve[x_ > n_Integer, x_] := {Interval[n+1, Infinity]}
  2141.  
  2142. IneqSolve[x_ >= n_Integer, x_] := {Interval[n, Infinity]}
  2143.  
  2144. IneqSolve[a_, x_] := a /; Length[Union[PatternList[a, 
  2145.                      _Symbol?(Context[#] =!= "System`"&) ]]] > 1
  2146.  
  2147. IneqSolve[a_ && b_, x_] := 
  2148.     Block[{aa = IneqSolve[a, x], 
  2149.            bb = IneqSolve[b, x],
  2150.            merge, accum, begin, nn},
  2151.         If[FreeQ[aa, IneqSolve] && FreeQ[bb, IneqSolve],
  2152.         aa = Append[aa, {}];
  2153.         bb = Append[bb, {}];
  2154.         merge = Sort[Join[
  2155.             Thread[{Join @@ (aa /. Interval -> List),
  2156.                     Join @@ (aa /. Interval[_,_] -> {-1,1}) }],
  2157.             Thread[{Join @@ (bb /. Interval -> List),
  2158.                     Join @@ (bb /. Interval[_,_] -> {-1,1}) }] ],
  2159.                     LEQ];
  2160.         accum = Thread[{First /@ merge, 
  2161.                         FoldList[Plus, First[Last /@ merge],
  2162.                 Drop[Last /@ merge, 1]]}];
  2163.         begin = Join @@ Append[Position[accum, {_,-2}], {}];
  2164.         Return[(Interval @@ # & /@ Thread[{begin, begin + 1}]) /. 
  2165.                              nn_Integer :> First[accum[[nn]]]],
  2166.         Return[aa && bb] ]]
  2167.         
  2168. IneqSolve[!a_, x_] := 
  2169.     Block[{aa = IneqSolve[a, x], merge, c, d},
  2170.         If[FreeQ[aa, IneqSolve],
  2171.             aa = Append[aa, {}];
  2172.             merge = Join @@ (aa /. Interval[c_, d_] -> {c-1, d+1});
  2173.             If[merge === {},
  2174.                 merge = {-Infinity, Infinity},
  2175.                 If[First[merge] === -Infinity, 
  2176.                     merge = Rest[merge],
  2177.                     PrependTo[merge, -Infinity] ];
  2178.                 If[Last[merge] === Infinity, 
  2179.                     merge = Drop[merge,-1],
  2180.                     AppendTo[merge, Infinity] ] ];
  2181.             Return[Interval @@ # & /@ Partition[merge, 2]],
  2182.           Return[!aa]] ]
  2183.         
  2184. IneqSolve[a_ || b_, x_] := 
  2185.     Block[{aa = IneqSolve[!(!a && !b), x]},
  2186.         If[FreeQ[aa, IneqSolve],
  2187.             Return[aa],
  2188.           Return[IneqSolve[a, x] || IneqSolve[b, x]]] ]
  2189.  
  2190. IneqSolve[Inequality[a_, b_, c_, d_, e__], x_] := 
  2191.     IneqSolve[And @@ (Inequality @@ # & /@ 
  2192.         Partition[{a,b,c,d,e},3,2]), x]
  2193.  
  2194.  
  2195. IneqSolve[h_[a_, b_, c__], x_] := 
  2196.     IneqSolve[And @@ (h @@ # & /@ 
  2197.                 Partition[{a,b,c},2,1]), x] /;
  2198.     MemberQ[{Less, LessEqual, Greater, GreaterEqual, Equal, 
  2199.              Unequal}, h]
  2200.             
  2201. IneqSolve[ineq: h_[a_, b_], x_] :=
  2202.     Block[{lhs = Together[a - b],
  2203.            num, den, zeroList, poleList, points, merge, c, d},
  2204.            
  2205.        {num, den} = {Numerator[lhs], Denominator[lhs]};
  2206.         Off[Roots::neq];
  2207.         zeroList = If[FreeQ[num, x], {},
  2208.             #[[1,2]]& /@ {ToRules[NRoots[num == 0, x]]} ];
  2209.         poleList = If[FreeQ[den, x], {},
  2210.             #[[1,2]]& /@ {ToRules[NRoots[den == 0, x]]} ];
  2211.         On[Roots::neq];
  2212.         zeroList = Select[zeroList, (Im[#] == 0)&] // Union;
  2213.         poleList = Select[poleList, (Im[#] == 0)&] // Union;
  2214.         points = Union[zeroList, poleList];
  2215.         merge = Partition[Append[Prepend[points,-Infinity],Infinity], 
  2216.             2, 1];
  2217.         merge = Select[merge, (ineq /. x :> 
  2218.                   Which[# === {-Infinity, Infinity}, 0,
  2219.                         #[[1]] === -Infinity, #[[2]] - 1,
  2220.                         #[[2]] ===  Infinity, #[[1]] + 1,
  2221.                         True, (Plus @@ #)/2 ] )&];
  2222.         merge = (Interval @@ # & /@ merge) /. 
  2223.             Interval[c_, d_] :> Interval[Ceiling[c], Floor[d]];
  2224.         merge = merge /. {Ceiling[-Infinity] -> -Infinity,
  2225.                           Floor[Infinity] -> Infinity};
  2226.         merge = Select[merge, #[[1]] <= #[[2]]&];
  2227.         zeroList = Sort[Join[Floor /@ zeroList, Ceiling /@ zeroList]];
  2228.         zeroList = Select[zeroList, ((den /. x->#) != 0)&];
  2229.         zeroList = Select[zeroList, (ineq /. x->#)&];
  2230.         merge = Sort[Join[merge, Thread[Interval[zeroList, zeroList]]],
  2231.             LEQ];
  2232.         merge = merge /. 
  2233.             Interval[c_?((NumberQ[#] && ((den /. x:>#) == 0 ||
  2234.                 !ineq/.x->#))&), d_] :> Interval[c+1, d];
  2235.         merge = merge /. 
  2236.             Interval[c_, d_?((NumberQ[#] && ((den /. x:>#) == 0 ||
  2237.                 !ineq/.x->#))&)] :> Interval[c, d-1];
  2238.         merge = Select[merge, #[[1]] <= #[[2]]&];
  2239.         args = List @@ (Join @@ Append[merge, Interval[]]);
  2240.         If[args === {}, {},
  2241.            {first, last} = {First[args], Last[args]};
  2242.             Interval @@ # & /@ Partition[Append[Prepend[Join @@ 
  2243.                 Append[Select[Partition[Drop[Rest[args], -1], 2], 
  2244.                 (#[[2]] - #[[1]] > 1)&], {}], first], last], 2] ] ]  /; 
  2245.         MemberQ[
  2246.            {Less, LessEqual, Greater, GreaterEqual, Equal, Unequal}, h]
  2247.                 
  2248.  
  2249. ISolve[expr_] :=
  2250.     Block[{l = UserSymbols[expr], n, temp, p, nn},
  2251.         If[Length[l] == 1,
  2252.             n = First[l];
  2253.             temp = IneqSolve[expr /. n -> nn, nn] /. nn -> n;
  2254.             temp = temp /. List[p___Interval] :> Or[p];
  2255.             temp = temp /.
  2256.               {Interval[-Infinity, Infinity] -> True,
  2257.                Interval[-Infinity, b_] -> n <= b,
  2258.                Interval[a_, Infinity] -> n >= a,
  2259.                Interval[a_, a_] -> n == a,
  2260.                Interval[a_, b_] -> a <= n <= b};
  2261.             temp = temp /. IneqSolve[a_, n_] -> a;
  2262.             Return[temp],
  2263.           Return[expr]] ]
  2264.         
  2265. LEQ[h_[a_,b_], h_[c_,d_]] := (a < c) || (a === c && b <= d)
  2266.  
  2267.  
  2268. (* INFO *)
  2269.  
  2270. Info[_?NumberQ] = True
  2271. Info[_Symbol] = True
  2272. Info[_[a___]] := And @@ (Info /@ {a})
  2273.  
  2274.  
  2275. (* LEADING COEF *)
  2276.  
  2277. LeadingCoef[poly_, x_] := If[poly === 0, 0,
  2278.     Coefficient[poly, x, Exponent[poly, x]] ]
  2279.     
  2280.  
  2281. (* MAKE LIST *)
  2282.  
  2283. MakeList[expr_, head_:List] := 
  2284.     If[Head[expr] === head, List @@ expr, List @ expr]
  2285.     
  2286.     
  2287. (* MAKE TRINOMIAL and LIST COMPLEMENT *)
  2288.  
  2289. MakeTrinomial[a: Binomial[b_, c_] Binomial[d_, e_]] :=
  2290.     Block[{l1 = {b, d}, 
  2291.            l2 = {c, b - c, e, d - e}, l},
  2292.         l = Intersection[l1, l2];
  2293.         If[Length[l] != 1, Return[a]];
  2294.         l1 = ListComplement[l1, l];
  2295.         l2 = ListComplement[l2, l];
  2296.         If[{Plus @@ l2} != l1, Return[a]];
  2297.         Return[Sort[Multinomial @@ l2]] ]
  2298.                     
  2299. MakeTrinomial[a_] := a
  2300.  
  2301.  
  2302. ListComplement[l1_List, l2_List] :=
  2303.     Block[{l = l1},
  2304.         If[MemberQ[l, #], 
  2305.             l = Drop[l, {FirstPos[l, #]}] ]& /@ l2;
  2306.          Return[l] ]
  2307.  
  2308.  
  2309.  
  2310. (* PARSING ROUTINE *)
  2311.  
  2312. Parse[list1_, list2_, n_] :=
  2313.     Block[{l1 = MakeList[ReleaseHold @@ (Hold[list1] /. Condition -> Pair)],
  2314.            unknowns = MakeList[list2],
  2315.            eqns, conds, recur, extra, range, startValues, aa, kk,
  2316.            lo, hi, check},
  2317.          
  2318.         unknowns = Head /@ unknowns;
  2319.         
  2320.         check = If[Head[#] === Pair, #[[1]], #]& /@ l1;
  2321.         If[!MatchQ[check, {__Equal}], 
  2322.             Message[Parse::eqn, check];
  2323.             Return[$Failed] ];
  2324.             
  2325.         eqns =  Select[l1, !FreeQ[#, n]&];
  2326.         conds = Select[l1,  FreeQ[#, n]&];
  2327.         
  2328.             (* make  n >= 0  the default range of validity *)
  2329.         
  2330.         eqns = If[Head[#] =!= Pair, Pair[#, n >= 0], #]& /@ eqns; 
  2331.         
  2332.             (* solve inequalities, and separate recurrences from 
  2333.                initial conditions *)
  2334.                
  2335.         eqns = Pair[#[[1]], IneqSolve[#[[2]], n]]& /@ eqns;
  2336.         recur = Select[eqns, !FreeQ[#, Infinity]&];
  2337.         extra = Union[Select[eqns,  FreeQ[#, Infinity]&],
  2338.                       Pair[#[[1]], Drop[#[[2]], -1]]& /@ recur];
  2339.         recur = Pair[#[[1]], Last[#[[2]]][[1]]]& /@ recur;
  2340.         extra = Select[extra, #[[2]] =!= {} &];
  2341.         extra = Union @@ (Thread /@ extra);
  2342.         conds = Union[conds, Union @@ (Table[#[[1]], {n, #[[2,1]],
  2343.             #[[2,2]]}]& /@ extra) ];
  2344.             
  2345.             (* shift n so that all recurrences will be valid
  2346.                for  n >= 0  *)
  2347.  
  2348.         recur = (#[[1]] /. (n -> n + #[[2]]) )& /@ recur;
  2349.         
  2350.             (* determine starting values of n for all sequences *)
  2351.                
  2352.        {recur, conds} = {recur, conds} //. 
  2353.             Sum[a_, {k_, m_}] :> Sum[a, {k, 1, m}];
  2354.        {recur, conds} = {recur, conds} /.
  2355.             {Literal[Even[a_]] :> SymbolicMod[a, 2] == 0,
  2356.              Literal[Odd[a_]]  :> SymbolicMod[a, 2] == 1,
  2357.              Mod -> SymbolicMod};
  2358.         range = {conds, recur /. (Sum[aa_, {kk_, lo_, hi_}] :>
  2359.             {aa /. kk -> lo, aa /. kk -> hi} )};
  2360.         startValues = Min[First /@ (PatternList[range, #[_]] /. 
  2361.             n -> 0)]& /@ unknowns;
  2362.             
  2363.         Return[{recur, conds, unknowns, startValues}] ]
  2364.             
  2365.  
  2366. (* PATTERN LIST *)
  2367.  
  2368. PatternList[expr_?(FreeQ[#, Rule]&), pattern_] := 
  2369.     Apply[Part, Prepend[#, expr]]& /@ Position[expr, pattern]
  2370.  
  2371. PList[expr_, pattern_] := 
  2372.     Block[{xpr = expr /. ((a_ -> b_) -> {a,b})},
  2373.         Apply[Part, Prepend[#, xpr]]& /@ Position[xpr, pattern]]
  2374.         
  2375.  
  2376. (* POLE MULTIPLICITY *)
  2377.  
  2378. PoleMultiplicity[f_, {z_, a_}] :=
  2379.     Block[{pom = SafeSeries[f /. z -> z + a, {z, 0, 0}]},
  2380.         If[!FreeL[pom, {$Failed, Series}], Return[Infinity],
  2381.            pom = Normal[pom];
  2382.            If[pom == 0, Return[0]];
  2383.            Return[Max[0, Exponent[Expand[pom /. z -> 1/z], z]]] ]]
  2384.          
  2385.  
  2386. (* RESET *)
  2387.  
  2388. Reset[recur_, conds_, unknowns_, startValues_, s0_] :=
  2389.     Block[{start},
  2390.        (Evaluate[start /@ unknowns]) = startValues;
  2391.         Return[{recur, conds} /. 
  2392.             aa_?(MemberQ[unknowns,#]&)[kk_] :> 
  2393.             aa[kk + s0 - start[aa]] ]]
  2394.             
  2395.             
  2396. (* SAFE FIRST, SAFE SERIES *)
  2397.  
  2398. SafeFirst[l_] :=
  2399.     If[Length[l] == 0, Null, First[l]]
  2400.     
  2401. SafeSeries[f_, {z_, a_, n_Integer}] :=
  2402.     Block[{poles = z /. Solve[Denominator[Together[f]] == 0, z,
  2403.                    InverseFunctions -> True],
  2404.       ord, tem},
  2405.  
  2406.         ord = Max[n, Count[poles, a]];
  2407.         tem = Check[Series[f, {z, a, ord}], $Failed, Series::esss];
  2408.         If[tem =!= $Failed, tem = Series[tem, {z, a, n}]];
  2409.         Return[tem]
  2410.     ]
  2411.                 
  2412.  
  2413. (* SEPARATE PRODUCT *)
  2414.  
  2415. SeparateProduct[p_Times, predicate_] :=
  2416.    {Select[p, predicate], Select[p, !predicate[#]&]}
  2417.    
  2418.  
  2419. (* USER SYMBOLS *)
  2420.  
  2421. UserSymbols[expr_] :=
  2422.     Block[{h = If[Length[expr] == 0, expr, Head[expr]]},
  2423.          Which[!FreeQ[h, 
  2424.                 _Symbol?(Context[#] =!= "System`"&)],
  2425.                   {expr},
  2426.                Length[expr] == 0,
  2427.                   {},
  2428.                True,
  2429.                   Union @@ (UserSymbols /@ (List @@ expr))]]
  2430.                   
  2431.         
  2432.  
  2433. (* TTQ, WHEN *)
  2434.  
  2435. TTQ[True] := True
  2436. TTQ[a_] := Entails[Info[a], a]
  2437.  
  2438. When[cond_, a_, b_] :=
  2439.     Which[TrueQ[Entails[Info[cond],  cond]], a,
  2440.           TrueQ[Entails[Info[cond], !cond]], b,
  2441.           True, If[cond, a, b]
  2442.      ]
  2443.  
  2444.  
  2445. (* EVEN, ODD *)
  2446.  
  2447. Even[n_]:=EvenQ[n] /; !(SameQ[Head[n],Symbol] || MatchQ[n,K[_Integer]])  
  2448. Odd[n_]:=OddQ[n] /; !(SameQ[Head[n],Symbol] || MatchQ[n,K[_Integer]])
  2449.           
  2450.  
  2451.  
  2452. End[]
  2453.  
  2454. (* Protect[RSolve, PowerSum, ExponentialPowerSum, GeneratingFunction,
  2455. ExponentialGeneratingFunction, Gf, EGf, HSolve, HypergeometricF, K,
  2456. SeriesTerm, SymbolicMod, RealQ, ReP, ImP, ConjugateQ, SimplifySum, 
  2457. SimplifyTrig, SimplifyComplexE1, SimplifyComplexE2, SimplifyComplex2,
  2458. SimplifyComplex3, SimplifyComplex4, SimplifyComplexN, FactorialSimplify,
  2459. PartialFractions, ArgPi, Entails, FirstPos, FreeL, Info, ISolve, 
  2460. LeadingCoef, ListComplement, MakeList, MakeTrinomial, PatternList,
  2461. PoleMultiplicity, Reset, SafeFirst, SafeSeries, TTQ, UserSymbols, When,
  2462. Methods, MethodGF, MethodEGF, PrecisionHS, PrecisionST, UseApart, UseMod,
  2463. MakeReal, SpecialFunctions] *)
  2464.  
  2465. EndPackage[]
  2466.