template for general power, power of NRPerm and CyclePerm

mat.cc

@@ -28,6 +28,7 @@
 #include <errno.h>
 #include <unistd.h>
 #include <math.h>
+#include "nonclass.h"
 
 
 namespace LA {
@@ -3163,12 +3164,13 @@ for(int i=0; i<n; ++i)
 }
 
 template<typename T>
-NRMat<T>::NRMat(const NRPerm<int> &p, const bool direction)
+NRMat<T>::NRMat(const NRPerm<int> &p, const bool direction, const bool parity)
 {
 int n=p.size();
 resize(n,n);
 clear();
-axpy((T)1,p,direction);
+T alpha= parity? p.parity():1;
+axpy(alpha,p,direction);
 }
 
 
@@ -3358,6 +3360,22 @@ copyonwrite();
 
 
 
+template<>
+NRMat<double> NRMat<double>::inverse()
+{
+NRMat<double> tmp(*this);
+return calcinverse(tmp);
+}
+
+template<>
+NRMat<std::complex<double> > NRMat<std::complex<double> >::inverse()
+{
+NRMat<std::complex<double> >  tmp(*this);
+return calcinverse(tmp);
+}       
+
+
+
 
 
 

mat.h

@@ -22,6 +22,7 @@
 #define _LA_MAT_H_
 #include "la_traits.h"
 
+
 namespace LA {
 
 //forward declaration
@@ -129,7 +130,7 @@ public:
 	void scale_row(const int i, const T f); //in place
 	void scale_col(const int i, const T f); //in place
 	void axpy(const T alpha, const NRPerm<int> &p, const bool direction);
-	explicit NRMat(const NRPerm<int> &p, const bool direction); //permutation matrix
+	explicit NRMat(const NRPerm<int> &p, const bool direction, const bool parity=false); //permutation matrix
 
 
 	/***************************************************************************//**
@@ -155,6 +156,13 @@ public:
 
 	//! assign scalar value to the diagonal elements
 	NRMat & operator=(const T &a);
+
+	//! unit matrix for general power
+	void identity() {*this = (T)1;}
+
+	//! inverse matrix
+	NRMat inverse();
+
 	//! add scalar value to the diagonal elements
 	NRMat & operator+=(const T &a);
 	//! subtract scalar value to the diagonal elements

nonclass.h

@@ -25,6 +25,48 @@
 namespace LA {
 
 //MISC
+
+
+//positive power for a general class
+template<typename T>
+T positive_power(const T &x, const int n)
+{
+if(n<=0) laerror("zero or negative n in positive_power");
+int i=n;
+if(i==1) return x;
+T y,z;
+z=x;
+while(!(i&1))
+        {
+        z = z*z;
+        i >>= 1;
+        }
+y=z;
+while((i >>= 1)/*!=0*/)
+                {
+                z = z*z;
+                if(i&1) y = y*z;
+                }
+return y;
+}
+
+
+//general integer power for a class providing identity() and inverse()
+template<typename T>
+T power(const T &x, const int n)
+{
+int i=n;
+if(i==0) {T r(x); r.identity(); return r;}
+T z;
+if(i>0) z=x;
+else {z=x.inverse(); i= -i;}
+if(i==1) return z;
+return positive_power(z,i);
+}
+
+
+
+
 template <class T>
 const NRSMat<T> twoside_transform(const NRSMat<T> &S, const NRMat<T> &C, bool transp=0) //calculate C^dagger S C
 {
@@ -138,7 +180,7 @@ extern void cholesky(NRMat<std::complex<double> > &a, bool upper=1);
 
 //inverse by means of linear solve, preserving rhs intact
 template<typename T>
-const NRMat<T> inverse(NRMat<T> a, T *det=0)
+const NRMat<T> calcinverse(NRMat<T> a, T *det=NULL)
 {
 #ifdef DEBUG
 	if(a.nrows()!=a.ncols()) laerror("inverse() for non-square matrix");

permutation.cc

@@ -569,6 +569,7 @@ for(typename std::list<NRVec_from1<T> >::iterator l=cyclelist.begin(); l!=cyclel
 }
 
 
+
 template <typename T>
 bool CyclePerm<T>::is_valid() const
 {
@@ -641,6 +642,23 @@ NRPerm<T> rr=pp*qq;
 return CyclePerm<T>(rr);
 }
 
+template <typename T>
+void CyclePerm<T>::simplify(bool keep1)
+{
+int j=0;
+for(int i=1; i<=this->size(); ++i)
+	{
+	int il= (*this)[i].size();
+	if(keep1 && il>0 || il>1 )
+		{
+		++j;
+		if(j!=i) (*this)[j] = (*this)[i]; //keep this
+		}
+	}
+this->resize(j,true);
+}
+
+
 template <typename T>
 int CyclePerm<T>::parity() const
 {
@@ -788,6 +806,57 @@ return s;
 }
 
 
+template <typename T>
+CyclePerm<T> CyclePerm<T>::pow(const int n, const bool keep1) const
+{
+#ifdef DEBUG
+        if(!this->is_valid()) laerror("power of invalid permutation");
+#endif
+CyclePerm<T> r;
+if(n==0) {r.identity(); return r;}
+
+//probably negative n will work automatically using the modulo approach
+std::list<NRVec_from1<T> > cyclelist;
+T currentcycle=0;
+
+//go through all our cycles and compute power of each cycle separately
+for(int i=1; i<=this->size(); ++i)
+	{
+	int cyclesize = (*this)[i].size();
+	if(cyclesize>0 && keep1 || cyclesize>1)
+		{
+		int modulo = n%cyclesize;
+		if(modulo<0) modulo += cyclesize;
+		if(modulo==0) 
+			{
+			if(keep1) //keep cycles of length 1 to keep info about the size of the permutation
+				{
+				for(int j=1; j<=cyclesize; ++j) {NRVec_from1<T> c(1); c[1] = (*this)[i][j]; ++currentcycle; cyclelist.push_back(c);}
+				}
+			}
+		else //the nontrivial case
+			{
+			int nsplit=gcd(modulo,cyclesize);
+			int splitsize=cyclesize/nsplit;
+			for(int j=0; j<nsplit; ++j) //loop over split cycles
+				{
+				NRVec_from1<T> c(splitsize);
+				for(int k=1; k<=splitsize; ++k) //fill in the cycle
+					{
+					c[k] = (*this)[i][((k-1)*modulo)%cyclesize + 1 + j];
+					}
+				++currentcycle; cyclelist.push_back(c);	
+				}
+			}
+		}
+	}
+
+//convert list to NRVec
+r.resize(currentcycle);
+int i=1;
+for(typename std::list<NRVec_from1<T> >::iterator l=cyclelist.begin(); l!=cyclelist.end(); ++l) r[i++] = *l;
+return r;
+}
 
 ///////////////////////////////////////////////////////
 

permutation.h

@@ -23,6 +23,7 @@
 #include "la_traits.h"
 #include "vec.h"
 #include "polynomial.h"
+#include "nonclass.h"
 
 typedef   unsigned long long PERM_RANK_TYPE;
 
@@ -30,6 +31,7 @@ typedef   unsigned long long PERM_RANK_TYPE;
 
 namespace LA {
 
+
 //forward declaration
 template <typename T> class CyclePerm;
 template <typename T> class Partition;
@@ -65,6 +67,7 @@ public:
 	NRVec_from1<T> inversions(const int type, PERM_RANK_TYPE *prank=NULL) const; //inversion tables
 	explicit NRPerm(const int type, const NRVec_from1<T> &inversions); //compute permutation from inversions
 	explicit NRPerm(const int n, const PERM_RANK_TYPE rank); //compute permutation from its rank
+	NRPerm pow(const int n) const {return power(*this,n);};
 };
 
 extern PERM_RANK_TYPE factorial(const int n);
@@ -90,6 +93,10 @@ public:
 	void readfrom(const std::string &line);
 	CyclePerm operator*(const CyclePerm q) const; //q is rhs and applied first, this applied second
 	PERM_RANK_TYPE order() const; //lcm of cycle lengths
+	bool operator==(const CyclePerm &rhs) const {return NRPerm<T>(*this) == NRPerm<T>(rhs);}; //cycle representation is not unique, cannot inherit operator== from NRVec
+	void simplify(bool keep1=false); //remove cycles of size 0 or 1
+	CyclePerm pow_simple(const int n) const {return CyclePerm(NRPerm<T>(*this).pow(n));}; //do not call power() with our operator*
+	CyclePerm pow(const int n, const bool keep1=false) const; //a more efficient algorithm
 };
 
 
@@ -144,7 +151,7 @@ public:
 	bool is_valid() const {return this->size() == this->sum();}
 	explicit CompressedPartition(const Partition<T> &rhs) : NRVec_from1<T>(rhs.size()) {this->clear(); for(int i=1; i<=rhs.size(); ++i) if(!rhs[i]) break; else (*this)[rhs[i]]++; }
 	PERM_RANK_TYPE Sn_class_size() const; 
-	int parity() const; //of a permutation with given cycle lengths
+	int parity() const; //of a permutation with given cycle lengths, returns +/- 1
 
 };
 
@@ -166,7 +173,7 @@ public:
 	PERM_RANK_TYPE Sn_irrep_dim() const;
 	PERM_RANK_TYPE Un_irrep_dim(const int n) const;
 	PERM_RANK_TYPE generate_all(void (*callback)(const Partition<T>&), int nparts=0); //nparts <0 means at most to -nparts
-	int parity() const; //of a permutation with given cycle lengths
+	int parity() const; //of a permutation with given cycle lengths, returns +/- 1
 
 };
 

polynomial.cc

@@ -20,6 +20,7 @@
 #include <stdio.h>
 #include <string.h>
 #include <math.h>
+#include "nonclass.h"
 
 
 namespace LA {
@@ -224,24 +225,9 @@ return poly_gcd(big,small,thr,d);
 template <typename T>
 Polynomial<T> Polynomial<T>::pow(const int n) const
 {
-int i=n;
+if(n<0) laerror("negative exponent in polynomial::pow");
-if(i<0) laerror("negative exponent in polynomial::pow");
+if(n==0) {Polynomial<T> r(0); r[0]=(T)1; return r;}
-if(i==0) {Polynomial<T> r(0); r[0]=(T)1; return r;}
+return positive_power(*this,n);
-if(i==1) return *this;
-Polynomial<T>y,z;
-z= *this;
-while(!(i&1))
-        {
-        z = z*z;
-        i >>= 1;
-        }
-y=z;
-while((i >>= 1)/*!=0*/)
-                {
-                z = z*z;
-                if(i&1) y = y*z;
-                }
-return y;
 }
 
 template <typename T>

quaternion.h

@@ -70,6 +70,8 @@ public:
 
 	//compiler generates default copy constructor and assignment operator
 	
+	void identity() {q[0]=(T)1; q[1]=q[2]=q[3]=0;};
+	
 	//formal indexing
 	inline const T operator[](const int i) const {return q[i];};
 	inline T& operator[](const int i) {return q[i];};

smat.cc

@@ -858,6 +858,12 @@ void NRSMat<T>::storesubmatrix(const NRVec<int> &selection, const NRSMat &rhs)
                 }
 }
 
+template<>
+NRSMat<double>  NRSMat<double>::inverse() {return NRSMat<double>(NRMat<double>(*this).inverse());}
+
+template<>
+NRSMat<std::complex<double> >  NRSMat<std::complex<double> >::inverse() {return NRSMat<std::complex<double> >(NRMat<std::complex<double> >(*this).inverse());}
+
 
 
 

smat.h

@@ -93,6 +93,13 @@ public:
 	//! assign scalar value to diagonal elements
 	NRSMat & operator=(const T &a);
 
+	//! unit matrix for general power
+        void identity() {*this = (T)1;}
+
+	//! inverse matrix
+        NRSMat inverse();
+
+
 	//! permute matrix elements
         const NRSMat permuted(const NRPerm<int> &p, const bool inverse=false) const;
 
@@ -148,7 +155,7 @@ public:
 
 	inline int nrows() const;
 	inline int ncols() const;
-	inline size_t size() const;
+	inline size_t size() const; //NOTE it is the size of data n*(n+1)/2, not nrows 
 
 	inline bool transp(const int i, const int j) const {return i>j;} //this can be used for compact storage of matrices, which are actually not symmetric, but one triangle of them is redundant
 	const typename LA_traits<T>::normtype norm(const T scalar = (T)0) const;

sparsemat.h

@@ -82,6 +82,7 @@ public:
 	explicit SparseMat(const NRSMat<T> &rhs); //construct from a dense symmetric one
 	SparseMat & operator=(const SparseMat &rhs);
 	SparseMat & operator=(const T &a);          //assign a to diagonal
+	void identity() {*this = (T)1;};
     	SparseMat & operator+=(const T &a);         //assign a to diagonal
 	SparseMat & operator-=(const T &a);         //assign a to diagonal
         SparseMat & operator*=(const T &a);         //multiply by a scalar

t.cc

@@ -899,7 +899,7 @@ a(2,2)=-4;
 a(0,2)=1;
 cout <<a;
 double d;
-NRMat<double> c=inverse(a,&d);
+NRMat<double> c=calcinverse(a,&d);
 cout <<a<<c;
 }
 
@@ -1442,7 +1442,7 @@ for(int i=0; i<n; ++i) Y(i,i) -= tmp2(i,i);
 
 cout <<"Y matrix \n"<< Y;
 
-NRMat<double> vri = inverse(vr);
+NRMat<double> vri = calcinverse(vr);
 
 NRMat<double> numX = vri * vrd;
 cout <<" numerical X matrix \n"<< numX;
@@ -2087,15 +2087,26 @@ for(i=1; i<=n; ++i) v[i]=10.*i;
 cout <<v.permuted(p);
 }
 
-if(0)
+if(1)
 {
+int n;
 CyclePerm<int> c;
-cin>>c;
+cin>>c >>n;
+c.simplify();
 cout<<c<<endl;
 NRPerm<int> p(c);
 cout <<p;
 CyclePerm<int> cc(p);
 cout <<cc<<endl;
+NRPerm<int> pp(cc);
+cout <<(p == pp)<<endl;
+cout <<(c == cc)<<endl;
+CyclePerm<int> cc1,cc2;
+cc1=cc.pow_simple(n);
+cc2=cc.pow(n);
+cout<<cc<<endl;
+cout<<cc2<<endl;
+cout <<(cc1 == cc2)<<endl;
 }
 
 if(0)
@@ -2703,7 +2714,7 @@ diagonalize(d,e);
 //now try similarity transformed H
 NRMat<double> a(n,n);
 a.randomize(1);
-NRMat<double> ainv = inverse(a);
+NRMat<double> ainv = calcinverse(a);
 NRMat<double> atr = a.transpose();
 NRMat<double> ainvtr = ainv.transpose();
 cout <<"error of inverse = "<<(a*ainv).norm(1.)<<endl;
@@ -2748,7 +2759,7 @@ if(0)
 laerror("test exception");
 }
 
-if(1)
+if(0)
 {
 NRSMat<char> adj;
 cin >>adj;

vecmat3.h

@@ -114,6 +114,8 @@ public:
 	Mat3(void) {};
 	Mat3(const T (&a)[3][3]) {memcpy(q,a,3*3*sizeof(T));}
 	Mat3(const T x) {memset(q,0,9*sizeof(T)); q[0][0]=q[1][1]=q[2][2]=x;}; //scalar matrix
+	Mat3& operator=(const T &x) {memset(q,0,9*sizeof(T)); q[0][0]=q[1][1]=q[2][2]=x; return *this;}; //scalar matrix
+	void indentity() {*this = (T)1;};
 	Mat3(const T* x) {memcpy(q,x,9*sizeof(T));} 
 	Mat3(const T x00, const T x01,const T x02,const T x10,const T x11,const T x12,const T x20,const T x21,const T x22) 
 		{q[0][0]=x00; q[0][1]=x01; q[0][2]=x02; q[1][0]=x10; q[1][1]=x11; q[1][2]=x12; q[2][0]=x20; q[2][1]=x21; q[2][2]=x22;};