diff --git a/mat.h b/mat.h index 3c2086a..aab82fd 100644 --- a/mat.h +++ b/mat.h @@ -139,7 +139,7 @@ public: void axpy(const T alpha, const NRPerm &p, const bool direction); explicit NRMat(const NRPerm &p, const bool direction, const bool parity=false); //permutation matrix explicit NRMat(const WeightPermutation &p, const bool direction); - explicit NRMat(const PermutationAlgebra &p, const bool direction, const int nforce=0); //note that one cannot represent e.g. young projectors in this way, since the representation of S(n) by permutation matrices is reducible just to two irreps [n] and [n-1,1] + explicit NRMat(const PermutationAlgebra &p, const bool direction, const int nforce=0); //note that one cannot represent e.g. young projectors in this way, since the representation of S(n) by permutation matrices is reducible just to two irreps [n] and [n-1,1] since the charater of that RR = number of cycles of length=1 /***************************************************************************//** diff --git a/permutation.cc b/permutation.cc index 79b3f0e..8df6e95 100644 --- a/permutation.cc +++ b/permutation.cc @@ -22,6 +22,7 @@ #include #include #include "qsort.h" +#include "bitvector.h" namespace LA { @@ -320,6 +321,8 @@ return ret; } + + template bool PermutationAlgebra::operator==(PermutationAlgebra &rhs) { @@ -426,6 +429,23 @@ if(callback) (*callback)(*this); return np; } +template +PermutationAlgebra NRPerm::list_all_lex() +{ +PERM_RANK_TYPE number = factorial(this->size()); +PermutationAlgebra ret(number); +PERM_RANK_TYPE np=0; +this->identity(); +do{ +ret[np].perm = *this; ret[np].perm.copyonwrite(); +ret[np].weight=0; +++np; +}while(this->next()); +return ret; +} + + + template static T _n2; @@ -765,6 +785,7 @@ for(int i=0; isize(); ++i) return res; } + template PermutationAlgebra PermutationAlgebra::operator*(const PermutationAlgebra &rhs) const { @@ -777,6 +798,7 @@ res.simplify(); return res; } + template PermutationAlgebra PermutationAlgebra::operator+(const PermutationAlgebra &rhs) const { @@ -974,11 +996,12 @@ return r; template -CompressedPartition CyclePerm::cycles(const T n) const +CompressedPartition CyclePerm::cycles(T n) const { #ifdef DEBUG if(!this->is_valid()) laerror("operation with an invalid cycleperm"); #endif +if(n==0) n=max(); CompressedPartition r(n); r.clear(); T ncycles=this->size(); for(T i=1; i<=ncycles; ++i) @@ -2040,6 +2063,48 @@ return r; } +template +NRMat Multable(T n) +{ +NRPerm p(n); +PermutationAlgebra all = p.list_all_lex(); +NRMat r(all.size(),all.size()); +for(PERM_RANK_TYPE i=0; i tmp = all[i].perm * all[j].perm; + r(i,j) = tmp.rank(); + } +//consistency checks +#ifdef DEBUG +bitvector occ(all.size()); +for(PERM_RANK_TYPE i=0; i +NRMat RegularRepresentation(const PermutationAlgebra &a, const NRMat &mtable) +{ +NRMat r(mtable.nrows(),mtable.ncols()); +r.clear(); +for(int i=0; i PermutationAlgebra general_antisymmetrizer(const NRVec > &groups, int restriction_type, bool inverted) { @@ -2072,6 +2137,7 @@ template class Partition; \ template class YoungTableaux; \ template class Sn_characters; \ template class CycleIndex; \ +template NRMat Multable(T n); \ template PermutationAlgebra general_antisymmetrizer(const NRVec > &groups, int, bool); \ template std::istream & operator>>(std::istream &s, CyclePerm &x); \ template std::ostream & operator<<(std::ostream &s, const CyclePerm &x); \ @@ -2086,6 +2152,7 @@ template class WeightPermutation; \ template class PermutationAlgebra; \ template std::istream & operator>>(std::istream &s, WeightPermutation &x); \ template std::ostream & operator<<(std::ostream &s, const WeightPermutation &x); \ +template NRMat RegularRepresentation(const PermutationAlgebra &a, const NRMat &mtable); \ INSTANTIZE(int) diff --git a/permutation.h b/permutation.h index 48af5d1..a90fb6b 100644 --- a/permutation.h +++ b/permutation.h @@ -60,6 +60,7 @@ public: void identity(); bool is_valid() const; //is it really a permutation bool is_identity() const; + CompressedPartition cycles() const {return CyclePerm(*this).cycles(size());}; NRPerm inverse() const; NRPerm reverse() const; //backward order NRPerm operator&(const NRPerm &rhs) const; //concatenate the permutations this,rhs, renumbering rhs (not commutative) @@ -72,6 +73,7 @@ public: bool next(); //generate next permutation in lex order PERM_RANK_TYPE generate_all(void (*callback)(const NRPerm&), int parity_select=0); //Algorithm L from Knuth's vol.4, efficient but not in lex order! PermutationAlgebra list_all(int parity_select=0); + PermutationAlgebra list_all_lex(); PERM_RANK_TYPE generate_all_multi(void (*callback)(const NRPerm&)); //Algorithm L2 from Knuth's vol.4, for a multiset (repeated numbers, not really permutations) PERM_RANK_TYPE generate_all2(void (*callback)(const NRPerm&)); //recursive method, also not lexicographic PERM_RANK_TYPE generate_all_lex(void (*callback)(const NRPerm&)); //generate in lex order using next() @@ -126,6 +128,7 @@ public: bool operator<(const WeightPermutation &rhs) const {return this->perm < rhs.perm;}; bool operator>=(const WeightPermutation &rhs) const {return !(*this < rhs);}; bool operator<=(const WeightPermutation &rhs) const {return !(*this > rhs);}; + WeightPermutation & operator=(const WeightPermutation &rhs) {weight=rhs.weight; perm=rhs.perm; return *this;}; }; @@ -145,6 +148,7 @@ public: 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);}; + static void clear(WeightPermutation *v, int nn) {for(int i=0; isize(); ++i) {T mm= (*this)[i].max(); if(mm>m) m=mm;} return m;} - CompressedPartition cycles(const T n) const; + CompressedPartition cycles(T n = 0) const; void readfrom(const std::string &line); CyclePerm operator*(const CyclePerm &q) const; //q is rhs and applied first, this applied second NRPerm operator*(const NRPerm &r) const; @@ -390,6 +394,12 @@ else for(int i=1; i<=n; ++i) r[p[i]-1] = v[i-1]; return r; } +template +NRMat Multable(T n); + +template +NRMat RegularRepresentation(const PermutationAlgebra &a, const NRMat &mtable); + template PermutationAlgebra general_antisymmetrizer(const NRVec > &groups, int restriction_type=0, bool inverted=false); diff --git a/t.cc b/t.cc index f101b09..04ba9e1 100644 --- a/t.cc +++ b/t.cc @@ -88,10 +88,12 @@ cout< Snmtable; static int unitary_n; static PERM_RANK_TYPE space_dim; static NRVec > allyoung; static NRVec >allyoungmat; +static NRVec >allyoungregular; static NRVec allyoung_irrep; int current_irrep; int allyoung_index; @@ -103,7 +105,9 @@ if(!y.is_standard()) laerror("internal error in young"); allyoung[allyoung_index] = y.young_operator(); cout <<"Young "< r=allyoung[i]*allyoung[j]; NRMat rm(r,false,n); NRMat rm2 = allyoungmat[i]*allyoungmat[j]; - cout <<"Product of Young "< rreg=RegularRepresentation(r,Snmtable); + NRMat rreg2=allyoungregular[i]*allyoungregular[j]; + cout <<"Product of Young "<