continue implementing permutations

mat.h

@@ -26,6 +26,7 @@ namespace LA {
 
 //forward declaration
 template<typename T> class NRPerm;
+template<typename T> class CyclePerm;
 
 
 /***************************************************************************//**

permutation.cc

@@ -117,11 +117,173 @@ for(T i=2;i<=n;i++) for(T j=1;j<i;j++) if((*this)[j]>(*this)[i]) count++;
 return (count&1)? -1:1;
 }
 
+template <typename T>
+NRPerm<T>::NRPerm(const CyclePerm<T> &rhs, int n)
+: NRVec_from1<T>(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 <typename T>
+CyclePerm<T>:: CyclePerm(const NRPerm<T> &p)
+{
+#ifdef DEBUG
+        if(!p.is_valid()) laerror("invalid permutation");
+#endif
+T n=p.size();
+NRVec_from1<T> used(0,n),tmp(n);
+T firstunused=1;
+T currentcycle=0;
+std::list<NRVec<T> > 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<T> 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 <typename T>
+bool CyclePerm<T>::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 <typename T>
+bool CyclePerm<T>::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 <typename T>
+CyclePerm<T> CyclePerm<T>::inverse() const
+{
+#ifdef DEBUG
+if(!this->is_valid()) laerror("operation with an invalid cycleperm");
+#endif
+
+CyclePerm<T> 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 <typename T>
+int CyclePerm<T>::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 <typename T>
+Partition<T> CyclePerm<T>::cycles(const T n) const
+{
+#ifdef DEBUG
+if(!this->is_valid()) laerror("operation with an invalid cycleperm");
+#endif
+Partition<T> 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<int>;
+template class CyclePerm<int>;
+template class Partition<int>;
 }

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 <typename T> class NRVec_from1;
+template <typename T> class CyclePerm;
+template <typename T> class Partition;
 
 template <typename T>
 class NRPerm : public NRVec_from1<T> {
@@ -40,6 +41,7 @@ public:
 	NRPerm(const NRVec_from1<T> &rhs): NRVec_from1<T>(rhs) {};
 	NRPerm(const T &a, const int n): NRVec_from1<T>(a, n) {};
         NRPerm(const T *a, const int n): NRVec_from1<T>(a, n) {};
+	NRPerm(const CyclePerm<T> &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 <typename T>
+class CyclePerm : public NRVec_from1<NRVec_from1<T> > {
+public:
+	CyclePerm() :  NRVec_from1<NRVec_from1<T> >() {};
+	CyclePerm(const NRPerm<T> &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<T> 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 <typename T>
+class Partition : public NRVec_from1<T> {
+public:
+	Partition(): NRVec_from1<T>() {};
+	Partition(const int n) : NRVec_from1<T>(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

vec.h

@@ -31,6 +31,7 @@ template <typename T> void lawritemat(FILE *file, const T *a, int r, int c,
 		const char *form0, int nodim, int modulo, int issym);
 
 template <typename T> class NRPerm;
+template <typename T> 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<nn; i++){ sum += v[i];}
+		T sum(v[0]);
+		for(register int i=1; i<nn; i++){ sum += v[i];}
 		return sum;
 	};
 
+	//! compute the product of the vector elements
+        inline const T prod() const {
+                T prod(v[0]);
+                for(register int i=1; i<nn; i++){ prod *= v[i];}
+                return prod;
+        };
+
 	//! permute vector elements
 	const NRVec permute(const NRPerm<int> &p) const;