(*=========================================================================*)
(*                                                                         *)
(*                             R  C  O  M  P                               *)
(*                                                                         *)
(* Authors: Istvan Nemes and Marko Petkovsek                               *)
(*                                                                         *)
(* Domain: doubly infinite sequences satisfying a linear recurrence        *)
(*         with constant coefficients                                      *)
(*                                                                         *)
(* CHANGE LOG:                                                             *)
(*                                                                         *)
(*   Aug. 4, 2000: upgraded to Mma 4.0 (MP)                                *)
(*   Jan. 9, 1995: Homogenize removed; PartialSum corrected;               *)
(*                 Tabulate expanded (MP)                                  *) 
(*   May 22, 1994: CyclicToGfn, CyclicToRec added (IN)                     *)
(*   Oct.13, 1993: errors noticed at Cornell corrected (MP)                *)
(*                                                                         *)
(*                                                                         *)
(*=========================================================================*)


BeginPackage["RComp`"]


(* Usage messages *)

Convolve::usage = "Convolve[s1, s2] computes the convolution of sequences
s1 and s2. Here s1 and s2 must have head Rec or Gfn. Convolve[a1, a2, n]
does the same for explicitly given sequences a1 and a2."

Define::usage = "Define[symb, s] defines symb[n] to be the n-th term of
the sequence s which can be given either by recurrence or by generating
function."

ExponentialPolynomialQ::usage = "ExponentialPolynomialQ[a, n] determines
if a is an exponential polynomial of indeterminate n."

FindRec::usage = "FindRec[list] returns the C-finite sequence of minimal
order whose initial segment agrees with the given list of values, provided
that this sequence is unique."

Gfn::usage = "Gfn[f, z] represents the sequence of Taylor coefficients of
f as a function of z. It is assumed that f is a rational function of z."

Hadamard::usage = "Hadamard[s1, s2] returns the Hadamard product of sequences 
s1 and s2. Here s1 and s2 must have head Rec or Gfn. Hadamard[f, g, z] returns 
a rational function h such that the coefficients of the Laurent series 
representing h at z=0 are products of the corresponding coefficients for 
rational functions f and g."

Hankel::usage = "Hankel[{a0, a1, .., ak}] produces the Hankel matrix for a
segment of a sequence."

Interlace::usage = "Interlace[s1, s2, .., sk] computes the sequence
a[n] = am[Floor[n/k]] where m = Mod[n,k]+1, and am[n] is the sequence 
represented by sm (the si's are either recurrences or generating functions, 
but not mixed). Interlace[a1, a2, .., ak, n] does the same for explicitly 
given sequences a1[n], a2[n], .., ak[n]."

ListOfFunctions::usage = "ListOfFunctions[] returns the list of all functions
defined in the package."

Minimize::usage = "Minimize is an alias for SimplifySeq."

CyclicToGfn::usage = "CyclicToGfn[{{c11,c12,..,c1k}, {c21,c22,..,c2k}, .., 
{cm1,cm2,..,cmk}}, {a1,a2,..ak}, z] computes the generating function 
(with indeterminate z) of the sequence a[n] defined by the recurrence 
a[n] = cj1 a[n-1] + cj2 a[n-2] + .. + cjk a[n-k] where j is congruent to n-k+1 
(mod m), and by the initial conditions a[0] = a1, a[1] = a2, .., a[k-1] = ak."

CyclicToRec::usage = "CyclicToRec[{{c11,c12,..,c1k}, {c21,c22,..,c2k}, ..,
{cm1,cm2,..,cmk}}, {a1,a2,..ak}] computes the recurrence representation of the
sequence a[n] defined by the recurrence a[n] = cj1 a[n-1] + cj2 a[n-2] + .. 
+ cjk a[n-k] where j is congruent to n-k+1 (mod m), and by the initial 
conditions a[0] = a1, a[1] = a2, .., a[k-1] = ak."

Multisect::usage = "Multisect[s, m, k] computes the recurrence representing 
the sequence a[m n + k] where a[n] is the sequence represented by s. It is 
assumed that s has head Rec or Gfn."

PartialSum::usage = "PartialSum[s] (or PartialSum[s, n] for the explicit
representation) returns the sequence of partial sums of sequence s."

RationalQ::usage = "RationalQ[f, z] returns True if f is a rational function
of z."

Rec::usage = "Rec[{c1,..,cd}, {a1,..,ad}] represents the sequence a[n] defined
by the recurrence a[n] = c1 a[n-1] + ... + cd a[n-d] valid for all integer n, 
and by the initial conditions a[0] = a1, a[1] = a2, ..., a[d-1] = ad. 
Rec[{c1,..,cd}, {a1,..,ad}, s] represents the sequence a[n] defined likewise 
by the inhomogeneous recurrence a[n] = c1 a[n-1] + ... + cd a[n-d] + b[n] 
where b[n] is the sequence represented by s. It is assumed that s has head 
Rec or Gfn or is of the form {b[n], n}. Rec with 3 arguments is automatically 
converted to Rec with 2 arguments (i.e., homogenized)."

Shift::usage = "Shift[s, k] returns the sequence s shifted by k. Here s must
have head Rec or Gfn. Shift[a, n, k] does the same for the explicitly given
sequence a indexed by n. In both forms, if k is omitted it is assumed to be 1."

SimplifySeq::usage = "SimplifySeq[s] simplifies sequence s given by recurrence
or by generating function. SimplifySeq[a, n] tries to simplify the 
sequence given explicitly by a and indexed by n."

StrikeOut::usage = "StrikeOut[s, m, k] (or StrikeOut[s, n, m, k] for the
explicit representation) returns the sequence obtained from s by striking out
the terms whose indices are congruent to k mod m. If k is not given it is
assumed to be 0."

Tabulate::usage = "Tabulate[s, it] returns the table of terms of sequence s
given by recurrence, generating function, or explicitly, as specified by 
iterator it."

Term::usage = "Term[s, k] returns the k-th term of sequence s defined
by recurrence or by generating function. Here k must be an integer. 
Term[a, n, k] does the same for explicitly given sequence a indexed by n."

ToExpPoly::usage = "ToExpPoly[s, n] converts sequence s given by recurrence,
generating function, or explicitly, to an exponential polynomial in n."

ToGfn::usage = "ToGfn[a, n, z] converts explicit sequence a indexed by n
to the corresponding (minimal) generating function with indeterminate z. 
ToGfn[s, z] does the same for sequences given by recurrence or by generating
function."

ToRec::usage = "ToRec[a, n] converts explicit sequence a indexed by n
to the corresponding (minimal) recurrence. ToRec[s] does the same for sequences
given by generating function or by recurrence."

Verify::usage = "Verify[a == b, n] returns True if a[n] == b[n] for all
integer n. Verify[a == b] does the same when both a and b are given by
recurrences, or when they are both given by generating functions."

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

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

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

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

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

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

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

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




Begin["`Private`"]

(** ATTRIBUTES **)

Attributes[GeneratingFunction] = {HoldFirst}
Attributes[Parse] = {HoldFirst}


(** OPTIONS **)

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




(** MESSAGES **)


ST::badopt = "`` -> `` is not a valid option."
PF::badopt = "`` -> `` is not a valid option."
Parse::eqn = "`` is not an equation or a system of equations."




(** REDEFINITIONS OF SYSTEM FUNCTIONS **)

Unprotect[Binomial, Power]
    Binomial[n_, n_] := When[n >= 0, 1, 0]
    0^0 = 1
Protect[Binomial, Power];                                           
                               

Unprotect[IntegerQ];
   IntegerQ[a_ b_] := IntegerQ[a] && IntegerQ[b];
   IntegerQ[a_ + b_] := IntegerQ[a] && IntegerQ[b];
   IntegerQ[a_^n_] := IntegerQ[a] && IntegerQ[n] && TTQ[n >= 0];
Protect[IntegerQ];


(** NEW FUNCTIONS **)

CyclicToGfn[reclist_, initlist_, z_] :=
  Module[{m = Length[reclist], in = Length[initlist], d = 
          Max[Length /@ reclist], vars, a, n, laglist, j, k, args, eqns, 
          initargs, gf},
    vars = a[#][n] & /@ Range[0, m - 1];
    laglist = Join @@ Table[a[j][n - k], {k, 0, Ceiling[d/m]}, {j, m-1, 0, -1}];
    args = MapIndexed[Take[Drop[laglist, m - #2[[1]] + 1], Length[#1]] &, 
      reclist];
    eqns = 
      Thread[vars == Inner[Dot, myhead @@ reclist, myhead @@ args, List]];
    eqns = MapThread[Condition, {eqns, n >= Max[Cases[#, a[_][x_] -> n-x,
	  Infinity]] & /@ eqns}];
    initargs = Take[Join @@ Table[a[j][k], {k, 0, Ceiling[in/m] - 1}, 
	  {j, 0, m - 1}], in];
    eqns = Join[eqns, Thread[initargs == initlist]];
    gf = {eqns, vars, n, z};
    gf = GeneratingFunction @@ gf;
	If[gf === {}, Return[$Failed]];
    gf = gf[[1]];
    Return[Gfn[Factor[Sum[z^(j - 1) (gf[[j]] /. z -> z^m), {j, m}]], z]]
  ];

CyclicToRec[ rcl_, init_]:= ToRec[ CyclicToGfn[ rcl, init, z] ];


Convolve[Gfn[f1_, z1_], Gfn[f2_, z2_]] := 
   Gfn[Factor[f1 (f2 /. z2 -> z1)], z1];

Convolve[r1_Rec, r2_Rec] := 
   Module[{z},
      ToRec[Gfn[First[ToGfn[r1, z]] First[ToGfn[r2, z]], z]]
   ];

Convolve[ex1_, ex2_, n_] := 
   Module[{z},
      ToExpPoly[Convolve[ToGfn[ex1, n, z], ToGfn[ex2, n, z]], n]
   ];


Define[a_Symbol, rec:Rec[{},{}]] :=
   Module[{z},
      Clear[a];
      DefineRec[a, rec];
      DefineGfn[a, ToGfn[rec, z]]
   ];

Define[a_Symbol, rec_Rec] :=
   Module[{z},
      If[Last[First[rec]] === 0, Return[$Failed]];
      Clear[a];
      DefineRec[a, rec];
      DefineGfn[a, ToGfn[rec, z]]
   ];

Define[a_Symbol, f_Gfn] :=
   (Clear[a];
    DefineRec[a, ToRec[f]];
    DefineGfn[a, f]
   );

DefineGfn[a_Symbol, Gfn[f_, x_]] :=
   Module[{z}, Gf[a][z_] = f /. x -> z;];

DefineRec[a_Symbol, Rec[coef_, init_]] :=
   Module[{z, d = Length[coef]},
      Evaluate[Table[a[k], {k,0,d-1}]] = init;
      a[n_Integer?Negative] := a[n] = (a[n+d] - 
	 Drop[coef, -1] . Table[a[k], {k, n+d-1, n+1, -1}]) / Last[coef];
      a[n_Integer] := a[n] = coef . Table[a[k], {k, n-1, n-d, -1}]
   ]; 

ExponentialPolynomialQ[a_, n_Symbol] :=
   Module[{z},
      RationalQ[First[ToGfn[a, n, z]], z]
   ];

FindRec[l_List] := 
   Module[{m = Hankel[Drop[l, -1]], c, d, r},
      d = Length[m];
      c := LinearSolve[m, Take[l, {d+1, 2d}]] /. _LinearSolve -> $Failed;
      If[c === $Failed, Return[$Failed]];
      If[Together[Append[c, -1] . Take[l, -(d+1)]] =!= 0, Return[$Failed]];
      r = Minimize[Rec[Reverse[c], Take[l, d]]];
      If[Head[r] =!= Rec, Return[$Failed]];
      If[(Length[First[r]] =!= 0) && (Together[Last[First[r]]] === 0), 
	 Return[$Failed], Return[r]] ];


Gfn /: Gfn[f1_, z1_] + Gfn[f2_, z2_] := 
	   Gfn[Factor[f1 + (f2 /. z2 -> z1)], z1];

Gfn /: Gfn[f1_, z1_] Gfn[f2_, z2_] := 
	   Gfn[Hadamard[f1, (f2 /. z2 -> z1), z1], z1];

Gfn /: (d_ /; FreeQ[d, Gfn]) Gfn[f_, z_Symbol]  := 
	   Gfn[Factor[d f], z];

Gfn /: Gfn[a__]^n_Integer?Positive := Gfn[a]^(n-1) Gfn[a];


Hadamard[f_, g_, z_Symbol] :=
   Module[{f1 = Factor[f], g1 = Factor[g], p, q, r, d, c, a, y, dq},
     ({p, q} = Denominator /@ {f1, g1};
      dq = Exponent[q, z];
      r = Resultant[p, z^dq (q /. z -> y/z), z];
      c = CoefficientList[r, y];
      c = - Rest[c] / First[c];
      d = Length[c];
      a = Times @@ (Pad[CoefficientList[Series[#, {z,0,d-1}], z], 
	  d]& /@ {f,g});
      Return[First[ToGfn[Rec[c, a], z]]]) /; 
				      And @@ (RationalQ[#, z]& /@ {f,g}) ];

Hadamard[Gfn[f1_, z1_], Gfn[f2_, z2_]] :=
    Gfn[Hadamard[f1, (f2 /. z2 -> z1), z1], z1];

Hadamard[r1_Rec, r2_Rec] :=
    Module[{z},
	ToRec[Hadamard[ToGfn[r1, z], ToGfn[r2, z]]]
    ];

Hankel[a_List] := 
    Module[{n = Ceiling[Length[a]/2]}, 
	Take[Partition[a, n, 1], n]
    ];

Rec[c_List, a_List, inh : Rec[cr_List, ar_List]] :=
       Module[{in, ch, d = Length[c], lh, rhs, f, n, x},
	  in = Shift[inh, d];          
	  ch = -Drop[Reverse[CoefficientList[
	     (Prepend[-cr, 1] . Table[x^i, {i, Length[cr], 0, -1}])
	     (Prepend[-c, 1] . Table[x^i, {i, d, 0, -1}]), 
	      x]], 1];
	  lh = Length[ch];
	  rhs = Tabulate[in, {n, 0, lh - d - 1}];
	  f = (Append[#1, Take[#1, -d] . Reverse[c] + #2]) &;
	  Return[Minimize[Rec[ch, Fold[f, a, rhs]]]] ];

Rec[c_List, a_List, inh_Gfn] := Rec[c, a, ToRec[inh]];

Rec[c_List, a_List, {exp_, n_}] := Rec[c, a, ToRec[exp, n]];

Interlace[a__Gfn] := 
    Module[{alist = {a}, k, z},
	k = Length[alist];
	zlist = Last /@ alist;
	z = First[zlist];
	flist = (First /@ alist) /. Thread[zlist -> z];
	Gfn[Factor[Table[z^j, {j,0,k-1}] . (flist /. z -> z^k)], z] ];


Interlace[a__Rec] := 
    ToRec[Interlace @@ (ToGfn[#, z]& /@ {a})];

Interlace[a__, n_Symbol] := 
    ToExpPoly[Interlace @@ (ToGfn[#, n, z]& /@ {a}), n];


ListOfFunctions[] = Sort[{RationalQ, Hadamard, ToGfn, 
   ToExpPoly, ToRec, Convolve, Multisect, PartialSum, Hankel, Verify, Define, 
   Rec, Interlace, ListOfFunctions, Tabulate, Term, FindRec, Shift,
   SimplifySeq, Minimize, StrikeOut, ExponentialPolynomialQ, CyclicToGfn,
   CyclicToRec}];

Minimize[r_Rec] :=
   Module[{z},
      Return[ToRec[ToGfn[r, z]]]
   ];
     
Minimize[Gfn[f_, z_]] :=
   Gfn[Factor[f], z];
     
Minimize[s__] := SimplifySeq[s];

Multisect[r_Rec, m_, k_] := 
   Module[{a, n},
      Define[a, r];
      Return[ToRec[a[m n + k], n]]
   ];

Multisect[g_Gfn, m_, k_] :=
    ToGfn[Multisect[ToRec[g], m, k], Last[g]];

Pad[l_List, d_Integer?(# >= 0 &)] :=

   (* Pad l to length d with zeros from the right *)

   Module[{k = Length[l]},
      If[d <= k, l,
	  Join[l, Table[0, {d-k}]] ]
   ];

PartialSum[r_Rec] := Convolve[r, Rec[{1},{1}]];

PartialSum[Gfn[f_, z_]] := Gfn[Factor[f/(1-z)], z];

PartialSum[expr_, n_] := ToExpPoly[PartialSum[ToRec[expr, n]], n];


RationalQ[f_, z_Symbol] := 
      And @@ (PolynomialQ[#, z]& /@ {Numerator[f], Denominator[f]});


(* the zero sequence is normalized to Rec[{},{}] *)

Rec /: Rec[_, {(0)...}] := Rec[{},{}];

Rec /: Rec[_, {}] := Rec[{},{}]; (* This should be subsumed by the rule
				    above but is not - a bug in ... ? *)

(* automatic homogeneization *)

Rec[c_List, a_List, inh : Rec[cr_List, ar_List]] :=
       Module[{in, ch, d = Length[c], lh, rhs, f, n, x},
	  in = Shift[inh, d];          
	  ch = -Drop[Reverse[CoefficientList[
	     (Prepend[-cr, 1] . Table[x^i, {i, Length[cr], 0, -1}])
	     (Prepend[-c, 1] . Table[x^i, {i, d, 0, -1}]), 
	      x]], 1];
	  lh = Length[ch];
	  rhs = Tabulate[in, {n, 0, lh - d - 1}];
	  f = (Append[#1, Take[#1, -d] . Reverse[c] + #2]) &;
	  Return[Minimize[Rec[ch, Fold[f, a, rhs]]]] ];

Rec[c_List, a_List, inh_Gfn] := Rec[c, a, ToRec[inh]];

Rec[c_List, a_List, {exp_, n_}] := Rec[c, a, ToRec[exp, n]];


(* addition and multiplication of Rec objects *)

Rec /: Rec[a__] + Rec[b__] := 
       Module[{z},
	 ToRec[Gfn[Together[
	    First[ToGfn[Rec[a], z]] + First[ToGfn[Rec[b], z]]], z]] ];

Rec /: Rec[a__]  Rec[b__] := 
       Module[{z},
	 ToRec[Gfn[Hadamard[
	    First[ToGfn[Rec[a], z]], First[ToGfn[Rec[b], z]], z], z]] ];

Rec /: (d_ /; FreeQ[d, Rec]) Rec[a_List, b_List]  := 
       Rec[a, d#& /@ b];

Rec /: Rec[a__]^n_Integer?Positive := Rec[a]^(n-1) Rec[a];


Shift[r_Rec, k_Integer:1] := Multisect[r, 1, k];

Shift[f_Gfn, k_Integer:1] := ToGfn[Shift[ToRec[f], k], Last[f]];

Shift[a_, n_Symbol, k_Integer:1] := Simplify[a /. n -> n + k];


SimplifySeq[expr_, n_Symbol] :=
   Module[{z},
      Return[ToExpPoly[ToGfn[expr, n, z], n]]
   ];

SimplifySeq[s_] := Minimize[s];

StrikeOut[s_, m_Integer, k_Integer:0] := 
    Interlace @@ 
       (Multisect[s, m, #]& /@ Complement[Range[0, m-1], {Mod[k, m]}]);

StrikeOut[s_, n_Symbol, m_Integer, k_Integer:0] :=
    ToExpPoly[StrikeOut[ToRec[s, n], m, k], n];


Tabulate[s_, it_] := 
    Module[{n, iter},    
	iter = Switch[Length[it],
		   1,  {n, 0, it[[1]] - 1},
		   2,  Switch[Head[it[[1]]], 
			   Symbol,  {n, 1, it[[2]]},
			   Integer, {n, it[[1]], it[[2]]}],
		   3,  {n, it[[2]], it[[3]]}];
	Table[Term[s, n], Evaluate[iter]] /; 
				 MemberQ[{Rec, Gfn}, Head[s]]];

Tabulate[s_, it_] := Table[s, it];


Term[r_Rec, k_Integer] :=
    Module[{a},
	Define[a, r];
	Return[a[k]] ];

Term[f_Gfn, k_Integer] := Term[ToRec[f], k];

Term[a_, n_Symbol, k_Integer] := a /. n -> k;


ToExpPoly[Gfn[f_, z_], n_Symbol] :=
    SeriesTerm[f, {z,0,n}] /. If[n >= _, a_, _] -> a;


ToExpPoly[r_Rec, n_Symbol] :=
   Module[{z},
      Return[ToExpPoly[ToGfn[r, z], n]]
   ];

ToExpPoly[ToRec[s_], n_Symbol] := ToExpPoly[s, n];

ToExpPoly[expr_, n_] := SimplifySeq[expr, n];


ToGfn[Rec[c_List, a_List], z_Symbol] := 
   Module[{i, j, d, zz, z1, r, den, fc},
     (d = Length[c];
      zz = Table[z^i, {i, d}];
      z1 = If[zz =!= {}, Drop[zz, -1], zz];
      r = Factor[((a.Table[z^i - Take[c, d-i-1].Take[z1, -(d-i-1)], 
	  {i, 0, d-1}])) / (1 - c.zz)];
      Return[Gfn[r, z]]) /; Length[a] == Length[c] ];

ToGfn[a_, n_Symbol, z_Symbol] :=
   Gfn[Simplify[PowerSum[a, {z, n, 0}]], z];

ToGfn[g_] := Minimize[g];


ToRec[Gfn[f_ /; f =!= 0, z_Symbol]] :=
   Module[{dp = Exponent[num = Numerator[Factor[f]], z], 
	   dq = Exponent[den = Denominator[Factor[f]], z], rec, d, init},
    ( rec = CoefficientList[den, z];
      rec = - Rest[rec] / First[rec];
      d = Length[rec];
      init = Pad[CoefficientList[Series[f, {z, 0, d-1}], z], d];
      Return[Rec[rec, init]] ) /; dp < dq && Together[den /. z -> 0] =!= 0
   ];


ToRec[Gfn[0, z_Symbol]] := Rec[{},{}];

ToRec[a_, n_] :=
   Module[{z},
      ToRec[ToGfn[a, n, z]]
   ];

ToRec[r_Rec] := Minimize[r];

ToRec[$Failed] := $Failed;


Verify[lhs_ == rhs_, n_Symbol] :=
   Module[{z},
      Return[Together[First[ToGfn[lhs - rhs, n, z]]] === 0]
   ];

Verify[lhs_ == rhs_] := MemberQ[{{},0}, First[Minimize[lhs - rhs]]];

           
(* ROOTS OF UNITY *)

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



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

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

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


(* SOLVING ROUTINE *)

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

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

        init = Select[conds, Function[xx, MatchQ[xx,
            a_?(MemberQ[unknowns, #]&)[_?NumberQ] == _?NumberQ]]];
        subst = Solve[init, Union @@ (PatternList[init, #[_]]& /@ 
            unknowns)];
        If[subst === {}, Return[$Failed], eqns = eqns /. First[subst]];
            
(* remember the unknown functions *)
               
        functions = G[#][z]& /@ unknowns;
        fnHeads = G /@ unknowns;
  
(* solve the equations using Solve, DSolve, or HypergeometricSolve *)
	
        If[FreeQ[eqns, Derivative],
        
            temp = Solve[Thread[eqns == 0], functions],
            
    (* before using differential equation solvers, integrate
       equations of the form  lhs' == rhs'  *)
       
            times = Integrable[#, fnHeads, z]& /@ eqns;
            eqns = PreIntegrate[#[[1]], fnHeads, {z, #[[2]]}]& /@
                Thread[List[eqns, times]];
            eqns = Table[Sum[C[ii, jj] z^(jj - 1), {jj, times[[ii]]}] +
                eqns[[ii]], {ii, Length[eqns]} ];
            eqns = eqns /. Derivative[kk_][ff_] :> D[ff, {z, kk}];
            eqns = Thread[eqns == 0];
            
    (* use HSolve then DSolve *)
            
	    temp = HSolve[eqns, functions, z, opts];
(*            temp = temp /. C -> First[unknowns]; *)
            If[!FreeL[temp, {HSolve, $Failed}],
                temp = DSolve[eqns, functions, z] 
	    ]
	];
        If[temp === {}, Return[{}]];
        If[!FreeL[temp, {Solve, DSolve, $Failed}], Return[$Failed]];
        If[!FreeL[functions /. temp, fnHeads], Return[$Failed]];
        temp = temp //. SimplifyExpLog;

(* find all initial terms occurring in conditions, equations,
   or solutions *)
               
        initials = Union @@ (PatternList[{conds, eqns, 
            functions /. temp}, #[_]]& /@ unknowns);
            
            (* sort initials by falling indices *)
            
        initials = Join @@ (Append[Cases[initials, _[#]]& /@
            Reverse[Union[PatternList[initials, _?NumberQ]]], {}]);
        pinitials = Select[initials, (#[[1]] >= 0)&];
            
            (* How many series terms will be needed to determine 
               the initials? *)
               
       (order[#] = Max[Append[First /@ PatternList[initials, 
            #[_]], 0 ]])& /@ unknowns;
            
(* determine values of initials and eliminate constants
   of integration *)
               
        constants = PatternList[functions /. temp, C[__]] // Union;

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

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

Taylor[f_, z_, a_, n_] :=
    Block[{tmp, kk, normal = Normal[SafeSeries[f, {z, a, n}]]},
	tmp = If[!PolynomialQ[Expand[normal], z],
		{normal},
		CoefficientList[normal, z]
		];
       If[Length[tmp] > n + 1, tmp = Take[tmp, n + 1]];
       Return[Join[tmp, Table[0, 
           {kk, n + 1 - Length[tmp]} ]]]]
           
CList[l_List, x_] := CList[#, x]& /@ l

CList[0, x_] := {0}

CList[a_, x_] := CoefficientList[a, x] /; a =!= 0
        
        
(* INTEGRABLE and PREINTEGRATE *)        
                
Integrable[a_Plus, b___] := Min[Integrable[#, b]& /@ (List @@ a)]

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

Integrable[a_ b_, functions_, z_] := 
    Integrable[b, functions, z] /; FreeQ[a, z]
    
Integrable[a_, functions_, z_] := 
    FirstArg[a] /; FirstHead[a] === Derivative
                           
Integrable[a___] := 0


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

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

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

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


(* FIRSTHEAD, FIRSTARG and HEADCOUNT *)

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



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

HeadCount[a_] := Which[Length[a] == 0,
                           0,
                       True,
                           HeadCount[Head[a]] + 1 ]
                           
                           
                           
(** POWER SUM **)

(* MULTISECTION *)

PowerSum[a_[k_Integer?Positive n_ + h_.], z_, n_] :=
   Module[{f, ff, m, d, cp, i, j, u, order},
      f = PowerSum[a[n], {z, n}];
      If[RationalQ[f, z],
	 ff = Denominator[f];
	 d = Exponent[ff, z];
	 cp = Resultant[ff /. z -> u, z - u^k, u];
	 rec = CoefficientList[cp, z];
	 rec = - Rest[rec] / First[rec];
	 order = Length[rec];
	 init = Table[a[k i + h], {i, 0, order - 1}];
	 Return[First[ToGfn[Rec[rec, init], z]]] ]
   ];

(* Equations (5.56) and (5.57) from GKP *)

PowerSum[Binomial[n_ + a_., k_], z_, n_] :=
   Module[{nn},
      z^(k-a)/((1-z)^(k+1)) - When[a > k, 1, 0] Sum[Binomial[nn, k] z^(nn-a), 
	  {nn,0,a-1}]] /;  IntegerQ[a] && IntegerQ[k] && TTQ[k >= 0];


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

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

PowerSum[x: Literal[Times[(_^n_)..]], z_, n_] :=
    Module[{b = First /@ x},
	Factor[1/(1 - b z)] /; FreeQ[b, n]];


(* CONVOLUTION *)

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


    (* lists and equations *)
    
PowerSum[expr_List, rest__] :=
	PowerSum[#, rest]& /@ expr

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


    (* revert to inner syntax *)
     
PowerSum[expr_, {z_Symbol, n_, n0_:0}] :=
    Block[{e = expr //. {Sum[a_, {k_, m_}] :> Sum[a, {k, 1, m}],
                         Literal[Even[a_]] :> SymbolicMod[a, 2] == 0,
                         Literal[Odd[a_]]  :> SymbolicMod[a, 2] == 1,
                         Mod -> SymbolicMod} },
	PowerSum[e, z, n, n0] /. 
            Derivative[kk_][ff_] :> D[ff, {z, kk}]
    ]


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

    (* the geometric series *)
   
PowerSum[1, z_, n_] :=
	1/(1 - z)

    (* linearity *)
   
PowerSum[c_ expr_., z_, n_] := 
    c PowerSum[expr, z, n] 	/; FreeL[c, {n, Blank}]
                                                          
PowerSum[c_.(expr1_ + expr2_), rest__] := 
        PowerSum[c expr1, rest] + PowerSum[c expr2, rest] 
	
                    
PowerSum[Sum[a_, {k_, lo_, hi_}] b_., z_, n_] :=
    Sum[PowerSum[Expand[a b], z, n],
	{k, lo, hi}]	 /;      FreeQ[{lo, hi}, n] && FreeQ[b, k]
    
                       
    (* powers *)
    
PowerSum[c_^((a_ + b_)d_) e_., rest__] := 
    PowerSum[c^(a d + b d) e, rest]

PowerSum[c_^(a_ + b_) e_., z_, n_] := 
    c^a PowerSum[c^b e, z, n]	/; FreeL[c, {n, Blank}] && FreeL[a, {n, Blank}]
                        
PowerSum[(a_ + b_)^c_Integer?Positive d_., rest__] := 
    PowerSum[Expand[(a + b)^c d], rest]
	
                                          
                                          
    (* exponentials *)
    
PowerSum[e_. c_^((a_. n_ + b_.)d_.), z_, n_] :=
    (c^(b d) PowerSum[e, z, n] /. z -> c^(a d) z)	/;
			( (FreeL[#, {n, Blank}]& /@ And[a, b, c, d]) &&
			  FreeL[e, {z, Blank}]
			)
                    
                    
    (* binomials *)
    
PowerSum[Binomial[a_, n_ + k_.], z_, n_] :=
    z^(-k) *
	((1 + z)^a - Sum[Binomial[a, m] z^m, {m, 0, k - 1}])	/;
			FreeQ[a, n] && IntegerQ[k]

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


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

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

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

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



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


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

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

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



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



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

PowerSum[(n_ + a_)^m_Integer?Negative, z_, n_] :=
    Block[{j},
        z^(-a) PolyLog[-m, z] - Sum[z^(j-a) j^m, {j, 1, a-1}]
    ] /;	IntegerQ[a] && TTQ[a > 0]
        
                           
    (* convolutions *)
    
PowerSum[Sum[expr_, {k_, lo_, hi_}], z_, n_] :=
    Block[{nfree, nfull, s, aa = Expand[lo], bb = Expand[hi - n]},
         ({nfree, nfull} = If[Head[expr] === Times, 
           			{Select[expr,  FreeQ[#, n]&],
					Select[expr, !FreeQ[#, n]&]},
	        		If[FreeQ[expr, n],
					{expr, 1},
					{1, expr}
				]
			];
        Sum[z^k nfree PowerSum[Expand[nfull /. n -> s + k], z, s, -k],
		 {k, aa, bb - 1}] +
        PowerSum[nfree PowerSum[Expand[nfull /. n -> s + k], z, s, -bb],
                  z, k, Max[aa, bb] ]
    	)	 /; 	
                        FreeQ[aa, n] && FreeQ[bb, n] 
    ]


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

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


    (* PowerSum under Series *)
    
PowerSum /: 
    Series[PowerSum[a_, z_, n_], {z_, 0, m_}] :=
        Block[{k}, 
        Sum[(a /. n -> k) z^k, {k, 0, m}] + O[z]^(m+1)]
    
PowerSum[a_ b_, z_, n_] :=
   Module[{f, g},
      f = PowerSum[a, {z, n}];
      g = PowerSum[b, {z, n}];
      If[RationalQ[f, z] && RationalQ[g, z],
	 Return[Hadamard[f, g, z]] ]
   ];

PowerSum[a_^k_Integer?Positive, z_, n_] :=
   Module[{f, g},
      f = PowerSum[a, {z, n}];
      g = PowerSum[a^(k-1), {z, n}];
      If[RationalQ[f, z],
	 Return[Hadamard[f, g, z]] ]
   ];


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

PowerSum[a_[n_ + d_?Negative], z_, n_] :=
    Module[{j, k},
       (PowerSum[a[k], z, k] + Sum[a[k] z^k, {k, d, -1}])/z^d // Factor
    ] /; IntegerQ[d];

PowerSum[a_[n_ + d_.], z_, n_] :=
    Module[{j, k},
       (PowerSum[a[k], z, k] - Sum[a[k] z^k, {k, 0, d - 1}])/z^d // Factor
    ] /; IntegerQ[d];
	    

    

(* GENERATING FUNCTION *)

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


(* SERIES TERM *)

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

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

    (* finite n *)
        
SeriesTerm[f_, z_Symbol, 0] := 
	SeriesTerm[z f, z, 1] /; FreeQ[f, HypergeometricF]

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

    (* linearity *)
    
SeriesTerm[f_Plus, rest__] :=
	SeriesTerm[#, rest]& /@ f

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

    (* SeriesTerm is inverse of Gf and of PowerSum *)
    
SeriesTerm[Gf[a_][z_], z_Symbol, n_] := a[n]

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

    (* SeriesTerm is close to the inverse of EGf and of 
       ExponentialPowerSum *)
    
SeriesTerm[EGf[a_][z_], z_Symbol, n_] := a[n] / n!

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



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

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


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

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


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

        
     (* exponentials *)
     
SeriesTerm[E^(k_.z_ + m_.), z_Symbol, n_] :=
	E^m k^n/n!  /;	FreeQ[k, z] && FreeQ[m, z]
                                
SeriesTerm[a_^b_, z_Symbol, n_] :=
	SeriesTerm[E^(b Log[a]), z, n] /;	a =!= E && FreeQ[a, z]
                                    
    (* logarithms *)
    
SeriesTerm[Log[(a_ + b_.z_)^k_.], z_Symbol, n_] :=
    When[n == 0, Release[k Log[a]], Release[-k(-b/a)^n / n]]  /;
                   FreeQ[k, z] && FreeQ[a, z] && FreeQ[b, z]
                   
SeriesTerm[Log[a_Times], rest___] :=
    Plus @@ (SeriesTerm[Log[#], rest]& /@ a) 
    
                                    
    (* powers and constants *)
    
SeriesTerm[z_^k_.f_, z_Symbol, n_] :=
	SeriesTerm[f, z, n-k]  /; IntegerQ[k] && FreeQ[k, z] && !FreeQ[f, z]                                                         
SeriesTerm[z_^k_., z_Symbol, n_] :=
	When[n == k, 1, 0]  /; IntegerQ[k] && FreeQ[k, z]                                                        
SeriesTerm[f_, z_Symbol, n_] :=
	f When[n == 0, 1, 0]  /; FreeQ[f, z]


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

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

SeriesRat[f_. g_^k_Integer?Negative, z_Symbol, n_] :=
  
    Block[{ff = Expand[f], gg = Expand[g], poleList, poleSet,
           mp, denom = Expand[g^(-k)], i, deg, mod, modList,
           displace, temp, l, pr = precision, inverseRoot, kk},
         deg = Exponent[gg, z];
         If[useMod, 
             l = CoefficientList[ff, z];
             mod = Exponent[gg, z, GCD];
             modList = Union[Mod[#, mod]& /@ (Select[Range[Length[l]],
                 (l[[#]] =!= 0)&] - 1)];
             If[Length[modList] == 1,
                 displace = First[modList],
                {displace, mod} = {0, 1}
	     ],
             {displace, mod} = {0, 1}
	 ];     
         {ff, gg, denom} = {Expand[ff/z^displace], gg, denom} /.
             z -> z^(1/mod);
         deg = deg/mod;
         Which[precision === Automatic && deg > 2, 
                   pr = Indeterminate,
               precision === Automatic && deg <= 2, 
                   pr = Infinity
	 ];
         If[pr == Infinity && !FreeQ[Roots[gg == 0, z], Roots],
             pr = Indeterminate];
         Which[
             pr === Indeterminate, 
             poleList = Cancel[#[[1,2]]& /@ {ToRules[
                            NRoots[gg == 0, z] ]}];
             poleSet = Union[poleList];
             mp = -k Count[poleList, #]& /@ poleSet;
             poleUse = poleSet;
             temp = SeriesTerm[Plus @@
			Table[
			  SeriesDivide[
			    Horner[ff, {z, poleSet[[i]], mp[[i]]}],
	     		    Drop[Horner[denom, {z, poleSet[[i]], 2 mp[[i]]}],
	     		      mp[[i]]
			    ],
			    mp[[i]]
			  ] .
                          Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
                          {i, Length[poleSet]}
			],
			z, Expand[(n-displace)/mod]
		    ],
             
             pr =!= Infinity, 
                 poleList = Cancel[#[[1,2]]& /@ {ToRules[
                            NRoots[gg == 0, z, pr] ]}];
                 poleSet = Union[poleList];
                 mp = -k Count[poleList, #]& /@ poleSet;
                 poleUse = poleSet;
                 temp = SeriesTerm[Plus @@
			  Table[
			    SeriesDivide[
			      Horner[ff, {z, poleSet[[i]], mp[[i]]}],
	     		      Drop[Horner[denom, {z, poleSet[[i]], 2 mp[[i]]}],
	     			mp[[i]]
			      ],
			      mp[[i]]
			    ] .
             		    Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
             		    {i, Length[poleSet]}
			  ],
			  z, Expand[(n-displace)/mod]
			],
             
             pr === Infinity && deg <= 2, 
                 poleList = Cancel[#[[1,2]]& /@ {ToRules[
                            Roots[gg == 0, z] ]}];
                 poleSet = Union[poleList];
                 mp = -k Count[poleList, #]& /@ poleSet;
                 inverseRoot = Expand[(1 - gg/(gg /. z -> 0))/z];
                 poleUse = ArgPi[Apart[
                 	inverseRoot /. z -> #]]^(-1)& /@ poleSet;
                 temp = SeriesTerm[Plus @@
			  Table[
			    SeriesDivide[
			      Horner[ff, {z, poleSet[[i]], mp[[i]]}],
	     		      Drop[Horner[denom, {z, poleSet[[i]], 2 mp[[i]]}],
	     		        mp[[i]]
			      ],
			      mp[[i]]
			    ] .
             		    Table[(z - poleUse[[i]])^(-j), {j, mp[[i]], 1, -1}],
             		    {i, Length[poleSet]}
			  ],
			  z, Expand[(n-displace)/mod]
			],
             
             pr === Infinity && deg > 2, 
                 poleList = Cancel[#[[1,2]]& /@ {ToRules[
                            Roots[CoefficientList[gg, z].
                            Table[z^kk, {kk, deg, 0, -1}] == 0, z] ]}];
                 poleSet = Union[poleList];
                 mp = -k Count[poleList, #]& /@ poleSet;
                 poleUse = Complex[Together[ReP[#]//.SimplifyCubic],
                           Together[ImP[#]//.SimplifyCubic]]& /@ poleSet;
                 poleUse = poleUse /. Complex[a_, 0] -> a;
                 If[FreeQ[poleUse, Complex],
			makeReal = False
		 ];
                 temp = SeriesTerm[Plus @@
			  Table[
			    SeriesDivide[
			      Horner[ff /. z -> z / poleUse[[i]],
				{z, 1, mp[[i]]}],
                              Drop[Horner[denom /. z -> z / poleUse[[i]],
				{z, 1, 2 mp[[i]]}],
			        mp[[i]]
			      ],
			      mp[[i]]
			    ] .
             		    Table[(z poleUse[[i]] - 1)^(-j),
				{j, mp[[i]], 1, -1}],
             		    {i, Length[poleSet]}
			  ],
			  z, Expand[(n-displace)/mod]
			]
              ];

         temp = temp //. a_ c_If + b_ c_ -> (a + b) c;
	If[makeReal,
	  If[pr === Infinity,
		temp = temp /. SimplifyComplexE1];
	  If[pr === Infinity && deg == 2,
		temp = temp /. SimplifyComplex2];
	  If[pr === Infinity && deg == 3,
		temp = temp /. SimplifyComplex3];
	  If[pr === Infinity && deg == 4,
		temp = temp //. SimplifyComplex4];
	  If[pr === Infinity,
		temp = temp /. SimplifyComplexE2];
	  If[pr =!= Infinity,
		temp = temp //. SimplifyComplexN]
	];
        temp = temp /. Complex[a_, b_] :> a + b I;
        Return[If[Evaluate[SymbolicMod[n, mod] == Mod[displace, mod]], 
                Evaluate[temp], 0]]
        ] /; PolynomialQ[Expand[f], z] && PolynomialQ[Expand[g], z] &&
		!FreeQ[g, z]
         

    (* special functions *)

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

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

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

    

    (* products with a rational function *)
    
SeriesTerm[a_ b_, z_Symbol, n_] :=
    Block[{temp, deg, i, m, bn = SeriesTerm[b, z, m], sd = sumDepth},
         Increment[sumDepth];
         Off[General::itervar];
         If[PolynomialQ[Expand[a], z],
             deg = Exponent[a, z];
             temp = CoefficientList[a, z] . 
                 Table[bn /. m -> i, {i, n, n - deg, -1}] /. If -> When,
             temp = Sum[SeriesTerm[Apart[a, z], z, n - K[sd]] bn /. m -> K[sd], 
               		Evaluate[{K[sd], -PoleMultiplicity[b, {z, 0}],
              		n + PoleMultiplicity[a, {z, 0}]}] ]
         ];
         On[General::itervar];
         Return[temp] ] /; SeparateProduct[a b, (PolynomialQ[Expand[#], z] || 
                           PolynomialQ[Expand[#^(-1)], z])& ] == {a, b}
        
    (* general products *)
    
SeriesTerm[a_ b_, z_Symbol, n_] := 
    Block[{temp, sd = sumDepth},
         Increment[sumDepth];
         Off[General::itervar];
         temp = Sum[SeriesTerm[a, z, n - K[sd]] * 
             SeriesTerm[b, z, K[sd]],
            Evaluate[{K[sd], -PoleMultiplicity[b, {z, 0}],
                   n + PoleMultiplicity[a, {z, 0}]}] ];
         On[General::itervar];
         Return[temp] ] /; SeparateProduct[a b, (PolynomialQ[Expand[#], z] || 
                           PolynomialQ[Expand[#^(-1)], z])& ] == {1, a b}
        
SeriesTerm[a_^k_Integer?Positive, z_Symbol, n_] := 
    Block[{temp, sd = sumDepth},
         Increment[sumDepth];
         Off[General::itervar];
         temp = Sum[SeriesTerm[a, z, n - K[sd]] * 
             SeriesTerm[a^(k - 1), z, K[sd]],
            Evaluate[{K[sd], -PoleMultiplicity[a^(k - 1), {z, 0}],
                   n + PoleMultiplicity[a, {z, 0}]}] ];
         On[General::itervar];
         Return[temp] ] /; !FreeQ[a, z]


   (* trinomials *)
   
SeriesTerm[(a_ + b_.z_ + c_.z_^2)^alpha_, z_Symbol, n_] := 
    Block[{temp},
         Off[General::itervar];
         temp = Sum[Binomial[alpha, K[sumDepth]] Binomial[K[sumDepth],
             n - K[sumDepth]] a^(alpha - K[sumDepth]) b^(2 K[sumDepth] -
             n) c^(n - K[sumDepth]), Evaluate[{K[sumDepth], 0, n}]];
         On[General::itervar];
         ++sumDepth;
         Return[temp] ] /; FreeQ[{a,b,c}, z]
         
         
   (* derivatives *)
   
SeriesTerm[Derivative[k_][f_][z_], z_Symbol, n_] :=
    Pochhammer[n + 1, k] SeriesTerm[f[z], z, n + k]
      


(* HORNER and SERIES DIVIDE *)

Horner[p_, {z_, alpha_, m_}] :=
    Block[{qh = Expand[p], ih, kh, nh, ah},
        {nh, ah} = {Exponent[qh, z], CoefficientList[qh, z]};
        Do[ah[[ih]] = alpha ah[[ih + 1]] + ah[[ih]], {kh, 1, Min[nh, m]},
            {ih, nh, kh, -1}];
	If[m > nh + 1,
        	Do[AppendTo[ah, 0], {ih, m - nh - 1}]
	];
        Return[Take[ah, m]]
    ]
        
        
SeriesDivide[a_, b_, n_] :=
    Block[{is, cs = Range[n], ks},
        Do[cs[[is]] = (a[[is]] - Sum[cs[[ks]] b[[is - ks + 1]], 
            {ks, 1, is - 1}])/b[[1]], 
            {is, 1, n}];
        Return[cs] ]
        
        

(* SYMBOLIC MOD *)

SymbolicMod[n_, 1] = 0

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

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

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

     
(* PARTIAL FRACTIONS *)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

ReP[x_^Rational[p_,q_]] := 
  Module[{u = x^Quotient[p,q], v = x^(Mod[p,q]/q)},
    ReP[u] ReP[v] - ImP[u] ImP[v]
  ]	/; 	AbsV[p] >= AbsV[q]

ImP[x_^Rational[p_,q_]] :=
  Module[{u = x^Quotient[p,q], v = x^(Mod[p,q]/q)},
    ReP[u] ImP[v] + ImP[u] ReP[v]
  ]	/; 	AbsV[p] >= AbsV[q]
    


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

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

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

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

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

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



Angle[z_] := 
    Block[{angle =
	     Block[{x = ReP[z], y = ImP[z]},
                 Which[
		   TrueQ[Expand[x] == 0], Pi/2 Sgn[y],
                   TrueQ[Expand[y] == 0], Pi/2 (1 - Sgn[x]),
                   True, ArcTan[y/x] - Pi/2 (Sgn[x]-1) Sgn[y]
	         ]
	     ]},
	  angle /; FreeQ[angle,Sgn]
    ] /;  FreeQ[{ReP[z],ImP[z]},ReP] && FreeQ[{ReP[z],ImP[z]},ImP]
                  
Sgn[0] = 0
Sgn[a_?Positive] = 1
Sgn[a_?Negative] = -1
Sgn[a_ b_] := Sgn[a] Sgn[b]
Sgn[a_^b_] := Sgn[a]^b
Sgn[Sqrt[a_]] = 1
Sgn[a_ + b_] := Sgn[Chop[N[a + b]]]

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

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

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

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


SimplifyCubic = Join[SimplifyTrig, SimplifySqrt]
    

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

SimplifySum = {
    Sum[0, {__}] -> 0,
    Sum[1, {k_, a_, b_}] :> When[a <= b, b - a + 1, 0],
    Sum[a_. b_, {k_, c__}] :> b Sum[a, {k, c}] /; 
        FreeQ[b, k],
    Sum[If[k_ == m_, a_, b_] c_., {k_, m_, n_}] :>
       ((a c) /. k -> m) When[n >= m, 1, 0] + 
        Sum[b c, {k, m + 1, n}],
    Sum[(If[-k_ + n_ >= 0, a_, b_] c_. + d_.) e_., {k_, m_, n_}] :>
         Sum[Evaluate[Expand[a c e]], {k, m, n}] + Expand[d e],
    Sum[If[k_ >= m_, a_, b_] c_., {k_, m_, n_}] :>
         Sum[Evaluate[Expand[a c]], {k, m, n}] }
                  
SimplifyComplexE1 = {
    a_.I^c_ + d_.(-I)^c_ :> ArgPi[a] E^Expand[I Pi c/2] +
        ArgPi[d] E^Expand[-I Pi c/2] /;
	  (ConjugateQ[a, d] && FreeQ[{ArgPi[a],ArgPi[d]},ArgPi]),
    a_.E^b_ + d_.E^e_ :> 2 a Cos[ImP[b]] /; 
       (ImP[a] == a - d == ReP[b] == PolynomialMod[b + e, 2 Pi I] == 0),     
    a_.E^b_ + d_.E^e_ :> -2 ImP[a] Sin[ImP[b]] /;  
       (ReP[a] == a + d == ReP[b] == PolynomialMod[b + e, 2 Pi I] == 0),     
    a_.E^b_ + d_.E^e_ :> ArgPi[a] E^b + ArgPi[d] E^e /;
       FreeQ[{ArgPi[a],ArgPi[d]},ArgPi]
		    }


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

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

        
SimplifyComplexN = {
    a_.b_Complex^c_ + d_.e_Complex^c_ :> 
        2 Abs[a]Abs[b]^c Cos[Arg[a] + c Arg[b]] /; 
        (ConjugateQ[a, d] && ConjugateQ[b, e] && ImP[c] == 0)}
        
SimplifyCos[x_] :=
    Block[{x1 = (x /. Cos[a_] :> Cos[Expand[a]]) //. SimplifyTrig,
           x2 =  x //. SimplifyTrig },
        If[LeafCount[x1] <= LeafCount[x2], x1, x2 ]]
        
SimplifyFloat = {1. -> 1, -1. -> -1, 0. -> 0}


(* FACTORIAL SIMPLIFY *)

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


DeFact[expr_] :=
 
    Block[{kin, j, sub, diff, t, tt, arg = {}, ex = expr},
    
        While[!FreeQ[ex, Factorial],
            t = First[Position[ex, Factorial[_]]];
            tt = First[Part @@ Prepend[t, ex]];
            kin = Scan[If[IntegerQ[Expand[# - tt]],
                Return[Fact[#]] ]&, arg];
            If[kin === Null,
                AppendTo[arg, tt]; sub = Fact[tt],
                diff = Expand[First[kin] - tt];
                If[diff > 0,
                    sub = kin / Product[tt + j, {j, diff}],
                    sub = kin Product[tt - j + 1, {j, -diff}] ] ];
            ex = MapAt[sub&, ex, {t}] ];
        Factor[ex] /. Fact -> Factorial ]
        
        
PostSimplify[expr_] :=
    (* DEBUG; added "_" to n! below... blank was missing I think! *)
    Block[{ex = expr /. n_! :> Expand[n]!, i},

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

(* ARGPI *)

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

(* ENTAILS *)

Entails[True, a_] := a;

Entails[a_, a_] = True

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

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

(* FIRST POS and FREE L *)

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

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


(* ISOLVE, INEQSOLVE, LEQ *)

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

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

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

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

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

IneqSolve[a_ && b_, x_] := 
    Block[{aa = IneqSolve[a, x], 
           bb = IneqSolve[b, x],
           merge, accum, begin, nn},
        If[FreeQ[aa, IneqSolve] && FreeQ[bb, IneqSolve],
        aa = Append[aa, {}];
        bb = Append[bb, {}];
        merge = Sort[Join[
            Thread[{Join @@ (aa /. MyInterval -> List),
                    Join @@ (aa /. MyInterval[_,_] -> {-1,1}) }],
            Thread[{Join @@ (bb /. MyInterval -> List),
                    Join @@ (bb /. MyInterval[_,_] -> {-1,1}) }] ],
                    LEQ];
        accum = Thread[{First /@ merge, 
                        FoldList[Plus, First[Last /@ merge],
				Drop[Last /@ merge, 1]]}];
        begin = Join @@ Append[Position[accum, {_,-2}], {}];
        Return[(MyInterval @@ # & /@ Thread[{begin, begin + 1}]) /. 
                             nn_Integer :> First[accum[[nn]]]],
        Return[aa && bb] ]]
        
IneqSolve[!a_, x_] := 
    Block[{aa = IneqSolve[a, x], merge, c, d},
        If[FreeQ[aa, IneqSolve],
            aa = Append[aa, {}];
            merge = Join @@ (aa /. MyInterval[c_, d_] -> {c-1, d+1});
            If[merge === {},
                merge = {-Infinity, Infinity},
                If[First[merge] === -Infinity, 
                    merge = Rest[merge],
                    PrependTo[merge, -Infinity] ];
                If[Last[merge] === Infinity, 
                    merge = Drop[merge,-1],
                    AppendTo[merge, Infinity] ] ];
            Return[MyInterval @@ # & /@ Partition[merge, 2]],
          Return[!aa]] ]
        
IneqSolve[a_ || b_, x_] := 
    Block[{aa = IneqSolve[!(!a && !b), x]},
        If[FreeQ[aa, IneqSolve],
            Return[aa],
          Return[IneqSolve[a, x] || IneqSolve[b, x]]] ]

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


IneqSolve[h_[a_, b_, c__], x_] := 
    IneqSolve[And @@ (h @@ # & /@ 
                Partition[{a,b,c},2,1]), x] /;
    MemberQ[{Less, LessEqual, Greater, GreaterEqual, Equal, 
             Unequal}, h]
            
IneqSolve[ineq: h_[a_, b_], x_] :=
    Block[{lhs = Together[a - b],
           num, den, zeroList, poleList, points, merge, c, d},
           
       {num, den} = {Numerator[lhs], Denominator[lhs]};
        Off[Roots::neq];
        zeroList = If[FreeQ[num, x], {},
            #[[1,2]]& /@ {ToRules[NRoots[num == 0, x]]} ];
        poleList = If[FreeQ[den, x], {},
            #[[1,2]]& /@ {ToRules[NRoots[den == 0, x]]} ];
        On[Roots::neq];
        zeroList = Select[zeroList, (Im[#] == 0)&] // Union;
        poleList = Select[poleList, (Im[#] == 0)&] // Union;
        points = Union[zeroList, poleList];
        merge = Partition[Append[Prepend[points,-Infinity],Infinity], 
            2, 1];
        merge = Select[merge, (ineq /. x :> 
                  Which[# === {-Infinity, Infinity}, 0,
                        #[[1]] === -Infinity, #[[2]] - 1,
                        #[[2]] ===  Infinity, #[[1]] + 1,
                        True, (Plus @@ #)/2 ] )&];
        merge = (MyInterval @@ # & /@ merge) /. 
            MyInterval[c_, d_] :> MyInterval[Ceiling[c], Floor[d]];
        merge = merge /. {Ceiling[-Infinity] -> -Infinity,
                          Floor[Infinity] -> Infinity};
        merge = Select[merge, #[[1]] <= #[[2]]&];
        zeroList = Sort[Join[Floor /@ zeroList, Ceiling /@ zeroList]];
        zeroList = Select[zeroList, ((den /. x->#) != 0)&];
        zeroList = Select[zeroList, (ineq /. x->#)&];
        merge = Sort[Join[merge, Thread[MyInterval[zeroList, zeroList]]],
            LEQ];
        merge = merge /. 
            MyInterval[c_?((NumberQ[#] && ((den /. x:>#) == 0 ||
                !ineq/.x->#))&), d_] :> MyInterval[c+1, d];
        merge = merge /. 
            MyInterval[c_, d_?((NumberQ[#] && ((den /. x:>#) == 0 ||
                !ineq/.x->#))&)] :> MyInterval[c, d-1];
        merge = Select[merge, #[[1]] <= #[[2]]&];
        args = List @@ (Join @@ Append[merge, MyInterval[]]);
        If[args === {}, {},
           {first, last} = {First[args], Last[args]};
            MyInterval @@ # & /@ Partition[Append[Prepend[Join @@ 
                Append[Select[Partition[Drop[Rest[args], -1], 2], 
                (#[[2]] - #[[1]] > 1)&], {}], first], last], 2] ] ]  /; 
        MemberQ[
           {Less, LessEqual, Greater, GreaterEqual, Equal, Unequal}, h]
                
 
ISolve[expr_] :=
    Block[{l = UserSymbols[expr], n, temp, p, nn},
        If[Length[l] == 1,
            n = First[l];
            temp = IneqSolve[expr /. n -> nn, nn] /. nn -> n;
            temp = temp /. List[p___MyInterval] :> Or[p];
            temp = temp /.
              {MyInterval[-Infinity, Infinity] -> True,
               MyInterval[-Infinity, b_] -> n <= b,
               MyInterval[a_, Infinity] -> n >= a,
               MyInterval[a_, a_] -> n == a,
               MyInterval[a_, b_] -> a <= n <= b};
            temp = temp /. IneqSolve[a_, n_] -> a;
            Return[temp],
          Return[expr]] ]
        
LEQ[h_[a_,b_], h_[c_,d_]] := (a < c) || (a === c && b <= d)


(* INFO *)

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


(* LEADING COEF *)

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

(* MAKE LIST *)

MakeList[expr_, head_:List] := 
    If[Head[expr] === head, List @@ expr, List @ expr]
    
    
(* MAKE TRINOMIAL and LIST COMPLEMENT *)

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


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



(* PARSING ROUTINE *)

Parse[list1_, list2_, n_] :=
    Module[{l1 = MakeList[ReleaseHold @@ (Hold[list1] /. Condition -> Pair)],
           unknowns = MakeList[list2],
           eqns, conds, recur, extra, range, startValues, aa, kk,
           lo, hi, check},
         
        unknowns = Head /@ unknowns;
        check = If[Head[#] === Pair, #[[1]], #]& /@ l1;
        If[!MatchQ[check, {__Equal}], 
            Message[Parse::eqn, check];
            Return[$Failed] ];
            
        eqns =  Select[l1, !FreeQ[#, n]&];
        conds = Select[l1,  FreeQ[#, n]&];
        
            (* make  n >= 0  the default range of validity *)
        
        eqns = If[Head[#] =!= Pair, Pair[#, n >= 0], #]& /@ eqns; 
        
            (* solve inequalities, and separate recurrences from 
               initial conditions *)
               
        eqns = Pair[#[[1]], IneqSolve[#[[2]], n]]& /@ eqns;
        recur = Select[eqns, !FreeQ[#, Infinity]&];
        extra = Union[Select[eqns,  FreeQ[#, Infinity]&],
                      Pair[#[[1]], Drop[#[[2]], -1]]& /@ recur];
        recur = Pair[#[[1]], Last[#[[2]]][[1]]]& /@ recur;
        extra = Select[extra, #[[2]] =!= {} &];
        extra = Union @@ (Thread /@ extra);
        conds = Union[conds, Union @@ (Table[#[[1]], {n, #[[2,1]],
            #[[2,2]]}]& /@ extra) ];
            
            (* shift n so that all recurrences will be valid
               for  n >= 0  *)

        recur = (#[[1]] /. (n -> n + #[[2]]) )& /@ recur;
        
            (* determine starting values of n for all sequences *)
               
       {recur, conds} = {recur, conds} //. 
            Sum[a_, {k_, m_}] :> Sum[a, {k, 1, m}];
       {recur, conds} = {recur, conds} /.
            {Literal[Even[a_]] :> SymbolicMod[a, 2] == 0,
             Literal[Odd[a_]]  :> SymbolicMod[a, 2] == 1,
             Mod -> SymbolicMod};
        range = {conds, recur /. (Sum[aa_, {kk_, lo_, hi_}] :>
            {aa /. kk -> lo, aa /. kk -> hi} )};
        startValues = Min[First /@ (PatternList[range, #[_]] /. 
            n -> 0)]& /@ unknowns;
            
        Return[{recur, conds, unknowns, startValues}] ]
            

(* PATTERN LIST *)

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

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

(* POLE MULTIPLICITY *)

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

(* RESET *)

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

SafeFirst[l_] :=
    If[Length[l] == 0, Null, First[l]]
    
SafeSeries[f_, {z_, a_, n_Integer}] :=
    Block[{poles = z /. Solve[Denominator[Together[f]] == 0, z,
           		InverseFunctions -> True],
	  ord, tem},

        ord = Max[n, Count[poles, a]];
        tem = Check[Series[f, {z, a, ord}], $Failed, Series::esss];
        If[tem =!= $Failed, tem = Series[tem, {z, a, n}]];
        Return[tem]
    ]
                

(* SEPARATE PRODUCT *)

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

(* USER SYMBOLS *)

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

(* TTQ, WHEN *)

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

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


(* EVEN, ODD *)

Even[n_]:=EvenQ[n] /; !SameQ[Head[n],Symbol]
Odd[n_]:=OddQ[n] /; !SameQ[Head[n],Symbol]
          


End[]

EndPackage[];

