regular representation of permutations implemented

This commit is contained in:
Jiri Pittner 2024-01-18 23:45:42 +01:00
parent 2cb5258cd0
commit 13c23fb85d
4 changed files with 103 additions and 8 deletions

2
mat.h
View File

@ -139,7 +139,7 @@ public:
void axpy(const T alpha, const NRPerm<int> &p, const bool direction); void axpy(const T alpha, const NRPerm<int> &p, const bool direction);
explicit NRMat(const NRPerm<int> &p, const bool direction, const bool parity=false); //permutation matrix explicit NRMat(const NRPerm<int> &p, const bool direction, const bool parity=false); //permutation matrix
explicit NRMat(const WeightPermutation<int,T> &p, const bool direction); explicit NRMat(const WeightPermutation<int,T> &p, const bool direction);
explicit NRMat(const PermutationAlgebra<int,T> &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<int,T> &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
/***************************************************************************//** /***************************************************************************//**

View File

@ -22,6 +22,7 @@
#include <string.h> #include <string.h>
#include <list> #include <list>
#include "qsort.h" #include "qsort.h"
#include "bitvector.h"
namespace LA { namespace LA {
@ -320,6 +321,8 @@ return ret;
} }
template <typename T, typename R> template <typename T, typename R>
bool PermutationAlgebra<T,R>::operator==(PermutationAlgebra<T,R> &rhs) bool PermutationAlgebra<T,R>::operator==(PermutationAlgebra<T,R> &rhs)
{ {
@ -426,6 +429,23 @@ if(callback) (*callback)(*this);
return np; return np;
} }
template <typename T>
PermutationAlgebra<T,T> NRPerm<T>::list_all_lex()
{
PERM_RANK_TYPE number = factorial(this->size());
PermutationAlgebra<T,T> 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 <typename T> template <typename T>
static T _n2; static T _n2;
@ -765,6 +785,7 @@ for(int i=0; i<this->size(); ++i)
return res; return res;
} }
template <typename T, typename R> template <typename T, typename R>
PermutationAlgebra<T,R> PermutationAlgebra<T,R>::operator*(const PermutationAlgebra<T,R> &rhs) const PermutationAlgebra<T,R> PermutationAlgebra<T,R>::operator*(const PermutationAlgebra<T,R> &rhs) const
{ {
@ -777,6 +798,7 @@ res.simplify();
return res; return res;
} }
template <typename T, typename R> template <typename T, typename R>
PermutationAlgebra<T,R> PermutationAlgebra<T,R>::operator+(const PermutationAlgebra<T,R> &rhs) const PermutationAlgebra<T,R> PermutationAlgebra<T,R>::operator+(const PermutationAlgebra<T,R> &rhs) const
{ {
@ -974,11 +996,12 @@ return r;
template <typename T> template <typename T>
CompressedPartition<T> CyclePerm<T>::cycles(const T n) const CompressedPartition<T> CyclePerm<T>::cycles(T n) const
{ {
#ifdef DEBUG #ifdef DEBUG
if(!this->is_valid()) laerror("operation with an invalid cycleperm"); if(!this->is_valid()) laerror("operation with an invalid cycleperm");
#endif #endif
if(n==0) n=max();
CompressedPartition<T> r(n); r.clear(); CompressedPartition<T> r(n); r.clear();
T ncycles=this->size(); T ncycles=this->size();
for(T i=1; i<=ncycles; ++i) for(T i=1; i<=ncycles; ++i)
@ -2040,6 +2063,48 @@ return r;
} }
template<typename T>
NRMat<PERM_RANK_TYPE> Multable(T n)
{
NRPerm<T> p(n);
PermutationAlgebra<T,T> all = p.list_all_lex();
NRMat<PERM_RANK_TYPE> r(all.size(),all.size());
for(PERM_RANK_TYPE i=0; i<all.size(); ++i) r[0][i] = r[i][0]=i; //identity
for(PERM_RANK_TYPE i=1; i<all.size(); ++i)
for(PERM_RANK_TYPE j=1; j<all.size(); ++j)
{
NRPerm<T> 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<all.size(); ++i)
{
occ.clear();
for(PERM_RANK_TYPE j=0; j<all.size(); ++j) {if(occ[r(i,j)]) laerror("inconsistency in Multable"); occ.set(r(i,j));}
occ.clear();
for(PERM_RANK_TYPE j=0; j<all.size(); ++j) {if(occ[r(j,i)]) laerror("inconsistency in Multable"); occ.set(r(j,i));}
}
#endif
return r;
}
template<typename T, typename R>
NRMat<R> RegularRepresentation(const PermutationAlgebra<T,R> &a, const NRMat<PERM_RANK_TYPE> &mtable)
{
NRMat<R> r(mtable.nrows(),mtable.ncols());
r.clear();
for(int i=0; i<a.size(); ++i)
{
PERM_RANK_TYPE rx=a[i].perm.rank();
for(PERM_RANK_TYPE j=0; j<mtable.nrows();++j) r(mtable(rx,j),j) += a[i].weight;
}
return r;
}
template<typename T> template<typename T>
PermutationAlgebra<T,T> general_antisymmetrizer(const NRVec<NRVec_from1<T> > &groups, int restriction_type, bool inverted) PermutationAlgebra<T,T> general_antisymmetrizer(const NRVec<NRVec_from1<T> > &groups, int restriction_type, bool inverted)
{ {
@ -2072,6 +2137,7 @@ template class Partition<T>; \
template class YoungTableaux<T>; \ template class YoungTableaux<T>; \
template class Sn_characters<T>; \ template class Sn_characters<T>; \
template class CycleIndex<T>; \ template class CycleIndex<T>; \
template NRMat<PERM_RANK_TYPE> Multable(T n); \
template PermutationAlgebra<T,T> general_antisymmetrizer(const NRVec<NRVec_from1<T> > &groups, int, bool); \ template PermutationAlgebra<T,T> general_antisymmetrizer(const NRVec<NRVec_from1<T> > &groups, int, bool); \
template std::istream & operator>>(std::istream &s, CyclePerm<T> &x); \ template std::istream & operator>>(std::istream &s, CyclePerm<T> &x); \
template std::ostream & operator<<(std::ostream &s, const CyclePerm<T> &x); \ template std::ostream & operator<<(std::ostream &s, const CyclePerm<T> &x); \
@ -2086,6 +2152,7 @@ template class WeightPermutation<T,R>; \
template class PermutationAlgebra<T,R>; \ template class PermutationAlgebra<T,R>; \
template std::istream & operator>>(std::istream &s, WeightPermutation<T,R> &x); \ template std::istream & operator>>(std::istream &s, WeightPermutation<T,R> &x); \
template std::ostream & operator<<(std::ostream &s, const WeightPermutation<T,R> &x); \ template std::ostream & operator<<(std::ostream &s, const WeightPermutation<T,R> &x); \
template NRMat<R> RegularRepresentation(const PermutationAlgebra<T,R> &a, const NRMat<PERM_RANK_TYPE> &mtable); \
INSTANTIZE(int) INSTANTIZE(int)

View File

@ -60,6 +60,7 @@ public:
void identity(); void identity();
bool is_valid() const; //is it really a permutation bool is_valid() const; //is it really a permutation
bool is_identity() const; bool is_identity() const;
CompressedPartition<T> cycles() const {return CyclePerm<T>(*this).cycles(size());};
NRPerm inverse() const; NRPerm inverse() const;
NRPerm reverse() const; //backward order NRPerm reverse() const; //backward order
NRPerm operator&(const NRPerm &rhs) const; //concatenate the permutations this,rhs, renumbering rhs (not commutative) 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 bool next(); //generate next permutation in lex order
PERM_RANK_TYPE generate_all(void (*callback)(const NRPerm<T>&), int parity_select=0); //Algorithm L from Knuth's vol.4, efficient but not in lex order! PERM_RANK_TYPE generate_all(void (*callback)(const NRPerm<T>&), int parity_select=0); //Algorithm L from Knuth's vol.4, efficient but not in lex order!
PermutationAlgebra<T,T> list_all(int parity_select=0); PermutationAlgebra<T,T> list_all(int parity_select=0);
PermutationAlgebra<T,T> list_all_lex();
PERM_RANK_TYPE generate_all_multi(void (*callback)(const NRPerm<T>&)); //Algorithm L2 from Knuth's vol.4, for a multiset (repeated numbers, not really permutations) PERM_RANK_TYPE generate_all_multi(void (*callback)(const NRPerm<T>&)); //Algorithm L2 from Knuth's vol.4, for a multiset (repeated numbers, not really permutations)
PERM_RANK_TYPE generate_all2(void (*callback)(const NRPerm<T>&)); //recursive method, also not lexicographic PERM_RANK_TYPE generate_all2(void (*callback)(const NRPerm<T>&)); //recursive method, also not lexicographic
PERM_RANK_TYPE generate_all_lex(void (*callback)(const NRPerm<T>&)); //generate in lex order using next() PERM_RANK_TYPE generate_all_lex(void (*callback)(const NRPerm<T>&)); //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->perm < rhs.perm;};
bool operator>=(const WeightPermutation &rhs) const {return !(*this < rhs);}; bool operator>=(const WeightPermutation &rhs) const {return !(*this < rhs);};
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<T,R>& x){return x.weight;}; static R coefficient(const WeightPermutation<T,R>& x){return x.weight;};
static R& coefficient(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);}; static typename LA_traits<R>::normtype abscoefficient(const WeightPermutation<T,R>& x){return LA_traits<R>::abs2(x.weight);};
static void clear(WeightPermutation<T,R> *v, int nn) {for(int i=0; i<nn; ++i) {v[i].weight=0; v[i].perm.clear();}}
}; };
@ -211,7 +215,7 @@ public:
CyclePerm inverse() const; //reverse all cycles CyclePerm inverse() const; //reverse all cycles
int parity() const; //negative if having odd number of even-length cycles int parity() const; //negative if having odd number of even-length cycles
T max() const {T m=0; for(int i=1; i<=this->size(); ++i) {T mm= (*this)[i].max(); if(mm>m) m=mm;} return m;} T max() const {T m=0; for(int i=1; i<=this->size(); ++i) {T mm= (*this)[i].max(); if(mm>m) m=mm;} return m;}
CompressedPartition<T> cycles(const T n) const; CompressedPartition<T> cycles(T n = 0) const;
void readfrom(const std::string &line); void readfrom(const std::string &line);
CyclePerm operator*(const CyclePerm &q) const; //q is rhs and applied first, this applied second CyclePerm operator*(const CyclePerm &q) const; //q is rhs and applied first, this applied second
NRPerm<T> operator*(const NRPerm<T> &r) const; NRPerm<T> operator*(const NRPerm<T> &r) const;
@ -390,6 +394,12 @@ else for(int i=1; i<=n; ++i) r[p[i]-1] = v[i-1];
return r; return r;
} }
template<typename T>
NRMat<PERM_RANK_TYPE> Multable(T n);
template<typename T, typename R>
NRMat<R> RegularRepresentation(const PermutationAlgebra<T,R> &a, const NRMat<PERM_RANK_TYPE> &mtable);
template<typename T> template<typename T>
PermutationAlgebra<T,T> general_antisymmetrizer(const NRVec<NRVec_from1<T> > &groups, int restriction_type=0, bool inverted=false); PermutationAlgebra<T,T> general_antisymmetrizer(const NRVec<NRVec_from1<T> > &groups, int restriction_type=0, bool inverted=false);

28
t.cc
View File

@ -88,10 +88,12 @@ cout<<p;
} }
static NRMat<PERM_RANK_TYPE> Snmtable;
static int unitary_n; static int unitary_n;
static PERM_RANK_TYPE space_dim; static PERM_RANK_TYPE space_dim;
static NRVec<PermutationAlgebra<int,int> > allyoung; static NRVec<PermutationAlgebra<int,int> > allyoung;
static NRVec<NRMat<int> >allyoungmat; static NRVec<NRMat<int> >allyoungmat;
static NRVec<NRMat<int> >allyoungregular;
static NRVec<int> allyoung_irrep; static NRVec<int> allyoung_irrep;
int current_irrep; int current_irrep;
int allyoung_index; int allyoung_index;
@ -103,7 +105,9 @@ if(!y.is_standard()) laerror("internal error in young");
allyoung[allyoung_index] = y.young_operator(); allyoung[allyoung_index] = y.young_operator();
cout <<"Young "<<allyoung_index<<" (irrep "<<current_irrep<<") = "<<allyoung[allyoung_index]<<endl; cout <<"Young "<<allyoung_index<<" (irrep "<<current_irrep<<") = "<<allyoung[allyoung_index]<<endl;
allyoungmat[allyoung_index] = NRMat<int>(allyoung[allyoung_index],false); allyoungmat[allyoung_index] = NRMat<int>(allyoung[allyoung_index],false);
cout <<"Matrix representation = "<<allyoungmat[allyoung_index]; allyoungregular[allyoung_index] = RegularRepresentation(allyoung[allyoung_index],Snmtable);
//cout <<"Matrix representation = "<<allyoungmat[allyoung_index];
cout <<"Regular representation = "<<allyoungregular[allyoung_index];
allyoung_irrep[allyoung_index]=current_irrep; allyoung_irrep[allyoung_index]=current_irrep;
allyoung_index++; allyoung_index++;
} }
@ -2246,9 +2250,14 @@ cin >>n >>unitary_n;
Sn_characters<int> Sn(n); Sn_characters<int> Sn(n);
cout <<Sn; cout <<Sn;
if(!Sn.is_valid()) laerror("internal error in Sn character calculation"); if(!Sn.is_valid()) laerror("internal error in Sn character calculation");
Snmtable = Multable(n);
cout <<"Multiplication table = "<<Snmtable<<endl;
cout <<"allyoung.resize "<<Sn.sumirrepdims()<<endl; cout <<"allyoung.resize "<<Sn.sumirrepdims()<<endl;
allyoung.resize(Sn.sumirrepdims()); allyoung.resize(Sn.sumirrepdims());
allyoungmat.resize(Sn.sumirrepdims()); allyoungmat.resize(Sn.sumirrepdims());
allyoungregular.resize(Sn.sumirrepdims());
allyoung_irrep.resize(Sn.sumirrepdims()); allyoung_irrep.resize(Sn.sumirrepdims());
allyoung_index=0; allyoung_index=0;
@ -2264,14 +2273,23 @@ for(int i=0; i<allyoung.size(); ++i)
{ {
for(int j=0; j<allyoung.size(); ++j) for(int j=0; j<allyoung.size(); ++j)
{ {
//cout <<"Young "<<i<<" "<<allyoung[i]<<endl; // cout <<"Young "<<i<<" "<<allyoung[i]<<endl;
//cout <<"Young "<<j<<" "<<allyoung[j]<<endl; // cout <<"Young "<<j<<" "<<allyoung[j]<<endl;
PermutationAlgebra<int,int> r=allyoung[i]*allyoung[j]; PermutationAlgebra<int,int> r=allyoung[i]*allyoung[j];
NRMat<int> rm(r,false,n); NRMat<int> rm(r,false,n);
NRMat<int> rm2 = allyoungmat[i]*allyoungmat[j]; NRMat<int> rm2 = allyoungmat[i]*allyoungmat[j];
cout <<"Product of Young "<<i<<" and "<<j<<" = "<<r<<"\n"<<"matrix "<<rm<<endl; NRMat<int> rreg=RegularRepresentation(r,Snmtable);
if(rm!=rm2) laerror("internal error in matrix representation of permutationalgebra"); NRMat<int> rreg2=allyoungregular[i]*allyoungregular[j];
cout <<"Product of Young "<<i<<" and "<<j<<" = "<<r<<endl;
//cout<<"matrix "<<rm<<endl;
if(i!=j && !r.is_zero()) cout <<"NONORTHOGONAL Young operators found "<<i<< " "<<j<<" (irreps "<<allyoung_irrep[i]<<" "<<allyoung_irrep[j]<<")\n"; if(i!=j && !r.is_zero()) cout <<"NONORTHOGONAL Young operators found "<<i<< " "<<j<<" (irreps "<<allyoung_irrep[i]<<" "<<allyoung_irrep[j]<<")\n";
if(rreg!=rreg2)
{
cout <<"Representation of product = "<<rreg;
cout <<"Product of representations = "<<rreg2;
laerror("internal error in multiplication of permutationalgebra");
}
if(rm!=rm2) laerror("internal error in matrix representation of permutationalgebra");
if(allyoung_irrep[i]!=allyoung_irrep[j] && !r.is_zero()) laerror("internal error in PermutationAlgebra"); if(allyoung_irrep[i]!=allyoung_irrep[j] && !r.is_zero()) laerror("internal error in PermutationAlgebra");
} }
} }