From 3288e51fba93e9bf0c301541e6857f96dd322cda Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Sat, 26 Jun 2021 22:41:40 +0200 Subject: [PATCH] some more functions in polynomials and permutations --- permutation.cc | 77 ++++++++++++++++++++++++++++++++++++++++++++++++++ permutation.h | 16 +++++++++++ polynomial.cc | 66 +++++++++++++++++++++++++++++++++++++++++-- polynomial.h | 17 +++++++++-- t.cc | 22 +++++++++++++-- 5 files changed, 192 insertions(+), 6 deletions(-) diff --git a/permutation.cc b/permutation.cc index b3be6a2..ee5cf7f 100644 --- a/permutation.cc +++ b/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 @@ -1599,6 +1641,40 @@ for(int i=1; i<=c.chi.nrows(); ++i) return s; } + +template +bool CycleIndex::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 +Polynomial CycleIndex::substitute(const Polynomial &p, PERM_RANK_TYPE *denom) const +{ +Polynomial r(0); +r[0]=(T)0; +*denom =0; + +for(int i=1; i<=classes.size(); ++i) + { + Polynomial 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; \ template class Partition; \ template class YoungTableaux; \ template class Sn_characters; \ +template class CycleIndex; \ template std::istream & operator>>(std::istream &s, CyclePerm &x); \ template std::ostream & operator<<(std::ostream &s, const CyclePerm &x); \ template std::ostream & operator<<(std::ostream &s, const CompressedPartition &x); \ diff --git a/permutation.h b/permutation.h index 3ab787d..9427835 100644 --- a/permutation.h +++ b/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 chi; //characters bool is_valid() const; //check internal consistency }; +template class Polynomial; //forward declaration + +template +class CycleIndex { +public: +NRVec_from1 > classes; +NRVec_from1 classsizes; + + CycleIndex(const Sn_characters &rhs): classes(rhs.classes),classsizes(rhs.classsizes) {}; + bool is_valid() const; //check internal consistency + Polynomial substitute(const Polynomial &p, PERM_RANK_TYPE *denom) const; + +}; + template extern std::ostream & operator<<(std::ostream &s, const Sn_characters &c); diff --git a/polynomial.cc b/polynomial.cc index 80c4852..7d531d5 100644 --- a/polynomial.cc +++ b/polynomial.cc @@ -90,6 +90,16 @@ laerror("roots not implemented for integer polynomials"); return NRVec::complextype>(1); } +template<> +NRVec::complextype> Polynomial::roots() const +{ +laerror("roots not implemented for integer polynomials"); +return NRVec::complextype>(1); +} + + + + template NRVec::complextype> Polynomial::roots() const { @@ -107,7 +117,7 @@ NRVec rr(degree()); int nr=0; for(int i=0; i +Polynomial svd_gcd(const Polynomial &p, const Polynomial &q, const typename LA_traits::normtype thr) +{ +laerror("SVD gcd only for floating point numbers"); +return p; +} + + + template Polynomial svd_gcd(const Polynomial &p, const Polynomial &q, const typename LA_traits::normtype thr) { @@ -202,12 +221,54 @@ int d=big.degree()+small.degree()-rank; return poly_gcd(big,small,thr,d); } +template +Polynomial Polynomial::pow(const int n) const +{ +int i=n; +if(i<0) laerror("negative exponent in polynomial::pow"); +if(i==0) {Polynomial r(0); r[0]=(T)1; return r;} +if(i==1) return *this; +Polynomialy,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 +Polynomial Polynomial::powx(const int i) const +{ +if(i<0) laerror("negative exponent in polynomial::pow"); +if(i==0) {Polynomial r(0); r[0]= this->sum(); return r;} +if(i==1) return *this; +Polynomial r(i*this->degree()); +r.clear(); +for(int j=0; j<=this->degree(); ++j) r[j*i] = (*this)[j]; +return r; +} + +template +void Polynomial::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; +template class Polynomial; template class Polynomial; template class Polynomial >; @@ -220,6 +281,7 @@ template Polynomial svd_gcd(const Polynomial &p, const Polynomial &q, c INSTANTIZE(int) +INSTANTIZE(unsigned int) INSTANTIZE(double) INSTANTIZE(std::complex) diff --git a/polynomial.h b/polynomial.h index f06eebc..73e9678 100644 --- a/polynomial.h +++ b/polynomial.h @@ -26,6 +26,13 @@ namespace LA { +template +inline typename LA_traits::normtype MYABS(const T &x) {return abs(x);} + +template <> +inline unsigned int MYABS(const unsigned int &x) {return x;} + + template class Polynomial : public NRVec { 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::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::complextype> roots() const; //implemented for complex and double only NRVec realroots(const typename LA_traits::normtype thr) const; T newton(const T x0, const typename LA_traits::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 C value(const Polynomial &p, const C &x) @@ -199,6 +210,8 @@ return p; } + + template extern Polynomial lagrange_interpolation(const NRVec &x, const NRVec &y); diff --git a/t.cc b/t.cc index d262c0f..f142500 100644 --- a/t.cc +++ b/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<>n ; Sn_characters Sn(n); cout < p(1); +p[0]=p[1]=1; +CycleIndex c(Sn); +PERM_RANK_TYPE d; +Polynomial 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 @@ Polynomialpp=q.derivative(2); cout<<"test deriv. "<<(pp-p).norm()<>n; @@ -2268,5 +2274,17 @@ cout <<"test gcd "<<(rr-rrr).norm()<>n; +Polynomial p; +p.binomial(n); +Polynomial q(1); +q[0]=q[1]=1.; +Polynomial qq=q.pow(n); +cout <<"test binom "<<(p-qq).norm()<