From 8bc2179f3b8ad8ac81a8fe516bcc84743befedc0 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Wed, 17 Jan 2024 22:53:59 +0100
Subject: [PATCH] permutation - implemented general antisymmetrizer

---
 permutation.cc | 40 ++++++++++++++++++++++++++++++----------
 permutation.h  | 12 +++++++-----
 t.cc           | 26 ++++++++++++++++++++++++++
 3 files changed, 63 insertions(+), 15 deletions(-)

diff --git a/permutation.cc b/permutation.cc
index 0efdd46..79b3f0e 100644
--- a/permutation.cc
+++ b/permutation.cc
@@ -440,8 +440,6 @@ static NRPerm<T> *_perm2;
 template <typename T>
 static const T *pclasses2;
 
-template <typename T>
-static T *_pperm;
 
 static int sameclass;
 
@@ -454,7 +452,7 @@ if(n<=_n2<T>)
         T c=pclasses2<T>[n];
         for(i=1;i<=_n2<T>;i++)
                 {
-                if(_pperm<T>[i]>=1) goto skipme; /*place already occupied*/
+                if((*_perm2<T>)[i]>=1) goto skipme; /*place already occupied*/
                 if (sameclass)
                         {
                         if(sameclass==1 && c!=pclasses2<T>[i]) goto skipme;
@@ -463,14 +461,14 @@ if(n<=_n2<T>)
                                 {
                                 T j,k;
                                 for(j=i+1; j<=_n2<T>; j++)
-                                        if((k=_pperm<T>[j])>=1)
+                                        if((k=(*_perm2<T>)[j])>=1)
                                                  if(/* automatically fulfilled k<n &&*/ c==pclasses2<T>[k]) goto skipme;
                                 }
                         }
                 /*put it there and next number*/
-                _pperm<T>[i]=n;
+                (*_perm2<T>)[i]=n;
                 permg2(n+1);
-                _pperm<T>[i]=0;
+                (*_perm2<T>)[i]=0;
                 skipme:;
                 }
         }
@@ -493,7 +491,6 @@ _n2<T> = this->size();
 _callback2<T> =callback;
 _sumperm2=0;
 _perm2<T> = this;
-_pperm<T> = &(*this)[1]-1;
 pclasses2<T> = &classes[1]-1;
 sameclass=restriction_type;
 permg2<T>(1);
@@ -505,22 +502,25 @@ static PermutationAlgebra<T,T>  *list_restricted_return;
 
 static PERM_RANK_TYPE list_restricted_index;
 
+static bool list_restricted_invert;
+
 template <typename T> 
 static void list_restricted_callback(const NRPerm<T> &p)
 {
 (*list_restricted_return<T>)[list_restricted_index].weight=p.parity();
-(*list_restricted_return<T>)[list_restricted_index].perm=p;
-(*list_restricted_return<T>)[list_restricted_index].copyonwrite();
+(*list_restricted_return<T>)[list_restricted_index].perm= list_restricted_invert?p.inverse():p;
+(*list_restricted_return<T>)[list_restricted_index].perm.copyonwrite();
 ++list_restricted_index;
 }
 
 template <typename T>
-PermutationAlgebra<T,T>  NRPerm<T>::list_restricted(const NRVec_from1<T> &classes, int restriction_type)
+PermutationAlgebra<T,T>  NRPerm<T>::list_restricted(const NRVec_from1<T> &classes, int restriction_type, bool invert)
 {
 PERM_RANK_TYPE number = this->generate_restricted(NULL,classes,restriction_type);
 PermutationAlgebra<T,T>  ret(number);
 list_restricted_return<T> = &ret;
 list_restricted_index=0;
+list_restricted_invert=invert;
 generate_restricted(list_restricted_callback,classes,restriction_type);
 return ret;
 }
@@ -2040,6 +2040,25 @@ return r;
 }
 
 
+template<typename T>
+PermutationAlgebra<T,T> general_antisymmetrizer(const NRVec<NRVec_from1<T> > &groups, int restriction_type, bool inverted)
+{
+PermutationAlgebra<T,T> r;
+int ngroups=groups.size();
+if(ngroups==0) return r;
+NRVec<PermutationAlgebra<T,T> > lists(ngroups);
+for(int i=0; i<ngroups; ++i)
+	{
+	int ni = groups[i].size();
+	NRPerm<T> tmp(ni);
+	lists[i] = tmp.list_restricted(groups[i],restriction_type,inverted);
+	}
+//cross-product the lists
+r=lists[0];
+for(int i=1; i<ngroups; ++i) r= r&lists[i];
+return r;
+}
+
 
 /***************************************************************************//**
  * forced instantization in the corresponding object file
@@ -2053,6 +2072,7 @@ template class Partition<T>; \
 template class YoungTableaux<T>; \
 template class Sn_characters<T>; \
 template class CycleIndex<T>; \
+template PermutationAlgebra<T,T> general_antisymmetrizer(const NRVec<NRVec_from1<T> > &groups, int, bool); \
 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); \
diff --git a/permutation.h b/permutation.h
index 596c1d3..a43c5f3 100644
--- a/permutation.h
+++ b/permutation.h
@@ -76,7 +76,7 @@ public:
 	PERM_RANK_TYPE generate_all2(void (*callback)(const NRPerm<T>&)); //recursive method, also not lexicographic
 	PERM_RANK_TYPE generate_all_lex(void (*callback)(const NRPerm<T>&)); //generate in lex order using next()
 	PERM_RANK_TYPE generate_restricted(void (*callback)(const NRPerm<T>&), const NRVec_from1<T> &classes, int restriction_type=0);
-	PermutationAlgebra<T,T>  list_restricted(const NRVec_from1<T> &classes, int restriction_type=0); //weight is set to parity (antisymmetrizer) by default
+	PermutationAlgebra<T,T>  list_restricted(const NRVec_from1<T> &classes, int restriction_type=0, bool inverted=false); //weight is set to parity (antisymmetrizer) by default
 	PERM_RANK_TYPE rank() const; //counted from 0 to n!-1
 	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
@@ -94,8 +94,8 @@ NRVec_from1<X> applypermutation(const NRPerm<T> &p, const NRVec_from1<X> &set, b
 if(p.size()!=set.size()) laerror("size mismatch in applypermutation");
 #endif
 NRVec_from1<X> 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[p[i]] = set[i];
+else for(int i=1; i<=p.size(); ++i) r[i] = set[p[i]];
 return r;
 }
 
@@ -190,8 +190,6 @@ public:
 };
 
 
-//TODO@@@ permutationalgebra  for Kucharski style antisymmetrizers
-
 
 extern PERM_RANK_TYPE factorial(const int n);
 extern PERM_RANK_TYPE binom(int n, int k);
@@ -392,5 +390,9 @@ return r;
 }
 
 
+template<typename T>
+PermutationAlgebra<T,T> general_antisymmetrizer(const NRVec<NRVec_from1<T> > &groups, int restriction_type=0, bool inverted=false);
+
+
 }//namespace
 #endif
diff --git a/t.cc b/t.cc
index d706401..5f6cf38 100644
--- a/t.cc
+++ b/t.cc
@@ -3083,8 +3083,34 @@ if(!dif.is_zero()) laerror("error in gf 2^n sqrt");
 }
 
 
+if(0)
+{
+int n;
+cin>>n;
+NRPerm<int> p(n);
+PermutationAlgebra<int,int> a=p.list_all();
+for(int i=0; i<a.size(); ++i)
+	if(a[i].weight!=a[i].perm.parity()) laerror("internal error in parity");
+}
+
 if(1)
 {
+//this is kucharskiP antisymmetrizer just without parsing the input
+NRVec<NRVec_from1<int> > groups;
+cin >> groups;
+int ntot=0;
+for(int i=0; i<groups.size(); ++i) ntot+=groups[i].size();
+NRVec_from1<string> indices(ntot);
+cin >>indices;
+
+PermutationAlgebra<int,int> a=general_antisymmetrizer(groups,-2,true);
+for(int i=0; i<a.size(); ++i)
+	{
+	cout << (a[i].weight==1?'+':'-');
+	NRVec_from1<string> pindices=applypermutation(a[i].perm,indices,false);
+	for(int j=1; j<=pindices.size(); ++j) cout <<pindices[j];
+	cout <<endl;
+	}
 }
 
 }