(*^
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,...,xn],{x1,...},{x2,...},...,{xn,...}] is the
sum over x1, x2, É, 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,...,xn],{x1,x2,...,xn},
{{a11,a12,...,a1n},{a21,a22,...,a2n},...,{aJ1,aJ2,...,aJn}},
{b1,b2,...,bJ}]\nis
the n-dimensional integral of f[...] 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 + ... + ajn xn <= bj,	1<=j<=J.\nInput form is
ipoly[f,x,{c1<=c2,c3<=c4,...}] 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.";
	
convert::usage="convert[expr,f,{v1,v2,...}] converts variables
number v1, v2, ... 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[...],{q1,q2,...}] where qi are 
expressions linear in the integration variables xi, 
defines the new value of each variable to be xi->qi,
i=1,2,... .";

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

rotate::usage="rotate[ipoly[...],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
//.Sequence[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//.Sequence[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/.sequence[rule,g] and expr//.Sequence[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 and signs using TestValue.  Define extensions to 
inequalities, Positive, Negative, NonNegative, and Sign.


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


Extensions to Conjugate[É] 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[y]},
		Left, 
		Above
]
		
(* There is some trouble with mixing Derivatives and conjugates
Conjugate/:D[Conjugate[x_],y_Symbol]:=
	Conjugate[D[x,y]];*)

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[É] is defined near the definitions of ipolyexpand and sumexpand


Built in definitions are put after the following six definitions.

Commented out in 1995

Big/:TestValue[Big]=1000000;
Attributes[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__}]:=
	Module[
		{cef},
		Integrate[Exp[cef],{x,xx}]/;
		f=!=(cef=Collect[Expand[f],x])
	];
Integrate[x__]:=
	Apply[Integrate,{x}/.Big->Infinity]/;
	Not[FreeQ[{x},Big]];
Integrate[f_. delta[a_. x_ +b_.], {x_,-Infinity,Infinity}]:=
	Sign[a] (f/.x->-b/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];
DownValues[Integrate]=RotateLeft[#,Length[#]-6]&[
	DownValues[Integrate]
];
Protect[Integrate];


Define various functions


sum[f[x1,x2,É,xn],{x1,É},{x2,É},É,{xn,É}] is the sum over x1, x2, É, 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,É,xn],{x1,x2,É,xn},{{a11,a12,É,a1n},{a21,a22,É,a2n},É,{aJ1,aJ2,É,
aJn}},{b1,b2,É,bJ}] is the n-dimensional integral of f[É] 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 + É + ajn xn <= bj,	1<=j<=J. 
 InputForm is ipoly[f,x,{c1<=c2,c3<=c4,É}] 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_]:=
	Module[
		{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,É}] converts variables number v1, v2, É 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]:=
	Module[
		{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]:=
	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[
					Sign[a[[#,y]]]==1,
					AppendTo[upper,#],
					AppendTo[lower,#],
					Message[
						convert::unknown,
						Sign[TestValue[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[
					Sign[a[[#,1]]]==1,
					AppendTo[upper,#],
					AppendTo[lower,#],
					Message[
						convert::unknown,
						Sign[TestValue[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[É],r] transforms the coordinate vector as x=r.x.


*)
rotate::zero="the matrix `` has zero determinant.";
rotate[ipoly[f_,x_,a_,b_],r_?MatrixQ]:=
Module[
	{rinverse},
	rinverse=Inverse[r];
	ipoly[
		(Sign[#] #)&[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[É],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[
		Inner[LessEqual,a.x,b,List]
	];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,a0,b0,
					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[[#]]];a0=a[[#]];
				b0=b[[#]];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,a0[[#]],b0[[#]]]&[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[#[[1,c]]==0,c++];#/Abs[#[[1,c]]])&,
			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::nosol];
		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::nosol];
		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[É],{q1,q2,É}] where qi are expressions linear in the variables 
defining the new value of each variable, xi=qi, i=1,2,É .


*)
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[É] and CnumberQ[É].


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_?ScalarQ bqb_,cqc_]:>
		NonCommutativeMultiply[aqa,Dot[bqb,cqc]],
	Dot[cqc_,aqa_?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_?ScalarQ b_.]:>a trace[b]
}];
(*


*)
timesexpand=Dispatch[{
	NonCommutativeMultiply[aqa_?CnumberQ bqb_,cqc_]:>
		aqa NonCommutativeMultiply[bqb,cqc],
	NonCommutativeMultiply[aqa_,bqb_?CnumberQ cqc_]:>
		bqb NonCommutativeMultiply[aqa,cqc],
	NonCommutativeMultiply[ala_,blb_]:>
		blb ala/;Head[blb]=!=
			NonCommutativeMultiply&&CnumberQ[blb],
	NonCommutativeMultiply[ala_,blb_]:>
		blb ala/;Head[ala]=!=
			NonCommutativeMultiply&&CnumberQ[ala],
	NonCommutativeMultiply[ala_,blb_+clc_]:>
		NonCommutativeMultiply[ala,blb]+
		NonCommutativeMultiply[ala,clc],
	NonCommutativeMultiply[ala_+blb_,clc_]:>
		NonCommutativeMultiply[ala,clc]+
		NonCommutativeMultiply[blb,clc],
	NormalOrder[a_?CnumberQ b_]:>
			a NormalOrder[b],
	NormalOrder[a_]:>
		a/;
		(FreeQ[a,NonCommutativeMultiply]&&
			FreeQ[a,Dot])||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[# aa,z,zz],
			aa #
		]&[If[
			IntegerQ[TestValue[(b-a)/c]]&&
			TestValue[(b-a)/c+1]>0,
			(b-a)/c+1,
			Max[Floor[(b-a)/c+1],0]
		]])/;FreeQ[{aa,z},x],
	sum[
		a_. delta[n_,x1_],
		z1___,{x1_,b1_,b2_,b3_:1},z2___
	]:>Block[
		{nn,bb1,bb2,bb3},
		If[
			bb1<=nn<=bb2&&Mod[nn-bb1,bb3]==0,
			If[Length[{z1,z2}]==0,a,sum[a,z1,z2]]/.
				x1->n,
			0
		]/;NumberQ[bb1=TestValue[b1]]&&
		NumberQ[bb2=TestValue[b2]]&&
		NumberQ[bb3=TestValue[b3]]&&
		NumberQ[nn=TestValue[n]]
	],
	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];
(*The integration variable x[[n]] must be substituted
out AFTER checkpoly is applied for the case where
the substition gives an error and checkpoly gives zero.*)
			 checkpoly[ipoly[
				Sign[c[[n]]] f/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],
				(Sign[c] f/c)/.x->Expand[-d/c],
				0,
				Message[
					ipolyexpand::unknown,
					Inner[
						Map[TestValue,#1<=#2]&,
						a.{Expand[-d/c]},
						b,
						And
					]
				];
				(Sign[c] f/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 to give new variable names to variables when 
combining sum and ipoly.


*)
variables[varx_List,vary_List]:=
	Module[
		{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}
			]
		]
	}];	
(*


ipolyfix[f,x,a,b]  rescales a and chooses an ordering for the rows of a based 
on the elements of a.  It also reorders and scales b a consistant manner.


*)
ipolysort[a_,b_]:=
	Reverse[Transpose[Sort[
		Map[
			(#/Apply[
				Min,
				Map[
					(# Sign[#])&,
					Select[Last[#],(#!=0)&]
				]
			])&,
			Transpose[{b,a}]
		]
	]]];
ipolyfactor=Dispatch[{
	ipoly[f1_,x1_,a1_,b1_]+ipoly[f2_,x2_,a2_,b2_]:>
		Module[
			{aa1,bb1,aa2,bb2},
			 If[
			 	{aa1,bb1}==={aa2,bb2},
			 	ipoly[
			 		(f1/.Inner[(#1:>#2)&,x1,x2,List])+
			 			f2,
			 			x2,aa2,bb2
			 	],
			 	ipoly[f1,x1,aa1,bb1]+
			 	ipoly[f2,x2,aa2,bb2]
			 ]/;
			 ({aa1,bb1}=ipolysort[a1,b1];
			 {aa2,bb2}=ipolysort[a2,b2];
			 {a1,b1}=!={aa1,bb1}||
			 {a2,b2}=!={aa2,bb2}||
			 {aa1,bb1}=={aa2,bb2})
			 
		],
	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[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[É] 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[];
(*


^*)