From 13c23fb85d15ebaefe0649961bf314e923a8d626 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Thu, 18 Jan 2024 23:45:42 +0100
Subject: [PATCH] regular representation of permutations implemented

---
 mat.h          |  2 +-
 permutation.cc | 69 +++++++++++++++++++++++++++++++++++++++++++++++++-
 permutation.h  | 12 ++++++++-
 t.cc           | 28 ++++++++++++++++----
 4 files changed, 103 insertions(+), 8 deletions(-)

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<int> &p, const bool direction);
 	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 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
 
 
 	/***************************************************************************//**
diff --git a/permutation.cc b/permutation.cc
index 79b3f0e..8df6e95 100644
--- a/permutation.cc
+++ b/permutation.cc
@@ -22,6 +22,7 @@
 #include <string.h>
 #include <list>
 #include "qsort.h"
+#include "bitvector.h"
 
 
 namespace LA {
@@ -320,6 +321,8 @@ return ret;
 }
 
 
+
+
 template <typename T, typename R>
 bool PermutationAlgebra<T,R>::operator==(PermutationAlgebra<T,R> &rhs)
 {
@@ -426,6 +429,23 @@ if(callback) (*callback)(*this);
 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>
 static T _n2;
 
@@ -765,6 +785,7 @@ for(int i=0; i<this->size(); ++i)
 return res;
 }
 
+
 template <typename T, typename R>
 PermutationAlgebra<T,R> PermutationAlgebra<T,R>::operator*(const PermutationAlgebra<T,R> &rhs) const
 {
@@ -777,6 +798,7 @@ res.simplify();
 return res;
 }
 
+
 template <typename T, typename R>
 PermutationAlgebra<T,R> PermutationAlgebra<T,R>::operator+(const PermutationAlgebra<T,R> &rhs) const
 {
@@ -974,11 +996,12 @@ return r;
 
 
 template <typename T>
-CompressedPartition<T> CyclePerm<T>::cycles(const T n) const
+CompressedPartition<T> CyclePerm<T>::cycles(T n) const
 {
 #ifdef DEBUG
 if(!this->is_valid()) laerror("operation with an invalid cycleperm");
 #endif
+if(n==0) n=max();
 CompressedPartition<T> r(n); r.clear();
 T ncycles=this->size();
 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>
 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 Sn_characters<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 std::istream & operator>>(std::istream &s, 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 std::istream & operator>>(std::istream &s, 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)
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<T> cycles() const {return CyclePerm<T>(*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<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_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_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()
@@ -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<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 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
 	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;}
-	CompressedPartition<T> cycles(const T n) const;
+	CompressedPartition<T> 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<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;
 }
 
+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>
 PermutationAlgebra<T,T> general_antisymmetrizer(const NRVec<NRVec_from1<T> > &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<<p;
 }
 
 
+static NRMat<PERM_RANK_TYPE> Snmtable;
 static int unitary_n;
 static PERM_RANK_TYPE space_dim;
 static NRVec<PermutationAlgebra<int,int> > allyoung;
 static NRVec<NRMat<int> >allyoungmat;
+static NRVec<NRMat<int> >allyoungregular;
 static NRVec<int> 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 "<<allyoung_index<<" (irrep "<<current_irrep<<") = "<<allyoung[allyoung_index]<<endl;
 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_index++;
 }
@@ -2246,9 +2250,14 @@ cin >>n >>unitary_n;
 Sn_characters<int> Sn(n);
 cout <<Sn;
 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;
 allyoung.resize(Sn.sumirrepdims());
 allyoungmat.resize(Sn.sumirrepdims());
+allyoungregular.resize(Sn.sumirrepdims());
 allyoung_irrep.resize(Sn.sumirrepdims());
 allyoung_index=0;
 
@@ -2264,14 +2273,23 @@ for(int i=0; i<allyoung.size(); ++i)
 	{
 	for(int j=0; j<allyoung.size(); ++j)
 		{
-		//cout <<"Young "<<i<<" "<<allyoung[i]<<endl;
-		//cout <<"Young "<<j<<" "<<allyoung[j]<<endl;
+		//	cout <<"Young "<<i<<" "<<allyoung[i]<<endl;
+		//	cout <<"Young "<<j<<" "<<allyoung[j]<<endl;
 		PermutationAlgebra<int,int> r=allyoung[i]*allyoung[j];
 		NRMat<int> rm(r,false,n);
 		NRMat<int> rm2 = allyoungmat[i]*allyoungmat[j];
-		cout <<"Product of Young "<<i<<" and "<<j<<" = "<<r<<"\n"<<"matrix "<<rm<<endl;
-		if(rm!=rm2) laerror("internal error in matrix representation of permutationalgebra");
+		NRMat<int> rreg=RegularRepresentation(r,Snmtable);
+		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(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");
 		}
 	}