(* ::Package:: *)

(* ::Text:: *)
(*(* tensalgv3.m *)*)
(*(* Version 3. Started 2016.12.11. finished 2017.8.10.*)*)
(*(* V.1: Date: February 26, 2013; v.1.1, May 17, 2014 *)*)
(**)
(*(*Author Rolf Sulanke *)*)
(**)
(*(* This package contains definitions and procedures for n-dimensional linear and*)
(*tensor algebra. The functions and modules defined have been developed together *)
(*with the notebook vectensalgv2.nb, which can  serve as the basis for  the applications *)
(*of tensor algebra and as a tool for learning tensor algebra. The package contains the *)
(*functions and modules of my package vectorcalc.m. Version 2 contains some new *)
(*functions  and many corrections, see "Changes".   The notebook "Affine Geometry"  *)
(*(Affine Geom.nb) contains selected applications of the newly in v.3 introduced functions*)*)
(**)
(*(* Changes:*)
(*2013.03.23. SetAtributes[tensorToArrays,Listable];*)
(*2013.03.24. linform revised. lf introduced as short form of linform.*)
(*2013.03.25. covec introduced.*)
(*2013,04.01, Wedge (== wedge) introduced.*)
(*2014.01.25. plus introduced for the addition of functions. The zero-function*)
(*defined as a pure function 0&.*)
(*2014.02.19, vector and covector become synonyms of vec resp. covec.*)
(*2014.03.5. Function mathf: Definition changed. *)
(*2014.03.11. Functions vecTotensor,CovecTotensor introduced.*)
(*2014.03.12. Functions wedgeform and wedgevector corrected.*)
(*Function ArrayTotensor corrected. *)
(*2014.03.12. Functions formbasis, vecbasis corrected.*)
(*2014.03.12. epsilon tensor: Functions eps, epst introduced.*)
(*2014.05.17. extprod introduced.*)
(*2015.03.08. tensorToArray. Exception k=h=0 added.*)
(*2015.03.09. tensor. Head[a[i,j]]== List excluded.*)
(*2015.03.23. Infix operators \[Wedge],  \[CircleTimes] for tensors defined.*)
(*2015.05.21. ArrayTootensor korrigiert.*)
(*2015.05.23. strictcoords: Antisymmetry control added.*)
(*2015.05.26. wedge: Controls anmd operation optimized.*)
(*2015.05.27. tensorproduct accepts now also a list of tensors as argument.*)
(*2015.06.05. New Module for plus; sum deleted, replaces by plus.*)
(*2015.06.09. equi corrected and simplified.cosymtensor*)
(*2015.06.09. ArrayTotensor corrected and simplified.*)
(*2015.06.24. cova,covalt,contra,contralt changed; nuow circa 10x faster.*)
(*2015.06.28. formbasis ersetzt durch wstb; vecbasis ersetzt durch vstb.*)
(*2015.06.29. Introduced symform,altform,cosymtensor,contrasymtensor,coalttensor,contraalttensor.*)
(*2015.07.01. equi revised, gives detailed results now.*)
(*2015.07.09. Infix operator \[CircleTimes] as tensorproduct defined at the end of the package.*)
(*2015.07.26.  strictcoords extended to antisymmetric Arrays.*)
(*2016.03.01. Functions one and tpl added.*)
(*2016.12.16. Old obsolete function vector deleted.*)
(*	    New:point,vector, afftr, cgeo, cbary introduced. *)
(*2018.01.02. Definition of type changed.*)
(*	    *)*)
(**)
(*(* Keywords: vector objects, null vector, basis vectors, random vector,  *)
(*matrix, random matrix, rank, transposed matrix, tensor,  symmetric, *)
(*antisymmetric, covariant, contravariant, mixed type, symmetrize, alternate, *)
(*tensor product, contraction, basis tensors, tensor coordinates, wedge product,*)
(*exterior product, strict coordinates, Array, point, geometric center, barycenter,*)
(* affine transformation. *)*)
(**)
(*(*  Most of the concepts depend on the dimension dim, which must be defined in*)
(*the  Global Context! *)*)
(**)
(*(* Free software. Copyright.  This and the accompanying packages and notebooks and packages are  public. *)
(*Authors who intend to publish a changed or completed version of them should do it under their *)
(*own names, and under the condition that they cite the original package or notebook with the *)
(*Internet address or other source where they got it. I am not responsible for errors or damages *)
(*originated by the use of the procedures contained in my notebooks or packages; anybody who *)
(*applies them should test carefully whether they are appropriate for his purposes.*)*)
(**)
(*(*Warning. The package tensalgv3.m can not substitute the package euvec.m in *)
(*the item Spheres, MathSource, Enhancements/Geometry 0210-092. *)*)


Off[ToExpression::notstr]; 


(* ::Section:: *)
(*The usages *)


BeginPackage["tensalgv3`"]


(* ::Subsubsection::Closed:: *)
(*Functions already defined in the obsolete package vectorcalc.m *)


  (* dim::usage="dim is the dimension of the vector space under consideration; 
 it must be set within the Global Context." *)

stb::usage="stb[i, d:dim] is the i-th standard basis vector, 1 <= i <= d."

vec::usage="vec[x, n:dim] is the n-tuple of the coordinates x[i] of the vector vec[x].";

null::usage="null[d:dim] is the null vector of length d.";

tpl::usage = "tpl[n,k] is the List of k-tuples of the set {1,2,...,n}." ;

dotnorme::usage="dotnorme[b] yields the unit vector with respect to the Dot
  product in direction of b, b!=o.";

projpt::usage = "projpt[vx] is the affine point in the projective model of the 
affine space corresponding to the vector vx; if the first coordinate of vx 
equals 0, projpt returns ComplexInfinity."

projpt3D::usage = "projpt3D[a] is the Graphics primitive Point of the 
corresponding projective point corresponding to the 4-dimensional vector a."


randomv::usage = "randomv[n:dim,range:{-1,1}, type:Real] is a random n-vector 
with values of type in range; type and range like in Random[].
(New version: order of arguments has changed!)"

smoothing::usage = "smoothing[a, b] is a listable function returning zero if 
Abs[a]< b: 10^-10, and a else; for b = 1 a is returned: smoothing[a, 1] = a."

outzero::usage = "outzero[b_List,dd] is the sublist of the vectors b[[i]] with 
Abs[b[[i]].b[[i]]] > dd; the default value of dd is 10^(-20)."

basis::usage = "basis[v] yields a sublist of v being a basis of the span of 
  the vectors v[[i]] of v; v must be a matrix."

noprops::usage = "noprops[c] deletes rows c[[j]] of the matrix c which are 
proportional to an element c[[i]] with i<j, and returns the matrix 
originating this way."

matrix::usage = "matrix[b,n:dim] is a square matrix of order n with elements 
b[i,j]."

nullmatrix::usage = "nullmatrix[d] is the square nullmatrix of order d; the 
default value of d is dim."

trace::usage = "trace[a] yields the trace of the matrix a."

blf::usage = "blf[a, d:dim] is a symmetric matrix of order d with elements
a[i,j], often used as matrix of a bilinear form."

randommatrix::usage = "randommatrix[k,n:dim,range:{-1,1},type_:Real]is a matrix with k 
rows and n columns, with random elements of type, within range. (New version: 
order of arguments has changed!)"

rank::usage="rank[m] calculates the rank of the rectangular matrix m of a not
   too large order (< 10000).rank is obsolete now; use the built-in MatrixRank instead."

trp::usage = "trp[i,j,n:dim] is the unit matrix of order n where the i-th and 
the j-th row are transposed." 

renorm::usage = 
  "renorm[b] for a given matrix b  is the matrix with  normed vectors 
b[[i]]/|b[[i]]|, where |x| denotes the usual Euclidean norm of the vector x. 
Before norming the null vectors of b are deleted.";

cr::usage = "cr[x1,x2,x3,x4] is the cross ratio of the elements x1,...,x4 who are supposed 
to belong to a division algebra (a field,in particular).";


(* ::Subsubsection::Closed:: *)
(* New functions *)


andtest::usage = "If v is a list of Booleans (i.e. values False or True), andtest[v] returns True,
iff all the members of the list v are True.";

ortest::usage = "If v is a list of Booleans (i.e. values False or True), ortest[v] returns True,
iff at least one of  the members of the list v is True.";

lts::usage = "If h is a List,lts[h] transforms h into the corresponding Sequence of arguments 
and else it returns Null with an error message.";

cc::usage = "The double test function cc[i,j][k,l]  yields 1 iff i==k 
and j==l, and 0 else.";

vector::usage="vector[x, n:dim] is as vec[x, n:dim]  prepended with the number 0 modelling the 
translation with vec[x] in the affine point space.";

point::usage="point[x, n:dim] is as vec[x, n:dim]  prepended with the number 1 modelling the 
points with coordinates vec[x] in the affine point space.";

cgeo::usage = "cgeo[p1,p2,...,pk] is the affine center of the points p1,...,pk.";

cbary::usage = "cbary[p1,..,pk][m1,...,mk]yields the barycenter of the mass
 points (p1,m1),...,(pk,mk) if m1+...+mk \[NotEqual] 0, and an error message else.";

afftr::usage = "afftr[a,b,d:dim] is the matrix of the affine transformation generated 
by the linear isomorphism defined by the matrix a and the translation defined 
by the vector b.";

covector::usage="covector[a, n:dim] is the same as covec, it is the n-tuple of the 
coordinates a[i] of covector[a] with respect to the dual standard basis; they are transformed 
covariantly under a linear transformation of the basis.";

covec::usage = "covec[a,n:dim] is the List of the coordinates of a covector with respect to the 
dual standard basis."

linform::usage = "linform[a][vlist] = lf[a][v] is the linear function on the vector space 
of dimension d= Length[v] corresponding to covector[a].";

lf::usage = "lf[a][v] is the short form of linform[a][v]."

tensor::usage = "tt:=tensor[a,k,h,d:dim] defines a tensor of  type (k,h) with coordinate function a 
over a d-dimensional vector space V on a field. The coordinates are obtained by 
tt[{i1,...,ik},{j1,...jh}] as a[{i1,...,ik},{j1,...jh}]. they can be numbers, functions, symbols, but not Lists. 
For a tensor the index arguments {i...}, {j...} are assumed to be related to a given basis of 
the vector space V, the coordinates transform k times covariantly and h times contravariantly 
with respect to a basis transformation in V. These transformation rules are additional assumptions 
not implemented in the Mathematica object tensor. Thus it is possible to use the function tensor formally for 
objects depending on the basis having another transformation behavior, e.g. Christoffel symbols in differential geometry.";

type::usage = "type[tt] yields the coordinate function a, the type {k, h}, and the dimension d, if 
tt = tensor[a,k,h,d] is a tensor, and a failure message else.";

tstb::usage = "tstb[{i...},{j...},d:dim] is the standard basis tensor in the space of tensors of 
type (k, h) corresponding to the k-Tupel {i...} and the h-tuple {j...}, on the vector space of 
finite dimension d with default value d=dim over a field; the dimension dim of the underlying 
vector space must be defined globally."

zero::usage = "zero[i___] = 0 is the general null function.";

one::usage = "one[i___] = 1 is the constant function 1.";


times::usage = "times is an auxiliary function needed for the definition of the Product of tensors
with scalars."

symform::usage = "symform[a,p,d:dim] creates an Array whose elements
 a[[\!\(\*SubscriptBox[\(i\), \(1,  ... , \)]\)\!\(\*SubscriptBox[\(i\), \(p\)]\)]], \!\(\*SubscriptBox[\(i\), \(j\)]\)\[Element]Range[d],
are symmetric in these indices. Such an Array may serve as the coordinates of a symmetric multilinear   
form depending on p vectors of a d-dimensional vector space or as the Array of coordinates of 
a symmetric covariant or contravariant tensor of type {a,p,0,d}, resp. {a,0,p,d}.";

altform::usage = "altform[a,p,d:dim] creates an Array whose elements 
 a[[\!\(\*SubscriptBox[\(i\), \(1,  ... , \)]\)\!\(\*SubscriptBox[\(i\), \(p\)]\)]], \!\(\*SubscriptBox[\(i\), \(j\)]\) \[Element] Range[d],  are antisymmetric in these indices. Such an Array may serve as the coordinates of an antisymmetric 
multilinear form depending on p vectors of a d-dimensional vector space or as the Array of coordinates 
of an antisymmetric covariant or contravariant tensor of type {a,p,0,d}, resp. {a,0,p,d}.";

nulltensor::usage = "nulltensor[k,h,d:dim] defines the nulltensor of type (k,h) and dimension d.";

coords::usage = "coords[tt] yields the family of coordinates of the tensor tt as a matrix depending 
on the covariant and the contravariant index tuples .";

strictcoords::usage ="strictcoords[ta] yields the list of the coordinates of a covariant 
or contravariant antisymmetric tensor (or Array) ta of degree k over a d-dimensional vector space 
corresponding to the subsets of k elements of the set {1,...,d} called the strict coordinates of ta.";
cosymtensor::usage = "cosymtensor[a,p,d]creates a symmetric covariant tensor of type {a,p,0,d:dim}.";

contrasymtensor::usage = "contrasymtensor[a,p,d]creates a symmetric contravariant tensor of 
type {a,0,p,d:dim}.";

coalttensor::usage = "coalttensor[a,p,d] creates an antiymmetric covariant tensor of 
type {a,p,0,d:dim}.";

contraalttensor::usage = "contraalttensor[a,p,d] creates a symmetric contravariant tensor 
of type {a,0,p,d:dim}.";

equi::usage = "equi[tt,ts] returns True iff the corresponding types {k,h}, the dimensions d, and the 
coordinates of the tensors tt, ts are equal, and False if this is definitely not the case; 
it prints Investigate! and returns Null, if the decision is not found. If equi[tt,ts] returns True then 
ts equals tt since the corresponding coordinates coincide.";

tp::usage = "For two tensors ta, tb the function tp[ta,tp] is the coordinate function of the 
product tensor of ta,tb.";

tensorproduct::usage ="Given tensors ta,... of type (ka,ha), ..., the function tensorproduct[ta,...]  
(also tensorproduct[{ta,...}]), ta\[CircleTimes]tb, yields the tensor product of the tensors ta,.... 
The tensors ta,...must be related to the same vector space of given dimension d with 
default value dim."; 

contraction::usage = "contraction[tt,a,b] yields the contraction of the mixed tensor tt 
corresponding to the a-th covariant and the b-th contravariant slot of the indices. If slots is a
suitable list of slot pairs, then contraction[tt, slots] is the contracrion of tt over all 
these pairs."

covpermute::usage = "covpermute[a,perm][i, j] permutes the covariant indices of index list i with 
the permutation perm and leaves the kontravariant index list j unchanged.";

contrapermute::usage = "contrapermute[a,perm][i, j] permutes the contravariant indices of 
index list j with the permutation perm and leaves the covariant index list i unchanged.";

covs::usage = "covs[a][ilist, jlist] is the sum of the expressions a[Permute[i,pr],j]
over all permutations pr of ilist, jlist remains fixed.";

contras::usage = "contras[a][ilist, jlist] is the sum of the expressions a[i,Permute[j,pr]]  
over all permutations pr of jlist, ilist remains fixed.";

covsym::usage = "covsym[tt] yields the symmetrization of the tensor tt in the covariant indices." ;

contrasym::usage = "contrasym[tt] yields the symmetrization of the tensor tt in the 
contravariant indices.";

cova::usage = "cova[a,d:dim][ilist, jlist] is the alternating sum of the expressions 
a[Permute[i,pr],j]over all permutations pr of ilist, jlist remains fixed.";

contra::usage = "contra[a,d:dim][ilist, jlist] is the alternating sum of the expressions
 a[i, Permute[j,pr]] over all permutations pr of jlist, ilist remains fixed.";

covalt::usage = "covalt[tt] yields the alternation of the tensor tt in the covariant indices." ;

contralt::usage = "contralt[tt] yields the alternation of the tensor tt in the contravariant indices." ;

wedgeform::usage = "wedgeform[w1,w2,...] is the wedge product of the covectors w1,w2,... 
entered by their coordinate tuples or by a matrix whose rows are the coordinate tuples of the 
covectors.";

wedgevector::usage = "wedgevector[v1,v2,...] is the wedge product of the vectors v1,v2,... 
entered by their coordinate tuples or by a matrix whose rows are the coordinate tuples of the vectors.";

transpose::usage = "transpose[ta] is the tensor with coordinate function b[i,j] = a[j,i]; 
the type is (b,(h,k),d) if type[ta]=(a,(k,h),d).";

wedge::usage = "wedge[t1,...,tk] (or wedge[{t1,...,tk}]) is the antisymmetric product of covariant 
(resp. contravariant) tensors t1,...,tk. Besides of tensors the arguments ti may be numeric quantities 
serving as scalar factors. For antisymmetric tensors t1,...,tk wedge defines the product of the 
exterior algebra. Non-antisymmetric tensors will be alternated before the multiplication. 
In place of wedge the system symbol Wedge may be used; what permits to utilize the symbol \[Wedge], 
 which can be entered as  Esc\!\(\*
StyleBox[\"\[ThinSpace]\", \"InlineCode\"]\)\!\(\*
StyleBox[\"^\", \"InlineCode\"]\)\!\(\*
StyleBox[\"\[ThinSpace]\", \"InlineCode\"]\)Esc.";

Wedge::usage = "Wedge is defined to coincide with wedge; it permits the usage of the symbol \[Wedge]."

tensorbasis::usage = "tensorbasis[bb,k,h,d:dim] yields a formal tensor basis {bb[i,j]} in the space of 
tensors of type (k,h) over a vector space of dimension d with default value d = dim.
Setting bb =tstb yields the standard tensor basis.";

wtstb::usage = "wtstb[i,d:dim] is the antisymmetric standard basis form corresponding to the 
subset i \[Element] Subsets[Range[d]].";

vtstb::usage = "vtstb[i,d:dim]is the antisymmetric standard basis vector corresponding to the 
subset i \[Element] Subsets[Range[d]].";

ptind::usage = "ptind[\!\(\*SubsuperscriptBox[\(a\), \(j\), \(k\)]\)] = \!\(\*SubsuperscriptBox[\(a\), \(k\), \(j\)]\) transposes the index sets of a tensor expression.";

tradf::usage = "tradf[a[ilist,jlist]] returns the traditional form \!\(\*SuperscriptBox[SubscriptBox[\(a\), \(i\)], \(j\)]\) 
of the tensor coordinates. The function tradf gives Supersubscript output only for arguments of the shape just indicated, other arguments are returned unchanged.";

mathff::usage = "mathff[x] returns the Input form  a[i,j] of the tensor coordinates if 
x = \!\(\*SuperscriptBox[SubscriptBox[\(a\), \(i\)], \(j\)]\), 1 if x is empty, and x else.";

tensorToArray::usage = "If ta is a tensor of type (k,h) over a vector space of dimension d, 
tensorToArray[ta] yields the Array of TensorRank rk=Max[{1,k+h}] and TensorDimensions Table[d,{rk}] whose 
elements are the coordinates of ta.";

ArrayTotensor ::usage = "If ar is an  array of rank n and dimensions {d,d,...d} 
then ArrayTotensor[ar,k:0] generates a tensor of typ (k,n-k) over a d-dimensional vector  space
with coordinates a[(i1,...,ik},{i(k+1),...,in}] = ar[i1,...,ik,...,in];
the default value of k is zero."; 

plus::usage = "plus[f,..][x,...] (or plus[{f,g,...}][x,...]) is a function yielding the sum of 
the functions f,..., where the arguments f,... stand for the names of  functions or for  numbers 
being considered as the corresponding constant functions. All the functions must 
be defined in the same domain and must have values in the same set in which 
a commutative addition + is defined."; 

vecTotensor::usage = "vecTotensor[x] is the tensor of type (0,1) whose coordinates are 
the coordinates of x.";

covecTotensor::usage = "covecTotensor[x] is the tensor of type (1,0) whose coordinates 
are the coordinates of x.";

eps::usage = "If a, b are lists of positive Integers between 1 and d of equal length k, 
then eps[a,b,d:dim] is the corresponding coordinate of the antisymmetric epsilon tensor of 
order k and dimension d; if b is a permutation of a, then 
eps[a,b]=Signature[a]*Signature[b], and zero else.";

epst::usage = "epst[k,d:dim] is the epsilon tensor of type(k,k) in dimension d with coordinate
 function eps; it is antisymmetric with respect to its covariant and its contravariant indices." ;

extprod::usage = "extprod[w1,w2,...] yields the exterior product of antisymmetric forms or 
multivectors represented by their Arrays of coordinates. It avoids the concept of tensor."; 



(* ::Section:: *)
(*The Definitions*)


(* ::Subsubsection::Closed:: *)
(*Functions once defined in the package vectorcalc.m *)


Begin["`Private`"]

(*The dimension dim must be set within the Global Context after loading the package.*)

dimm:= Global`dim;

(* Argumentwise operations of functions *)

Attributes[plus] ={Flat,NumericFunction,OneIdentity,Orderless};
plus[f___][x__]:= Module[{ff,n,i,hh={}},ff=Flatten[{f}];n=Length[ff];
For[i=1,i<n+1,i++,If[NumericQ[N[ff[[i]]]],AppendTo[hh,ff[[i]]],AppendTo[hh,ff[[i]][x]]]];
Sum[hh[[i]],{i,1,n}]];

(* The next command admits the multiplication of a function f with a numeric quantity.*)

 Unprotect[Times];Times[h_?NumericQ,f_][x___]:=Times[h,f[x]];
Times[h_?NumericQ,f_][x___][i__]:=Times[h,f[x][i]];
Protect[Times];


(*A function which makes small entities to become zero.*)
  
  smoothing[a_, b_:10^(-10)]:=
  Module[{ c= a}, If[b == 1, Return[c]]; If[Abs[c]< b, c=0]; c];
Attributes[smoothing]=Listable;  

(*The standard basis.*)

stb[i_Integer,d_:dimm]:=IdentityMatrix[d][[i]]

(* A general vector *) 

vec[x_,n_:dimm]:=Array[x,n];

(* The null vector *)

null[n_:dimm]:=Table[0,{n}];

(* Norming a vector *)

dotnorme[b_]:=If[ b.b==0, Return[b], Return[b/Sqrt[b.b]],b/Sqrt[b.b]]

projpt[vx_] := Module[{v = vx.vx},
    If[Chop[v] == 0, Return["nopt"]];
    If[Chop[vx[[1]]]== 0, Return[ComplexInfinity]];
  vx/vx[[1]]]
    

projpt3D[a_] := Point[projpt[a]]

(* Random vectors *)

randomv[n_:dimm,range_:{-1,1}, type_:Real] := Table[Random[type,range],{i, n}]

(* Vector sequences *)

outzero[b_List, dd_:10^(-20)]:=
  Module[{v={},l=Length[b]},
    For[i=1,i<=l,i++,
      If[Not[Chop[b[[i]].b[[i]],dd]==0],AppendTo[v,b[[i]]], ,
        AppendTo[v,b[[i]]]]];Return[v]]

basis[v_]:=
  Module[{r,l,w,i,fdrop,new},
    If[Not[MatrixQ[v]],Print["Argument is not a matrix, no span!"];Return[]];
    r=rank[v];
    If[r==0,Print["The span is the null space!"];Return[]];
    w=v;
    fdrop[i_]:=Drop[w,{l-i}];
    Label[new];
    l=Length[w];
    If[r==l,Return[w]];
    For[i=0,i<=l-1,i++,
      If[rank[fdrop[i]]==r,w=fdrop[i];Goto[new]]]];

(* The function wedg is not exported; it is an auxiliary function for noprops only. *)

wedg[a_,b_]:= 
  Module[{d=Length[a]},Table[a[[i]]*b[[j]]-a[[j]]*b[[i]],{i,1,d},{j,1,d}]];

noprops[c___]:= Module[{m=1,cc, d,ttt,i,j},
    If[(Length[{c}]>1)||Not[MatrixQ[c]], Print["incorrect argument!"];Return[]];
    If[Length[c]<=1,Return[c]];
    d=Length[c[[1]]];
    cc=outzero[c];
    If[Length[cc]<=1,Return[cc]];
    Label[ttt];
    For[j=m+1,j<= Length[cc],j++, 
      If[Chop[wedg[cc[[m]],cc[[j]]],10^(-20)]== nullmatrix[d],
        cc[[j]]=null[d]]];
    cc=outzero[cc];m = m+1;
    If[Length[cc] <= m, Return[cc],Goto[ttt]]];

(* Matrices *)

matrix[b_,n_:dimm]:= Array[b,{n,n}];

nullmatrix[d_:dimm]:= Table[0,{i,1,d},{j,1,d}];

rank[m_] := MatrixRank[m];

trace[a_?MatrixQ] := Sum[a[[i,i]],{i,dimm}];
               
blf[a_,d_:dimm]:= Module[{b = matrix[a, d]},
    For[i = 1,i<d,i++,For[j=i+1,j<=d,j++,a[j,i]=a[i,j]]];
    b]

randommatrix[k_,n_:dimm,range_:{-1,1}, type_:Real] := Table[randomv[n,range,type], {k}]

trp[i_,j_,n_:dimm]:=Module[{m=IdentityMatrix[n]}, m[[i]]=stb[j,n];
    m[[j]]=stb[i,n];
    m]


(* The renorming procedure (useful for badly conditioned matrices) *)

renorm[b_?MatrixQ]:= Module[{bo = outzero[b],i },
	Table[dotnorme[bo[[i]]], {i, Length[bo]}]];

cr[x1_,x2_,x3_,x4_]:= (x1-x3)*((x1-x4)^(-1))*((x2-x3)*(x2-x4)^(-1))^(-1);


(* ::Subsubsection::Closed:: *)
(* New functions *)


zero[HoldPattern[___]] := 0;

one[x___] := 1


Unprotect[Null]; Null[___] := Null;  Protect[Null];

(* Test functions *)

andtest[v_List]:=Module[{i= 1,ww=True},
While[i<= Length[v],ww=And[v[[i]],ww];i++];ww];

ortest[v_List]:=Module[{i= 1,ww=False},
While[i<= Length[v],ww=Or[v[[i]],ww];i++];ww];

cc[i_, j_][k_, l_]:= If[(i==k)&& (j ==l),1,0,0];

lts[h__]:= If[Head[h]=== List,h/.List->Sequence, Print[h," is not a List!"];Null];

(* vector, covectors and corresponding linear forms. *)

vector[a_,d_:dimm] := Prepend[vec[a,d],0];

point[a_,d_:dimm] := Prepend[vec[a,d],1];

cgeo[p__] := Module[{m = Length[{p}], ft,i,k},
ft= {p};k=Length[ft[[1]]];  For [i =1, i<m+1,i++,
If[(ft[[i]][[1]]!= 1), 
Print["Invalid argument ",ft[[i]]," !"]; Return[]]];
 For [i =1, i<m+1,i++,If[(Length[ft[[i]]]!=k), 
Print["Dimensions of points differ!"]; Return[]]];
Sum[ft[[i]],{i,m}]/m];

cbary[p__][m__]:=  Module[{mm = Length[{p}], ft,i,M,k},
ft= {p};
k=Length[ft[[1]]];  For [i =1, i<mm+1,i++,
If[(ft[[i]][[1]]!= 1), 
Print["Invalid argument ",ft[[i]]," !"]; Return[]]];
 For [i =1, i<mm+1,i++,
If[(Length[ft[[i]]]!=k), 
Print["Dimensionsof points differ!"]; Return[]]];
If[Length[{m}]!= mm, Print["Masses incorrectly defined!"];Return[]];
M=Sum[{m}[[i]],{i,mm}];
If[M==0,
 Print["Sum[",{m}, " equals zero!"];Return[]];
Sum[{m}[[i]]*ft[[i]],{i,mm}]/M];

afftr[a_,b_,d_:dimm]:= Module[{aa,i,j},
If[(Head[b]==List)&&Dimensions[b]== {d},For[i=1,i<d+1,i++,aa[i_,0]:= b[[i]]],
Print["Invalid argment ",b," !"];Return[]];
If[Not[(Head[a]==List)&&(Dimensions[a]== {d,d})],
Print["Invalid argment ",a," !"];Return[]];For[i=1,i<d+1,i++,For[j=1,j<d+1,j++,aa[i,j]= a[[i,j]]]];
aa[0,0]=1;
For[j=1,j<d+1,j++,aa[0,j]=0];
Array[h,{d+1,d+1},{0,0}]/.h[i_,j_]->aa[i,j]];

covector[a_,n_:dimm]:=  Array[a,n];

covec[a_,n_:dimm]:= Array[a,n];

linform[a_][b_]:=  If[Head[b]==List,covec[a,Length[b]].b,
Print["Argument ", b, " incorrect!"];Return[Null],
Print["Argument ", b, " incorrect!"];Return[Null]];

lf[a_][b_]:= linform[a][b];

(* Multilinearforms *)

symform[a_,p_,d_:dimm]:=Module[{rr,i},
If[Not[IntegerQ[p]]||p<0,Print["Wrong p = ",p];Return[]];
If[p==0, Return[a]]; rr=SymmetrizedArray[{i__}->a[i],Table[d,{p}],{Symmetric[Range[p]]}];
Normal[rr]];

altform[a_,p_,d_:dimm]:=Module[{rr,i},
If[Not[IntegerQ[p]]||p<0,Print["Wrong p = ",p];Return[]];
If[p== 0, Return[a]];
rr=SymmetrizedArray[{i__}->a[i],Table[d,{p}],{Antisymmetric[Range[p]]}];
Normal[rr]];

(*tensors *)

tpl[n_,k_] := Tuples[Range[n],k];

tensor[a_, k_Integer, h_Integer, d_:dimm][i_,j_] := Module[{},
If[((Head[a]===List)||(Head[a[i,j]]=== List)),Print["Head[",a[i,j],"]= ",Head[a[i,j]], ": Invalid argument!"];Return[]];
If[MemberQ[Tuples[Range[d],k],i]&&MemberQ[Tuples[Range[d],h],j], 
Return[a[i,j]], Print["Incorrect indices!"];Return[], 
Print["Incorrect indices?"];Return[]]]/. s_[{},{}]-> s;


tensor[a_, k_Integer, h_Integer, d_:dimm][i_] := 
Print["A tensor needs two arguments being Lists!"];Return[];

tensor[a_, k_Integer, h_Integer, d_:dimm][] := 
Print["A tensor needs two arguments being Lists!"];Return[];

(* The nulltensor *)

nulltensor[k_Integer, h_Integer,d_:dimm] := tensor[zero,k,h,d];
(* nulltensor[k_Integer, h_Integer] := tensor[zero,k,h]; *)


Unprotect[Plus]; tensor[a_,k_,l_]+tensor[b_,k_,l_]:= tensor[plus[a,b],k,l]; 
tensor[a_,k_,l_]+tensor[b_,k_,l_,d_]:= If[ d== dimm,tensor[a,k,l]+tensor[b,k,l],
Print["Unequal dimensions!"]; Return[]]; 
tensor[a_,k_,l_,d_:dimm]+tensor[b_,k_,l_,d_:dimm]:= tensor[plus[a,b],k,l,d]; 
Protect[Plus];

Unprotect[Times]; -tensor[b_,k_,l_]:= tensor[-b,k,l];  
-tensor[b_,k_,l_,d_:dimm]:= tensor[-b,k,l,d]; Protect[Times];
Unprotect[Plus];tensor[a_,k_,l_]-tensor[b_,k_,l_]:= tensor[plus[a,-b],k,l];
tensor[a_,k_,l_,d_:dimm]-tensor[b_,k_,l_,d_:dimm]:= tensor[plus[a,-b],k,l,d];
Protect[Plus];


type[tt_] := Module[{typ, ty},
typ[tensor[a_,k_,h_,d_:dimm]] := {a,{k,h},d};
If[(Head[tt] === tensor), ty=typ[tt],  Print["Argument is not a tensor!"];ty=Null,
Print["Argument is not a tensor!"];ty=Null];ty];

(* type[tt_] := Module[{typ, a,k,h,d,ty},
typ[tensor[a_,k_,h_,d_:dimm]] := {a,{k,h},d};
If[(Head[tt] === tensor), ty=typ[tt],  Print["Argument is not a tensor!"];ty=Null,
Print["Argument is not a tensor!"];ty=Null];ty];
*)

(* Scalar times tensor *)
Attributes[times]= {Listable,NumericFunction,OneIdentity};
times[aa_, tensor[b_,k_,l_,d_:dimm]] := tensor[times[aa,b],k,l,d];
times[aa_, tensor[b_,k_,l_]] := tensor[times[aa,b],k,l];
times[xx_,aa_][i___, j___]:=Times[xx,aa[i, j]];

times/: times[0,aa_] := 0&;
times/: times[x_,0&] := 0&;


times/: times[-a_,tt_] := -times[a,tt];
times/: times[a_,-tt_] := -times[a,tt];
times/: times[a_ + b_,tt_] := times[a,tt] + times[b,tt];
times/: times[a_ ,tt_+ ts_] := times[a,tt] + times[a,ts];

Unprotect[Times];Times/: Times[x_, tensor[a_,k_,l_,d_:dimm]]:= tensor[times[x,a],k,l,d];
Times/: Times[x_, tensor[a_,k_,l_]]:=tensor[times[x,a],k,l];Protect[Times]; 




(* Standard tensor basis *)

tstb[i_,j_,d_:dimm] := tensor[cc[i,j],Length[i],Length[j],d];

SetAttributes[coords, Listable];
coords[tt_] := Module[{ti,k,h,dd,rk,rh, i,j},
If[Not[Head[tt]=== tensor],Print["The argument must be a tensor!"];Return[]];
ti = type[tt]; k = ti[[2,1]]; h =ti[[2,2]]; dd =ti[[3]];
rk = Tuples[Range[dd],k];rh = Tuples[Range[dd],h];
Table[ tt[rk[[i]],rh[[j]]],{i,1,dd^k},{j,1,dd^h}]];

cosymtensor[a_,p_,d_:dimm] :=If[p===0,Return[a],ArrayTotensor[symform[a,p,d],p]];

contrasymtensor[a_,p_,d_:dimm] := If[p===0,Return[a],ArrayTotensor[symform[a,p,d]]];

coalttensor[a_,p_,d_:dimm] :=If[p== 0,Return[tensor[a,0,0,d]], ArrayTotensor[altform[a,p,d],p]];

contraalttensor[a_,p_,d_:dimm] := If[p== 0,Return[tensor[a,0,0,d]],ArrayTotensor[altform[a,p,d]]];

strictcoords[ts_] := Module[{aa,ttp,taa, k,di,rdk,st,i},
If[Not[(Head[ts]=== tensor)||(Head[ts]=== List)],
Print["invalid argument!"];Return[]];
If[Head[ts]=== tensor,
ttp = type[ts];aa = ttp[[1]];k= ttp[[2,1]] +ttp[[2,2]];di=ttp[[3]];
If[k==0, Return[{aa[{},{}]}]];
If[Not[ttp[[2,1]]*ttp[[2,2]]==0], Print[ts," is a mixed tensor!"];Return[]]];
If[Head[ts]=== tensor,taa =tensorToArray[ts],taa=ts,
Print["Valid argument?"];Return[]];
If[Head[ts]=== List, k=TensorRank[ts];di=TensorDimensions[ts][[1]]];
If[Not[TensorDimensions[taa]== Table[di,{k}]], Print["Invalid argument!!"];Return[]];
If[k<=1, Return[taa]];
If[Not[(TensorSymmetry[taa] === Antisymmetric[Range[k]])
||(TensorSymmetry[taa] === ZeroSymmetric[{}])],
Print["Tensor ",ts," is not antisymmetric!"];Return[]];
If[(k>di)||( aa==zero),Return[{0}]];
rdk=Subsets[Range[di],{k}];
If[(Head[ts]=== tensor)&&(ttp[[2,1]]==0), Return[Table[aa[{},rdk[[i]]],{i,Binomial[di,k]}]]];
If[(Head[ts]=== tensor)&&(ttp[[2,2]]==0), Return[Table[aa[rdk[[i]],{}],{i,Binomial[di,k]}]]];
Table[taa[[rdk[[i]]/.List -> Sequence]],{i,Binomial[di,k]}]]


equi[tt_,ts_ ]:= Module[{kt,ks,n,i,ttt,tts},
If [Not[(Head[tt]===  tensor)&&(Head[ts]===  tensor)], 
Print["Arguments must be tensors!"];Return[]];
kt =type[tt]; ks = type[ts];
If[Not[(kt[[2]]===  ks[[2]])&&(kt[[3]]===  ks[[3]])],Return[False]];
ttt=Flatten[coords[tt]];tts=Flatten[coords[ts]];
n=Length[ttt];
i=1;While[i<n+1, If[Not[Simplify[ttt[[i]]==tts[[i]]]],
Print[ttt[[i]]," = ",tts[[i]], " is "];Return[False],i++,
Print["Investigate ",ttt[[i]], " = ", tts[[i]]," ?"];Return[]]];
True];



tensorbasis[bb_,k_,h_,d_:dimm] :=Module[
{rk = Tuples[Range[d],k],rh = Tuples[Range[d],h],i,j},
Flatten[Table[ bb[rk[[i]],rh[[j]]],{i,1,d^k},{j,1,d^h}]]];


(* Tensor product *)

tp[ta_tensor,tb_tensor][i_List, j_List]:=Module[{ka,kb,ha,hb},
ka =type[ta][[2,1]];ha =type[ta][[2,2]];
kb =type[tb][[2,1]];hb =type[tb][[2,2]];
ta[Part[i,1;;ka],Part[j,1;;ha]]*tb[Part[i,ka+1;;ka+kb],Part[j,ha+1;;ha+hb]]];

tenpro[ta_,tb_] := Module[
{ka,kb,ha,hb,dd},
If[Not[(Head[ta]===  tensor)&&(Head[tb]===  tensor)], Print["Incorrect arguments!"];Return[],
ka=type[ta][[2,1]];ha=type[ta][[2,2]];dd=type[ta][[3]];
kb=type[tb][[2,1]];hb=type[tb][[2,2]];
If[ta == nulltensor[ka,ha,dd], Return[nulltensor[ka+kb,ha+hb,dd]]];
If[tb == nulltensor[kb,hb,dd], Return[nulltensor[ka+kb,ha+hb,dd]]];
If[(ka+kb==0)&&(ha+hb==0), Return[tensor[type[ta][[1]]*type[tb][[1]],0,0,dd]]];
tensor[tp[ta,tb],ka+kb,ha+hb,dd]]
];

Attributes[tensorproduct] ={Flat, HoldAll};
(* tensorproduct[ta___]:= tensorproduct [{ta}]; *)
tensorproduct[ta___]:=Module[{ttt, tt,k,i,dd,ee},
ttt=Flatten[{ta}]; k = Length[ttt];
If[k==0, Return[1]];
ee=0;For[i=1,i<k+1,i++,
If[Not[Head[ttt[[i]]]=== tensor],ee=ee+1]];
If[ee>0, Print["Incorrect arguments!"];Return[]];
tt = ttt[[1]]; dd = type[tt][[3]]; 
If[k== 1, Return[tt]];
ee=0; For[i=1,i<k+1,i++,If[Not[type[ttt[[i]]][[3]]=== dd],ee=ee+1]];
If[ee>0,Print["Different dimensions!"];Return[]];
i=2;While[i<k+1,tt =tenpro[tt,ttt[[i]]];i++];
tt];

xx_\[CircleTimes]yy_ := tensorproduct[xx,yy];



(* contraction *)

tract[tt_,i_,j_][x_,y_] := Module[{ty,h,k,d,aa, bb,cov,cont,x1,y1,jj,xx,yy},
If[Head[tt] === tensor, ty = type[tt], Return[type[tt]],Return[type[tt]]];
aa = ty[[1]];h=ty[[2,1]]; k= ty[[2,2]]; d=ty[[3]];
cov =Tuples[Range[d],h-1];cont =Tuples[Range[d],k-1];
If[Not[MemberQ[cov,x]&&MemberQ[cont,y]],Print["Incorrect Arguments!"];Return[]];
x1[jj_] :=Join[x[[1;;i-1]],{jj},x[[i;;h-1]]];
y1[jj_] :=Join[y[[1;;j-1]],{jj},y[[j;;k-1]]];
bb[xx_,yy_]:=Sum[aa[x1[jj],y1[jj]],{jj,1,d}];
bb[x,y]];

contractio[tt_,i_,j_] := Module[{ty},If[Head[tt] === tensor, 
ty = type[tt], Return[type[tt]],Return[type[tt]]]; ca = tract[tt,i,j]; 
tensor[ca, type[tt][[2,1]]-1,type[tt][[2,2]]-1, type[tt][[3]]]];

contraction[tt_, slots_]:=Module[{ty, k,h,aa,d,
nn, eps = 0,err,vv,  vi,vj,i,w,m,j,ww,ta,sa},
If[Head[tt]=== tensor, ty =type[tt],Print["First argument must be a tensor"];Return[]];
aa=ty[[1]];h =ty[[2,1]];k =ty[[2,2]];d =ty[[3]];
If[MatrixQ[slots], nn= Length[slots],eps =1;Goto[err]];
If[Not[Length[slots[[1]]] == 2],eps =1;Goto[err]];
vv =Reverse[ Sort[slots]];
vi =Table[vv[[i,1]],{i,nn}];vj= Table[vv[[i,2]],{i,nn}];
If[Not[andtest[Map[IntegerQ,Join[vi,vj]]]],eps =1;Goto[err]];
If[(Min[vi]<1)||(Min[vj]<1),eps =1;Goto[err]];
If[(Max[vi]>h)||(Max[vj]>k),eps =1;Goto[err]];
If[(Length[Union[vi]]< nn)||(Length[Union[vj]]< nn),eps =1;Goto[err]];
Label[err];
If[eps== 1, Print["slots:Arguments incorrect!"];Return[]];
w[i_,1] := vv[[i,1]];w[1,2] = vv[[1,2]];
For[ i =2,i<= nn,i++,
For[m=0;j=1,j<i,j++,If[vv[[j,2]]<vv[[i,2]],m =m+1]];
w[i,2]=vv[[i,2]]-m];ww=Array[w,{nn,2}];
ta = tt;i = 1;While[i<nn+1,sa = contractio[ta,w[i,1],w[i,2]]; ta =sa;i++];
sa];

contraction[tt_,i_Integer,j_Integer] := contraction[tt,{{i,j}}]


(* Symmetry properties of tensors *)

covpermute[a_,pe_List][i_List, j_List]:= a[Permute[i,pe], j];

contrapermute[a_,pe_List][i_List, j_List]:= a[i, Permute[j,pe]]; 

covs[aa_][ii_, jj_]:=If[(Length[ii]==0), aa[ii, jj],
Sum[covpermute[aa,Permutations[Range[Length[ii]]][[l]]][ii, jj],{l,Length[ii]!}]];
covs/: covs[a_ +b_] := covs[a] + covs[b];
covs/:  covs[times[c_,a_]] := times[c,covs[a]];

contras[aa_][ii_, jj_]:=If[(Length[jj]==0), aa[ii, jj],
Sum[contrapermute[aa,Permutations[Range[Length[jj]]][[l]]][ii, jj],{l,Length[jj]!}]];
contras/: contras[a_ +b_] := contras[a] + contras[b];
contras/: contras[times[c_,a_]] := times[c,contras[a]];

covsym[tt_]:= Module[{ttp, k,h,dd,aa},
If[Not[Head[tt] ===  tensor], Print["Argument must be a tensor!"];Return[]];
ttp = type[tt]; k =ttp[[2,1]];h =ttp[[2,2]];dd =ttp[[3]]; aa = ttp[[1]];
If[(k==0), Return[tt],Return[tensor[covs[aa],k,h,dd]]]];

contrasym[tt_]:= Module[{ttp, k,h,dd,aa},
If[Not[Head[tt] ===  tensor], Print["Argument must be a tensor!"];Return[]];
ttp = type[tt]; k =ttp[[2,1]];h =ttp[[2,2]];dd =ttp[[3]];aa = ttp[[1]];
If[(h==0), tt,tensor[contras[aa],k,h,dd]]];


cova[aa_,d_:dimm][ii_,jj_]:=If[(Length[ii]==0), aa[ii,jj],
If[MemberQ[Subsets[Range[d],{Length[ii]}],Sort[ii]],
Return[Sum[Signature[Permutations[Range[Length[ii]]][[l]]]*covpermute[aa,
Permutations[Range[Length[ii]]][[l]]][ii,jj],{l,Length[ii]!}]],0]];
cova/: cova[a_ + b_,d_:dimm] := cova[a,d] + cova[b,d];
cova/: cova[times[c_,a_],d_:dimm] := times[c,cova[a,d]];

contra[aa_,d_:dimm][ii_,jj_]:=If[(Length[jj]==0), aa[ii,jj],
If[MemberQ[Subsets[Range[d],{Length[jj]}],Sort[jj]],
Return[Sum[Signature[Permutations[Range[Length[jj]]][[l]]]*contrapermute[aa,Permutations[Range[Length[jj]]][[l]]][ii,jj],{l,Length[jj]!}]],0]];
contra/: contra[a_ +b_,d_:dimm] := contra[a,d] + contra[b,d];
contra/: contra[times[c_,a_],d_:dimm] := times[c,contra[a,d]];


covalt[tt_]:= Module[{ttp, k,h,dd,aa},
If[Not[Head[tt] ===  tensor], Print["Argument must be a tensor!"];Return[]];
ttp = type[tt]; k =ttp[[2,1]];h =ttp[[2,2]]; dd =ttp[[3]];
 aa = ttp[[1]];
If[k>dd, Return[nulltensor[k,h,dd]]];
If[(k<=1), tt,tensor[cova[aa,dd],k,h,dd]]];
covalt/: covalt[times[c_,a_]] := times[c,covalt[a]];

contralt[tt_]:= Module[{ttp, k,h,dd,aa},
If[Not[Head[tt] ===  tensor], Print["Argument must be a tensor!"];Return[]];
ttp = type[tt]; k =ttp[[2,1]];h =ttp[[2,2]];dd =ttp[[3]];
 aa = ttp[[1]];
If[h>dd,Return[nulltensor[k,h,dd]]];
If[(h<= 1), tt,tensor[contra[aa,dd],k,h,dd]]];  
contralt/: contralt[times[c_,a_]] := times[c,contralt[a]];


(* Wedge product *)
Attributes[extprod] = {Flat,OneIdentity};
extprod[w___] := Normal[TensorWedge[w]];

wedgeform[x___] := Module[{xx,rk,h,j,k,tv,tvv,talt,d},
If[Length[{x}]==0,Return[1]];
If[(Length[{x}]==1) && Not[ListQ[x]],Return[x]];If[(Length[{x}]==1) && VectorQ[x],
Return[covecTotensor[x]]];
If[(Length[{x}]==1) && MatrixQ[x],xx=x];
If[(Length[{x}]>1),xx={x}];
If[Not[MatrixQ[xx]], Print["Invalid Argument!"];Return[]];
rk = MatrixRank[xx];h=Length[xx];d = Length[xx[[1]]];
If[rk<h, Return[nulltensor[h,0,d]]];
tv[j_] := covecTotensor[xx[[j]]];tvv := Table[tv[j],{j,h}];
talt :=covalt[tensorproduct[tvv/.List-> Sequence]];Return[talt]];

transpose[tt_]:= Module[{ttp,aa,bb,i,j},
If[Not[Head[tt] === tensor], Return[tt]];
ttp = type[tt];aa=ttp[[1]];
If[aa === zero, Return[nulltensor[ttp[[2,2]], ttp[[2,1]],ttp[[3]]]]];
 bb[i_,j_ ] := aa[j,i];
tensor[bb,ttp[[2,2]], ttp[[2,1]],ttp[[3]]]];

wedgevector[y___] := transpose[wedgeform[y]];

Attributes[wedge]= {Flat, OneIdentity};
wedge[x___] := Module[{ii,k,kk,j,ko,dd, kc,CON,num,tlist},
ii=Flatten[{x}];k = Length[ii];
If[k== 0,Return[1]];
For[tlist={};num=1;j=1,j<k+1,j++,
If[Not[NumericQ[ii[[j]]]||(Head[ii[[j]]]=== tensor)],
Print["Invalid Argument ", ii[[j]]," !"];Return[],
If[NumericQ[ii[[j]]],num =num*ii[[j]],
AppendTo[tlist,ii[[j]]]]]];
 k = Length[tlist];
ko = Sum[type[tlist[[j]]][[2,1]],{j,k}];
dd=type[tlist[[1]]][[3]];
kc = Sum[type[tlist[[j]]][[2,2]],{j,k}];
If[ko>0&&kc>0,Print["Invalid: Mixed type arguments!"];Return[]];
For[j=1,j<k+1,j++,If[Not[type[tlist[[j]]][[3]]===dd],
Print["Arguments have varying dimensions!"];Return[]]];
If[ko==0, Goto[CON]];
kk =Product[1/((type[tlist[[j]]][[2,1]])!),{j,k}];
If[ko>dd, Return[nulltensor[ko,0,dd]],
Return[num*kk*covalt[tensorproduct[tlist]]]];
Label[CON];
kk =Product[1/((type[tlist[[j]]][[2,2]])!),{j,k}];
If[kc>dd, Return[nulltensor[0,kc,dd]],
num*kk*contralt[tensorproduct[tlist]]]];

Wedge[x___] := wedge[x];

wtstb[i_,d_:dimm]:= Module[{p,pp,ppi,j},
If[Head[i]===List,p = Length[i],Print["Invalid argument!"];Return[]];
If[p==1, Return[tstb[i,{},d]]]; 
If[(p<=d),pp = Subsets[Range[d],{p}],Print["Invalid argument!"];Return[]];
If[MemberQ[pp,i],ppi =Permutations[Range[p]],Print["Invalid argument!"];Return[]];
Sum[Signature[ppi[[j]]]*tstb[Permute[i,ppi[[j]]],{},d],{j,p!}]];

vtstb[i_,d_:dimm]:= Module[{p,pp,ppi,j},
If[Head[i]===List,p = Length[i],Print["Invalid argument!"];Return[]];
If[p==1, Return[tstb[{},i,d]]]; 
If[(p<=d),pp = Subsets[Range[d],{p}],Print["Invalid argument!"];Return[]];
If[MemberQ[pp,i],ppi =Permutations[Range[p]],Print["Invalid argument!"];Return[]];
Sum[Signature[ppi[[j]]]*tstb[{},Permute[i,ppi[[j]]],d],{j,p!}]];


(* Traditional form *)

trad[a_[{i___},{j___}]]:=Module[{ii = ToString[{i}],jj = ToString[{j}]},
Superscript[Subscript[a,ii],jj]]

 SetAttributes[tradf, Listable]; 
tradf[x___]:=Module[{xx={x},zz,i,j,a},
If[Length[xx] ==  0,Return[1]];
If[(Head[x]===Symbol)|| NumberQ[x], Return[x]];
zz=Cases[xx,a_[{i___},{j___}]];
If[Length[zz]==0, Return[x]];
If[Length[zz]== 1, Return[\!\(\*
TagBox[
RowBox[{"trad", "[", 
RowBox[{"zz", "[", 
RowBox[{"[", "1", "]"}], "]"}], "]"}],
FullForm]\)],
Print["Invalid argument!"];Return[],
Print["Invalid argument!"];Return[]]]

tradf/: tradf[x_+y_] := tradf[x]+tradf[y]; 

tradf/: tradf[z_*x__]:= tradf[z]*tradf[x];

tradf/: tradf[z_/x__]:= tradf[z]/tradf[x];

(* SetAttributes[mathf, Listable]; *)

mathf[Superscript[Subscript[a_,i_],j_]] :=a[ToExpression[i,InputForm],ToExpression[j,InputForm]];

SetAttributes[mathff, Listable]; 
mathff[x___]:= Module[{n=Length[{x}],xx},
If[n==0, Return[1]];
If[Not[MatchQ[x,Superscript[Subscript[a_, i_],j_]]], Return[x], xx = mathf[x],xx=mathf[x]];
xx]

mathff/: mathff[x_+y_] := mathff[x]+mathff[y]; 

mathff/: mathff[x_*y_]:= mathff[x]*mathff[y];

mathff/: mathff[z_*x__]:= mathff[z]*mathff[x];

mathff/: mathff[z_/x__]:= mathff[z]/mathff[x];

ptind[Superscript[Subscript[a_, i_],j_]] := Superscript[Subscript[a, j],i]; 



(* tensor to Array *)
SetAttributes[tensorToArray,Listable];
tensorToArray[tt_]:= Module[{ti,aa,k,h,dd,bb,i},
If[Not[Head[tt]=== tensor],Print["The argument must be a tensor!"];Return[]];
ti = type[tt];aa = ti[[1]]; k = ti[[2,1]]; h =ti[[2,2]]; dd =ti[[3]];
If[k==0&&h==0, Return[{coords[tt][[1,1]]}]];
Off[Part::take];
 bb[i___]:= aa[{i}[[1;;k]],{i}[[k+1;;k+h]]];
Array[bb,Table[dd,{k+h}]]];

vecTotensor[x_] := Module[{k,xx,i},
If[Not[VectorQ[x]], Print["Argument ", x, " isn't a vector!"]; Return[]];
k= Length[x];xx[{},{i_} ]:= x[[i]];tensor[xx,0,1,k]];

covecTotensor[x_] := Module[{k,xx,i},
If[Not[VectorQ[x]], Print["Argument ", x, " isn't a covector!"]; Return[]];
k= Length[x];xx[{i_},{} ]:= x[[i]];tensor[xx,1,0,k]];

ArrayTotensor[a_,k_:0] := Module[{rg,dd ,tt,ii,jj,b},
If[ArrayQ[a],rg = TensorRank[a],
Print["Argument must be an Array!"]; Return[],
Print["Argument must be an Array!"]; Return[]];
If[(k<0)||(k>rg)||Not[IntegerQ[k]], Print[k ," bad argument!"]; Return[]];
dd=Length[a[[1]]];If[dd==0,b[{},{}] =a[[1]];Return[tensor[b,0,0,1]]];
If[VectorQ[a]&&(k== 1),Return[covecTotensor[a]]];
If[VectorQ[a]&&(k== 0),Return[vecTotensor[a]]];
If[Not[TensorDimensions[a] ==Table[dd,{rg}]], 
Print["Array does't fit!"];Return[]];
tt = tensor[b,k,rg-k,dd];
b[ii_,jj_] := a[[Join[ii,jj]/.List-> Sequence]];
tt];

(* The epsilon tensor *)

eps[a_List,b_List,d_:dimm]:= Module[{k=Length[a],m=Length[b],rg,ee=0},
If[((k>d)||(m>d)||(k !=m)), Print["False arguments!"];Return[]];
rg = Tuples[Range[d],k];
If[Not[MemberQ[rg,Sort[a]]]||Not[MemberQ[rg,Sort[b]]], Print["False arguments!"];Return[]];
If[Sort[a] == Sort[b],ee =Signature[a]*Signature[b]];ee];

epst[k_,d_:dimm]:=Module[{tt},
If[(k<0)||Not[IntegerQ[k]],Print["incorrect arguments!"];Return[]];
If[k<=d,tt= tensor[eps,k,k,d],tt=nulltensor[k,k,d]];tt];

<<Notation`;
Notation[ParsedBoxWrapper[RowBox[{"a_", "\[CircleTimes]", "b_"}]] \[DoubleLongLeftRightArrow] ParsedBoxWrapper[RowBox[{"tensorproduct", "[", RowBox[{"a_", ",", "b_"}], "]"}]]];
(*Notation[ParsedBoxWrapper[RowBox[{"a_", "\[Wedge]", "b_"}]] \[DoubleLongLeftRightArrow] ParsedBoxWrapper[RowBox[{"extprod", "[", RowBox[{"a_", ",", "b_"}], "]"}]]]; *)
Notation[ParsedBoxWrapper[RowBox[{"a_", "\[Wedge]", "b_"}]] \[DoubleLongLeftRightArrow] ParsedBoxWrapper[RowBox[{"wedge", "[", RowBox[{"a_", ",", "b_"}], "]"}]]];




(* ::Section:: *)
(*End*)


End[]

Protect["`*"];

EndPackage[];




