working on permutation

permutation.cc

@@ -163,12 +163,19 @@ return q.inverse()*((*this)*q);
 }
 
 
+//NOTE: for larger permutations it might be more efficient to use cycle decomposition but where exactly the breakeven is?
+//at the moment we do it the trivial way for small permutations and using sort for larger ones
 template <typename T>
 int NRPerm<T>::parity() const
 {
 if(!this->is_valid()) return 0;
 T n=this->size();
 if(n==1) return 1;
+if(n>=10) //???
+	{
+	NRPerm<T> tmp(*this);
+	return (tmp.sort()&1) ? -1:1;
+	}
 T count=0;
 for(T i=2;i<=n;i++) for(T j=1;j<i;j++) if((*this)[j]>(*this)[i]) count++;
 return (count&1)? -1:1;
@@ -313,6 +320,17 @@ return ret;
 }
 
 
+template <typename T, typename R>
+bool PermutationAlgebra<T,R>::operator==(PermutationAlgebra<T,R> &rhs)
+{
+simplify();
+rhs.simplify();
+if(size()!=rhs.size()) return false;
+for(int i=0; i<size(); ++i) if((*this)[i].weight!=rhs[i].weight || (*this)[i].perm!=rhs[i].perm) return false;
+return true;
+}
+
+
 
 //Algorithm L2 from Knuth's volume 4 for multiset permutations
 //input must be initialized with (repeated) numbers in any order
@@ -755,6 +773,7 @@ PermutationAlgebra<T,R> res(this->size()*rhs.size());
 for(int i=0; i<this->size(); ++i)
         for(int j=0; j<rhs.size(); ++j)
                 res[i*rhs.size()+j] = (*this)[i]*rhs[j];
+res.simplify();
 return res;
 }
 
@@ -762,17 +781,12 @@ template <typename T, typename R>
 PermutationAlgebra<T,R> PermutationAlgebra<T,R>::operator+(const PermutationAlgebra<T,R> &rhs) const
 {
 if(this->size()>0 && rhs.size()>0 && (*this)[0].perm.size()!=rhs[0].perm.size()) laerror("permutation sizes do not match");
-return this->concat(rhs);
+PermutationAlgebra<T,R> res=this->concat(rhs);
+res.simplify();
+return res;
 }
 
 
-template <typename T, typename R>
-void PermutationAlgebra<T,R>::simplify()
-{
-copyonwrite();
-sort();
-//@@@@implement merging weights of identical permutations and resize like we did in qubithamiltonian with other stuff
-}
 
 
 

permutation.h

@@ -107,6 +107,9 @@ public:
 	NRPerm<T> perm;
 	
 	int size() const {return perm.size();};
+	bool is_zero() const {return weight==0;}
+	bool is_scaledidentity() const {return perm.is_identity();}
+	bool is_identity() const {return weight==1 && is_scaledidentity();}
 	bool is_plaindata() const {return false;};
 	WeightPermutation() : weight(0) {};
 	WeightPermutation(const R w, const NRPerm<T> &p) : weight(w), perm(p) {};
@@ -116,7 +119,6 @@ public:
 	WeightPermutation  operator|(const WeightPermutation &rhs) const {return WeightPermutation(weight*rhs.weight,perm|rhs.perm);};
 	WeightPermutation  operator*(const WeightPermutation &rhs) const {return WeightPermutation(weight*rhs.weight,perm*rhs.perm);};
 	WeightPermutation operator*(const R &x) const {return WeightPermutation(weight*x,perm); }
-//@@@operator+ and - yielding Permutationalgebra
 
 	bool operator==(const WeightPermutation &rhs) const {return this->perm == rhs.perm;}; //NOTE for sorting, compares only the permutation not the weight!
 	bool operator>(const WeightPermutation &rhs) const {return  this->perm > rhs.perm;};
@@ -135,9 +137,13 @@ public:
         static bool is_plaindata() {return false;};
         static void copyonwrite(WeightPermutation<T,R>& x) {x.perm.copyonwrite();};
 	typedef typename LA_traits<R>::normtype  normtype;
+	typedef R coefficienttype;
+	typedef NRPerm<T> elementtype;
 	static inline bool smaller(const WeightPermutation<T,R>& x, const WeightPermutation<T,R>& y) {return x.perm<y.perm;};
 	static inline bool bigger(const WeightPermutation<T,R>& x, const WeightPermutation<T,R>& y) {return x.perm>y.perm;};
-	
+	static R coefficient(const WeightPermutation<T,R>& x){return x.weight;};
+	static R& coefficient(WeightPermutation<T,R>& x) {return x.weight;};
+	static typename LA_traits<R>::normtype abscoefficient(const WeightPermutation<T,R>& x){return LA_traits<R>::abs2(x.weight);};
 };
 
 
@@ -168,7 +174,7 @@ public:
 	PermutationAlgebra(const NRVec<WeightPermutation<T,R> > &x) :  NRVec<WeightPermutation<T,R> >(x) {};
 	int size() const {return NRVec<WeightPermutation<T,R> >::size();};
 	void copyonwrite() {NRVec<WeightPermutation<T,R> >::copyonwrite();};
-	void sort(int direction = 0, int from = 0, int to = -1, int *permut = NULL, bool stable=false) {NRVec<WeightPermutation<T,R> >::sort(direction,from, to,permut,stable);};
+	int sort(int direction = 0, int from = 0, int to = -1, int *permut = NULL, bool stable=false) {return NRVec<WeightPermutation<T,R> >::sort(direction,from, to,permut,stable);};
 
 	PermutationAlgebra operator&(const PermutationAlgebra &rhs) const;
 	PermutationAlgebra operator|(const PermutationAlgebra &rhs) const;
@@ -176,10 +182,13 @@ public:
 	PermutationAlgebra operator+(const PermutationAlgebra &rhs) const;
 	PermutationAlgebra &operator*=(const R &x) {this->copyonwrite(); for(int i=1; i<=size(); ++i) (*this)[i].weight *= x; return *this;};
 	PermutationAlgebra operator*(const R &x) const {PermutationAlgebra r(*this); return r*=x;};
-	void simplify(); 
+	void simplify(const typename LA_traits<R>::normtype thr=0) {NRVec_simplify(*this,thr);};
+	bool operator==(PermutationAlgebra &rhs); //do NOT inherit from NRVec, as the underlying one ignores weights for the simplification; also we have to simplify before comparison
+	bool is_zero() const {return size()==0;}; //assume it was simplified
+	bool is_scaled_identity() const {return size()==1 && (*this)[0]. is_scaledidentity();}; //assume it was simplified
+	bool is_identity() const {return size()==1 && (*this)[0]. is_identity();}; //assume it was simplified
 };
 
-//TODO@@@ also simplification after operation
 
 //TODO@@@ permutationalgebra  for Kucharski style antisymmetrizers
 
@@ -343,6 +352,8 @@ NRMat_from1<T> chi; //characters
 
 	Sn_characters(const int n0); //compute the table
 	bool is_valid() const; //check internal consistency
+	T irrepdim(T i) const {return chi(i,1);};
+	T sumirrepdims() const {T s=0; for(T i=1; i<=chi.nrows(); ++i) s+=irrepdim(i); return s;};
 };
 
 template <typename T>  class Polynomial; //forward declaration

t.cc

@@ -90,28 +90,38 @@ cout<<p;
 
 static int unitary_n;
 static PERM_RANK_TYPE space_dim;
+static NRVec<PermutationAlgebra<int,int> > allyoung;
+static NRVec<int> allyoung_irrep;
+int current_irrep;
+int allyoung_index;
 
 void yprintme(const YoungTableaux<int>&y)
 {
 cout <<y;
 if(!y.is_standard()) laerror("internal error in young");
-cout << y.young_operator();
+allyoung[allyoung_index] = y.young_operator();
+cout <<"Young "<<allyoung_index<<" (irrep "<<current_irrep<<") = "<<allyoung[allyoung_index]<<endl;
+allyoung_irrep[allyoung_index]=current_irrep;
+allyoung_index++;
 }
 
 void pprintme(const Partition<int> &p)
 {
 CompressedPartition<int> pc(p);
+++current_irrep;
 cout<<'['<<pc<<"]\n";
 
+
+PERM_RANK_TYPE snd=p.Sn_irrep_dim();
+cout<<"Sn IR dim "<<snd<<endl;
+
 YoungTableaux<int> y(p);
 PERM_RANK_TYPE dim=y.generate_all_standard(yprintme);
 
 Partition<int> q=p.adjoint();
-PERM_RANK_TYPE snd=p.Sn_irrep_dim();
-cout<<"IR dim "<<snd<<endl;
 if(dim!=snd) laerror("inconsistency in standard tableaux generation");
 PERM_RANK_TYPE und=p.Un_irrep_dim(unitary_n);
-cout<<"U("<<unitary_n<<") ir dim "<<und<<endl;
+cout<<"U("<<unitary_n<<") IR dim "<<und<<endl;
 space_dim += und*snd;
 CompressedPartition<int> qc(q);
 cout <<"("<<qc<<')';
@@ -2211,7 +2221,7 @@ int tot=p.generate_all_multi(printme0);
 cout <<"generated "<<tot<<endl;
 }
 
-if(1)
+if(0)
 {
 int n; cin >>n;
 NRPerm<int> p(n); 
@@ -2230,14 +2240,36 @@ if(0)
 {
 int n;;
 cin >>n >>unitary_n;
+Sn_characters<int> Sn(n);
+cout <<Sn;
+if(!Sn.is_valid()) laerror("internal error in Sn character calculation");
+cout <<"allyoung.resize "<<Sn.sumirrepdims()<<endl;
+allyoung.resize(Sn.sumirrepdims());
+allyoung_irrep.resize(Sn.sumirrepdims());
+allyoung_index=0;
+
+
 Partition<int> p(n);
 space_dim=0;
 int tot=p.generate_all(pprintme,0);
 cout <<"partitions generated "<<tot<<endl;
 if(tot!=partitions(n)) laerror("internal error in partition generation or enumerations");
 if(space_dim!=longpow(unitary_n,n)) {cout<<space_dim<<" "<<ipow(unitary_n,n)<<endl;laerror("integer overflow or internal error in space dimensions");}
+
+for(int i=0; i<allyoung.size(); ++i)
+	for(int j=0; j<allyoung.size(); ++j)
+		{
+		PermutationAlgebra<int,int> r=allyoung[i]*allyoung[j];
+		//cout <<"Young "<<i<<" "<<allyoung[i]<<endl;
+		//cout <<"Young "<<j<<" "<<allyoung[j]<<endl;
+		cout <<"Product of Young "<<i<<" and "<<j<<" = "<<r<<"\n";
+		if(i!=j && !r.is_zero()) cout <<"NONORTHOGONAL Young operators found "<<i<< " "<<j<<" (irreps "<<allyoung_irrep[i]<<" "<<allyoung_irrep[j]<<")\n";
+		if(allyoung_irrep[i]!=allyoung_irrep[j] && !r.is_zero()) laerror("internal error in PermutationAlgebra");
+		}
 }
 
+
+
 if(0)
 {
 int n;
@@ -3051,4 +3083,8 @@ if(!dif.is_zero()) laerror("error in gf 2^n sqrt");
 }
 
 
+if(1)
+{
+}
+
 }

vec.h

@@ -1986,6 +1986,51 @@ for(typename std::list<T>::const_iterator i=l.begin(); i!=l.end(); ++i) (*this)[
 }
 
 
+//general simplification template for a NRVec of a class consisting from a coefficient and an element
+//the class T must have traits for sorting, normtype, elementtype, coefficienttype, coefficient, abscoefficient, and operator== which ignores the coefficient and uses just the element
+//it is not a member function to avoid the need of extra traits when this is not needed
+
+template<typename T>
+void NRVec_simplify(NRVec<T> &x, const typename LA_traits<T>::normtype thr=0, bool alwayskeepfirst=false)
+{
+if(x.size()==0) return;
+
+x.copyonwrite();
+
+//first sort to identify identical terms
+x.sort();
+
+//the following operations conserve the sorting, so no need to reset the issorted flag
+
+int newsize=1;
+//add factors of identical elements
+for(int next=1; next<x.size(); ++next)
+	{
+	if(x[next] == x[newsize-1])
+		{
+		LA_traits<T>::coefficient(x[newsize-1]) += LA_traits<T>::coefficient(x[next]);
+		}
+	else
+		{
+		if(next!=newsize) x[newsize] = x[next];
+		++newsize;
+		}
+	}
+
+//now skip terms with zero coefficients
+int newsize2=0;
+for(int next=0; next<newsize; ++next)
+	{
+	if(next==0 && alwayskeepfirst || !(LA_traits<T>::coefficient(x[next])==0 || LA_traits<T>::abscoefficient(x[next]) <thr))
+		{
+		if(next!=newsize2) x[newsize2] = x[next];
+		++newsize2;
+		}
+	}
+
+x.resize(newsize2,true);
+}
+
 
 
 }//namespace