From 63b9e8c92ea76dfab0383d80da204a5f14117ca1 Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Wed, 28 Oct 2026 21:53:40 +0100 Subject: [PATCH] bugfixes in permutation --- permutation.cc | 33 ++++++++++++++++++++------------- permutation.h | 12 +++++++----- t.cc | 2 +- 3 files changed, 28 insertions(+), 19 deletions(-) diff --git a/permutation.cc b/permutation.cc index 82d2ca1..cd95b25 100644 --- a/permutation.cc +++ b/permutation.cc @@ -110,6 +110,7 @@ return smallest; template void NRPerm::inverseme() { +this->copyonwrite(); for(int i=1; i<=size(); ++i) if(cycle_leader(i)==i) reverse_cycle(i); } @@ -365,7 +366,7 @@ while (i < j) { return true; } -//Algorithm L from Knuth's volume 4 +//Algorithm P (plain changes) from Knuth's volume 4; parity of the permutations alternates template PERM_RANK_TYPE NRPerm::generate_all(void (*callback)(const NRPerm&), int select) { @@ -381,7 +382,7 @@ PERM_RANK_TYPE sumperm=0; p2: sumperm++; - if(callback) {if(!select || (select&1) == (sumperm&1)) {(*callback)(*this); this->copyonwrite();}} + if(callback) {if(!select || (select>0) == (sumperm&1)) {(*callback)(*this); this->copyonwrite();}} //odd sumperm = even permutation j=n; s=0; p4: q=c[j]+d[j]; @@ -409,11 +410,12 @@ template static PermutationAlgebra *list_all_return; static PERM_RANK_TYPE list_all_index; +static int list_all_parity; template static void list_all_callback(const NRPerm &p) { -(*list_all_return)[list_all_index].weight= (list_all_index&1)?-1:1; +(*list_all_return)[list_all_index].weight= list_all_parity? list_all_parity : ((list_all_index&1)?-1:1); //weight = parity (*list_all_return)[list_all_index].perm=p; (*list_all_return)[list_all_index].perm.copyonwrite(); ++list_all_index; @@ -427,6 +429,7 @@ if(parity_sel) number/=2; PermutationAlgebra ret(number); list_all_return = &ret; list_all_index=0; +list_all_parity=parity_sel>0?1:(parity_sel<0?-1:0); this->generate_all(list_all_callback,parity_sel); return ret; } @@ -460,7 +463,7 @@ while(1) T t; sumperm++; - if(callback) (*callback)(*this); + if(callback) {(*callback)(*this); this->copyonwrite();} j=this->size()-1; while(j>0 && (*this)[j]>=(*this)[j+1]) j--; if(j==0) break; /*finished*/ @@ -509,7 +512,7 @@ if(n<= _n) else { _sumperm++; - if(_callback) (*_callback)(*_perm); + if(_callback) {(*_callback)(*_perm); _perm->copyonwrite();} } } @@ -549,7 +552,7 @@ PERM_RANK_TYPE np=0; this->identity(); do{ ret[np].perm = *this; ret[np].perm.copyonwrite(); -ret[np].weight=0; +ret[np].weight=1; //unnormalized symmetrizer ++np; }while(this->next()); return ret; @@ -731,18 +734,19 @@ for(j=1;j<=n;++j) /*number of elements left from j bigger >j*/ } break; default: laerror("illegal type in inversions"); +} if(prank) { if(type!=3) laerror("rank can be computed from inversions only for type 3"); - l=1; + *prank=0; + PERM_RANK_TYPE f=1; for(j=1;j CyclePerm::operator*(const CyclePerm &q) const { int m=this->max(); int mm=q.max(); -if(mm>m) mm=m; +if(mm>m) m=mm; NRPerm qq(q,m); NRPerm pp(*this,m); NRPerm rr=pp*qq; @@ -1045,6 +1049,7 @@ return tt*r; template void CyclePerm::simplify(bool keep1) { +this->copyonwrite(); int j=0; for(int i=1; i<=this->size(); ++i) { @@ -1365,7 +1370,7 @@ template void partgen(int remain, int pos) { int hi,lo; -if(remain==0) {++partitioncount; if(_pcallback) (*_pcallback)(*partition); return;} +if(remain==0) {++partitioncount; if(_pcallback) {(*_pcallback)(*partition); partition->copyonwrite();} return;} if(partitiontyp) lo=(remain+partitiontypabs-pos)/(partitiontypabs-pos+1); else lo=1; hi=remain; if(partitiontyp>0) hi -= partitiontypabs-pos; @@ -1764,7 +1769,7 @@ if(ncyc> _ncycles) { if(_callyoung) { - if(_young_callback) _young_callback(*_tchi); + if(_young_callback) {_young_callback(*_tchi); _tchi->copyonwrite();} _character ++; } else _character += _tchi->character_contribution(_ncycles); @@ -2208,11 +2213,13 @@ for(int i=1; i template template //this works void NRPerm::testik(U u) { } +#endif /***************************************************************************//** * forced instantization in the corresponding object file diff --git a/permutation.h b/permutation.h index f4a4cd4..17603f1 100644 --- a/permutation.h +++ b/permutation.h @@ -64,7 +64,7 @@ public: bool is_identity() const; CompressedPartition cycles() const {return CyclePerm(*this).cycles(size());}; NRPerm inverse() const; - void reverse_cycle(int i); //aux for inverseme + void reverse_cycle(int i); //aux for inverseme, no copyonwrite (caller's responsibility) int cycle_leader(int i) const; //aux for inverseme void inverseme(); //in-place NRPerm reverse() const; //backward order @@ -81,10 +81,12 @@ public: int parity() const; //returns +/- 1 void randomize(void); //uniformly random by Fisher-Yates shuffle bool next(); //generate next permutation in lex order - PERM_RANK_TYPE generate_all(void (*callback)(const NRPerm&), int parity_select=0); //Algorithm L from Knuth's vol.4, efficient but not in lex order! + PERM_RANK_TYPE generate_all(void (*callback)(const NRPerm&), int parity_select=0); //Algorithm P (plain changes) from Knuth's vol.4, efficient but not in lex order! parity_select +1 even, -1 odd permutations only PermutationAlgebra list_all(int parity_select=0); PermutationAlgebra list_all_lex(); +#if 0 //reminder how a member template of a class template is defined in .cc template void testik(U u); //demostrate that in .cc one has to declare the nested templates separately +#endif PERM_RANK_TYPE generate_all_multi(void (*callback)(const NRPerm&)); //Algorithm L2 from Knuth's vol.4, for a multiset (repeated numbers, not really permutations) PERM_RANK_TYPE generate_all2(void (*callback)(const NRPerm&)); //recursive method, also not lexicographic PERM_RANK_TYPE generate_all_lex(void (*callback)(const NRPerm&)); //generate in lex order using next() @@ -107,8 +109,8 @@ NRVec_from1 applypermutation(const NRPerm &p, const NRVec_from1 &set, b if(p.size()!=set.size()) laerror("size mismatch in applypermutation"); #endif NRVec_from1 r(set.size()); -if(inverse) for(int i=1; i<=p.size(); ++i) r[p[i]] = set[i]; -else for(int i=1; i<=p.size(); ++i) r[i] = set[p[i]]; +if(inverse) for(int i=1; i<=p.size(); ++i) r[i] = set[p[i]]; //same convention as NRVec::permuted +else for(int i=1; i<=p.size(); ++i) r[p[i]] = set[i]; return r; } @@ -208,7 +210,7 @@ public: PermutationAlgebra operator|(const PermutationAlgebra &rhs) const; //each term with each PermutationAlgebra operator*(const PermutationAlgebra &rhs) const; //each term with each PermutationAlgebra operator+(const PermutationAlgebra &rhs) const; - PermutationAlgebra &operator*=(const R &x) {this->copyonwrite(); for(int i=1; i<=size(); ++i) (*this)[i].weight *= x; return *this;}; + PermutationAlgebra &operator*=(const R &x) {this->copyonwrite(); for(int i=0; i::normtype thr=0) {NRVec_simplify(*this,thr);}; bool operator==(PermutationAlgebra &rhs); //do NOT inherit from NRVec, as the underlying one ignores weights for the simplification; also we have to simplify before comparison diff --git a/t.cc b/t.cc index eb93a4b..d75b9d1 100644 --- a/t.cc +++ b/t.cc @@ -3157,7 +3157,7 @@ cout <<"PA = "< pindices=applypermutation(a[i].perm,indices,false); + NRVec_from1 pindices=applypermutation(a[i].perm,indices,true); for(int j=1; j<=pindices.size(); ++j) cout <