Functions

Brett van de Sande
September 1991

This package defines extensions to ReplaceAll, ReplaceRepeated, Dot, **, inequalities, Sign and Integrate using Big.  It also defines new functions NormalOrder, ScalarQ, CnumberQ, TestValue, trace, sum and ipoly.  It defines associated functions used in manipulating them:  delta, convert, showpoly, checkpoly, transform, translate, rotate, blockdiagonal, LinearQ, and coefficients.  It defines lists of simplification rules:  Dotexpand, timesexpand, sumexpand, ipolyexpand, sumfactor, and ipolyfactor.

BeginPackage["Functions`"];

Define usage of various functions

sum::usage="sum[f[x1,x2,I,xn],{x1,I},{x2,I},I,{xn,I}] is the
sum over x1, x2, I, xn in that order; the iteration form is
the same as that of Sum.  The order of summation is oppisite
that of Sum.";

ipoly::usage = "ipoly[f[x1,x2,I,xn],{x1,x2,I,xn},{{a11,a12,I,a1n},{a21,a22,I,a2n},I,{aJ1,aJ2,I,aJn}},{b1,b2,I,bJ}]\nis
the n-dimensional integral of f[I] over a finite volume bounded
by an n-dimensional convex polytope P.  P is defined to be
all points which satisfy the J inequalities:\naj1 x1 + aj2
x2 + I + ajn xn <= bj,	1<=j<=J.\nInput form is
ipoly[f,x,{c1<=c2,c3<=c4,I}] where ci's are linear in x.";

NormalOrder::usage="NormalOrder[a**b] (or a.b) is the normal 
order of a**b.  Conjugate must be explicit.  It is expanded by
timesexpand.";

Big::usage = "Big has a TestValue of 1,000,000. 
Integrate recognizes Big as a limit for producing
delta functions.";
	
convert::usage="convert[expr,f,{v1,v2,I}] converts variables
number v1, v2, I of sum or ipoly to function f using a
syntax compatable with f=Sum or f=Integrate.  expr is
linear in sum or ipoly; convert[expr,f] converts all
variables.";

showpoly::usage="showpoly[expr] Prints the polytope for each 
ipoly in expr.";

transform::usage="transform[ipoly[I],{q1,q2,I}] where qi are 
expressions linear in the integration variables xi, 
defines the new value of each variable to be xi->qi,
i=1,2,I .";

translate::usage="translate[ipoly[I],t] transforms the 
coordinate vector as x->x+t.";

rotate::usage="rotate[ipoly[I],r] transforms the coordinate 
vector as x->r.x.";

Dotexpand::usage="Dotexpand is a List of rules which
expands Dot. It also expands trace.";

timesexpand::usage="timesexpand is a List of rules which
expands NonCommutativeMultiply.  It also expands 
NormalOrder.";

sumexpand::usage="sumexpand is a List of rules which
attempts to move sum to the inside of an expression,
and remove delta functions.";

ipolyexpand::usage="ipolyexpand is a List of rules which
attempts to move ipoly to the inside of an expression,
and remove delta functions.";

sumfactor::usage="sumfactor is a List of rules which attempts
to move sum to the outside of an expression, and combines 
various products of sum.  It is often faster to use expr
//.(sumfactor,sum).";

checkpoly::usage="checkpoly[ipoly[f,x,a,b]] checks that the 
matrix equation a.x<=b defines a valid convex polytope.  It
can take factorial time for complicated poytopes.  It
returns zero or a valid polytope with unneeded hyperplanes 
removed.";

TestValue::usage="TestValue[x] is used to to give a
dummy numerical value to x for use in evaluating
inequalities.  One defines TestValue for various
symbols used in inequalities.";

ipolyfactor::usage="ipolyfactor is a List of rules which 
attempts to move ipoly to the outside of an expression, and 
combines various products of ipoly.  It is often faster to
use expr//.(sumfactor,sum).";

delta::usage="delta[x] is the Dirac delta-function.  
delta[x,y] is the Kroneker delta function.";

coefficients::usage="coefficients[expr,x] where expr is 
linear in the variables in the List x, gives a List of 
each coefficient and the constant term.";

trace::usage="trace[x] takes the trace of matrix x.  It is
 expanded by Dotexpand.";

LinearQ::usage="LinearQ[expr,x] checks if expr is linear 
in the variables in the List x.";

CnumberQ::usage="CnumberQ[x] returns True if x is a
c-number and False otherwise.  One defines CnumberQ to 
be False for non-commuting objects";

ScalarQ::usage="ScalarQ[x] returns True if x is not
a matrix or vector and False otherwise.  One defines
ScalarQ to be False for nonscalar objects.";

Begin["Functions`internal`"];

Extensions to built-in functions

Extensions to ReplaceAll and ReplaceRepeated

expr/.(rule,g) and expr//.(rule,g) are extensions of ReplaceAll and ReplaceRepeated which apply the rule or list of rules only to subexpressions which contain the object g (g may also be a list of objects) as a subexpression.  This is completely anologous to UpSet and is useful to increase speed.

Unprotect[ReplaceAll,ReplaceRepeated];
ReplaceRepeated[e_,r_,x_]:=
	FixedPoint[ReplaceAll[#,r,x]&,e];
ReplaceAll[e_,r_,x_?AtomQ]:=ReplaceAll[e,r,{x}];
ReplaceAll[e_,r_,{x__?AtomQ}]:=
	MapAt[
		Replace[#,r]&,
		e,
		Union[
			Map[
				If[Last[#]==0,Drop[#,-2],Drop[#,-1]]&,
				Select[
					Apply[Join,Map[Position[e,#]&,{x}]],
					(#=!={}&&#=!={0})&
				]
			]
		]
	];
Protect[ReplaceAll,ReplaceRepeated];

Handle inequalities using TestValue.  Define extensions to Positive and Negative and Sign.

Unprotect[LessEqual,GreaterEqual,Less,Greater,Sign];
Sign[x_]:=
	Block[
		{xx},
		Sign[xx]/;
		Not[NumberQ[x]]&&NumberQ[xx=TestValue[x]]
	];
LessEqual[x__]:=
	Block[
		{xx},
		Apply[LessEqual,xx]/;
		Not[Apply[And,Map[NumberQ,{x}]]]&&
		(Apply[And,Map[NumberQ,xx=Map[TestValue,{x}]]])
	];
Less[x__]:=
	Block[
		{xx},
		Apply[Less,xx]/;
		Not[Apply[And,Map[NumberQ,{x}]]]&&
		(Apply[And,Map[NumberQ,xx=Map[TestValue,{x}]]])
	];
GreaterEqual[x__]:=
	Block[
		{xx},
		Apply[GreaterEqual,xx]/;
		Not[Apply[And,Map[NumberQ,{x}]]]&&
		(Apply[And,Map[NumberQ,xx=Map[TestValue,{x}]]])
	];
Greater[x__]:=
	Block[
		{xx},
		Apply[Greater,xx]/;
		Not[Apply[And,Map[NumberQ,{x}]]]&&
		(Apply[And,Map[NumberQ,xx=Map[TestValue,{x}]]])
	];
 Protect[LessEqual,GreaterEqual,Less,Greater,Sign];

Extensions Conjugate[I] which takes the Hermetian conjugate.

Unprotect[Conjugate];
Format[Conjugate[x_?AtomQ]]:=
	ColumnForm[
		{Array[
			"_"&,{StringLength[ToString[x]]},
			1,StringJoin],x},
		Left,
		Above
];
Format[Conjugate[x_[y___]]]:=
	ColumnForm[
		{Array[
			"_"&,{StringLength[ToString[x]]},
			1,StringJoin],x},
		Left,
		Above
][y];
(* There is some trouble with mixing
 Derivatives and conjugates
Conjugate/:D[Conjugate[x_],y_Symbol]:=
	Conjugate[D[x,y]];*)
Conjugate[E]=E;
Conjugate[Pi]=Pi;
Conjugate[Times[x_,y_]]=Times[Conjugate[x],Conjugate[y]];
Conjugate[Plus[x_,y_]]=Plus[Conjugate[x],Conjugate[y]];
Conjugate[Power[x_,y_]]=Power[Conjugate[x],Conjugate[y]];
Conjugate[Dot[x_,y_]]=Dot[Conjugate[y],Conjugate[x]];
Conjugate[NonCommutativeMultiply[x_,y_]]=
		NonCommutativeMultiply[Conjugate[y],Conjugate[x]];
Conjugate[sum[x_,y__]]:=
	sum[
		Conjugate[x]/.Map[(Conjugate[#[[1]]]:>#[[1]])&,{y}],
		y
	];
Conjugate[ipoly[f_,x_,a_,b_]]:=
	ipoly[
		Conjugate[f]/.Map[(Conjugate[#]:>#)&,x],
		x,
		a,
		b
	];
Conjugate[delta[x__]]=delta[x];
Conjugate[Conjugate[x_]]:=x;
Protect[Conjugate];

Dot

Unprotect[Dot];
Dot[0,_]=0;
Dot[_,0]=0;
Protect[Dot];

Extensions to **.

Unprotect[NonCommutativeMultiply];
Attributes[NonCommutativeMultiply]={Flat,OneIdentity};
Protect[NonCommutativeMultiply];

Additions to Integrate (for one integration variable), producing and integrating over delta functions.

constants[I] is defined near the definitions of ipolyexpand and sumexpand

Big/:TestValue[Big]=1000000;
Unprotect[Integrate];
Integrate[a_ + b_,x__]:=Integrate[a,x]+Integrate[b,x];
Integrate[a_Times ,{x_,xx__}]:=
	Block[
		{const=1,nonconst=1},
		const Integrate[nonconst,{x,xx}]/;
		constants[a,{x}]
	];
Integrate[Exp[f_],{x_,xx__}]:=
	Block[
		{cef=Collect[Expand[f],x]},
		Integrate[Exp[cef],{x,xx}]/;f=!=cef
	];
Integrate[x__]:=
	Apply[Integrate,{x}/.Big->Infinity]/;
	Not[FreeQ[{x},Big]];
Integrate[f_. delta[a_. x_ +b_.], {x_,-Infinity,Infinity}]:=
	(f/.x->-b/a)/If[Positive[a],a,-a]/;FreeQ[{a,b},x];
Integrate[Exp[p_ x_ + b_.],{x_,DirectedInfinity[-1],
	DirectedInfinity[1]}]:=
	Exp[b] 2 Pi delta[Expand[p/I]]/;
	FreeQ[{b,p},x]&&FreeQ[Expand[p/I],Complex];
Protect[Integrate];

Define various functions

sum[f[x1,x2,I,xn],{x1,I},{x2,I},I,{xn,I}] is the sum over x1, x2, I, xn in that order; the iteration form is the same as that of Sum.

sum[sum[f_,x__],y__]:=sum[f,x,y];
sum[f_List,x__]:=Map[sum[#,x]&,f];
sum[0,__]=0;
sum/: D[sum[x_,y__],z_]:=
	sum[D[x,z],y]/;FreeQ[{y},z];
sum/:Integrate[sum[x_,y__],z__]:=
	sum[Integrate[x,z],y]/;
	freelist[{y},Map[First,{z}]];

ipoly[f[x1,x2,I,xn],{x1,x2,I,xn},{{a11,a12,I,a1n},{a21,a22,I,a2n},I,{aJ1,aJ2,I,aJn}},{b1,b2,I,bJ}] is the n-dimensional integral of f[I] over a region bounded by an n-dimensional convex polytope P.  P is defined to be all points which satisfy the J inequalities:	aj1 x1 + aj2 x2 + I + ajn xn <= bj,	1<=j<=J.  InputForm is ipoly[f,x,{c1<=c2,c3<=c4,I}] where ci's are linear in x and the inequalities define the convex polytope.

The most general case is for b linear in y; it is not treated here.

ipoly[ipoly[f_,x_,a_,b_],y_,aa_,bb_]:=
	ipoly[
		f,
		Join[x,y],
		blockdiagonal[a,aa],
		Join[b,bb]
	]/;freelist[{a,b},y];
ipoly/:Integrate[ipoly[x_,y__],z__]:=
	ipoly[Integrate[x,z],y]/;
	freelist[{y},Map[First,{z}]];
ipoly[f_List,x__]:=Map[ipoly[#,x]&,f];
ipoly[0,__]=0;
(* output format *)
ipoly/:Format[Literal[ipoly[f_,x_,a_,b_]]]:=
	ipoly[f,x,".."];
ipoly/: D[ipoly[x_,y__],z_]:=
	ipoly[D[x,z],y]/;FreeQ[{y},z];
(* convert input format to internal format *)
ipoly[f_,x_,q_]:=
	Block[
		{xx=q/.{
			LessEqual[a_,b_]:>{Expand[a],Expand[b]},
			GreaterEqual[b_,a_]:>{Expand[a],Expand[b]}
		}},
		checkpoly[
			ipoly[
				f,
				x,
				Transpose[Drop[#,-1]],
				-#[[-1]]
			]
		]&[Transpose[
			Map[
				(coefficients[#[[1]],x]-
					coefficients[#[[2]],x])&,
				xx
			]
		]]/;
		Apply[
			And,
			Map[
				(Head[#]===List&&
				LinearQ[#[[1]],x]&&
				LinearQ[#[[2]],x])&,
				xx
			]
		]
	]

NormalOrder[a**b] (or a.b) is the normal order of a**b, Conjugate must be explicit.

Format[NormalOrder[a_]]:=
	SequenceForm[":",a,":"];
NormalOrder[NormalOrder[a_]]:=NormalOrder[a];
NormalOrder[a_+b_]:=
	NormalOrder[a]+NormalOrder[b];
NormalOrder[Literal[a_ ** b_]]:=
	NormalOrder[a]**b/;FreeQ[b,Conjugate];
NormalOrder[Literal[Conjugate[a_]**b_]]:=
	Conjugate[a]**NormalOrder[b];
NormalOrder[Literal[a_ ** b_ ** c_]]:=
	If[
		Head[b]===Conjugate,
		NormalOrder[NormalOrder[a**b]**c],
		NormalOrder[a**NormalOrder[b**c]]
	]/;Head[b]=!=NormalOrder

TestValue[x] gives a numerical value for the expression x for use in inequalities.  One defines a TestValue for various symbols used in inequalities.

TestValue[x_?NumberQ]=x;
TestValue[a_ +b_]=TestValue[a]+TestValue[b];
TestValue[a_ b_]=TestValue[a] TestValue[b];
TestValue[a_ ^b_]=TestValue[a]^TestValue[b];

Define trace which takes the trace a square matrix.

trace[x_?MatrixQ]:=Sum[x[[i,i]],{i,1,Length[x]}];
trace[a_+b_]=trace[a]+trace[b];
trace[x_?NumberQ]=x;

convert[expr,f,{v1,v2,I}] converts variables number v1, v2, I of sum or ipoly to function f using a syntax compatable with f=Sum or f=Integrate.  expr is linear in sum or ipoly; convert[expr,f] converts all variables.

convert[a_ + b_ ,c___]:=convert[a,c]+convert[b,c];
convert[a_ b_,c___]:=
	a convert[b,c]/;FreeQ[a,ipoly]&&FreeQ[a,sum];

convert[sum[f_,x__],g_]:=Apply[g[f,##]&,Reverse[{x}]];
convert[sum[f_,x__],g_,y_List]:=
	Block[
		{xin={x}[[y]],
		xout={x}[[Complement[Range[1,Length[{x}]],y]]]},
		If[
			xout=={},
			Apply[g[f,##]&,Reverse[xin]],
			Apply[
				sum[Apply[g[f,##]&,Reverse[xin]],##]&,
				xout
			]
		]/;freelist[xout,Map[First,xin]]
	]

convert[ipoly[f_,x_,a_,b_],g_]:=
	convert[ipoly[f,x,a,b],g,Range[1,Length[x]]];
convert[ipoly[f__],g_,y_List]:=
	Map[
		convert[convert[ipoly[f],g,{First[#]}],g,Rest[#]-1]&,
		Sort[y]
	]/;Length[y]>1;
convert::unknown="``, is unknown.";
convert[ipoly[f_,x_,a_,b_],g_,{y_}]:=
	Block[
		{lower={},upper={},neither={},
			xx,aa,al,bl,au,bu,an,bn},
		Map[
			If[
				a[[#,y]]===0,
				AppendTo[neither,#],
				If[
					Positive[a[[#,y]]],
					AppendTo[upper,#],
					AppendTo[lower,#],
					Message[
						convert::unknown,
						Positive[a[[#,y]]]
					]
				]
			]&,
			Range[1,Length[a]]
		];
		xx=Drop[x,{y,y}];
		aa=Map[Drop[#,{y,y}]&,a];
		au=Map[-aa[[#]]/a[[#,y]]&,upper];
		bu=Map[b[[#]]/a[[#,y]]&,upper];
		al=Map[-aa[[#]]/a[[#,y]]&,lower];
		bl=Map[b[[#]]/a[[#,y]]&,lower];
		an=Map[aa[[#]]&,neither];
		bn=Map[b[[#]]&,neither];
		Sum[
			checkpoly[ipoly[
				g[
					f,
					Release[{
						x[[y]],al[[l]].xx+bl[[l]],
						au[[u]].xx+bu[[u]]
					}]
				],
				xx,
				Join[
					an,
					Map[-#+au[[u]]&,Drop[au,{u,u}]],
					Map[#-al[[l]]&,Drop[al,{l,l}]],
					{-au[[u]]+al[[l]]}
				],
				Join[
					bn,
					Map[#-bu[[u]]&,Drop[bu,{u,u}]],
					Map[-#+bl[[l]]&,Drop[bl,{l,l}]],
					{bu[[u]]-bl[[l]]}
				]
			]],
			{l,1,Length[lower]},
			{u,1,Length[upper]}
		]
	]/;2<=Length[x];
convert[ipoly[f_,x_,a_,b_],g_,{1}]:=
	Block[
		{lower={},upper={},neither={},aa,al,bl,au,bu,an,bn},
		Map[
			If[
				a[[#,1]]===0,
				AppendTo[neither,#],
				If[
					Positive[a[[#,1]]],
					AppendTo[upper,#],
					AppendTo[lower,#],
					Message[
						convert::unknown,
						Positive[a[[#,1]]]
					]
				]
			]&,
			Range[1,Length[a]]
		];
		bu=Map[b[[#]]/a[[#,1]]&,upper];
		bl=Map[b[[#]]/a[[#,1]]&,lower];
		bn=Map[b[[#]]&,neither];
		If[
			Apply[And,Map[(0<=#)&,bn]],
			Sum[
				If[
					bl[[l]]<=bu[[u]],
					g[
						f,
						Release[{x[[1]],bl[[l]],bu[[u]]}]
					],
					0,
					Message[
						convert::unknown,
						Map[TestValue,bl[[l]]<=bu[[u]]]
					];0
				],
				{l,1,Length[lower]},
				{u,1,Length[upper]}
			],
			0,
			Message[
				convert::unknown,
				Apply[
					And,
					Map[(0<=TestValue[#])&,bn]
				]
			];0
		]
	]/;1==Length[x]

Define delta.

Define Kroneker delta function:

Attributes[delta]={Orderless};
delta[ajj_,ajj_]=1;
delta[all_,bll_]:=0/;all!=bll;
delta/: delta[a_,b_]^n_?Positive=delta[a,b];

rotate[ipoly[I],r] transforms the coordinate vector as x=r.x.

rotate::zero="the matrix `` has zero determinant.";
rotate[ipoly[f_,x_,a_,b_],r_?MatrixQ]:=
Block[
	{rinverse},
	rinverse=Inverse[r];
	ipoly[
		Abs[Det[rinverse]]*
		(f/.Inner[Rule,x,Expand[rinverse.x],List]),
		x,
		a.rinverse,
		b
	]
]/;If[
		Det[r]===0,
		Message[rotate::zero,r];False,
		True
	]&&freelist[r,x]

translate[ipoly[I],c] transforms the coordinate vector as x=x+c.

translate[ipoly[f_,x_,a_,b_],c_?VectorQ]:=
	ipoly[
		f/.Inner[Rule,x,Expand[x-c],List],
		x,
		a,
		b+a.c
	]/;freelist[c,x]

showpoly[expr] Prints the polytope for each ipoly in expr.

showpoly[expr_]:=
	expr/.ipoly[f_,x_,a_,b_]:>
	(Print[
		Block[
			{ff},
			Hold[
				Release[Inner[ff,Release[a.x],b,List]]
			]/.ff[aa_,bb_]->LessEqual[aa,bb]
		]
	];ipoly[showpoly[f],x,a,b])

checkpoly[ipoly[f,x,a,b]] checks that the matrix equation a.x<=b defines a valid convex polytope.

checkpoly uses a lot of subroutines:  divider,undouble, keepplane, polyQ, findvertex, and findrows.  It assumes there are no "boundaries at Infinity."

checkpoly::missing="Limits on variables `` are missing.";
checkpoly::unknown="Value of `` is unknown.";
checkpoly[a_ b_]:=a checkpoly[b]/;FreeQ[a,ipoly];
checkpoly[a_ +b_]:=checkpoly[a]+checkpoly[b];
checkpoly[ipoly[f_,x_,a_?MatrixQ,b_]]:=
	Block[
		{mat,ee,div,j=1,zrow={},n,
					nonzero={},matt,jj,nn,aa,bb},
(***** make sure that Testvalue gives numerical results **)
		(Map[
			If[
				Not[NumberQ[#]],
				Message[checkpoly::unknown,#]
			]&,
			Append[
				aa=Map[TestValue,a,{2}],
				bb=Map[TestValue,b]
			],
			{2}
		];
		mat=Map[(#==0)&,aa,{2}];
		While[
			j<=jj&&If[
				Apply[And,mat[[j]]],
				0<=bb[[j]],
				nonzero={nonzero,j};True
			],
			j++
		];
		Which[
			j<=jj,
			0,
			(matt=Transpose[mat=mat[[#]]];aa=aa[[#]];
				bb=bb[[#]];jj=Length[#])&[Flatten[nonzero]];
			Do[
				If[Count[matt[[n]],False]<2,zrow={zrow,n}],
				{n,1,nn}
			];zrow=!={},
			Message[
				checkpoly::missing,
				x[[Flatten[zrow]]]
			];ipoly[f,x,a,b],
(***This is the part that uses keepplane, etc.***)
			div=Apply[
				{#1[[undouble[aa[[##]],bb[[#1]]]]],#2}&,
				divider[mat],
				{1}
			];
			(*Print["div:  ",div];*)
			ee=Apply[
				(keepplane[aa[[##]],bb[[#1]]])&,
				div,
				{1}
			];MemberQ[ee,Null],
			0,
			True,
			ee=Map[div[[#,1,ee[[#]]]]&,Range[1,Length[div]]];
			(*Print["ee:  ",ee];*)
			ipoly[f,x,a[[#]],b[[#]]]&[Sort[Flatten[ee]]]
		])/;(({jj,nn}=Dimensions[a])=={Length[b],Length[x]})
	]

Takes matrix and returns a list of {list of rows in a group, list of columns in a group}.

divider[mat_]:=
	Block[
		{rr,ss,s=0,u,l,j,n},
		rr=Array[0&,jj];
		ss=Array[0&,nn];
		Do[
			Which[
				mat[[j,n]],
				Null,
				rr[[j]]==ss[[n]]==0,
				s++;rr[[j]]=ss[[n]]=s,
				rr[[j]]==0,
				rr[[j]]=ss[[n]],
				ss[[n]]==0,
				ss[[n]]=rr[[j]],
				rr[[j]]!=ss[[n]],
				l=Min[rr[[j]],ss[[n]]];
				u=Max[rr[[j]],ss[[n]]];s--;
				{rr,ss}=Map[
					Which[#<u,#,#>u,#-1,#==u,l]&,
					{rr,ss},
					{2}
				]
			],
			{j,1,jj},
			{n,1,nn}
		];
	(*Print[rr];Print[ss];Print[s];*)
		Table[
			{Select[Range[1,jj],(rr[[#]]==ii)&],
			Select[Range[1,nn],(ss[[#]]==ii)&]},
			{ii,1,s}
		]
	]

Removes all superimposed hyperplanes.  Returns a list of rows to keep.  It assumes all zero rows have been removed.

undouble[a_,b_]:=
	Block[
		{aa,which={1},k,j,c},
(***** append b to a and rescale a       *****)
		aa=Map[
			(c=1;While[#[[c]]==0,c++];#/Abs[#[[c]]])&,
			Transpose[Append[Transpose[a],b]]
		];
(***** find pairs of rows that are equal  *****)
		Do[
			k=1;While[k<j&&aa[[j]]=!=aa[[k]],k++];
			If[k==j,which={which,j}],
			{j,2,Length[a]}
		];
		Flatten[which]
	]

keepplane[a,b] returns a list of the faces which should be kept or Null if there is no nn-polytope.  A face is kept if a.x<=b restricted to that face is itself an nn-1 polytope.  It assumes that all zero rows and superimposed planes have been removed.

keepplane[a_,b_]:=
	Block[
		{j=1,nn,jj,facelist={},c,ta},
		(jj=Length[b];ta=Transpose[a];
		Off[LinearSolve::nsl];
		While[
			j<=jj,
			c=1;While[a[[j,c]]==0,c++];
			If[
				polyQ[
					Abs[a[[j,c]]]*
						Map[Drop[#,{c}]&,Drop[a,{j}]]-
					Sign[a[[j,c]]]*Outer[
						Times,#,Drop[a[[j]],{c}]],
					Abs[a[[j,c]]] Drop[b,{j}]-
					Sign[a[[j,c]]] b[[j]] (#)
				]&[Drop[ta[[c]],{j}]],
				facelist={facelist,j}
			];j++
		];
		On[LinearSolve::nsl];
		If[Length[#]>nn,#]&[Flatten[facelist]])/;
		(nn=Length[First[a]])>1
	]

keepplane[a_,b_]:=
	Block[
		{aa,j=1,jj,min=-Infinity,max=Infinity,jmin,jmax},
		(*Print["(1) a,b:  ",{a,b}];*)
		aa=Flatten[a];jj=Length[b];
		While[
			j<=jj&&
			If[
				aa[[j]]>0,
				(max<#||(jmax=j;min<(max=#)))&[
					b[[j]]/aa[[j]]],
				(min>#||(jmin=j;max>(min=#)))&[
					b[[j]]/aa[[j]]]
			],
			j++
		];
		If[
			(-Infinity<min<max<Infinity),
			{jmin,jmax}
		]
	]/;Length[First[a]]==1

Returns the result of whether a.x<=b describes an nn-polytope.  It uses checks that it has a vetex and, given that vertex, whether it has an (nn-1)-face not in that vertex.

polyQ[a_,b_]:=
	Block[
		{aa,bb,j=1,nn,jj,k=1,kl,nonzero={},c,taa,vertex},
		((*Print["a,b:  ",{a,b}];*)
		jj=Length[b];
(***** Remove rows with all zeros         *******)
		While[
			j<=jj&&
			If[
				Apply[And,Map[(#==0)&,a[[j]]]],
				0<=b[[j]],
				nonzero={nonzero,j};True
			],
			j++
		];j==jj+1&&
		(taa=Transpose[aa=a[[#]]];bb=b[[#]];
			nn<(kl=Length[#]))&[Flatten[nonzero]]&&
(***** Find a vertex                      ********)
		(vertex=findvertex[aa,bb])=!=Null&&
(***** find an n-1 face not in that vertex *******)
		(While[
			k<=kl&&
			(aa[[k]].vertex>=bb[[k]]||
			(c=1;While[aa[[k,c]]==0,c++];
			Not[polyQ[
				Abs[aa[[k,c]]]*
					Map[Drop[#,{c}]&,Drop[aa,{k}]]-
				Sign[aa[[k,c]]] Outer[
					Times,#,Drop[aa[[k]],{c}]],
				Abs[aa[[k,c]]] Drop[bb,{k}]-
				Sign[aa[[k,c]]] bb[[k]] #
			]&[Drop[taa[[c]],{k}]]])),
			k++
		];k<=kl))/;(nn=Length[First[a]])>1
	]

The one dimensional case is handled directly:

polyQ[a_,b_]:=
	Block[
		{aa,j=1,jj,min=-Infinity,max=Infinity},
		(*Print["(1) a,b:  ",{a,b}];*)
		aa=Flatten[a];jj=Length[b];
		While[
			j<=jj&&
			Which[
				aa[[j]]==0,
				0<=b[[j]],
				aa[[j]]>0,
				min<(max=Min[max,b[[j]]/aa[[j]]]),
				True,
				(min=Max[min,b[[j]]/aa[[j]]])<max
			],
			j++
		];j>jj&&(-Infinity<min<max<Infinity)
	]/;Length[First[a]]==1

Checks whether a.x<=b has any vertices.  First it finds an intersection of nn hyperplanes.  Then it moves inward, switching pairs of planes, until a vertex is found.  It returns either the coordinates of the vertex or Null if no vertex exists.  This assumes that any zero rows have been removed from a.

findvertex[a_,b_,v_:Null,w_:Null]:=
	Block[
		{jj,j=1,nn,n,vertex,v2,aa,bb,which},
		((*Print["a,b:  ",{a,b,w}];*)
(***** handle case where v and w are not given    *****)
		Which[
			w=!=Null,
			which=w;vertex=v,
			(which=findrows[a,jj,nn])=!=Null,
			vertex=LinearSolve[a[[which]],b[[which]]]
		];
(***** find any discrepencies with initial choice *****)
		If[
			which=!=Null,
			While[
				j<=jj&&(MemberQ[which,j]||
				a[[j]].vertex<=b[[j]]||(n=1;
				aa=Append[a[[which]],a[[j]]];
				bb=Append[b[[which]],b[[j]]];
(***** given a discrepency find which old plane to drop **)
				While[
					n<=nn&&(Not[
						VectorQ[
							v2=LinearSolve[
								Drop[aa,{n}],
								Drop[bb,{n}]
							]
						]
					]||aa[[n]].v2>=bb[[n]]),
					n++
				];False)),
				j++
			];
			Which[
				j>jj,
				vertex,
				n<=nn,
				findvertex[a,b,v2,Append[Drop[which,{n}],j]]
			]
		])/;(jj=Length[b])>(nn=Length[First[a]])
	]

Return a set of nn hyperplanes from the jj hyperplanes given which intersect at a point or Null if none exists.

findrows[a_,jj_,nn_]:=
	Block[
		{mat,j=1,n=1,which},
		mat=RowReduce[Transpose[a]];
		which=Array[0&,nn];
		While[
			j<=jj&&n<=nn,
			If[mat[[n,j]]!=0,which[[n]]=j;n++];j++
		];
		If[Not[MemberQ[which,0]],which]
	]

transform[ipoly[I],{q1,q2,I}] where qi are expressions linear in the variables defining the new value of each variable, xi=qi, i=1,2,I .

transform[a_ + b_,y_]:=transform[a,y]+transform[b,y];
transform[a_ b_,y_]:=
	a transform[b,y]/;Head[a]=!=Times&&FreeQ[a,ipoly];
transform[ipoly[f_,x_,a_,b_],y_List]:=
	translate[
		rotate[ipoly[f,x,a,b],Transpose[Drop[#,-1]]],
		#[[-1]]
	]&[
		Transpose[Map[coefficients[#,x]&,y]]
	]/;Apply[And,Map[LinearQ[#,x]&,y]]

Dotexpand, timesexpand, sumexpand, and ipolyexpand are Lists  of rules which Distribute the related functions over Plus, remove constants and other functions to the outside, and simplify delta functions.  The first two use ScalarQ[I] and CnumberQ[I].

freelist extends FreeQ to a List.  Constants divide a product into constants and nonconstants with respect to y;  The approach used here is MUCH faster for large expressions than the standard pattern recognition approach.

freelist[x_ ,{y_ ,yy__}]:=FreeQ[x,y]&&freelist[x,{yy}];
freelist[x_ , {y_}]:=FreeQ[x,y];
constants[x_,y_]:=
	(Map[
		If[
			freelist[#,y],
			const=const #,
			nonconst=nonconst #
		]&,
		Apply[List,x]
	];const=!=1)

Dotexpand=Dispatch[{
	Dot[aqa_?Global`ScalarQ bqb_,cqc_]:>
		NonCommutativeMultiply[aqa,Dot[bqb,cqc]],
	Dot[cqc_,aqa_?Global`ScalarQ bqb_]:>
		NonCommutativeMultiply[Dot[cqc,bqb],aqa],
	Dot[aja_+bjb_,cjc_]:>Dot[aja,cjc]+Dot[bjb,cjc],
	Dot[cjc_,aja_+bjb_]:>Dot[cjc,aja]+Dot[cjc,bjb],
	trace[a_?Global`ScalarQ b_.]:>a trace[b]
}];

timesexpand=Dispatch[{
	NonCommutativeMultiply[aqa_?Global`CnumberQ bqb_,cqc_]:>
		aqa NonCommutativeMultiply[bqb,cqc],
	NonCommutativeMultiply[aqa_,bqb_?Global`CnumberQ cqc_]:>
		bqb NonCommutativeMultiply[aqa,cqc],
	NonCommutativeMultiply[ala_,blb_]:>
		blb ala/;Head[blb]=!=
			NonCommutativeMultiply&&Global`CnumberQ[blb],
	NonCommutativeMultiply[ala_,blb_]:>
		blb ala/;Head[ala]=!=
			NonCommutativeMultiply&&Global`CnumberQ[ala],
	NonCommutativeMultiply[ala_,blb_+clc_]:>
		NonCommutativeMultiply[ala,blb]+
		NonCommutativeMultiply[ala,clc],
	NonCommutativeMultiply[ala_+blb_,clc_]:>
		NonCommutativeMultiply[ala,clc]+
		NonCommutativeMultiply[blb,clc],
	NormalOrder[a_?Global`CnumberQ b_]:>
			a NormalOrder[b],
	NormalOrder[a_]:>
		a/;
		(FreeQ[a,NonCommutativeMultiply]&&
			FreeQ[a,Dot])||Global`CnumberQ[a]
}];

Simplifies only the simplest delta-functions.

sumexpand=Dispatch[{
	sum[a_ + b_,x__]:>sum[a,x]+sum[b,x],
	sum[a_Times ,x__]:>
		Block[
			{const=1,nonconst=1},
			const sum[nonconst,x]/;
			constants[a,Map[First,{x}]]
		],
	sum[a_ ,x__]:>
		a sum[1,x]/;
		Head[a]=!=Times&&a=!=1&&
		freelist[a,Map[First,{x}]],
	sum[a_ **b_,x__]:>
		a**sum[b,x]/;
		Head[a]=!=NonCommutativeMultiply&&
		freelist[a,Map[First,{x}]],
	sum[b_ **a_,x__]:>
		sum[b,x]**a/;
		Head[a]=!=NonCommutativeMultiply&&
		freelist[a,Map[First,{x}]],
	sum[NormalOrder[a_],x__]:>NormalOrder[sum[a,x]],
	sum[Dot[a_ ,b_],x__]:>
		Dot[a,sum[b,x]]/;
		Head[a]=!=Dot&&freelist[a,Map[First,{x}]],
	sum[Dot[b_ ,a_],x__]:>
		Dot[sum[b,x],a]/;
		Head[a]=!=Dot&&freelist[a,Map[First,{x}]],
	sum[ipoly[f_,x_,a_,b_],y__]:>
		ipoly[sum[f,y],x,a,b]/;
		freelist[{a,b},Map[First,{y}]],
	sum[aa_,z___,{x_,a_,b_,c_:1},zz___]:>
		If[
			Length[{z,zz}]>=1,
			sum[Floor[Max[b-a+c,0]/c] aa,z,zz],
			aa Floor[Max[b-a+c,0]/c]
		]/;FreeQ[{aa,z},x],
	sum[
		a_. delta[n_?NumberQ,x1_],
		z1___,{x1_,b1_?NumberQ,b2_,b3_:1},z2___
	]:>If[
			b1<=n<=b2&&Mod[n-b1,b3]==0,
			If[Length[{z1,z2}]==0,a,sum[a,z1,z2]]/.
				x1->n,
			0
	]/;NumberQ[b3]&&NumberQ[TestValue[b2]],
	sum[
		a_. delta[x1_ ,x2_ ],
		z1___,{x1_,b1__},z2___,{x2_,b1__},z3___
	]:>(sum[a,z1,z2,{x2,b1},z3]/.x1->x2)
	}];

ipolyexpand::unknown="The value of `` is unknown.";
ipolyexpand=Dispatch[{
	ipoly[f_ + g_,x__]:>ipoly[f,x]+ipoly[g,x],
	ipoly[f_Times ,x_,a_,b_]:>
		Block[
			{const=1,nonconst=1},
			const ipoly[nonconst,x,a,b]/;
			constants[f,x]
		],
	ipoly[f_ ,x_,a_,b_]:>
		f ipoly[1,x,a,b]/;
		Head[f]=!=Times&&f=!=1&&freelist[f,x],
	ipoly[f_ **g_,x_,a_,b_]:>
		f**ipoly[g,x,a,b]/;
		Head[f]=!=NonCommutativeMultiply&&
		freelist[f,x],
	ipoly[g_ **f_,x_,a_,b_]:>
		ipoly[g,x,a,b]**f/;
		Head[f]=!=NonCommutativeMultiply&&
		freelist[f,x],
	ipoly[Dot[f_ ,g_],x_,a_,b_]:>
		Dot[f,ipoly[g,x,a,b]]/;
		Head[f]=!=Dot&&freelist[f,x],
	ipoly[NormalOrder[a_],x__]:>
		NormalOrder[ipoly[a,x]],
	ipoly[Dot[g_ ,f_],x_,a_,b_]:>
		Dot[ipoly[g,x,a,b],f]/;
		Head[f]=!=Dot&&freelist[f,x],
	ipoly[sum[f_,x__],y_,a_,b_]:>
		sum[ipoly[f,y,a,b],x]/;freelist[{x},y],
	ipoly[f_. delta[gg_],x_,a_,b_]:>
		Block[
			{c,d,n,xx,cc,alpha,g=Expand[gg]},
			(c=coefficients[g,x];
			 d=c[[-1]];c=Drop[c,-1];n=1;
			 While[c[[n]]==0&&n<Length[x],n++];
			 cc=Drop[c,{n,n}];xx=Drop[x,{n,n}];
			 alpha=Map[#[[n]]&,a];
			 checkpoly[ipoly[
				f/Abs[c[[n]]],
				xx,
				Map[Drop[#,{n,n}]&,a]-
					Outer[Times,alpha,cc]/c[[n]],
				b+alpha d/c[[n]]
			 ]]/.x[[n]]->Expand[-(d+cc.xx)/c[[n]]]
			)/;LinearQ[g,x]&&1<Length[x]
		],
	ipoly[f_. delta[gg_],{x_},a_,b_]:>
		Block[
			{c,d,g=Expand[gg]},
			({c,d}=coefficients[g,{x}];
			If[
				Inner[(#1<=#2)&,a.{Expand[-d/c]},b,And],
				f/Abs[c]/.x->Expand[-d/c],
				0,
				Message[
					ipolyexpand::unknown,
					Inner[
						Map[TestValue,#1<=#2]&,
						a.{Expand[-d/c]},
						b,
						And
					]
				];
				f/Abs[c]/.x->Expand[-d/c]
			])/;LinearQ[g,{x}]
		]
	}];

sumfactor and ipolyfactor are Lists  of rules which move constants and other functions to the inside of the related function, and combine various products of functions.

variables[a,b] is used internally.

variables[varx_List,vary_List]:=
	Block[
		{i,variablesa,variablesb,variablesc},
		variablesa=Intersection[varx,vary];
		variablesb=Union[varx,vary];
		For[
			i=1;variablesc=variablesa,
			Intersection[variablesc,variablesb]!={},
			++i,
			variablesc=Map[
				ToExpression[StringJoin[
					ToString[#],
					"is",
					ToString[i]
				]]&,
				variablesa
			]
		];
		Inner[(#1->#2)&,variablesa,variablesc,List]
	]

Local variables are assumed to be distinct between sum and ipoly

sumfactor=Dispatch[{
	a_ sum[b_,x__]:>sum[a b,x]/;FreeQ[a,sum],
	NormalOrder[sum[x_,y__]]:>sum[NormalOrder[x],y],
	a_ **sum[b_,x__]:>sum[a**b,x]/;FreeQ[a,sum],
	sum[b_,x__]**a_ :>sum[b**a,x]/;FreeQ[a,sum],
	Dot[a_ ,sum[b_,x__]]:>sum[Dot[a,b],x]/;FreeQ[a,sum],
	Dot[sum[b_,x__],a_]:>sum[Dot[b,a],x]/;FreeQ[a,sum],
	ipoly[sum[f_,y__],x_,a_,b_]:>
		sum[ipoly[f,x,a,b],y]/;freelist[{y},x],
	sum[a_,x__]+sum[b_,y__]:>
		sum[
			(a/.Inner[
				(#1:>#2)&,
				Map[First,{x}],
				Map[First,{y}],
				List
			])+b,
			y
		]/;Map[Rest,{x}]==Map[Rest,{y}],
	sum[a_,x__] sum[b_,y__]:>
		Apply[
			sum[
				(a/.variables[Map[First,{x}],
					Map[First,{y}]]) b,
				##
			]&,
			Join[
				{x}/.variables[Map[First,{x}],
					Map[First,{y}]],
				{y}
			]
		],
	sum[a_,x__]**sum[b_,y__]:>
		Apply[
			sum[
				(a/.variables[Map[First,{x}],
					Map[First,{y}]])**b,
				##
			]&,
			Join[
				{x}/.variables[Map[First,{x}],
					Map[First,{y}]],
				{y}
			]
		],
	Dot[sum[a_,x__],sum[b_,y__]]:>
		Apply[
			sum[
				Dot[a/.variables[Map[First,{x}],
					Map[First,{y}]],b],
				##
			]&,
			Join[
				{x}/.variables[Map[First,{x}],
					Map[First,{y}]],
				{y}
			]
		]
	}];

ipolyfactor=Dispatch[{
	ipoly[f_,x_,a_,b_] :> 
		Block[
			{ss=Transpose[
				Sort[Transpose[Append[Transpose[a],b]]]
			]},
			ipoly[f,x,Transpose[Drop[ss,-1]],ss[[-1]]]
		]/;Not[OrderedQ[a]],
	f_ ipoly[g_,x__]:>ipoly[f g,x]/;FreeQ[f,ipoly],
	NormalOrder[ipoly[g_,x__]]:>ipoly[NormalOrder[g],x],
	f_ **ipoly[g_,x__]:>ipoly[f**g,x]/;FreeQ[f,ipoly],
	ipoly[g_,x__]**f_ :>ipoly[g**f,x]/;FreeQ[f,ipoly],
	Dot[f_ ,ipoly[g_,x__]]:>
		ipoly[Dot[f,g],x]/;FreeQ[f,ipoly],
	Dot[ipoly[g_,x__],f_]:>
		ipoly[Dot[g,f],x]/;FreeQ[f,ipoly],
	sum[ipoly[f_,x_,a_,b_],y__]:>
		ipoly[sum[f,y],x,a,b]/;
		freelist[{a,b},Map[First,{y}]],
	ipoly[f_,x_,a_,b_]+ipoly[g_,y_,a_,b_]:>
		ipoly[(f/.Inner[(#1:>#2)&,x,y,List])+g,y,a,b],
	ipoly[f_,x_,a_,b_]+ipoly[g_,y_,c_,d_] :> 
		Block[
			{ss=Transpose[
				Sort[Transpose[Append[Transpose[a],b]]]
			]},
			ipoly[f,x,Transpose[Drop[ss,-1]],ss[[-1]]]+
			ipoly[g,y,c,d]
		]/;Not[OrderedQ[a]],
	ipoly[f_,x_,a_,b_] ipoly[ff_,xx_,aa_,bb_]:>
		ipoly[
			(f/.variables[x,xx]) ff,
			Join[x/.variables[x,xx],xx],
			blockdiagonal[a,aa],
			Join[b,bb]
		],
	ipoly[f_,x_,a_,b_]**ipoly[ff_,xx_,aa_,bb_]:>
		ipoly[
			(f/.variables[x,xx])**ff,
			Join[x/.variables[x,xx],xx],
			blockdiagonal[a,aa],
			Join[b,bb]
		],
	ipoly[f_,x_,a_,b_].ipoly[ff_,xx_,aa_,bb_]:>
		ipoly[
			(f/.variables[x,xx]).ff,
			Join[x/.variables[x,xx],xx],
			blockdiagonal[a,aa],
			Join[b,bb]
		]
	}];

blockdiagonal[a,b] combines matrices a and b in block diagonal form.

blockdiagonal[a_?MatrixQ,aa_?MatrixQ]:=
	Apply[
		Transpose[Join[
			Transpose[Join[
				a,
				Array[0&,{#3,#2}]
			]],
			Transpose[Join[
					Array[0&,{#1,#4}],
					aa
			]]
		]]&,
		Join[Dimensions[a],Dimensions[aa]]
	]

Define CnumberQ[x] to be False for noncommuting operators and ScalarQ[x] to be false for vectors and matrices.  Power[I] is assumed to be True for each (this is done for speed).

ScalarQ[x_,y__]:=ScalarQ[x]&&ScalarQ[y];
ScalarQ[delta[__]]=True;
ScalarQ[Power[_,_]]=True;
ScalarQ[_[x__]]:=ScalarQ[x];
ScalarQ[_?AtomQ]=True;

There is a bug in Mathematica reordering of the rules.  The order used here will give the correct results. (the problem comes from the first one)

CnumberQ[x_,y__]:=CnumberQ[x]&&CnumberQ[y];
CnumberQ[delta[__]]=True;
CnumberQ[Power[_,_]]=True;
CnumberQ[_[x__]]:=CnumberQ[x];
CnumberQ[_?AtomQ]=True;

LinearQ[expr,x] checks if expr is linear in the variables in the List x.  coefficients[expr,x] where expr is linear in the variables in the List x, gives a List of each coefficient and the constant term.

LinearQ[a_ + b_,x_]:=LinearQ[a,x]&&LinearQ[b,x];
LinearQ[a_ b_,x_List]:=LinearQ[b,x]/;freelist[a,x];
LinearQ[a_ ,x_List]:=True/;freelist[a,x];
LinearQ[a_,x_List]:=Count[x,a]==1

coefficient is local.

coefficients[a_,x_List]:=
	Append[
		Map[coefficient[a,#]&,x],
		a/.Map[(#->0)&,x]
	];
coefficient[a_+b_,x_]:=coefficient[a,x]+coefficient[b,x];
coefficient[a_,x_]:=0/;FreeQ[a,x];
coefficient[a_. x_,x_]:=a

End[];
EndPackage[];