Support for non-commuting fieldsv1.3

1 files changed, 53 insertions, 7 deletions

diff --git a/src/polynomial.c b/src/polynomial.c
index 639728a..d082b60 100644
--- a/src/polynomial.c
+++ b/src/polynomial.c
@@ -688,7 +688,7 @@ int polynomial_simplify(Polynomial* polynomial, Fields_Table fields){
   // sort monomials and factors
   for(i=0;i<(*polynomial).length;i++){
     sign=1;
-    monomial_sort((*polynomial).monomials[i],0,(*polynomial).monomials[i].length-1,fields,&sign);
+    monomial_sort((*polynomial).monomials[i],fields,&sign);
     number_Qprod_chain(quot(sign,1),(*polynomial).nums+i);
     int_array_sort((*polynomial).factors[i],0,(*polynomial).factors[i].length-1);
   }
@@ -787,7 +787,30 @@ int exchange_polynomial_terms(int i, int j, Polynomial* polynomial){
 }
 
 // sort a monomial (with sign coming from exchanging two Fermions)
-int monomial_sort(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign){
+// if the monomial contains noncommuting elements, put them at the beginning of the monomial
+int monomial_sort(Int_Array monomial, Fields_Table fields, int* sign){
+  int i,j;
+  int tmp;
+  // first index after noncommuting indices
+  int post_nc=0;
+
+  for(i=0;i<monomial.length;i++){
+    if(is_noncommuting(monomial.values[i], fields)){
+      tmp=monomial.values[i];
+      for(j=i;j>post_nc;j--){
+	monomial.values[j]=monomial.values[j-1];
+      }
+      monomial.values[post_nc]=tmp;
+      post_nc++;
+    }
+  }
+
+  monomial_sort_nonc(monomial, post_nc, monomial.length-1, fields, sign);
+
+  return(0);
+}
+// without noncommuting terms
+int monomial_sort_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign){
   int i;
   int index;
   // the pivot: middle of the monomial
@@ -812,8 +835,8 @@ int monomial_sort(Int_Array monomial, int begin, int end, Fields_Table fields, i
     exchange_monomial_terms(monomial, index, end, fields, sign);
 
     // recurse
-    monomial_sort(monomial, begin, index-1, fields, sign);
-    monomial_sort(monomial, index+1, end, fields, sign);
+    monomial_sort_nonc(monomial, begin, index-1, fields, sign);
+    monomial_sort_nonc(monomial, index+1, end, fields, sign);
   }
   return(0);
 }
@@ -872,7 +895,30 @@ int exchange_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table
 
 
 // sort a monomial by putting each group together
-int monomial_sort_groups(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign){
+// if the monomial contains noncommuting elements, put them at the beginning of the monomial
+int monomial_sort_groups(Int_Array monomial, Fields_Table fields, Groups groups, int* sign){
+  int i,j;
+  int tmp;
+  // first index after noncommuting indices
+  int post_nc=0;
+
+  for(i=0;i<monomial.length;i++){
+    if(is_noncommuting(monomial.values[i], fields)){
+      tmp=monomial.values[i];
+      for(j=post_nc;j<i;j++){
+	monomial.values[j+1]=monomial.values[j];
+      }
+      monomial.values[post_nc]=tmp;
+      post_nc++;
+    }
+  }
+
+  monomial_sort_groups_nonc(monomial, post_nc, monomial.length-1, fields, groups, sign);
+
+  return(0);
+}
+// without noncommuting terms
+int monomial_sort_groups_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign){
   int i;
   int index;
   // the pivot: middle of the monomial
@@ -897,8 +943,8 @@ int monomial_sort_groups(Int_Array monomial, int begin, int end, Fields_Table fi
     exchange_monomial_terms(monomial, index, end, fields, sign);
 
     // recurse
-    monomial_sort(monomial, begin, index-1, fields, sign);
-    monomial_sort(monomial, index+1, end, fields, sign);
+    monomial_sort_groups_nonc(monomial, begin, index-1, fields, groups, sign);
+    monomial_sort_groups_nonc(monomial, index+1, end, fields, groups, sign);
   }
   return(0);
 }