(* :Context: DiscreteMath`DiscreteStep` *)

(* :Title: DiscreteStep *)

(* :Author: E.C.Martin, John Novak *)

(* :Summary:
This package defines the DiscreteStep function of
an integer-valued variable.  The package extends Sum to accept
DiscreteStep.
*)

(* :Keywords: step *)

(* :Package Version: 1.0 *)

(* :Mathematica Version: 3.0 *)

(* :Copyright: Copyright 1994-1999, Wolfram Research, Inc.*)

(* :History:
    Version 1.0, June 1998. Based on V1.1 of KroneckerDelta.m by ECM.
*)

(* :Discussion:
*)

(* :Reference:
*)
		
(* :Warning:
*)

(* : Limitations:
      There is little or no support for Sum[expr, {i, imin, imax, di}]
        where expr includes DiscreteStep and di != 1.
      There is little or no support for Sum expressions involving
        DiscreteStep[c0 + c1 i] where Abs[c1] != 1.
      In general, there is limited support of filtering via
        summing an expression involving DiscreteStep.
*)


BeginPackage["DiscreteMath`DiscreteStep`"]

DiscreteStep::usage =
"DiscreteStep[n] is a function of an integer-valued variable n that is 1 for 
for n >= 0 and 0 for n < 0.  For integer-valued n1, n2, ...,
DiscreteStep[n1, n2, ...] is 1 for (n1 >= 0) && (n2 >= 0) && ... and 0 for
(n1 < 0) || (n2 < 0) || ... ."

SimplifyDiscreteStep::usage =
"SimplifyDiscreteStep[expr] returns a form of expr with fewer instances of
DiscreteStep."

Begin["`Private`"]

(* Note: negative and positive are used to determine whether certain
        assumptions made about parameters are valid. *)

positive[a_] := Module[{boolean},
        boolean /; (boolean = Positive[N[a]] || Negative[N[Simplify[-a]]];
                    boolean === True || boolean === False)]
positive[a_Plus] := With[{list = Map[positive, Apply[List,a]]},
                        Apply[And, list] /; FreeQ[list, positive]]
positive[a_Times] := Module[{l},
                        EvenQ[Length[Cases[l, 1]]] /;
                          (l = Map[If[TrueQ[negative[#]], 1,
                                   If[TrueQ[positive[#]], 0,
                                        Null]]&, Apply[List, a]];
                           FreeQ[l, Null])
                     ]
positive[Power[a_, b_]] := True /; (positive[a] || negative[Simplify[-a]]) &&
         FreeQ[N[b], Complex]
positive[Power[a_, Rational[n_, m_]]] := False /;
        NumberQ[n] && NumberQ[m] && FreeQ[N[{n, m}], Complex] &&
        (negative[a] || positive[Simplify[-a]])
  (* next rule assumes that a^(1/2) is positive, for a unspecified *)
positive[Power[a_, Rational[n_, m_]]] := True /;
        NumberQ[n] && NumberQ[m] && FreeQ[N[{n, m}], Complex] &&
        !TrueQ[negative[a] || positive[Simplify[-a]] || !FreeQ[N[a], Complex]]
positive[Power[a_, b_?EvenQ]] := True /; (negative[a] || positive[Simplify[-a]])
positive[a_] := True /; !NumberQ[N[a]] && Sign[a]===1

negative[a_] := Module[{boolean},
        boolean /; (boolean = Negative[N[a]] || Positive[N[Simplify[-a]]];
                    boolean === True || boolean === False)]
negative[a_Plus] := With[{list = Map[negative, Apply[List,a]]},
                        Apply[And, list] /; FreeQ[list, negative]]
negative[a_Times] := Module[{l = Apply[List,a]},
                        OddQ[Length[Cases[l, 1]]] /;
                          (l = Map[If[TrueQ[negative[#]], 1,
                                If[TrueQ[positive[#]], 0,
                                        Null]]&, Apply[List, a]];
                           FreeQ[l, Null])
                     ]
negative[Power[a_, b_]] := False /; (positive[a] || negative[Simplify[-a]]) &&
         FreeQ[N[b], Complex]
negative[Power[a_, Rational[n_, m_]]] := False /;
        NumberQ[n] && NumberQ[m] && FreeQ[N[{n, m}], Complex] &&
        (negative[a] || positive[Simplify[-a]])
negative[Power[a_, b_?OddQ]] := True /; (negative[a] || positive[Simplify[-a]])
negative[a_] := True /; !NumberQ[N[a]] && Sign[a]===-1


ConditionToIList[cond_, symbolList_] :=
 Module[{m, m0, m1},
  Switch[cond,
   GreaterEqual[n_ + n1_., n0_] /; ({m, m0, m1} = {n, n0, n1};
				 MemberQ[symbolList, n] && FreeQ[{n0, n1}, n]),
                        {m, m0-m1, Infinity},
   GreaterEqual[n0_, n_ + n1_.] /; ({m, m0, m1} = {n, n0, n1};
				 MemberQ[symbolList, n] && FreeQ[{n0, n1}, n]),
                        {m, -Infinity, m0-m1},
   Greater[n_ + n1_., n0_] /; ({m, m0, m1} = {n, n0, n1};
				MemberQ[symbolList, n] && FreeQ[{n0, n1}, n]),
                        {m, m0-m1+1, Infinity},
   Greater[n0_, n_ + n1_.] /; ({m, m0, m1} = {n, n0, n1};
				MemberQ[symbolList, n] && FreeQ[{n0, n1}, n]),
                        {m, -Infinity, m0-m1-1},
   LessEqual[n_ + n1_., n0_] /; ({m, m0, m1} = {n, n0, n1};
				 MemberQ[symbolList, n] && FreeQ[{n0, n1}, n]),
			{m, -Infinity, m0-m1},
   LessEqual[n0_, n_ + n1_.] /; ({m, m0, m1} = {n, n0, n1};
				MemberQ[symbolList, n] && FreeQ[{n0, n1}, n]),
			{m, m0-m1, Infinity},
   Less[n_ + n1_., n0_] /; ({m, m0, m1} = {n, n0, n1};
			MemberQ[symbolList, n] && FreeQ[{n0, n1}, n]),
			{m, -Infinity, m0-m1-1},
   Less[n0_, n_ + n1_.] /; ({m, m0, m1} = {n, n0, n1};
				MemberQ[symbolList, n] && FreeQ[{n0, n1}, n]),
			{m, m0-m1+1, Infinity},
   Equal[n_ + n1_., n0_] /; ({m, m0, m1} = {n, n0, n1};
				MemberQ[symbolList, n] && FreeQ[{n0, n1}, n]),
			{m, m0-m1, m0-m1},
   Equal[n0_, n_ + n1_.] /; ({m, m0, m1} = {n, n0, n1};
				MemberQ[symbolList, n] && FreeQ[{n0, n1}, n]),
			{m, m0-m1, m0-m1},
	_,
			$Failed
  ]
 ]

IListToDiscreteStep[ilist_, a_, b_] :=
 Module[{m, m0},
  Switch[ilist,
    {n_, -Infinity, n0_},
       (m = ilist[[1]]; m0 = ilist[[3]];
	b DiscreteStep[m-(m0+1)] + a DiscreteStep[-m+m0] ),
    {n_, n0_, Infinity},
       (m = ilist[[1]]; m0 = ilist[[2]];
	a DiscreteStep[m-m0] + b DiscreteStep[-m+(m0-1)] ),
    {n_, n0_, n0_},
       (m = ilist[[1]]; m0 = ilist[[2]];
	a DiscreteDelta[m-m0] +
		b(DiscreteStep[m-(m0+1)] + DiscreteStep[-m+(m0-1)])   ),
    _,
	$Failed
  ]
 ]

(***************************** DiscreteStep ********************************)

DiscreteStep[x__] := 0 /; Apply[Or,
	 Map[(  (IntegerQ[#] && # < 0) ||
		(# === DirectedInfinity[-1]) ||
		(NumberQ[N[#]] && !IntegerQ[#])  )&, {x}]]

(* Some nonnegative arguments. *)

DiscreteStep[x__] :=
  Module[{select, comp},
    (comp = Complement[{x}, select];
     If[comp === {}, 1, Apply[DiscreteStep, comp]]
    ) /; FreeQ[{x}, Complex] &&
	 Length[select = Select[{x},
		 ((IntegerQ[#] && # >= 0) || # === DirectedInfinity[1])&]] > 0
  ]

(* Products of DiscreteSteps. *)

DiscreteStep /: DiscreteStep[x_ + aa_.]*DiscreteStep[-x_ + bb_.] :=
  Module[{},
        0
  ] /; FreeQ[{aa, bb}, x] && -aa > bb


(************************ sums, simplifying summand *************************)

(* Note that Sum[g[n] v[n] + h[n] v[n], {n, n1, n2}] ->
	Sum[(g[n] + h[n]) v[n], {n, n1, n2}] is NOT built-in to Mathematica. *)
	
Unprotect[Sum]

  (* collect subexpressions before summing *)
  (* NOTE: need DSKDFreeQ[HoldPattern[f]] rather than DSKDFreeQ[f]
	because otherwise there are problems with expressions like
	Sum[g[[i]], {i, 3}] where g = {g1, g2, g3}.
  *)
Sum[f_, x__?IndexListQ, opts___?OptionQ] :=
  Module[{result},
    (
    Sum[Evaluate[result], x, opts]
    ) /; (If[!FreeQ[result = SimplifyDiscreteStep[f], DiscreteStep],
		result = f]; 
	 result = collect[result, {DiscreteStep[__], DiscreteDelta[__]}];
	 !SameQ[f, result])
  ] /; ( !DSKDFreeQ[HoldPattern[f]])

  (* Take factors not involving the summing variables outside of the sum;
	this is automatic with Integrate, but Integrate is not HoldAll. *)
  (* NOTE: need DSKDFreeQ[HoldPattern[f1 fx]] rather than DSKDFreeQ[{f1, fx}]
	because otherwise there are problems with expressions like
	Sum[k g[[i]], {i, 3}] where g = {g1, g2, g3}.
  *)
Sum[f1_ fx_ , x__?IndexListQ, opts___?OptionQ] :=
  Module[{symbolList, head, sum},
    (
    f1 Sum[Evaluate[fx], x, opts]
    ) /; (symbolList = Map[First, {x}];
	  FreeListQ[f1, symbolList] &&
	  NotFreeListQ[HoldPattern[fx], symbolList]) 
  ] /; ( !DSKDFreeQ[HoldPattern[f1 fx]])


(******************************* sums, linearity *****************************)

(* NOTE: this rule simplifies
	Sum[-DiscreteStep[-2 + x] + DiscreteStep[x], {x, -Infinity, Infinity}]
*)
Sum[k1_. DiscreteStep[x_ + b1_.] + k2_. DiscreteStep[x_ + b2_.],
	{x_, -Infinity, Infinity}, opts___?OptionQ] :=
   Module[{},
	If[b1 > b2,
	   Sum[k1, {x, -b1, -b2-1}, opts],
	   Sum[k2, {x, -b2, -b1-1}, opts]
	]
   ] /; k1+k2==0 && IntegerQ[b1] && IntegerQ[b2] && b1 != b2


  (* NOTE: need DSKDFreeQ[HoldPattern[sum1]] rather than DSKDFreeQ[sum1]
	because otherwise there are problems with expressions like
	Sum[f[[i]] + g[[i]], {i, 3}] where f = {f1, f2, f3} and
	 g = {g1, g2, g3}.
  *)
Sum[sum1_Plus, x__?((!OptionQ[#])&), opts___?OptionQ] :=
   Module[{inlist, successlist = {}, parsuccesslist = {}, failedlist = {},
           scan, i1, sum},
     (
     SimplifyDiscreteStep[ Apply[Plus, successlist] ] +
        (Apply[Plus, parsuccesslist] /. {sum :> Sum}) +
        Sum[ Evaluate[ Apply[Plus, failedlist] ], x, opts]
     ) /; (
	   inlist = Apply[List, sum1];
           scan = Scan[(
		        If[FreeQ[i1 = (Sum[#, x, opts] /. Sum:>sum), sum],
                           (* Sum evaluated completely *)
                           If[i1 === Indeterminate ||
			      Head[i1] === DirectedInfinity,
                              Return[$Failed] ];
                           (* Sum converged *)
                           AppendTo[successlist, i1],
                           (* Sum did not evaluate completely *)
                           If[Head[i1] === sum &&
                              Rest[Apply[List, i1]] =!= {x, opts},
                              (* Sum limits were trimmed so Sum
                                 partially succesful *)
                              AppendTo[parsuccesslist, i1],
                              AppendTo[failedlist, #] ]
                        ])&, inlist];
           (* evaluate if none of the sums fail to converge *)
           scan =!= $Failed)
   ] /; ( !DSKDFreeQ[HoldPattern[sum1]])

  (* constant multiplying DiscreteStep *)
Sum[a_ DiscreteStep[y_], x_?IndexListQ, opts___?OptionQ] :=
        a Sum[DiscreteStep[y], x, opts]  /; (  FreeQ[a, x[[1]]])

  (* constant and nonconstant multiplying DiscreteStep *)
Sum[a_ f_ DiscreteStep[y_], x__?IndexListQ, opts___?OptionQ] :=
       a Sum[Evaluate[f DiscreteStep[y]], x, opts]  /;
  With[{slist = Map[First, {x}]},
   Apply[And, Map[FreeQ[a, #]&, slist]] && Apply[Or, Map[!FreeQ[f, #]&, slist]]
  ]

  (* integral of DiscreteDelta or DiscreteStep terms *)
  (* NOTE: need DSKDFreeQ[HoldPattern[f]] rather than DSKDFreeQ[f]
	because otherwise there are problems with expressions like
	Sum[Integrate[f[[i]], t], {i, 3}] where f = {f1[t], f2[t], f3[t]}.
  *)
Sum[a_. Integrate[f_, i__] + b_., x_?IndexListQ, opts___?OptionQ] :=
        Module[{result = Module[{head, fF},
                                head = Head[F = Sum[Evaluate[a f], x, opts]];
                                If[head =!= Sum,
                                    Integrate[fF, i] + Integrate[b, x],
                                    $Failed]]},
                result /; result =!= $Failed
        ] /; (
		!DSKDFreeQ[HoldPattern[f]] &&
		FreeQ[{i}, x[[1]]] && !FreeQ[HoldPattern[f], x[[1]]]
	     )



(******************************* definite sums ******************************)

Sum[f_. If[cond_, a_, b_], z__?IndexListQ, opts___?OptionQ] :=
  Module[{symbolList, ilist, result},
    (
	Sum[Evaluate[f result], z, opts]
    ) /; ( symbolList = Map[First, {z}];
	   ilist = ConditionToIList[cond, symbolList];
	   (ilist =!= $Failed) &&
	   (result = IListToDiscreteStep[ilist, a, b];  result =!= $Failed)
         ) 
  ]

Sum[f_. Sign[n_ + n0_.], z__?IndexListQ, opts___?OptionQ] :=
  Module[{symbolList},
    (
	(* n==n0: 0
	   n > -n0 (n >= -n0+1, n-(-n0+1) >= 0): f
	   n < -n0 (n <= -n0-1, -n0-1-n >= 0): -f *)
	
	Sum[ Evaluate[f (DiscreteStep[n-(-n0+1)] - DiscreteStep[-n+(-n0-1)])],
		 z, opts]
    ) /; ( symbolList = Map[First, {z}];
	   MemberQ[symbolList, n] && FreeQ[n0, n])
  ] 

        (* Use DiscreteStep to trim sum limits. *)
Sum[f1_ DiscreteStep[y1__]^(n1_.), z__?IndexListQ, opts___?OptionQ] :=
  Module[{f, ylist, n,
	  scan, i, imin, imax, di, nimin, nimax, c0, c1, tau, result},
     (
       result
     ) /;
       (* Sum is HoldAll.  Simplify first argument of Sum. *)	
       MatchQ[f1 DiscreteStep[y1]^n1, f2_ DiscreteStep[y2__]^(n2_.) /;
		({f, ylist, n} = {f2, {y2}, n2};  IntegerQ[n2] && n2 > 0)] &&
       (scan = Scan[	
        (
	   Switch[Length[#], 2, {i, imin, imax, di} = {#[[1]], 1, #[[2]], 1},
			    3, {i, imin, imax, di} = Append[#, 1],
			    _, {i, imin, imax, di} = #];
	   {nimin, nimax} = N[{imin, imax}];
	   If[ (NumberQ[nimin] || Head[imin] === DirectedInfinity) &&
               (NumberQ[nimax] || Head[imax] === DirectedInfinity) &&	
	       !(NumberQ[nimin] && NumberQ[nimax]) &&
	       TrueQ[di == 1] &&  (* NOTE: relax this condition on di someday *)
               MemberQ[ylist, p_ /; PolynomialQ[p, i] && Exponent[p, i] == 1 &&
               		( {c0, c1} = CoefficientList[p, i]; (*   c0 + c1*i   *)
                          NumberQ[N[c0]] &&
		 	  TrueQ[Abs[c1] == 1] && (* NOTE: relax this condition
							 on c1 someday. *)
                          (tau = -c0/c1;
		           IntegerQ[tau] &&
			   If[N[c1] > 0, tau > nimin, tau < nimax])
                        )
               ],
	    		Return[True] (* return from Scan, the limits of this
					sum CAN be trimmed *)
	   ]
        )&, {z}];
	TrueQ[scan]
       ) &&
       Module[{newf, newz, i, imin, imax, nimin, nimax, c0, c1, nc1, tau, f4,
	       div, head},
         newf = ( SimplifyDiscreteStep[f *
		 	(DiscreteStep @@ ylist)^n] ) /. ToUniRules;
	 (* find newz representing the trimmed limits *)
         newz =
           Map[( Switch[Length[#], 2, {i, imin, imax} = {#[[1]], 1, #[[2]]},
		               3, {i, imin, imax} = #,
			       _, {i, imin, imax} = Take[#, 3]	];
		 {nimin, nimax} = N[{imin, imax}];
                 While[MatchQ[newf,
		   	 f3_ DiscreteStep[p_]^(m_.) /;
			 PolynomialQ[p, i] && Exponent[p, i] == 1 &&
			 ( {c0, c1} = CoefficientList[p, i];
                           NumberQ[N[c0]] &&
			   TrueQ[Abs[c1] == 1] ) && (* NOTE: relax this 
						    condition on c1 someday *)
                         ( {nc1, tau, f4} = {N[c1], -c0/c1, f3};
			   IntegerQ[tau] &&
                           If[nc1 > 0, tau > nimin, tau < nimax] )
		       ],
                         If[nc1 > 0,
                            If[nimax <= tau, newf = 0,
                                     newf = f4; imin = tau; nimin = N[imin] ],
                            If[nimin >= tau, newf = 0,
                                       newf = f4; imax = tau; nimax = N[imax] ]
			 ]
                 ];
	         {i, imin, imax}
	       )&,			{z}];
         div = (Head[Sum::div] === $Off);
         Off[Sum::div];
         If[newf === 0,
            result = 0,
	    newf = newf /. ToMultiRules;
            head = Head[result = Apply[Sum, Join[{newf}, newz, {opts}]]];
            If[ head === Indeterminate || head === DirectedInfinity,
                 If[!div, On[Sum::div]];  Return[False]
            ]
         ]; (* result OK *)
         If[!div, On[Sum::div]];
         True
       ] (* end Module *)

  ] /; IntegerQ[n1] && n1 > 0	(* end Module *)


        (* Eliminate DiscreteStep terms that are unity. *)
Sum[f1_. DiscreteStep[y1__]^(n1_.), z__?IndexListQ, opts___?OptionQ] :=
  Module[{f, ylist, n,
	  scan, x, a, b, arg, unit, newunit, symbolList},
    (
    newunit = (DiscreteStep @@ ylist)^n;
    Scan[({x, a, b} = #;
            If[(unit = (newunit //. {x -> a})) === (newunit //. {x -> b}),
               newunit = unit]
         )&, {z}];
    symbolList = Map[First, {z}];
    If[Apply[And, Map[FreeQ[newunit, #]&, symbolList]],
       newunit Sum[Evaluate[f], z, opts],
       Sum[Evaluate[f newunit], z, opts]
    ]
    ) /;
	(* Sum is HoldAll.  Simplify first argument of Sum. *)
	MatchQ[f1 DiscreteStep[y1]^n1, f2_ DiscreteStep[y2__]^(n2_.) /;
                ({f, ylist, n} = {f2, {y2}, n2};  IntegerQ[n2] && n2 > 0)] &&
        (scan = Scan[
         ( MatchQ[#, {x1_Symbol, a1_, b1_} /;
                ( {x, a, b} = {x1, a1, b1};
                  (IntegerQ[a] || Head[a] === DirectedInfinity) &&
                  (IntegerQ[b] || Head[b] === DirectedInfinity)
                )
           ] &&
	   (Length[arg = Select[ylist, !FreeQ[#, x]&]] == 1) &&
	   PolynomialQ[arg[[1]], x] && Exponent[arg[[1]], x] == 1 &&
	   Apply[And, Map[IntegerQ, CoefficientList[arg[[1]], x] ]] &&
           (  ((DiscreteStep @@ ylist)^n) /. {x -> a}  ) ===
	   (  ((DiscreteStep @@ ylist)^n) /. {x -> b}  )
         )&, {z}];
	TrueQ[scan])
  ]


        (* Expand if it appears that the DiscreteStep will allow the
                summation limits to be trimmed. *)
Sum[f_ (DiscreteStep[y__]^(n_.) + h_) + g_., z__?IndexListQ, opts___?OptionQ] :=
  Module[{x, a, b, c, d, t},
    (
    Sum[
	Evaluate[
	  Expand[f (DiscreteStep[y]^n + h), DiscreteStep[__]] + g
        ], z, opts]
    ) /;
       Apply[Or, Map[
        ( MatchQ[#, {x1_Symbol, a1_, b1_} /;
               ( {x, a, b} = {x1, a1, b1};
                 (IntegerQ[a] || Head[a] === DirectedInfinity) &&
                 (IntegerQ[b] || Head[b] === DirectedInfinity)
               )
          ] &&
          MemberQ[{y}, p_ /; PolynomialQ[p, x] && Exponent[p, x] == 1 &&
               ( {d, c} = CoefficientList[p, x];
                 t = -d/c; If[N[c] > 0, N[t] > N[a], N[t] < N[b]] )
          ]
        )&,
        {z}]]
  ]


(****************************** multiple sums *******************************)

        (* Nest multiple sums. *)
  (* NOTE: need DSKDFreeQ[HoldPattern[f]] rather than DSKDFreeQ[f]
	because otherwise there are problems with expressions like
	Sum[g[[i]] h[[j]], {i, 3}, {j, 3}] where g = {g1, g2, g3} and
	h = {h1, h2, h3}. *)
Sum[f_, iseq__, opts___?OptionQ] :=
  Module[{result, ilist, sum},
     (
     result = result /. {Sum :> sum};
     If[Head[result] === sum,
	ilist = {};
	While[Head[result] === sum,
	      ilist = Join[ilist, Rest[Apply[List, result]]];
              result = result[[1]] ];
        result = Apply[Sum, Prepend[ilist, result]]
     ];
     result /. {sum :> Sum}
     )  /;  Module[{summand = f, scan, oldsummand, arglist, head},
             scan = Scan[(
			  oldsummand = summand /. sum :> Sum;
			  arglist = Join[{oldsummand}, {#}, {opts}];
			  head = Head[summand =
				 (Apply[Sum, arglist] /. Sum :> sum)];
                          If[!(head =!= sum ||
		               summand[[1]] =!= oldsummand ||	
			       summand[[2]] =!= #),
                                Return[$Failed]
			  ])&,
                         Reverse[{iseq}]];
	     result = summand;
	     scan =!= $Failed] 
  ]  /;  (  !DSKDFreeQ[HoldPattern[f]] && Length[{iseq}] > 1)




(* NOTE:  Indefinite integrals are supported in DiracDelta.m, but there is
	no similar indefinite sum in Mathematica, so no support in
	KroneckerDelta.m. *)

(********** definite sums involving KroneckerDelta and DiscreteStep ********)

(* NOTE: the complexity of these rules is due to trying to express the
	filtering effect of DiscreteStep[x] (nonzero only for integer x >= 0)
	for abstract limits of Sum.  Not all possible cases are covered...
	some Sums that should evaluate (i.e., simplify to an expression
	with no DiscreteStep) do not because the rules are missing. *)


(* Sum of a univariate step function with linear argument. *)

Sum[f_. DiscreteStep[y_], x_?IndexListQ, opts___?OptionQ] :=
  Module[{i, di, k0, k1, n1, d1, k, ilist, arg, result, sum},
    (
	(* the following is possible if imin and imax are integer-valued *)
	If[Sign[di] == -1, {imin, imax, di} = {imax, imin, -di}];
	f1 = f /. {i->k};  f1 = f1 /. {k->d1 i};
	ilist = Which[
		 n1 == d1 == di == 1, (* di == 1 *)
			(
			arg = imax-Max[imin, -k0];
		       {i, Max[imin, -k0], imax}),
		n1 == d1 == 1,	     (* di is a positive integer > 1 *)
			(
			arg = imax-(imin + Max[0, di Ceiling[(-k0-imin)/di]]);
		       {i, imin + Max[0, di Ceiling[(-k0-imin)/di]], imax, di}),
		n1 == -1 && d1 == di == 1, (* di = 1 *)
			(
			arg = Min[imax, k0]-imin;
		       {i, imin, Min[imax, k0]}),
		n1 == -1 && d1 == 1, (* di is a positive integer > 1 *)
			(
			arg = Min[imax, k0]-imin;
		       {i, imin, Min[imax, k0], di}),
		positive[n1] && d1 == di == 1,
			(
			arg = imax-Max[imin, Ceiling[-k0/n1]];
		       {i, Max[imin, Ceiling[-k0/n1]], imax}),
		negative[n1] && d1 == di == 1,
			(
			arg = Min[imax, -k0/n1]-imin;
	               {i, imin, Min[imax, -k0/n1]}),
		positive[n1] && d1 != 1 && di == 1,
			(
			arg = imax/d1 - Max[Ceiling[imin/d1], Ceiling[-k0/n1]];
	               {i, Max[Ceiling[imin/d1], Ceiling[-k0/n1]], imax/d1}),
		negative[n1] && d1 != 1 && di == 1,  
			(
			arg = Min[imax/d1, -k0/n1] - Ceiling[imin/d1];
		       {i, Ceiling[imin/d1], Min[imax/d1, -k0/n1]})   ];
	result = DiscreteStep[arg] *
		Sum[Evaluate[f1], Evaluate[ilist], Evaluate[opts]];
	result = result /. {Sum :> sum};
        If[Head[result] === sum,
	   ilist = {};
	   While[Head[result] === sum,
	      ilist = Join[ilist, Rest[Apply[List, result]]];
              result = result[[1]] ];
           result = Apply[Sum, Prepend[ilist, result]]
        ];
        result /. {sum :> Sum}
    ) /; (Switch[Length[x], 2, {i, imin, imax, di} = {x[[1]], 1, x[[2]], 1},
			    3, {i, imin, imax, di} = Append[x, 1],
			    _, {i, imin, imax, di} = x];
	  PolynomialQ[y, i] && Exponent[y, i] == 1) &&
	 (* assume imin and imax are integer-valued *)
	 (!NumberQ[N[imin]] || IntegerQ[imin]) &&
	 (!NumberQ[N[imax]] || IntegerQ[imax]) &&
	 ({k0, k1} = CoefficientList[y, i];
	  {n1, d1} = {Numerator[k1], Denominator[k1]};
	  negative[n1] || positive[n1]) && (!NumberQ[N[k0]] || IntegerQ[k0]) &&
	 (IntegerQ[di] && Abs[di] > 0) && (Abs[di] == 1 || d1 == 1)
	 (* Note: the preceeding requirement that either Abs[di] == 1 
		  or d1 == 1 is not really required;  it is just an
		  attempt to limit the scope of this rule for now! *)
  ] /; !Apply[And, Map[NumberQ, N[Drop[x, 1]] ]] 

Sum[f_. DiscreteStep[y_], x_?IndexListQ, opts___?OptionQ] :=
  Module[{i, imin, imax, di, k0, k1, n1, d1, result},
    (
	result
    ) /; (Switch[Length[x], 2, {i, imin, imax, di} = {x[[1]], 1, x[[2]], 1},
			    3, {i, imin, imax, di} = Append[x, 1],
			    _, {i, imin, imax, di} = x];
	  PolynomialQ[y, i] && Exponent[y, i] == 1) &&
	 (* assume imin and imax are integer-valued *)
	 (!NumberQ[N[imin]] || IntegerQ[imin]) &&
	 (!NumberQ[N[imax]] || IntegerQ[imax]) &&
	 ({k0, k1} = CoefficientList[y, i];
	  {n1, d1} = {Numerator[k1], Denominator[k1]};
	  negative[n1] || positive[n1]) && (NumberQ[N[k0]] && !IntegerQ[k0]) &&
	 (IntegerQ[di] && Abs[di] > 0) && (Abs[di] == 1 || d1 == 1) &&
	 (* Note: the preceeding requirement that either Abs[di] == 1 
		  or d1 == 1 is not really required;  it is just an
		  attempt to limit the scope of this rule for now! *)
	 (result = sum1[f, {k0, k1}, {n1, d1}, {i, imin, imax, di}, opts]) =!=
		 $Failed
  ] /; !Apply[And, Map[NumberQ, N[Drop[x, 1]] ]] 

sum1[f_, {k0_, k1_}, {n1_, d1_}, {i_, imin_, imax_, di_}, opts___] :=
  Module[{n0, d0, q},
	(* the following is possible if imin and imax are integer-valued *)
	If[Sign[di] == -1, {imin, imax, di} = {imax, imin, -di}];
	If[ d1 == 1,
	    0, (* (k0 + k1 i) will never be integer-valued *)
	    (* d1 != 1 *)
	    {n0, d0} = {Numerator[k0], Denominator[k0]};
	    If[ d0 == 1,
	        0, (* a numerical k0 implies that that (k0 + k1 i) will never
			 be integer-valued *)
	       If[ Mod[d1, d0] != 0,
		   0, (* (n0/d0 + n1/d1 i) will only be integer-valued if
				d0 divides d1 *)
		   q = Quotient[d1, d0];
		   $Failed
	       ] (* If Mod[d1, d0] != 0 *)
	    ] (* If d0 == 1 *) 
	] (* If d1 == 1 *)
  ]

(* Sum of a univariate delta function with arbitrary argument.  Delta's
   located off the real line contribute 0. *)

Sum[f_. DiscreteDelta[y_], x_?IndexListQ, opts___?OptionQ] :=
  Module[{solution, select, infy, flist, i, imin, imax, di, singularitylist,
		which, j},
    (
	(* the following is possible if imin and imax are integer-valued *)
	If[Sign[di] == -1, {imin, imax, di} = {imax, imin, -di}];
 	singularitylist = i /. select;	
	pairs = Transpose[{singularitylist, flist}];
	result =
	 If[ di == 1,
	       Apply[Plus, Map[ (#[[2]] *
		 which[ imin === -Infinity && imax === Infinity, 1,
			imin === -Infinity, DiscreteStep[imax-#[[1]]],
			imax === Infinity, DiscreteStep[#[[1]]-imin],
			True, DiscreteStep[#[[1]]-imin] -
					DiscreteStep[#[[1]]-(imax+1)] ])&,
				 pairs]],
	  (* di != 1 *)
		(* NOTE: check the case of di != 1 *)
	       Apply[Plus, Map[
	         (j = (#[[1]]-imin)/di;
		  If[Evaluate[ !NumberQ[j] || IntegerQ[j] ],
		     Evaluate[ (#[[2]] *
		        which[ imin === -Infinity && imax === Infinity, 1,
			       imin === -Infinity, DiscreteStep[imax-#[[1]]],
			       imax === Infinity, DiscreteStep[#[[1]]-imin],
			       True, DiscreteStep[#[[1]]-imin] -
					DiscreteStep[#[[1]]-(imax+1)] ]) ],
		     0])&,
	         		pairs]]
	 ];
	result = result /. {which :> Which};
	SimplifyDiscreteStep[result]
    ) /; (Switch[Length[x], 2, {i, imin, imax, di} = {x[[1]], 1, x[[2]], 1},
                            3, {i, imin, imax, di} = Append[x, 1],
                            _, {i, imin, imax, di} = x];
          (* assume imin and imax are integer-valued *)
          (!NumberQ[N[imin]] || IntegerQ[imin]) &&
          (!NumberQ[N[imax]] || IntegerQ[imax])) &&
   (
    Head[solution = Solve[Map[(# == 0)&, {y}], i]] == List) &&
          solution =!= {{}} && If[FreeQ[y, x], True, solution =!= {}] &&
   (select = Select[solution, (!NumberQ[#] || IntegerQ[#])&];
    infy = (Head[Power::infy] === $Off);
    Off[Power::infy];  flist = f /. select;  If[!infy, On[Power::infy]];
    FreeQ[N[flist], DirectedInfinity] && FreeQ[flist, Indeterminate]) 
  ] /; !Apply[And, Map[NumberQ, N[Drop[x, 1]] ]]


(* Sum of a multivariate step function with linear arguments. *)

Sum[f_. DiscreteStep[yseq__]^(n_.), x_?IndexListQ, opts___?OptionQ] :=
  Module[{ylist, arg},
        (
      argc = Select[ylist, FreeQ[#, First[x]]&];
      If[argc === {}, 1, Apply[DiscreteStep, argc]] *
                Sum[Evaluate[f Apply[Times, Map[DiscreteStep, arg]]], x, opts]
        )
       /;  Length[ylist = {yseq}] > 1 &&
	   (arg = Select[ylist, (!FreeQ[#, First[x]])&];
            Length[arg] >= 1) &&
	    Apply[And, Map[PolynomialQ, arg]] &&
	    Apply[And, Map[(Exponent[#, First[x]] == 1)&, arg]]
  ] (* NOTE: /; FreeQ[f, Sum] *)


(* Sum of a multivariate delta function with arbitrary arguments. *)

Sum[f_. DiscreteDelta[yseq__], x_?IndexListQ, opts___?OptionQ] :=
  Module[{ylist, arg},
      (
      Apply[DiscreteDelta, Select[ylist, FreeQ[#, First[x]]&]] *
                Sum[Evaluate[f Apply[DiscreteDelta, arg]], x, opts]
      ) /;  Length[ylist = {yseq}] > 1 &&
           Length[arg = Select[ylist, (!FreeQ[#, First[x]])&]] == 1
  ] /; FreeQ[f, Sum]



Protect[Sum]

(********************************** Simplify ********************************)

SimplifyDiscreteStep[expr_] :=
   Module[{new = expr //. ToUniRules, originalLength = Length[cases[expr]],
		convertTomulti, sum, c, collected, newLength},
          convertTomulti = originalLength < Length[cases[new]];
          new = simplifyExpandedDS[ collect[new, {DiscreteStep[__]}] ];
          new = new //. Sum :> sum;
          new = new //.
          	{(c0_. + c1_. DiscreteStep[p__]^m_.)*
	   	(d0_. + d1_. DiscreteStep[q__]^n_.) :>
           	c0 d0 + c0 d1 DiscreteStep[q]^n + c1 d0 DiscreteStep[p]^m +
                c1 d1 DiscreteStep[p]^m DiscreteStep[q]^n};
          new = new //. sum :> Sum;
          new = new //. DiscreteStepRules;
          new = simplifyExpandedDS[new];
	  c = cases[new];
	  collected = Collect[new, Union[c]];
 	  If[Length[cases[ collected ]] < Length[c], new = collected];
          If[convertTomulti, new = new //. ToMultiRules];
          newLength = Length[cases[new]];
	  Which[newLength < originalLength, new,
		newLength == originalLength,
		  If[LeafCount[new] < LeafCount[expr], new, expr],
		True, expr]	
   ] /; !DSFreeQ[expr] && Head[expr] =!= Sum

SimplifyDiscreteStep[expr_] := expr


simplifyExpandedDS[a_?(DSFreeQ[#]&) +
        (b:(_?(DSFreeQ[#1]&) DiscreteStep[__] | DiscreteStep[__])..)] :=
        Simplify[a] + Plus[b]

simplifyExpandedDS[expr_] := expr

cases[expr_] := Join[Cases[{expr}, DiscreteStep[__], Infinity],
                Cases[{expr}, DiscreteDelta[__], Infinity]]



  (* reducing some DiscreteStep differences to KroneckerDelta *)
DiscreteStepRules = {
	an_. DiscreteStep[n_] + am_. DiscreteStep[m_] :>
		 an DiscreteDelta[n] /; n-m==1 && an+am==0,
	an_. DiscreteStep[n_] + am_. DiscreteStep[m_] + an_. DiscreteStep[o_] :>
          	an(DiscreteDelta[n] - DiscreteDelta[m]) /;
		n-m==1 && n-o==2 && 2an + am==0
		}

ToUniRules = {DiscreteStep[x__] :>
                Apply[Times, Map[DiscreteStep[#]&, {x}]] /; Length[{x}] > 1,
              DiscreteDelta[x__] :>
                Apply[Times, Map[DiscreteDelta, {x}]] /; Length[{x}] > 1}

ToMultiRules = {DiscreteStep[x__] * DiscreteStep[y__] :>
                        Apply[DiscreteStep, Join[{x}, {y}]],
                DiscreteDelta[x__] * DiscreteDelta[y__] :>
                        Apply[DiscreteDelta, Join[{x}, {y}]],
                DiscreteStep[x_]^n_Integer?Positive :>
                        Apply[DiscreteStep[#]&, Table[x, {n}]],
                DiscreteDelta[x_]^n_Integer?Positive :>
                        Apply[DiscreteDelta, Table[x, {n}]]}



(***************************** specialized Q fcns ****************************)

FreeListQ[expr_, l_List] := Apply[And, Map[FreeQ[expr, #]&, l]]
NotFreeListQ[expr_, l_List] := Apply[Or, Map[(!(FreeQ[expr, #]))&, l]]
numberQ[n_] := NumberQ[N[n]] || Head[n] === DirectedInfinity
DSFreeQ[expr_] := ( FreeQ[expr, DiscreteStep]) 
KDFreeQ[expr_] := ( FreeQ[expr, DiscreteDelta])
DSKDFreeQ[expr_] := ( DSFreeQ[expr] && KDFreeQ[expr])
LimitFreeQ[expr_] := FreeQ[expr, Limit]
IndexListQ[list_] :=
	 FreeQ[list, Complex] && VectorQ[list] && 
	 2 <= Length[list] <= 4 &&
	 Head[list[[1]]]===Symbol (* if this is true then OptionQ[list] is
					 False *)


(****************************** collect utility *****************************)

(* check for a pattern that doesn't need collecting *)
collect[c0_. + c1_. x_, pattern_List] :=
   (
    c0 + c1 x
   ) /; DSKDFreeQ[{c0, c1}] &&
        (MemberQ[pattern, DiscreteStep[__]] ||
         MemberQ[pattern, DiscreteDelta[__]]) &&
        (MatchQ[x, DiscreteStep[__] | DiscreteStep[__]^_Integer?Positive |
                  DiscreteDelta[__]] ||
         MatchQ[x, HoldPattern[Times[
          (DiscreteStep[__] | DiscreteStep[__]^_Integer?Positive |
           DiscreteDelta[__]).. ]]])

(* only Expand terms involving patterns in the specified list *)
collect[(expr_)^n_Integer?Positive, pattern_List] :=
  Module[{cases, varlist, output},
        cases = Apply[Union, Map[Cases[{expr}, #, Infinity]&, pattern]];
        varlist = Flatten[Map[If[PolynomialQ[expr, #], #, {}]&, cases]];
        output = Expand[(Collect[expr, varlist])^n, First[pattern] ];
        Scan[(output = Expand[output, #])&, Rest[pattern]];
        output
  ]

collect[expr_, pattern_List] :=
  Module[{cases, varlist, output},
        cases = Apply[Union, Map[Cases[{expr}, #, Infinity]&, pattern]];
        varlist = Flatten[Map[If[PolynomialQ[expr, #], #, {}]&, cases]];
        output = Expand[Collect[expr, varlist], First[pattern]];
        Scan[(output = Expand[output, #])&, Rest[pattern]];
        output
  ]


(* NOTE:  the section of rules "positiveIntervals" from
	   DiracDelta.m is not relevant here
	  (...will anyone want the DiscreteStep of a polynomial?) *)

(************************** miscellaneous definitions ************************)

DiscreteStep /: DiscreteStep[z__] *
  DiscreteDelta[x___, a_. y_Symbol + b_.] :=
      (
        (DiscreteStep[z] /. y -> -b/a) * DiscreteDelta[x, a y + b]
      )		/; !FreeQ[{z}, y] && NumberQ[N[b]] && TrueQ[Abs[a]==1]

DiscreteStep /: DiscreteStep[z__]^n_?IntegerQ :=
	DiscreteStep[z] /; n > 0

DiscreteStep /:
  DiscreteStep[a_. y_Symbol + b_.] DiscreteStep[a_. y_Symbol + c_.] :=
   (
	DiscreteStep[y + Min[b, c]]
   ) /; Abs[a] == 1 && FreeQ[{b, c}, y]

DiscreteStep /:
  DiscreteStep[z1___, a_. y_Symbol + b_., z2___, a_. y_Symbol + c_., z3___] :=
   (
	DiscreteStep[z1, y + Min[b, c], z2, z3]
   ) /; Abs[a] == 1 && FreeQ[{b, c}, y]


(*****************************************************************************)

End[]

EndPackage[]


