some more functions in polynomials and permutations

permutation.cc

@@ -471,6 +471,48 @@ return a[n];
 }
 
 
+#define MAXBINOM 100
+#define ibidxmaxeven(n) ((n-2)*(n-2)/4)
+#define ibidx(n,k) (k-2+(n&1?(n-3)*(n-3)/4:(n/2-1)*(n/2-2)))
+PERM_RANK_TYPE binom(int n, int k)
+{
+PERM_RANK_TYPE p,value;
+register int d;
+static PERM_RANK_TYPE ibitab[ibidxmaxeven(MAXBINOM)]= /*only nontrivial are stored,
+                                            zero initialization by compiler assumed*/
+{
+6,
+10,
+15,20,
+21,35,
+28,56,70
+};
+
+
+if(k>n||k<0) return(0);
+if(k>n/2) k=n-k;
+if(k==0) return(1);
+if(k==1) return(n);
+int ind=0;
+if(n<=MAXBINOM)
+        {
+        ind=ibidx(n,k);
+        if (ibitab[ind]) return ibitab[ind];
+        }
+/* nonrecurent method used anyway */
+d=n-k;
+p=1;
+for(;n>d;n--) p *= n;
+value=p/factorial(k);
+if(n<=MAXBINOM) ibitab[ind]=value;
+return value;
+}
+#undef ibidx
+#undef ibidxmaxeven
+#undef MAXBINOM
+
+
+
 ////////////////////////////////////////////////////////
 
 template <typename T>
@@ -1599,6 +1641,40 @@ for(int i=1; i<=c.chi.nrows(); ++i)
 return s;
 }
 
+
+template <typename T>
+bool CycleIndex<T>::is_valid() const
+{
+if(classes.size()!=classsizes.size()) return false;
+T n=classes[1].sum();
+for(int i=2; i<=classes.size(); ++i) 
+	{
+	if(classes[i].sum()!=n) return false;
+	}
+return true;
+}
+
+template <typename T>
+Polynomial<T> CycleIndex<T>::substitute(const Polynomial<T> &p, PERM_RANK_TYPE *denom) const
+{
+Polynomial<T> r(0);
+r[0]=(T)0;
+*denom =0;
+
+for(int i=1; i<=classes.size(); ++i)
+	{
+	Polynomial<T> term(0);
+	term[0]=(T)1;
+	for(int j=1; j<=classes[i].size(); ++j) if(classes[i][j]) term *= (p.powx(j)).pow((int)classes[i][j]);
+	r += term*classsizes[i];
+	*denom += classsizes[i];
+	}
+
+return r;
+}
+
+
+
 /***************************************************************************//**
  * forced instantization in the corresponding object file
  ******************************************************************************/
@@ -1610,6 +1686,7 @@ template class CompressedPartition<T>; \
 template class Partition<T>; \
 template class YoungTableaux<T>; \
 template class Sn_characters<T>; \
+template class CycleIndex<T>; \
 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); \

permutation.h

@@ -22,6 +22,7 @@
 
 #include "la_traits.h"
 #include "vec.h"
+#include "polynomial.h"
 
 typedef   unsigned long long PERM_RANK_TYPE;
 
@@ -65,6 +66,7 @@ public:
 };
 
 extern PERM_RANK_TYPE factorial(const int n);
+extern PERM_RANK_TYPE binom(int n, int k);
 extern PERM_RANK_TYPE longpow(PERM_RANK_TYPE x, int i);
 
 //permutations represented in the cycle format
@@ -212,6 +214,20 @@ NRMat_from1<T> chi; //characters
 	bool is_valid() const; //check internal consistency
 };
 
+template <typename T>  class Polynomial; //forward declaration
+
+template <typename T>
+class CycleIndex {
+public:
+NRVec_from1<CompressedPartition<T> > classes;
+NRVec_from1<PERM_RANK_TYPE> classsizes;
+
+	CycleIndex(const Sn_characters<T> &rhs): classes(rhs.classes),classsizes(rhs.classsizes) {};
+	bool is_valid() const; //check internal consistency
+	Polynomial<T> substitute(const Polynomial<T> &p, PERM_RANK_TYPE *denom) const;
+
+};
+
 template <typename T>
 extern std::ostream & operator<<(std::ostream &s, const Sn_characters<T> &c);
 

polynomial.cc

@@ -90,6 +90,16 @@ laerror("roots not implemented for integer polynomials");
 return NRVec<typename LA_traits<int>::complextype>(1);
 }
 
+template<>
+NRVec<typename LA_traits<unsigned int>::complextype> Polynomial<unsigned int>::roots() const
+{
+laerror("roots not implemented for integer polynomials");
+return NRVec<typename LA_traits<unsigned int>::complextype>(1);
+}
+
+
+
+
 template <typename T>
 NRVec<typename LA_traits<T>::complextype> Polynomial<T>::roots() const
 {
@@ -107,7 +117,7 @@ NRVec<T> rr(degree());
 int nr=0;
 for(int i=0; i<degree(); ++i)
 	{
-	if(abs(r[i].imag())<thr)
+	if(MYABS(r[i].imag())<thr)
 		{
 		rr[nr++] = r[i].real();
 		}
@@ -151,7 +161,7 @@ T x=x0;
 for(int i=0; i<maxit; ++i)
 	{
 	T v=value(*this,x);
-	if(abs(v)<thr) break;
+	if(MYABS(v)<thr) break;
 	x -= v/value(d,x);
 	}
 
@@ -185,6 +195,15 @@ laerror("SVD gcd only for floating point numbers");
 return p;
 }
 
+template <>
+Polynomial<unsigned int> svd_gcd(const Polynomial<unsigned int> &p, const Polynomial<unsigned int> &q, const typename LA_traits<unsigned int>::normtype thr)
+{
+laerror("SVD gcd only for floating point numbers");
+return p;
+}
+
+
+
 template <typename T>
 Polynomial<T> svd_gcd(const Polynomial<T> &p, const Polynomial<T> &q, const typename LA_traits<T>::normtype thr)
 {
@@ -202,12 +221,54 @@ int d=big.degree()+small.degree()-rank;
 return poly_gcd(big,small,thr,d);
 }
 
+template <typename T>
+Polynomial<T> Polynomial<T>::pow(const int n) const
+{
+int i=n;
+if(i<0) laerror("negative exponent in polynomial::pow");
+if(i==0) {Polynomial<T> r(0); r[0]=(T)1; return r;}
+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>
+Polynomial<T> Polynomial<T>::powx(const int i) const
+{
+if(i<0) laerror("negative exponent in polynomial::pow");
+if(i==0) {Polynomial<T> r(0); r[0]= this->sum(); return r;}
+if(i==1) return *this;
+Polynomial<T> r(i*this->degree());
+r.clear();
+for(int j=0; j<=this->degree(); ++j) r[j*i] = (*this)[j];
+return r;
+}
+
+template <typename T>
+void Polynomial<T>::binomial(const int n)
+{
+resize(n,false);
+for(int i=0; i<=n; ++i) (*this)[i] = (T) binom(n,i);
+}
 
 
 /***************************************************************************//**
  * forced instantization in the corresponding object file
  ******************************************************************************/
 template class Polynomial<int>;
+template class Polynomial<unsigned int>;
 template class Polynomial<double>;
 template class Polynomial<std::complex<double> >;
 
@@ -220,6 +281,7 @@ template Polynomial<T> svd_gcd(const Polynomial<T> &p, const Polynomial<T> &q, c
 
 
 INSTANTIZE(int)
+INSTANTIZE(unsigned int)
 INSTANTIZE(double)
 INSTANTIZE(std::complex<double>)
 

polynomial.h

@@ -26,6 +26,13 @@
 
 namespace LA {
 
+template <typename T>
+inline  typename LA_traits<T>::normtype  MYABS(const T &x) {return abs(x);}
+
+template <>
+inline unsigned int MYABS(const  unsigned int &x) {return x;}
+
+
 template <typename T>
 class Polynomial : public NRVec<T> {
 public:
@@ -72,12 +79,13 @@ public:
 		for(int i=0; i<=rhs.degree(); ++i) for(int j=0; j<=degree(); ++j) r[i+j] += rhs[i]*(*this)[j];
 		return r;
 		};
+	Polynomial& operator*=(const Polynomial &rhs) {*this = (*this)*rhs; return *this;};
 	void simplify(const typename LA_traits<T>::normtype thr=0)
 		{
 		NOT_GPU(*this); 
 		this->copyonwrite();
 		int n=degree();
-		while(n>0 && abs((*this)[n])<=thr) --n;
+		while(n>0 && MYABS((*this)[n])<=thr) --n;
 		resize(n,true);
 		};
 	void normalize() {if((*this)[degree()]==(T)0) laerror("zero coefficient at highest degree - simplify first"); *this /= (*this)[degree()];};
@@ -148,10 +156,13 @@ public:
 	NRVec<typename LA_traits<T>::complextype> roots() const; //implemented for complex<double> and double only
 	NRVec<T> realroots(const typename LA_traits<T>::normtype thr) const;
 	T newton(const T x0, const typename LA_traits<T>::normtype thr=1e-14, const int maxit=1000) const; //solve root from the guess
+	Polynomial pow(const int i) const; //integer power
+	Polynomial powx(const int i) const; //substitute x^i for x in the polynomial
+	void binomial(const int n); //(1+x)^n
 
 };
 
-//this is very general, can be used also for nesting polynomials
+//this is very general, can be used also for composition (nesting) of polynomials
 //for an alternative algorithm which minimizes number of multiplications cf. also matexp.h
 template <typename T, typename C>
 C value(const Polynomial<T> &p, const C &x)
@@ -199,6 +210,8 @@ return p;
 }
 
 
+
+
 template <typename T>
 extern Polynomial<T> lagrange_interpolation(const NRVec<T> &x, const NRVec<T> &y);
 

t.cc

@@ -2158,13 +2158,19 @@ if(tot!=partitions(n)) laerror("internal error in partition generation or enumer
 if(space_dim!=longpow(unitary_n,n)) {cout<<space_dim<<" "<<ipow(unitary_n,n)<<endl;laerror("integer overflow or internal error in space dimensions");}
 }
 
-if(0)
+if(1)
 {
 int n;
 cin >>n ;
 Sn_characters<int> Sn(n);
 cout <<Sn;
 if(!Sn.is_valid()) laerror("internal error in Sn character calculation");
+Polynomial<int> p(1);
+p[0]=p[1]=1;
+CycleIndex<int> c(Sn);
+PERM_RANK_TYPE d;
+Polynomial<int> pp=c.substitute(p,&d);
+for(int i=0; i<=p.size(); ++i) if(d!=pp[i]) laerror("error in cycle index");
 }
 
 if(0)
@@ -2242,7 +2248,7 @@ Polynomial<double>pp=q.derivative(2);
 cout<<"test deriv. "<<(pp-p).norm()<<endl;
 }
 
-if(1)
+if(0)
 {
 int n;
 cin >>n;
@@ -2268,5 +2274,17 @@ cout <<"test gcd "<<(rr-rrr).norm()<<endl;
 cout <<"test svdgcd "<<(rr-rrrr).norm()<<endl;
 }
 
+if(0)
+{
+int n;
+cin >>n;
+Polynomial<double> p;
+p.binomial(n);
+Polynomial<double> q(1);
+q[0]=q[1]=1.;
+Polynomial<double> qq=q.pow(n);
+cout <<"test binom "<<(p-qq).norm()<<endl;
+}
+
 
 }