From 83b9463334ac862e452cda919a6d81cb326ab675 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Fri, 14 May 2021 17:39:22 +0200 Subject: [PATCH] continue implementing permutations --- mat.h | 1 + permutation.cc | 162 +++++++++++++++++++++++++++++++++++++++++++++++++ permutation.h | 40 +++++++++++- vec.h | 12 +++- 4 files changed, 211 insertions(+), 4 deletions(-) diff --git a/mat.h b/mat.h index 34198d2..0b1eb41 100644 --- a/mat.h +++ b/mat.h @@ -26,6 +26,7 @@ namespace LA { //forward declaration template class NRPerm; +template class CyclePerm; /***************************************************************************//** diff --git a/permutation.cc b/permutation.cc index cee6b63..642673b 100644 --- a/permutation.cc +++ b/permutation.cc @@ -117,11 +117,173 @@ for(T i=2;i<=n;i++) for(T j=1;j(*this)[i]) count++; return (count&1)? -1:1; } +template +NRPerm::NRPerm(const CyclePerm &rhs, int n) +: NRVec_from1(n) +{ +#ifdef DEBUG + if(!rhs.is_valid()) laerror("invalid cycle permutation"); +#endif +identity(); +T ncycles=rhs.size(); +for(T j=1; j<=ncycles; ++j) + { + T length= rhs[j].size(); + for(T i=1; i<=length; ++i) (*this)[rhs[j][i]] = rhs[j][i%length+1]; + } +#ifdef DEBUG +if(!is_valid()) laerror("internal error in NRPerm constructor from CyclePerm"); +#endif +} + + +//////////////////////////////////////////////////////// + +template +CyclePerm:: CyclePerm(const NRPerm &p) +{ +#ifdef DEBUG + if(!p.is_valid()) laerror("invalid permutation"); +#endif +T n=p.size(); +NRVec_from1 used(0,n),tmp(n); +T firstunused=1; +T currentcycle=0; +std::list > cyclelist={}; +do + { + //find a cycle starting with first unused element + T cyclelength=0; + T next = firstunused; + do + { + ++cyclelength; + used[next]=1; + tmp[cyclelength] = next; + next = p[next]; + } + while(used[next]==0); + if(cyclelength>1) //nontrivial cycle + { + NRVec_from1 cycle(&tmp[1],cyclelength); + cyclelist.push_front(cycle); + ++currentcycle; + } + while(used[firstunused]) {++firstunused; if(firstunused>n) break;} //find next unused element + } +while(firstunused<=n); +//convert list to NRVec +this->resize(currentcycle); +T i=1; +for(auto l=cyclelist.begin(); l!=cyclelist.end(); ++l) (*this)[i++] = *l; +} + + +template +bool CyclePerm::is_valid() const +{ +for(T i=1; i<=this->size(); ++i) + { + T n=(*this)[i].size(); + if(n<=0) return false; + for(T j=1; j<=n; ++j) + { + T x=(*this)[i][j]; + if(x<=0) return false; + + //now check for illegal duplicity of numbers withis a cycle or across cycles + for(T ii=i; ii<=this->size(); ++ii) + { + T nn=(*this)[ii].size(); + for(T jj=1; jj<=nn; ++jj) + { + T xx=(*this)[ii][jj]; + if(x==xx) return false; + } + } + } + } +return true; +} + + +template +bool CyclePerm::is_identity() const +{ +#ifdef DEBUG +if(!this->is_valid()) laerror("operation with an invalid cycleperm"); +#endif +for(T i=1; i<=this->size(); ++i) if((*this)[i].size()>1) return false; //at least one nontrivial cycle +return true; +} + + +template +CyclePerm CyclePerm::inverse() const +{ +#ifdef DEBUG +if(!this->is_valid()) laerror("operation with an invalid cycleperm"); +#endif + +CyclePerm r; +T ncycles=this->size(); +r.resize(ncycles); +for(T i=1; i<=ncycles; ++i) + { + T length=(*this)[i].size(); + r[i].resize(length); + for(T j=1; j<=length; ++j) r[i][j] = (*this)[i][length-j+1]; + } +return r; +} + + +template +int CyclePerm::parity() const +{ +#ifdef DEBUG +if(!this->is_valid()) laerror("operation with an invalid cycleperm"); +#endif +int n_even_cycles=0; +T ncycles=this->size(); +for(T i=1; i<=ncycles; ++i) + { + T length=(*this)[i].size(); + if((length&1)==0) ++n_even_cycles; + } +return (n_even_cycles&1)?-1:1; +} + + +template +Partition CyclePerm::cycles(const T n) const +{ +#ifdef DEBUG +if(!this->is_valid()) laerror("operation with an invalid cycleperm"); +#endif +Partition r(n); r.clear(); +T ncycles=this->size(); +for(T i=1; i<=ncycles; ++i) + { + T length=(*this)[i].size(); + if(length<=0||length>n) laerror("unexpected cycle length in permutation"); + r[length]++; + } +//fill in trivial cycles of length one +r[1] = n - r.sum(); +if(r[1]<0) laerror("inconsistent cycle lengths in CyclePerm::cycles"); +return r; +} + + +/////////////////////////////////////////////////////// /***************************************************************************//** * forced instantization in the corresponding object file ******************************************************************************/ template class NRPerm; +template class CyclePerm; +template class Partition; } diff --git a/permutation.h b/permutation.h index 1839fe6..8c651f4 100644 --- a/permutation.h +++ b/permutation.h @@ -22,7 +22,6 @@ #include "la_traits.h" #include "vec.h" -#include "mat.h" //permutations are always numbered from 1; offset is employed when applied to vectors and matrices @@ -30,6 +29,8 @@ namespace LA { //forward declaration template class NRVec_from1; +template class CyclePerm; +template class Partition; template class NRPerm : public NRVec_from1 { @@ -40,6 +41,7 @@ public: NRPerm(const NRVec_from1 &rhs): NRVec_from1(rhs) {}; NRPerm(const T &a, const int n): NRVec_from1(a, n) {}; NRPerm(const T *a, const int n): NRVec_from1(a, n) {}; + NRPerm(const CyclePerm &rhs, int n); //specific operations void identity(); @@ -57,10 +59,44 @@ public: //@@@next permutation //@@@lex rank //@@@inversion tables - //@@@conversion to cycle structure and back }; +//permutations represented in the cycle format +template +class CyclePerm : public NRVec_from1 > { +public: + CyclePerm() : NRVec_from1 >() {}; + CyclePerm(const NRPerm &rhs); + + bool is_valid() const; //is it really a permutation + bool is_identity() const; //no cycles of length > 1 + CyclePerm inverse() const; //reverse all cycles + int parity() const; //negative if having odd number of even-length cycles + Partition cycles(const T n) const; + //@@@efficient algorithm for multiplication? + //@@@operator >> and << + //@@@operation in place on matrix and vector +}; + + +//partitions stored as #of 1s, #of 2s, etc. +template +class Partition : public NRVec_from1 { +public: + Partition(): NRVec_from1() {}; + Partition(const int n) : NRVec_from1(n) {}; + T sum() const {T s=0; for(T i=1; i<=this->size(); ++i) s += i*(*this)[i]; return s;} + T nparts() const {T s=0; for(T i=1; i<=this->size(); ++i) s += (*this)[i]; return s;} + T nclasses() const {T s=0; for(T i=1; i<=this->size(); ++i) if((*this)[i]) ++s; return s;} + bool is_valid() const {return this->size() == this->sum();} + +//@@@generate all partitions, +//@@@enumerator of partitions of n to r parts and total +//@@@adjoint partition, +//@@@ output as in the group character table +//@@@Sn character table +}; }//namespace diff --git a/vec.h b/vec.h index 09da04a..5c40609 100644 --- a/vec.h +++ b/vec.h @@ -31,6 +31,7 @@ template void lawritemat(FILE *file, const T *a, int r, int c, const char *form0, int nodim, int modulo, int issym); template class NRPerm; +template class CyclePerm; /***************************************************************************//** * static constants used in several cblas-routines @@ -257,11 +258,18 @@ public: //! compute the sum of the vector elements inline const T sum() const { - T sum(0); - for(register int i=0; i &p) const;