Support for non-commuting fieldsv1.3

1 files changed, 104 insertions, 21 deletions

diff --git a/src/fields.c b/src/fields.c
index 1b221f2..bfd1f39 100644
--- a/src/fields.c
+++ b/src/fields.c
@@ -32,6 +32,7 @@ int init_Fields_Table(Fields_Table* fields){
   init_Identities(&((*fields).ids), FIELDS_SIZE);
   init_Symbols(&((*fields).symbols), FIELDS_SIZE);
   init_Int_Array(&((*fields).fermions),FIELDS_SIZE);
+  init_Int_Array(&((*fields).noncommuting),FIELDS_SIZE);
   return(0);
 }
 int free_Fields_Table(Fields_Table fields){
@@ -41,6 +42,7 @@ int free_Fields_Table(Fields_Table fields){
   free_Identities(fields.ids);
   free_Symbols(fields.symbols);
   free_Int_Array(fields.fermions);
+  free_Int_Array(fields.noncommuting);
   return(0);
 }
 
@@ -73,6 +75,16 @@ int is_fermion(int index, Fields_Table fields){
   }
 }
 
+// check whether a field is non-commuting
+int is_noncommuting(int index, Fields_Table fields){
+  if(int_array_find(abs(index), fields.noncommuting)>=0){
+    return(1);
+  }
+  else{
+    return(0);
+  }
+}
+
 
 // ------------------ Identities --------------------
 
@@ -180,13 +192,16 @@ int identities_concat(Identities input, Identities* output){
 
 // resolve the identities
 // requires both the monomials in polynomial and the ids in fields to be sorted
+// IMPORTANT: the sorting must be such that noncommuting fields must come before the other fields
 int resolve_ids(Polynomial* polynomial, Fields_Table fields){
   int i,j,k,l;
   int sign;
   int fermion_count;
+  int first_field;
   int at_least_one;
   int security;
-  Int_Array monomial;
+  Int_Array pre_monomial;
+  Int_Array post_monomial;
   Number num;
   Number tmp_num;
 
@@ -207,29 +222,38 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields){
 
       // loop over ids
       for(j=0;j<fields.ids.length;j++){
+
 	// check whether the monomial matches the id
-	if(int_array_is_subarray_ordered(fields.ids.lhs[j],(*polynomial).monomials[i])==1){
-	  init_Int_Array(&monomial, (*polynomial).monomials[i].length);
-	  
-	  // remove lhs from monomial
-	  // sign from moving the fields out of the monomial
+	first_field=int_array_is_subarray_noncommuting(fields.ids.lhs[j],(*polynomial).monomials[i],fields);
+	if(first_field>=0){
+	  init_Int_Array(&pre_monomial, (*polynomial).monomials[i].length);
+	  init_Int_Array(&post_monomial, (*polynomial).monomials[i].length);
+
+	  // add whatever is before the first field to pre
+	  for(k=0;k<first_field;k++){
+	    int_array_append((*polynomial).monomials[i].values[k],&pre_monomial);
+	  }
+
+	  // find others and move them together
+	  // sign from moving the fields
 	  sign=1;
-	  // number of Fermions to remove from the monomial
+	  // number of Fermions to jump over
 	  fermion_count=0;
-	  for(k=0,l=0;k<(*polynomial).monomials[i].length;k++){
+	  for(l=1,k=first_field+1;k<(*polynomial).monomials[i].length;k++){
 	    // check whether the field is identical to the "current" one in the id
 	    // if l is too large, then keep the field
 	    if(l>=fields.ids.lhs[j].length || (*polynomial).monomials[i].values[k]!=fields.ids.lhs[j].values[l]){
-	      int_array_append((*polynomial).monomials[i].values[k],&monomial);
-	      // sign correction
-	      if(fermion_count % 2 ==1 && is_fermion((*polynomial).monomials[i].values[k], fields)){
-		sign*=-1;
+	      // add to post
+	      int_array_append((*polynomial).monomials[i].values[k],&post_monomial);
+	      // count Fermions to jump
+	      if(is_fermion((*polynomial).monomials[i].values[k],fields)){
+		fermion_count++;
 	      }
 	    }
 	    else{
-	      // increment fermion_count
-	      if(is_fermion(fields.ids.lhs[j].values[l],fields)){
-		fermion_count++;
+	      // sign correction
+	      if(is_fermion(fields.ids.lhs[j].values[l],fields) && fermion_count % 2 == 1){
+		sign*=-1;
 	      }
 	      // increment "current" field in the id
 	      l++;
@@ -240,30 +264,33 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields){
 	  // add extra monomials (if there are more than 1)
 	  for(k=1;k<fields.ids.rhs[j].length;k++){
 	    number_prod(num, fields.ids.rhs[j].nums[k], &tmp_num);
-	    polynomial_append(monomial, (*polynomial).factors[i], tmp_num, polynomial);
+	    polynomial_append(pre_monomial, (*polynomial).factors[i], tmp_num, polynomial);
 	    free_Number(tmp_num);
 	    int_array_concat(fields.ids.rhs[j].monomials[k],(*polynomial).monomials+(*polynomial).length-1);
+	    int_array_concat(post_monomial,(*polynomial).monomials+(*polynomial).length-1);
 	    // re-sort monomial
 	    sign=1;
-	    monomial_sort((*polynomial).monomials[(*polynomial).length-1],0,(*polynomial).monomials[(*polynomial).length-1].length-1,fields,&sign);
+	    monomial_sort((*polynomial).monomials[(*polynomial).length-1],fields,&sign);
 	    number_Qprod_chain(quot(sign,1),(*polynomial).nums+(*polynomial).length-1);
 	  } 
 	  // correct i-th monomial
 	  free_Number((*polynomial).nums[i]);
 	  (*polynomial).nums[i]=number_prod_ret(num,fields.ids.rhs[j].nums[0]);
 	  free_Int_Array((*polynomial).monomials[i]);
-	  (*polynomial).monomials[i]=monomial;
+	  (*polynomial).monomials[i]=pre_monomial;
 	  int_array_concat(fields.ids.rhs[j].monomials[0],(*polynomial).monomials+i);
+	  int_array_concat(post_monomial,(*polynomial).monomials+i);
+	  free_Int_Array(post_monomial);
 	  // re-sort monomial
 	  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);
 
 	  // free num
 	  free_Number(num);
 
-	  // repeat the step (in order to perform all of the replacements if several are necessary)
-	  j--;
+	  // repeat the replacement (in order to perform all of the replacements if several are necessary)
+	  j=0;
 	  if(at_least_one==0){
 	    at_least_one=1;
 	  }
@@ -275,6 +302,62 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields){
   return(0);
 }
 
+// check whether an array is a sub-array of another
+// requires noncommuting elements to be next to each other
+// other elements may be separated, but the order must be respected
+// returns the first index of the sub-array
+// IMPORTANT: the noncommuting elements must precede all others in input and in test_array
+int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fields_Table fields){
+  int i,j;
+  int matches=0;
+  int post_nc=0;
+  int match_nc;
+  int first=-1;
+
+  // bound noncommuting elements
+  while(is_noncommuting(input.values[post_nc], fields)==1){
+    post_nc++;
+  }
+
+  for(i=0,match_nc=0;i<test_array.length;i++){
+    if(test_array.values[i]==input.values[0]){
+      match_nc=1;
+    }
+    for(j=1;j<post_nc;j++){
+      if(test_array.values[i+j]!=input.values[j]){
+	match_nc=0;
+      }
+    }
+    if(match_nc==1){
+      first=i;
+      break;
+    }
+  }
+
+  if(first<0){
+    return(-1);
+  }
+
+  if(post_nc>0){
+    matches=post_nc;
+  }
+  else{
+    matches=1;
+  }
+
+  for(i=first+1;i<test_array.length && matches<input.length;i++){
+    if(input.values[matches]==test_array.values[i]){
+      matches++;
+    }
+  }
+  if(matches==input.length){
+    return(first);
+  }
+  else{
+    return(-1);
+  }
+}
+
 
 // ------------------ Symbols --------------------