(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



Off[General::"spell",General::"spell1",Remove::"rmnsm",UpSet::"write"];


Unprotect["Global`*"];ClearAll["Global`*"];Remove["Global`*"];Unprotect[$];ClearAll[Wedge,K,$];Protect[$];Unprotect[In,Out];Clear[In,Out];Protect[In,Out];$Line=0;$RecursionLimit=256;$IterationLimit=4096;


Forms[i_]:={};AllScalForms={};matForms[i_]:={};AllMatForms={};AllDifForms={};AllMatrices={};AllSymbols={};FormVars={_Wedge,_d};ZeroVars={_Wedge,_tr};HeadList={};ScalHeadList={};MatFormHeadList={};ZeroHeads={};AllMatHeads={};nodHeads={tr,Bar,Pattern,Blank,Condition,RuleDelayed,SeriesData};$EDCversion=389;


zeroQ[0]=True;zeroQ[x_List]:=And@@(zeroQ/@Union[Flatten[x]]);zeroQ[x_SeriesData]:=If[x[[3]]==={},True,False];zeroQ[x_]:=False;


nonZeroL[x_]:=Select[Union[Flatten[{x}]],!zeroQ[#]&]


indepTerms[x_,opt_:{}]:=If[NumericQ[x/._Wedge->1],reWrite[indep$cTerms[Factor[x],opt]],indep$cTerms[x,opt]]


indep$cTerms[x_,opt_:{}]:=redPlus[redTimes[nonZeroL[x],opt]];


redTimes[x_,opt_:{}]:=Block[{tmp$=Select[x,Head[#1]===Times &],rst$},rst$=Complement[x,tmp$]; Union[rst$,Map[numb2[#,opt]&,tmp$]]]


numb1[x_,opt_:{}]:=1/;NumericQ[x];numb1[x_^y_.,opt_:{}]:=1/;MemberQ[opt,x];numb1[x_,opt_:{}]:=x;numb2[x_,opt_:{}]:=Times@@Map[numb1[#,opt]&,Level[x,1]]


redPlus[x_]:=Block[{tmp$=Select[x,Head[#1]===Plus&],cmn$,rst$,i$,negi$},rst$=Complement[x,tmp$];cmn$=Intersection[tmp$,-tmp$];If[cmn$==={},Return[x]];i$=1;While[i$<Length[cmn$],negi$=Select[cmn$,#1===-cmn$[[i$]]&];If[negi$=!={},cmn$=Complement[cmn$,negi$]];i$=i$+1];Union[rst$,Complement[tmp$,-cmn$]]]


indepCoefs[y_]:=indepTerms[If[FormDegree[y]>0||MatQ[y],Flatten[(#1\[Wedge]$\[Wedge]$&)/@indepTerms[y]//.a_Wedge x_+(z_:Null):>{x,z}]/.Null->Sequence[],y]];


Gn$Cf0[x_List,y__]:=Gn$Cf0[#,y]&/@x;Gn$Cf0[x_SeriesData,y_List]:=Gn$Cf0[x,#]&/@y;Gn$Cf0[x_SeriesData,y__]:=x/.(x[[3]]->Map[Coefficient[#,y]&,x[[3]]]);Gn$Cf0[x_,y__]:=Coefficient[x,y];GenCoef[x_,y__]:=FacSimp[Gn$Cf0[x,y]];


SetAttributes[Bar,{Listable}];Bar[Blank]=Blank;Bar[Bar[x_]]:=x;Bar[Complex[x_,y_]]:=Complex[x,-y];Bar[x_Plus|x_Times|x_Wedge|x_Power|x_Rule|x_Equal]:=Bar/@x;Bar[x_?NumericQ]:=x;Bar[x_RuleDelayed]:=Rule[Bar[x[[1]]],Bar[x[[2]]]]/.Rule->RuleDelayed;Bar[x_SeriesData]:=x/.{First[x]->Bar[First[x]],x[[2]]->Bar[x[[2]]],x[[3]]->Bar/@x[[3]]};Bar[DirectedInfinity[x_]]:=DirectedInfinity[Bar[x]];Bar[HoldPattern[tr[x_]]]:=tr[Bar[x]];Bar[d[x_]]:=d[Bar[x]];Bar[x_Condition]:=x/.x[[1]]->Bar[x[[1]]];Bar[Derivative[x__][y_][z__]]:=If[Union[Bar[{z}]]===Union[{z}],Derivative[Sequence@@Coefficient[{x}.{z},Bar[{z}]]][Bar[y]][z],Derivative[x][Bar[y]][Sequence@@Bar[{z}]]];Bar[h_[y__]/;FreeQ[nodHeads,h]&&FreeQ[h,Derivative]]:=If[Union[Bar[{y}]]===Union[{y}],Bar[h][y],Bar[h][Sequence@@Bar[{y}]]]/;FreeQ[{Integer,Blank,BlankSequence,Pattern,Condition},Head[First[{y}]]];


Qvalue[x_List]:=Length[Select[OrderedQ/@Transpose[{Drop[x,-1],Rest[x]}],#===False&]]


AllCaseTable[x_List]:=With[{rvstr=Reverse[Map[Transpose,x]/.-z_:>z]},If[traceOption===2,Join[Table[{Qvalue[RotateLeft[x,i]],RotateLeft[x,i],i},{i,0,Length[x]-1}],Table[{Qvalue[RotateRight[rvstr,i]],RotateRight[rvstr,i],-i-1},{i,0,Length[x]-1}]],Table[{Qvalue[RotateLeft[x,i]],RotateLeft[x,i],i},{i,0,Length[x]-1}]]]


bestOrd[x_List]:=Block[{tmp$=Sort[AllCaseTable[x]]},tmp$=Select[tmp$,First[#1]===tmp$[[1,1]]&];If[Length[tmp$]>1,First[Sort[Select[tmp$,#1[[2]]===tmp$[[1,2]]&],Abs[Last[#1]]<Abs[Last[#2]]&]],First[tmp$]]];


splitFDprod[n_,FDx_List]:=(Plus@@Take[FDx,n])(Plus@@Drop[FDx,n]);allFDprods[FDx_List]:=Sum[FDx[[i]]*(Plus@@Drop[FDx,i]),{i,Length[FDx]-1}];


traceOption=2;tr[x_]:=reWrite[Plus@@Transpose[x,{1,1}]]/;MatrixQ[x];tr[x_Plus]:=tr/@x;tr[x_SeriesData]:=(x/.x[[3]]->tr/@x[[3]])+(x/.x[[3]]->{});tr[x__ y_]:=(x tr[y])/;!MatQ[x];tr[d[x_]]:=d[tr[x]];tr[x__\[Wedge]y__]:=(x\[Wedge]tr[Wedge[y]])/;ScalFormQ[x];tr[x_]:=(tr[x]=0)/;AntisymQ[x];tr[Bar[x_]]:=With[{evalHx=tr[x]},Bar[evalHx]/;Head[evalHx]=!=tr];tr[x_Wedge/;Union[Head/@Level[x,1]]==={Bar}]:=With[{evalHx=tr[Bar/@x]},Bar[evalHx]/;Head[evalHx]=!=tr];tr[x_Wedge]:=Block[{inpList=List@@x,resList,FDx},resList=bestOrd[inpList];FDx=FormDegree/@inpList;Which[ Last[resList]>0,(-1)^splitFDprod[ Last[resList],FDx]tr[Wedge@@resList[[2]]], Last[resList]===-1,(-1)^(allFDprods[FDx]+Length[Select[inpList,AntisymQ]])tr[Wedge@@resList[[2]]], Last[resList]<-1,(-1)^(allFDprods[FDx]+splitFDprod[Abs[1+ Last[resList]],FDx]+Length[Select[inpList,AntisymQ]])tr[Wedge@@resList[[2]]]]]/;traceOption>0&&bestOrd[List@@x][[2]]=!=List@@x;tr[x__\[Wedge]y__]:=(tr[x\[Wedge]y]=0)/;(SymQ[Wedge[x]]&&AntisymQ[Wedge[y]]||SymQ[Wedge[y]]&&AntisymQ[Wedge[x]]);tr[x_]:=x matDimension/;!MatQ[x];


toComponents[x_+y_,rul_]:=toComponents[x,rul]+toComponents[y,rul];
toComponents[x_.*y_,rul_]:=reWrite[If[MatQ[x*y],x*y,x*y*IdentityMatrix[If[IntegerQ[matDimension],matDimension,Length[Select[(#1[[2]]&)/@rul,MatrixQ][[1]]]]]]/.rul];
toComponents[x_\[Wedge]y_,rul_]:=reWrite[If[MatQ[x\[Wedge]y],x\[Wedge]y,x\[Wedge]y*IdentityMatrix[If[IntegerQ[matDimension],matDimension,Length[Select[(#1[[2]]&)/@rul,MatrixQ][[1]]]]]]/.rul];


FormDegree[x_Plus]:=FormDegree[First[x]];FormDegree[x_Times]:=Plus@@FormDegree/@List@@x;FormDegree[x_Wedge]:=Plus@@FormDegree/@List@@x;FormDegree[d[x_]]:=FormDegree[d[x]]=1+FormDegree[x];FormDegree[x_List]:=FormDegree[Last[Union[Flatten[x]]]];FormDegree[tr[x_]]:=FormDegree[x];FormDegree[Bar[x_]]:=FormDegree[Bar[x]]=FormDegree[x];FormDegree[x_SeriesData]:=If[x[[3]]==={},0,FormDegree[x[[3]]]];FormDegree[x_]:=0;


FormCoef[x_List,y_]:=(FormCoef[#,y]&)/@x;FormCoef[x_Plus,y_]:=(FormCoef[#1,y]&)/@x;FormCoef[x_.*y_,y_]:=x;FormCoef[x_,y_Wedge]:=Fold[FormCoef,x,Reverse[List@@y]];FormCoef[x_,y_]:=0/;FreeQ[x/.d[y]->d[$J$],y];FormCoef[x_,y_]:=(((x/.d[y]->d[$J$])-(x/.d[y]->d[$J$]/.y->0))/.y->1/.d[$J$]->d[y])/;EvenQ[FormDegree[y]];FormCoef[x_,y_]:=(((x/.d[y]->d[$J$])-(x/.d[y]->d[$J$]/.y->0))/.d[$J$]->d[y]/.{HoldPattern[Wedge[p___,y,q___]]:>(-1)^Length[Select[{q},OddQ[FormDegree[#]]&]]Wedge[p,q]})/;OddQ[FormDegree[y]];


DeclareMatrixForms[z__]:=Module[{x={{z}},xi,rxi,h,ht,k,min1,tmp=AllMatrices},While[Head[x[[1,1]]]===List,x=First[x]];Unprotect[Transpose];Do[xi=x[[i]];If[Length[xi]<3,Print["Missing symbol for transposed matrix"];Abort[]];If[xi[[3]]===xi[[2]]||xi[[3]]===-xi[[2]],rxi={xi[[2]]},rxi=Rest[xi]];matForms[First[xi]]=Union[matForms[First[xi]],rxi];
AllMatrices=Union[AllMatrices,rxi];h=Head[xi[[2]]];If[First[xi]>0,AllMatForms=Union[AllMatForms,rxi],If[h=!=Symbol,ZeroHeads=Union[{h},ZeroHeads];If[Head[xi[[3]]]=!=Times,ZeroHeads=Union[{Head[xi[[3]]]},ZeroHeads]]]];If[h===Symbol,FormDegree[xi[[2]]]=First[xi];BasicMatQ[xi[[2]]]=True;BasicScalFormQ[xi[[2]]]=False;Transpose[xi[[2]]]=xi[[3]];If[xi[[3]]===-xi[[2]],tr[xi[[2]]]=0,If[xi[[3]]=!=xi[[2]],FormDegree[xi[[3]]]=First[xi];BasicMatQ[xi[[3]]]=True;BasicScalFormQ[xi[[3]]]=False;Transpose[xi[[3]]]=xi[[2]];tr[xi[[3]]]=tr[xi[[2]]]]],FormDegree[_h]=First[xi];BasicMatQ[_h]=True;BasicScalFormQ[_h]=False;ht=Head[xi[[3]]];If[ht===Times,min1=-1;ht=h;tr[h[k__]]=0,min1=1;If[ht=!=h,FormDegree[_ht]=First[xi];BasicMatQ[_ht]=True;BasicScalFormQ[_ht]=False;tr[ht[k__]]=tr[h[k]]]];Transpose[h[k__]]=min1 ht[k];Transpose[ht[k__]]=min1 h[k]],{i,Length[x]}];Protect[Transpose];
AllDifForms=Union[AllScalForms,AllMatForms];AllSymbols=Union[AllDifForms,AllMatrices];MatFormHeadList=Complement[Union[Head/@AllMatForms,MatFormHeadList],{Symbol}];AllMatHeads=Union[ZeroHeads,MatFormHeadList];
HeadList=Union[Head/@AllSymbols];DifFormSymbols=Drop[AllDifForms,-Length[Union[ScalHeadList,MatFormHeadList]]];nodHeads=Union[nodHeads,ZeroHeads];k=Thread[Blank[Union[ScalHeadList,MatFormHeadList]]];FormVars=Flatten[{_Wedge,_d,Union[k,Bar[k]],Union[DifFormSymbols,Bar[DifFormSymbols]],_tr}];ZeroMatSymbols=Drop[matForms[0],-Length[ZeroHeads]];ZeroVars=Flatten[{_Wedge,Union[Thread[Blank[ZeroHeads]],Bar[Thread[Blank[ZeroHeads]]]],Union[ZeroMatSymbols,Bar[ZeroMatSymbols]],_tr}];If[tmp=!={},Print["OK"]]];


DeclareForms[z__]:=Module[{h,x={{z}},xi,rxi,k,tmp=AllScalForms},While[Head[x[[1,1]]]===List,x=First[x]];Do[xi=x[[i]];rxi=Rest[xi];Forms[First[xi]]=Union[Forms[First[xi]],rxi];
AllScalForms=Union[AllScalForms,rxi];
Do[h=Head[rxi[[j]]];If[h===Symbol,FormDegree[rxi[[j]]]=First[xi];BasicScalFormQ[rxi[[j]]]=True;BasicMatQ[rxi[[j]]]=False;,FormDegree[_h]=First[xi];BasicScalFormQ[_h]=True;BasicMatQ[_h]=False;],{j,Length[rxi]}],{i,Length[x]}];
AllDifForms=Union[AllScalForms,AllMatForms];ScalHeadList=Complement[Union[Head/@AllScalForms,ScalHeadList],{Symbol}];AllSymbols=Union[AllDifForms,AllMatrices];
HeadList=Union[Head/@AllSymbols];DifFormSymbols=Drop[AllDifForms,-Length[Union[ScalHeadList,MatFormHeadList]]];k=Thread[Blank[Union[ScalHeadList,MatFormHeadList]]];FormVars=Flatten[{_Wedge,_d,Union[k,Bar[k]],Union[DifFormSymbols,Bar[DifFormSymbols]],_tr}];If[tmp=!={},Print["OK"]]];


If[$VersionNumber<6.,st$pr6$[x__]:=StylePrint[x],st$pr6$[x__]:=Print[Style[x]]];


stPrint[x_,y_:Null]:=If[y===Null,st$pr6$[x,FontSize->16,FontSlant->"Italic",FontColor->RGBColor[1,0,1]],st$pr6$[x,FontSize->20,FontWeight->"Bold",FontColor->RGBColor[1,0,1]]];


NoDif[z__]:=(nodHeads=Union[nodHeads,Flatten[{z}]];If[z===interiorProduct||z===int$Prod||z===LieDcartan,Return[],Print["OK"]])


BasicMatQ[Bar[x_]]:=BasicMatQ[x];BasicMatQ[d[x_]]:=BasicMatQ[x];BasicMatQ[x_List]=True;BasicMatQ[x_]:=False;


MatQ[x_Times|x_Wedge|x_Plus]:=Or@@MatQ/@List@@x;MatQ[x_SeriesData]:=Or@@MatQ/@x[[3]];MatQ[x_]:=BasicMatQ[x];


BasicScalFormQ[HoldPattern[tr[x_]]]:=FormDegree[x]>0;BasicScalFormQ[Bar[x_]]:=BasicScalFormQ[x];BasicScalFormQ[d[x_]]:=!BasicMatQ[x];BasicScalFormQ[x_]:=False;


ScalFormQ[x_Times]:=Or@@ScalFormQ/@List@@x;ScalFormQ[x_Wedge|x_Plus]:=And@@ScalFormQ/@List@@x;ScalFormQ[x_]:=BasicScalFormQ[x];


Unprotect[Transpose];Transpose[x_Plus]:=Transpose/@x;Transpose[x_SeriesData]:=x/.x[[3]]->Transpose/@x[[3]];Transpose[x__ y_]:=x Transpose[y]/;!MatQ[x];Transpose[x__\[Wedge]y__]:=(x\[Wedge]Transpose[Wedge[y]])/;ScalFormQ[x];Transpose[x_Wedge]:=(-1)^allFDprods[FormDegree/@(List@@x)]Wedge@@Map[Transpose,Reverse[List@@x]];Transpose[d[x_]]:=d[Transpose[x]];HoldPattern[Transpose[Bar[x_]]]:=Bar[Transpose[x]];Transpose[x_]:=x/;!MatQ[x];Protect[Transpose];


SymQ[Bar[x_]]:=SymQ[x];SymQ[x_]:=Union[Flatten[{reWrite[Transpose[x]-x]}]]==={0};AntisymQ[Bar[x_]]:=AntisymQ[x];AntisymQ[x_]:=Union[Flatten[{reWrite[Transpose[x]+x]}]]==={0};


SetAttributes[d,{Listable}];d[x_Times|x_Wedge]:=d[First[x]]\[Wedge]Rest[x]+(-1)^FormDegree[First[x]]*First[x]\[Wedge]d[Rest[x]];d[x_?NumericQ|x_d]=0;d[matDimension]=0;d[Power[y_,n_]]:=n y^(n-1) d[y]+y^n Log[y]d[n];d[x_Plus]:=d/@x;HoldPattern[d[tr[x_]]]:=With[{evalHx=d[x]},tr[evalHx]/;Head[evalHx]=!=d];HoldPattern[d[Bar[x_]]]:=With[{evalHx=d[x]},Bar[evalHx]/;Head[evalHx]=!=d];d[x_Rule|x_Equal]:=reWrite[d/@x];d[x_SeriesData]:=(x/.x[[3]]->d[x[[3]]])+Wedge[d[First[x]],D[x,First[x]]];d[h_[y__]/;FreeQ[nodHeads,h]]:=Sum[(Derivative[Sequence@@RotateRight[Append[Table[0,{Length[{y}]-1}],1],i]][h][y])d[{y}[[i]]],{i,Length[{y}]}]/;FormDegree[h[y]]===0&&FreeQ[{Integer,Blank,Pattern,Condition},Head[First[{y}]]];


DeclareForms[1,e[_]];


NoDif[interiorProduct];interiorProduct[x_List,0,$z_:Null]=0;interiorProduct[x_List,y_,$z_:Null]:=0/;FormDegree[y]===0;interiorProduct[x_List,y_,$z_:Null]:=reWrite[Which[zeroQ[y/.e[_]->0],If[(y/.d[_]->0)=!=y,er$Wrn[6]];int$Prod[x,y],zeroQ[y/.d[_]->0&&FormDegree[y/.d[k_]:>k]===0],If[$z===Null,er$Wrn[4];Return[],dvarTOe=MapThread[Rule,{d/@$z,e/@Range[Length[$z]]}];eTOdvar=MapThread[Rule,{e/@Range[Length[$z]],d/@$z}];If[(y/.dvarTOe)=!=(y/.dvarTOe/.d[_]->0),er$Wrn[6]];int$Prod[x,y/.dvarTOe]/.eTOdvar],zeroQ[y/.d[_]->0]&&FreeQ[y/.d[k_]:>k,Wedge],Return[],True,tmp=Union[Head/@Select[Union[Flatten[y/.Plus->List/.Times->List/.Wedge->List]],FormDegree[#]===1&]];If[Length[tmp]===1,tmp=tmp[[1]];int$Prod[x,y/.tmp->e]/.e->tmp,If[Length[tmp/.d->Sequence[]]===1,er$Wrn[6];int$Prod[x,y/.tmp->e]/.e->tmp,er$Wrn[5]]]]]/;Head[y]=!=Symbol;


er$Wrn[4]:=stPrint["Must give ordered list of variables (coordinates) as 3rd argument"];


er$Wrn[5]:=stPrint["Error: more than one symbol used for 1-forms: "<>ToString[tmp]];


er$Wrn[6]:=stPrint["To complete the calculation, all exterior derivatives must be expressed in terms of the basis 1-forms!"];


NoDif[int$Prod];int$Prod[x_,0]=0;int$Prod[x_,y_]:=0/;zeroQ[x];int$Prod[x_,y_List]:=int$Prod[x,#]&/@y;int$Prod[x_,y_]:=0/;FormDegree[y]===0;int$Prod[x_,y_Plus]:=int$Prod[x,#]&/@y;int$Prod[x_,u_*y_]:=u*int$Prod[x,y]/;FormDegree[u]===0;int$Prod[x_,e[b_]]:=x[[b]];int$Prod[x_,Wedge[e[k_],p__]]:=x[[k]]*Wedge[p]-Wedge[e[k],int$Prod[x,Wedge[p]]];


NoDif[LieDcartan];LieDcartan[x_,y_,$z_:Null]:=reWrite[d[interiorProduct[x,y,$z]]+interiorProduct[x,d[y],$z]];


newSer$[x_SeriesData,k_]:=SeriesData[First[x],x[[2]],Flatten[Transpose[Prepend[Table[Table[0,{Length[x[[3]]]}],{k-1}],x[[3]]]]],k x[[4]],k x[[5]],k Last[x]]


Wedge[x_]:=x/;Length[{x}]<2&&Head[{x}[[1]]]=!=Pattern


Default[Wedge]:=1;SetAttributes[Wedge,{Flat,OneIdentity}];Wedge[0,y__]=0;Wedge[x__,0]=0;Wedge[x_SeriesData,y_SeriesData]:=Block[{x$,y$,s1$,s2$,re$,x3$,y3$},If[Last[x]===Last[y],x$=x;y$=y,x$=newSer$[x,LCM[Last[x],Last[y]]/Last[x]];y$=newSer$[y,LCM[Last[x],Last[y]]/Last[y]]];s1$=x$[[-3]]+y$[[-3]];s2$=Min[x$[[-2]]+y$[[-3]],x$[[-3]]+y$[[-2]]];If[Length[x$[[3]]]<x$[[-2]]-x$[[-3]],x3$=Join[x$[[3]],Table[0,{x$[[-2]]-x$[[-3]]-Length[x$[[3]]]}]],x3$=x$[[3]]];If[Length[y$[[3]]]<y$[[-2]]-y$[[-3]],y3$=Join[y$[[3]],Table[0,{y$[[-2]]-y$[[-3]]-Length[y$[[3]]]}]],y3$=y$[[3]]];Which[s2$===s1$,re$={},s2$-s1$===x$[[-2]]-x$[[-3]],re$=Sum[Map[Wedge[x3$[[i]],#]&,Join[Table[0,{i-1}],Drop[y3$,-i+1]]],{i,Length[x3$]}],True,re$=Sum[Map[Wedge[#,y3$[[i]]]&,Join[Table[0,{i-1}],Drop[x3$,-i+1]]],{i,Length[y3$]}]];SeriesData[First[x$],x$[[2]],re$,s1$,s2$,Last[x$]]]/;First[x]===First[y]&&x[[2]]===y[[2]];Wedge[y_,x_SeriesData]:=x/.x[[3]]->Map[Wedge[y,#]&,x[[3]]];Wedge[x_SeriesData,y_]:=x/.x[[3]]->Map[Wedge[#,y]&,x[[3]]];Wedge[x__,y_Plus]:=Plus@@Map[Wedge[x,#]&,List@@y];Wedge[x_Plus,y__]:=Plus@@Map[Wedge[#,y]&,List@@x];Wedge[z__,Times[x_,y_]]:=Times[x,Wedge[z,y]]/;NumericQ[x]||(FormDegree[x]===0&&!MatQ[x]);Wedge[Times[x_,y_],z__]:=Times[x,Wedge[y,z]]/;NumericQ[x]||(FormDegree[x]===0&&!MatQ[x]);Wedge[x_^n_.,y_]:=x^n*y/;FormDegree[x]===0&&!MatQ[x];Wedge[y_,x_^n_.]:=x^n*y/;FormDegree[x]===0&&!MatQ[x];Wedge[x_,y___,x_]:=0/;OddQ[FormDegree[x]]&&ScalFormQ[x];Wedge[x__,y__]:=(-1)^(FormDegree[x]*FormDegree[y])*y\[Wedge]x/;MatQ[x]&&ScalFormQ[y];Wedge[x__]:=Block[{doubL=Transpose[{FormDegree/@{x},{x}}]},Signature[Select[doubL,OddQ[First[#]]&]]Wedge@@Map[Last[#1]&,Sort[doubL]]]/;Union[BasicScalFormQ/@{x}]==={True}&&Map[Last[#1]&,Sort[Transpose[{FormDegree/@{x},{x}}]]]=!={x};Wedge[x_List,y_List]:=If[(FormDegree[x]===0&&!MatQ[Last[Union[Flatten[x]]]])||(FormDegree[y]===0&&!MatQ[Last[Union[Flatten[y]]]]),Dot[x,y],Inner[Wedge,x,y]];Wedge[x_,y_List]:=If[FormDegree[x]===0||(FormDegree[y]===0&&!MatQ[Last[Union[Flatten[y]]]]),x *y,Map[Wedge[x,#]&,y]]/;!MatQ[x];Wedge[x_List,y_]:=If[FormDegree[y]===0||(FormDegree[x]===0&&!MatQ[Last[Union[Flatten[x]]]]),x*y,Map[Wedge[#,y]&,x]]/;!MatQ[y];


TrigRules={x_.*Cos[h_]^2+x_.*Sin[h_]^2:>x,x_.*Cot[h_]^2+y_.*Csc[h_]^2:>y/;x+y===0,(Cot[h_]-Csc[h_]) (Cot[h_]+Csc[h_])->-1};HypTrigRules={x_.*Cosh[h_]^2+y_.*Sinh[h_]^2:>x/;x+y===0,x_.*Coth[h_]^2+y_.*Csch[h_]^2:>x/;x+y===0,(Cosh[x_]-Sinh[x_]) (Cosh[x_]+Sinh[x_])->1,(Coth[x_]-Csch[x_]) (Coth[x_]+Csch[x_])->1};


simpRules = {x_.*Cos[h_]^2+x_.*Sin[h_]^2:>x};varList={};


SetAttributes[reWrite,{Listable}];reWrite[0]=0;reWrite[x_Equal]:=Equal[reWrite[First[x]-Last[x]],0];reWrite[x_Rule]:=Rule[First[x],reWrite[Last[x]]];reWrite[x_RuleDelayed]:=Rule[First[x],reWrite[x[[2]]]]/.Rule->RuleDelayed;reWrite[x_SeriesData]:=(x/.x[[3]]->reWrite[x[[3]]]);reWrite[x_Condition]:=x/.x[[1]]->reWrite[x[[1]]];reWrite[x_]:=(Collect[x/.simpRules,FormVars,bestFacRul]/.simpRules)/;FormDegree[x]>0;reWrite[x_]:=(Collect[x/.simpRules,ZeroVars,bestFacRul]/.simpRules)/;MatQ[x];reWrite[x_]:=bestFacRul[x/.simpRules]/.simpRules;


FreeALLQ[x_,y_List]:=And@@(FreeQ[x,#]&/@y);


bestFacRul[x_]:=If[varList==={}||FreeALLQ[x,varList],minFacRul[x],varFacRul[x]];varFacRul[x_]:=Collect[Expand[x]/.simpRules,varList,Factor]/.simpRules;minFacRul[x_]:=Factor[Expand[x]/.simpRules]/.simpRules;


Unprotect[SeriesData];If[$VersionNumber<5,Times[x_SeriesData,0]^=0];SeriesData[x1_,x2_,x3_,x4_,x5_,x6_]:= Plus@@If[x2===Infinity,Map[First[#]/x1^((x4+Last[#]-1)/x6)&,Transpose[{x3,Range[Length[x3]]}]]+SeriesData[x1,x2,{},x4,x5,x6],   Map[First[#]*(x1-x2)^((x4+Last[#]-1)/x6)&,Transpose[{x3,Range[Length[x3]]}]]+SeriesData[x1,x2,{},x4,x5,x6]]/;!(FreeQ[x3,x1]||NumericQ[x1]);Protect[SeriesData];


subMatDR[m_,p_,q_]:=Transpose[Drop[Transpose[Drop[m,p]],q]];
subMatTK[m_,p_,q_]:=Transpose[Take[Transpose[Take[m,p]],q]];


FlatOuter[x_?MatrixQ,y_?MatrixQ]:=Partition[Flatten[Transpose[Outer[Times,x,y],{1,3,2,4}]],Dimensions[x][[2]]*Dimensions[y][[2]]];


FlatOuter[x__]:=Fold[FlatOuter,First[{x}],Rest[{x}]]/;Length[{x}]>2;


DirectSum[x_?MatrixQ,y_?MatrixQ]:=Block[{n$1=Dimensions[x],n$2=Dimensions[y],zr$},zr$[__]=0;Join[MapThread[Join,{x,Array[zr$,{n$1[[1]],n$2[[2]]}]}],MapThread[Join,{Array[zr$,{n$2[[1]],n$1[[2]]}],y}]]];


DirectSum[x__]:=Fold[DirectSum,First[{x}],Rest[{x}]]/;Length[{x}]>2;


DeclareMatrixForms[{0,$,$t}];


CHrule=CHrul;CHcoef[x_,0]=1;CHcoef[x_,1]:=tr[x];CHcoef[x_,n_]:=CHcoef[x,n]=reWrite[tr[CHeq[x,n]]/n];


CHeq[x_,n_]:=reWrite[CHcoef[x,n-1]\[Wedge] x+Sum[(-1)^(i-1)CHcoef[x,n-i]\[Wedge]Wedge@@Table[x,{i}],{i,2,n}]];


CHrul[n_,x_:$]:=Module[{cn=CHeq[x,n],cn0,rul,tmp,r1,r2},If[FormDegree[x]>0,Print[x," is not a 0-form matrix"];Return[]]; matDimension=n;rul=First[Solve[cn==cn0,Wedge@@Table[x,{n}]]]/.cn0->tr[cn]/n;rul=reWrite[rul];If[x===$,tmp=MapAt[HoldPattern,rul/.Rule->RuleDelayed,{1,1}];r1=tmp[[1,1]];r2=r1/.$->Pattern[$,Blank[]];tmp={Rule[r2,tmp[[1,2]]/;MatQ[$]&&FormDegree[$]===0&&matDimension===n]}/.Rule->RuleDelayed,tmp=rul]];


detTOtr[n_,A_:$]:=Module[{expr=CHrul[n][[1,2,1]]},(-1)^(1+matDimension) Plus@@Select[List@@expr,MatQ[#1]===False&]/.$->A];


SetAttributes[FaSi,{Listable}]; FaSi[0]=0;
FaSi[x_SeriesData]:=(x/.x[[3]]->FaSi[x[[3]]]);
FaSi[x_]:=(#1/.simpRules&)/@If[varList==={}||FreeALLQ[x,varList],Factor[x/.simpRules]/.simpRules,Collect[x/.simpRules,varList,Factor]/.simpRules];


SetAttributes[FacSimp,{Listable}];FacSimp[0]=0;FacSimp[x_Equal]:=Equal[FacSimp[First[x]-x[[2]]],0];FacSimp[x_Rule]:=Rule[First[x],FacSimp[x[[2]]]];FacSimp[x_SeriesData]:=(x/.x[[3]]->FacSimp[x[[3]]]);FacSimp[x_]:=(ReplaceAll[#,simpRules]&/@Collect[x/.simpRules,FormVars,FaSi])/;FormDegree[x]>0;FacSimp[x_]:=(#1/.simpRules&)/@If[varList==={}||FreeALLQ[x,varList],Factor[x/.simpRules]/.simpRules,Collect[x/.simpRules,varList,Factor]/.simpRules];


TensRnk[x_]:=If[Head[x]===List,Length[Dimensions[x]],0]


multiDot[x_,y_,indPr__]:=Module[{rnk1=TensRnk[x],rnk2=TensRnk[y],indL=Sort/@Transpose[{indPr}],ordPr=Sort[{indPr}]},If[!And@@MapThread[SameQ[Dimensions[x][[#1]],Dimensions[y][[#2]]]&,Transpose[ordPr]],er$Wrn[3];Return[]];If[Union[First[indL]]=!=First[indL]||indL[[1,-1]]>rnk1||Union[indL[[2]]]=!=indL[[2]]||Last[indL[[2]]]>rnk2,er$Wrn[2];Return[]];Which[Length[ordPr]===1,FacSimp[dotOne[x,y,indPr]],rnk1===rnk2&&rnk1===Length[ordPr],FacSimp[FacSimp[dotAll[x,y,ordPr]]],True,OKindCon[dotOne[x,y,First[ordPr]],newPrList[ReplaceAll[#,{a_,b_}:>{a,b+rnk1}]&/@ordPr]]]];


dotOne[x_,y_,{i_,j_}]:=Block[{rnk1=TensRnk[x],rnk2=TensRnk[y],tmp$},If[(FormDegree[x]>0||MatQ[Last[Union[Flatten[x]]]])&&(FormDegree[y]>0||MatQ[Last[Union[Flatten[y]]]]),tmp$=Wedge,tmp$=Dot];tmp$[Transpose[x,Join[Range[i-1],RotateRight[Range[i,rnk1]]]],Transpose[y,Join[RotateLeft[Range[j]],Range[j+1,rnk2]]]]];


dotAll[x_,y_,indPr_]:=Module[{s$,allInd,sumInd,tmp$},allInd=s$/@Range[Length[indPr]];allInd={allInd,allInd/.Map[Rule[#[[1]],#[[2]]]&,RotateLeft/@indPr]};sumInd=Transpose[{First[allInd],Dimensions[x]}];If[(FormDegree[x]>0||MatQ[Last[Union[Flatten[x]]]])&&(FormDegree[y]>0||MatQ[Last[Union[Flatten[y]]]]),tmp$=Wedge,tmp$=Times];Sum@@{Unevaluated[tmp$[x[[Sequence@@First[allInd]]],y[[Sequence@@allInd[[2]]]]]],Sequence@@sumInd}];


Contract[x_,indPr__]:=Module[{rnk1=TensRnk[x],indL=Flatten[{indPr}]},If[!And@@Map[SameQ[Dimensions[x][[#[[1]]]],Dimensions[x][[#[[2]]]]]&,{indPr}],er$Wrn[3];Return[]];If[Length[Union[indL]]<Length[indL]||Length[indL]>rnk1||Sort[indL][[-1]]>rnk1,er$Wrn[2];Return[],OKindCon[x,Sort[Sort/@{indPr}]]]];


newPrList[x_]:=Rest[x]-1/.k_:>k-1/;k>x[[1,2]]-1;


OKindCon[x_,indPr_]:=FacSimp[ConOne[x,First[indPr]]]/;Length[indPr]===1;OKindCon[x_,indPr_]:=OKindCon[ConOne[x,First[indPr]],newPrList[indPr]];


ConOne[x_,{i_,j_}]:=Block[{rnk1=TensRnk[x],tr1List,tr2List},tr1List=(Range[rnk1]/.{j->i})/.k_:>k-1/;k>j;tr2List=Join[Range[i-1],RotateRight[Range[i,rnk1-1]]];Apply[Plus,Transpose[Transpose[x,tr1List],tr2List],{rnk1-2}]];


er$Wrn[2]:=stPrint["Wrong tensor rank or number(s) of indices!"];er$Wrn[3]:=stPrint["Summed indices have unequal dimensions!"];


$PrePrint=If[ByteCount[#]>200000,Short[#,100],#]&;


Protect[CHeq,CHrul,CHrule,detTOtr,numb1,numb2,redPlus,redTimes,subMatTK,subMatDR,FlatOuter,DirectSum,indepTerms,indep$cTerms,indepCoefs,FormCoef,ConOne,dotAll,dotOne,OKindCon,FaSi,FacSimp,GenCoef,interiorProduct,int$Prod];


On[General::"spell",General::"spell1",UpSet::"write"];