continueing on partitions

permutation.cc

@@ -602,6 +602,23 @@ return (n_even_cycles&1)?-1:1;
 }
 
 
+template <typename T>
+PERM_RANK_TYPE  CyclePerm<T>::order() const
+{
+#ifdef DEBUG
+if(!this->is_valid()) laerror("operation with an invalid cycleperm");
+#endif
+PERM_RANK_TYPE r=1;
+T ncycles=this->size();
+for(T i=1; i<=ncycles; ++i)
+        {
+        T length=(*this)[i].size();
+	r=lcm(r,(PERM_RANK_TYPE)length);
+        }
+return r;
+}
+
+
 template <typename T>
 CompressedPartition<T> CyclePerm<T>::cycles(const T n) const
 {
@@ -771,6 +788,20 @@ for(int i=1;i<=n;++i)
 return r;
 }
 
+template <typename T>
+Partition<T>::Partition(const YoungTableaux<T> &x)
+{
+#ifdef DEBUG
+if(x.is_valid()) laerror("operation with an invalid tableaux");
+#endif
+int nparts=x.size();
+int n=0;
+for(int i=1; i<=nparts; ++i) n+= x[i].size();
+this->resize(n);
+this->clear();
+for(int i=1; i<=nparts; ++i) (*this)[i]=x[i].size();
+}
+
 PERM_RANK_TYPE ipow(PERM_RANK_TYPE x, int i)
 {
 if(i<0) return 0;
@@ -846,6 +877,58 @@ return s;
 
 
 
+/*see M. Aigner - Combinatorial Theory - number of partitions of number n*/
+#define MAXPART 260
+#define MIN(x,y) ((x)<(y)?(x):(y))
+/* here even few more columns more could be saved, and also the property part(n,k>=n/2)=part(n-k) could be used */
+#define partind(n,k) (((n)-3)*((n)-2)/2+(k)-2)
+#define getpart(n,k) (((k)>(n) || (k)<=0 )?0:(((k)==(n) || (k)==1)?1:(part[partind((n),(k))])))
+
+PERM_RANK_TYPE partitions(int n, int k)
+{
+static PERM_RANK_TYPE part[partind(MAXPART+1,1)]={1};
+static PERM_RANK_TYPE npart[MAXPART+1]={0,1,2,3};
+static int ntop=3;
+int i,j,l;
+PERM_RANK_TYPE s;
+
+if( n<=0 || k<-1 ) laerror("Partitions are available only for positive numbers");
+if(k== -1) /* total partitions */
+	{
+	if(n>ntop)(void)partitions(n,3); /* make sure that table is calculated */
+	return npart[n];
+	}
+
+if(n==k||k==1) return 1;
+if(k==0||k>n) return 0;
+if(k==2) return n/2;
+if(n>MAXPART && k>=n/2) return(partitions(n-k,-1)); 
+if( n>MAXPART) laerror("sorry, too big argument to partition enumerator");
+if(n<=ntop) return(getpart(n,k));
+
+/*calculate, if necessary, few next rows*/
+for(i=4;i<=n;i++)
+    {
+    npart[i]=2;
+    for(j=2;j<=i-1;j++)
+	{
+	s=0;
+	for(l=1; l<=MIN(i-j,j);l++) s+=getpart(i-j,l);
+	npart[i]+=(part[partind(i,j)]=s);
+	}
+    }
+
+ntop=n;
+return(part[partind(n,k)]);
+
+}
+
+#undef getpart
+#undef partind
+#undef MIN
+
+
+
 //////////////////////////////////////////////////////////////////////////////////////
 
 template <typename T>
@@ -860,6 +943,65 @@ this->resize(nlines);
 for(int i=1; i<=nlines; ++i) (*this)[i].resize(frame[i]);
 }
 
+template <typename T>
+bool YoungTableaux<T>::is_valid() const
+{
+int nrows=this->size();
+for(int i=1; i<=nrows; ++i)
+	{
+	if((*this)[i].size()<=0) return false;
+	if(i>1 && (*this)[i].size()> (*this)[i-1].size()) return false;
+	}
+return true; 
+}
+
+
+template <typename T>
+int YoungTableaux<T>::sum() const
+{
+#ifdef DEBUG
+if(is_valid()) laerror("invalid young frame");
+#endif
+int sum=0;
+int nrows=this->size();
+for(int i=1; i<=nrows; ++i) sum += (*this)[i].size();
+return sum;
+}
+
+template <typename T>
+bool YoungTableaux<T>::is_standard() const
+{
+#ifdef DEBUG
+if(!is_valid()) laerror("invalid young frame");
+#endif
+int nrows=this->size();
+//check numbers monotonous in rows and columns
+for(int i=1; i<=nrows; ++i)
+	{
+	for(int j=1; j<=(*this)[i].size(); ++j) 
+		{
+		if(j>1 && (*this)[i][j] < (*this)[i][j-1]) return false;
+		if(i>1 && (*this)[i][j] < (*this)[i-1][j]) return false;
+		}
+	}
+return true;
+}
+
+template <typename T>
+std::ostream & operator<<(std::ostream &s, const YoungTableaux<T> &x)
+{
+int nrows=x.size();
+for(int i=1; i<=nrows; ++i)
+        {
+	for(int j=1; j<=x[i].size(); ++j)
+		{
+		s.width(4);
+		s<<x[i][j]<<' ';
+		}
+	s<<std::endl;
+	}
+return s;
+}
 
 
 
@@ -876,6 +1018,7 @@ template class YoungTableaux<int>;
 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 CompressedPartition<T> &x); \
+template std::ostream & operator<<(std::ostream &s, const YoungTableaux<T> &x); \
 
 
 

permutation.h

@@ -30,10 +30,10 @@ typedef   unsigned long long PERM_RANK_TYPE;
 namespace LA {
 
 //forward declaration
-template <typename T> class NRVec_from1;
 template <typename T> class CyclePerm;
 template <typename T> class Partition;
 template <typename T> class CompressedPartition;
+template <typename T> class YoungTableaux;
 
 template <typename T>
 class NRPerm : public NRVec_from1<T> {
@@ -83,10 +83,40 @@ public:
 	CompressedPartition<T> cycles(const T n) const;
 	void readfrom(const std::string &line);
 	CyclePerm operator*(const CyclePerm q) const; //q is rhs and applied first, this applied second
-	//@@@order = lcm of cycle lengths
+	PERM_RANK_TYPE order() const; //lcm of cycle lengths
 };
 
 
+template <typename T>
+T gcd(T big, T small)
+{
+if(big==0)
+        {
+        if(small==0) laerror("bad arguments in gcd");
+        return small;
+        }
+if(small==0) return big;
+if(small==1||big==1) return 1;
+
+T help;
+if(small>big) {help=big; big=small; small=help;}
+do      {
+        help=small;
+        small= big%small;
+        big=help;
+        }
+while(small != 0);
+return big;
+}
+
+
+template <typename T> 
+inline T lcm(T a, T b)
+{
+return (a*b)/gcd(a,b);
+}
+
+
 template <typename T>
 std::istream & operator>>(std::istream &s, CyclePerm<T> &x);
 
@@ -123,12 +153,12 @@ public:
         T nparts() const {T s=0; for(int i=1; i<=this->size(); ++i) if((*this)[i]) ++s; return s;}
         bool is_valid() const {if(this->size() != this->sum()) return false; for(int i=2; i<=this->size(); ++i) if((*this)[i]>(*this)[i-1]) return false; return true; }
 	explicit Partition(const CompressedPartition<T>  &rhs) : NRVec_from1<T>(rhs.size()) {this->clear(); int ithru=0; for(int i=rhs.size(); i>=1; --i) for(int j=0; j<rhs[i]; ++j) (*this)[++ithru]=i; }
+	explicit Partition(const YoungTableaux<T> &x); //extract a partition  as a shape of Young tableaux
 	Partition adjoint() const; //also called conjugate partition
 	PERM_RANK_TYPE Sn_irrep_dim() const;
 	PERM_RANK_TYPE generate_all(void (*callback)(const Partition<T>&), int nparts=0); //nparts <0 means at most to -nparts
 
 
-//@@@yamanouchi symbol
 };
 
 
@@ -138,14 +168,27 @@ public:
         YoungTableaux() :  NRVec_from1<NRVec_from1<T> >() {};
         explicit YoungTableaux(const Partition<T> &frame);
 
-	bool is_valid() const; //@@@shape forms a partition
+	bool is_valid() const; //check whether its shape forms a partition
-        bool filled_correctly() const; //is it filled correctly@@@
+	int nrows() const {return this->size();}
+	int ncols() const {return (*this)[1].size();}
+        bool is_standard() const; //is it filled in standard way (possibly with repeated numbers)
+	//@@@???void clear(); //clear the numbers but keep shape (different than inherited clear())
+	int sum() const; //get back sum of the partition
+	NRVec_from1<T> yamanouchi() const; //@@@yamanouchi symbol
+		//@@@ ??>young operator as a linear comb of permutations - maybe a class for itself, i.e. element of the Sn group algebra? or maybe as a vector with index being the rank of the permutation(n! length) or as a sparsemat thereof or maybe a list???
+		//@@@???action of group algebra elements on vectors and matrices
 };
 
+template <typename T>
+std::ostream & operator<<(std::ostream &s, const YoungTableaux<T> &x);
 
 
-//@@@enumerator of partitions of n to r parts and total from expr
+extern PERM_RANK_TYPE partitions(int n, int k= -1); //enumerate partitions to k parts; k== -1 for total # of partitions
-//@@@Sn character table computation from young - young frame filling
+
+
+//@@@Sn character table computation from young - young frame filling - routine for one character of one irrep and another for generation of the whole group table (maybe class for a group table too)
+//
+//@@@@Un irrep dimensions from gelfand
 
 }//namespace
 #endif

t.cc

@@ -77,6 +77,18 @@ cout<<"IR dim "<<p.Sn_irrep_dim()<<endl;
 CompressedPartition qc(q);
 cout <<"("<<qc<<')';
 cout<<" Class size "<<qc.Sn_class_size()<<endl;
+
+int nn=p.size();
+YoungTableaux<int> y(p);
+for(int i=1; i<=y.nrows(); ++i)
+        {
+        for(int j=1; j<=y[i].size(); ++j) y[i][j]= nn--;
+	}
+cout <<y;
+YoungTableaux<int> yy(y);
+yy.clear();
+cout<<yy;
+cout <<y;
 }
 
 
@@ -2113,8 +2125,10 @@ cin >>n;
 Partition<int> p(n);
 int tot=p.generate_all(pprintme,0);
 cout <<"generated "<<tot<<endl;
+if(tot!=partitions(n)) laerror("internal error in partition generation or enumerations");
 }
 
 
 
+
 }