From 32f2a1abd58e7ee3504b38f9a890d28cf58a2178 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Wed, 17 Jan 2024 17:59:19 +0100 Subject: [PATCH] working on permutation --- permutation.cc | 30 ++++++++++++++++++++++-------- permutation.h | 21 ++++++++++++++++----- t.cc | 46 +++++++++++++++++++++++++++++++++++++++++----- vec.h | 45 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 124 insertions(+), 18 deletions(-) diff --git a/permutation.cc b/permutation.cc index 4cf3878..0efdd46 100644 --- a/permutation.cc +++ b/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 int NRPerm::parity() const { if(!this->is_valid()) return 0; T n=this->size(); if(n==1) return 1; +if(n>=10) //??? + { + NRPerm tmp(*this); + return (tmp.sort()&1) ? -1:1; + } T count=0; for(T i=2;i<=n;i++) for(T j=1;j(*this)[i]) count++; return (count&1)? -1:1; @@ -313,6 +320,17 @@ return ret; } +template +bool PermutationAlgebra::operator==(PermutationAlgebra &rhs) +{ +simplify(); +rhs.simplify(); +if(size()!=rhs.size()) return false; +for(int i=0; i res(this->size()*rhs.size()); for(int i=0; isize(); ++i) for(int j=0; j PermutationAlgebra PermutationAlgebra::operator+(const PermutationAlgebra &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 res=this->concat(rhs); +res.simplify(); +return res; } -template -void PermutationAlgebra::simplify() -{ -copyonwrite(); -sort(); -//@@@@implement merging weights of identical permutations and resize like we did in qubithamiltonian with other stuff -} diff --git a/permutation.h b/permutation.h index 9978d49..596c1d3 100644 --- a/permutation.h +++ b/permutation.h @@ -107,6 +107,9 @@ public: NRPerm 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 &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& x) {x.perm.copyonwrite();}; typedef typename LA_traits::normtype normtype; + typedef R coefficienttype; + typedef NRPerm elementtype; static inline bool smaller(const WeightPermutation& x, const WeightPermutation& y) {return x.perm& x, const WeightPermutation& y) {return x.perm>y.perm;}; - + static R coefficient(const WeightPermutation& x){return x.weight;}; + static R& coefficient(WeightPermutation& x) {return x.weight;}; + static typename LA_traits::normtype abscoefficient(const WeightPermutation& x){return LA_traits::abs2(x.weight);}; }; @@ -168,7 +174,7 @@ public: PermutationAlgebra(const NRVec > &x) : NRVec >(x) {}; int size() const {return NRVec >::size();}; void copyonwrite() {NRVec >::copyonwrite();}; - void sort(int direction = 0, int from = 0, int to = -1, int *permut = NULL, bool stable=false) {NRVec >::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 >::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::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 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 class Polynomial; //forward declaration diff --git a/t.cc b/t.cc index 31ef544..d706401 100644 --- a/t.cc +++ b/t.cc @@ -90,28 +90,38 @@ cout< > allyoung; +static NRVec allyoung_irrep; +int current_irrep; +int allyoung_index; void yprintme(const YoungTableaux&y) { cout < y(p); PERM_RANK_TYPE dim=y.generate_all_standard(yprintme); Partition q=p.adjoint(); -PERM_RANK_TYPE snd=p.Sn_irrep_dim(); -cout<<"IR dim "< qc(q); cout <<"("<>n; NRPerm p(n); @@ -2230,14 +2240,36 @@ if(0) { int n;; cin >>n >>unitary_n; +Sn_characters Sn(n); +cout < p(n); space_dim=0; int tot=p.generate_all(pprintme,0); cout <<"partitions generated "< r=allyoung[i]*allyoung[j]; + //cout <<"Young "<::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 +void NRVec_simplify(NRVec &x, const typename LA_traits::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::coefficient(x[newsize-1]) += LA_traits::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::coefficient(x[next])==0 || LA_traits::abscoefficient(x[next])