Support for non-commuting fieldsv1.3

10 files changed, 266 insertions, 119 deletions

diff --git a/src/definitions.cpp b/src/definitions.cpp
index 982f7a7..1884488 100644
--- a/src/definitions.cpp
+++ b/src/definitions.cpp
@@ -17,7 +17,7 @@ limitations under the License.
 #ifndef DEFINITIONS_GCC
 #define DEFINITIONS_GCC
 
-#define VERSION "1.2.1"
+#define VERSION "1.3"
 
 // number of entries in a configuration file
 #define ARG_COUNT 10
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 --------------------
 
diff --git a/src/fields.h b/src/fields.h
index 3795d92..b6f05c2 100644
--- a/src/fields.h
+++ b/src/fields.h
@@ -29,6 +29,8 @@ int free_Fields_Table(Fields_Table fields);
 int field_type(int index, Fields_Table fields);
 // check whether a field anticommutes
 int is_fermion(int index, Fields_Table fields);
+// check whether a field is non-commuting
+int is_noncommuting(int index, Fields_Table fields);
 
 
 // init
@@ -51,6 +53,8 @@ int identities_concat(Identities input, Identities* output);
 
 // resolve the identities
 int resolve_ids(Polynomial* polynomial, Fields_Table fields);
+// check whether an array is a sub-array of another, support for noncommuting elements
+int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fields_Table fields);
 
 
 // init
diff --git a/src/kondo.c b/src/kondo.c
index 39d74bb..79b7c56 100644
--- a/src/kondo.c
+++ b/src/kondo.c
@@ -29,16 +29,17 @@ limitations under the License.
 #include "definitions.cpp"
 #include "rational.h"
 
-// dimension of A, B and h (must be <10)
+// dimension of A, B, h and t (must be <10)
 #define KONDO_DIM 3
 // number of spin components
 #define KONDO_SPIN 2
 
-// offsets for indices of A, B and h
+// offsets for indices of A, B, h and t
 // order matters for symbols table
 #define KONDO_A_OFFSET 1
 #define KONDO_B_OFFSET 2
 #define KONDO_H_OFFSET 3
+#define KONDO_T_OFFSET 4
 
 // parsing modes (from parse_file.c)
 #define PP_NULL_MODE 0
@@ -193,13 +194,19 @@ int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fiel
 
   //  parameters
   char_array_append_str("h:",str_fields);
+  // h
   for(i=0;i<KONDO_DIM;i++){
-    char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_H_OFFSET));
+    char_array_snprintf(str_fields, "%d,", 10*(i+10*KONDO_H_OFFSET));
+  }
+  // t
+  for(i=0;i<KONDO_DIM;i++){
+    char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_T_OFFSET));
     if(i<KONDO_DIM-1){
-      char_array_append(',',str_fields);
+      char_array_append(',', str_fields);
     }
   }
-  char_array_append('\n',str_fields);
+  char_array_append('\n', str_fields);
+
 
   // declare Fermions
   char_array_append_str("f:",str_fields);
@@ -226,6 +233,16 @@ int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fiel
   }
   char_array_append('\n',str_fields);
 
+  // declare noncommuting
+  char_array_append_str("a:",str_fields);
+  for(i=0;i<KONDO_DIM;i++){
+    char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_T_OFFSET));
+    if(i<KONDO_DIM-1){
+      char_array_append(',',str_fields);
+    }
+  }
+  char_array_append('\n', str_fields);
+
   // parse fields table
   parse_input_fields(*str_fields, fields);
 
@@ -245,7 +262,7 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){
 
   // loop over box index
   for(i=1;i<=box_count;i++){
-    // loop over letters
+    // loop over letters (A and B)
     for(j=0;j<2;j++){
       // loop over space dimension
       for(k=0;k<KONDO_DIM;k++){
@@ -255,7 +272,7 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){
 	init_Char_Array(&tmp_str, 6);
 	char_array_snprintf(&tmp_str, "%c%d%d", letters[j], k, i);
 	// compute corresponding polynomial
-	kondo_resolve_ABh(tmp_str.str, &poly, *fields);
+	kondo_resolve_ABht(tmp_str.str, &poly, *fields);
 	free_Char_Array(tmp_str);
 	// write to output
 	polynomial_sprint(poly, str_symbols);
@@ -320,45 +337,6 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){
   return(0);
 }
 
-// generate Kondo symbols (older method: one symbol for each scalar product)
-int kondo_symbols_scalarprod(Char_Array* str_symbols, int box_count, Fields_Table* fields){
-  int i,j,k;
-  Char_Array tmp_str;
-  Polynomial poly;
-  char letters[3]={'A','B','h'};
-
-  init_Char_Array(str_symbols, STR_SIZE);
-  char_array_snprintf(str_symbols, "#!symbols\n");
-
-  // loop over box index
-  for(i=1;i<=box_count;i++){
-    // loop over letters
-    for(j=0;j<3;j++){
-      for(k=0;k<3;k++){
-	// write index
-	char_array_snprintf(str_symbols, "%d=", 100*(10*(KONDO_A_OFFSET+j)+KONDO_A_OFFSET+k)+i);
-	// write the name of the scalar product
-	init_Char_Array(&tmp_str, 6);
-	char_array_snprintf(&tmp_str, "%c%d.%c%d", letters[j], i, letters[k], i);
-	// compute corresponding polynomial
-	kondo_resolve_scalar_prod(tmp_str.str, &poly, *fields);
-	free_Char_Array(tmp_str);
-	// write to output
-	polynomial_sprint(poly, str_symbols);
-	free_Polynomial(poly);
-	// add ,
-	if(i<box_count || j<2 || k<2){
-	  char_array_snprintf(str_symbols,",\n");
-	}
-      }
-    }
-  }
-
-  parse_input_symbols(*str_symbols, fields);
-
-  return(0);
-}
-
 
 // generate Kondo groups (groups of independent variables)
 int kondo_groups(Char_Array* str_groups, int box_count){
@@ -395,12 +373,21 @@ int kondo_groups(Char_Array* str_groups, int box_count){
 
 // generate Kondo identities
 // h_3^2=1-h_1^2-h_2^2
+// and Pauli matrices
 int kondo_identities(Char_Array* str_identities){
   int i;
 
   init_Char_Array(str_identities,STR_SIZE);
   char_array_snprintf(str_identities, "#!identities\n");
 
+  // Pauli matrices
+  for(i=0;i<KONDO_DIM;i++){
+    char_array_snprintf(str_identities,"[f%d][f%d]=(1),\n",10*(10*KONDO_T_OFFSET+i),10*(10*KONDO_T_OFFSET+i));
+    char_array_snprintf(str_identities,"[f%d][f%d]=(s{-1})[f%d],\n",10*(10*KONDO_T_OFFSET+i),10*(10*KONDO_T_OFFSET+(i+1)%3),10*(10*KONDO_T_OFFSET+(i+2)%3));
+    char_array_snprintf(str_identities,"[f%d][f%d]=((-1)s{-1})[f%d],\n",10*(10*KONDO_T_OFFSET+(i+2)%3),10*(10*KONDO_T_OFFSET+(i+1)%3),10*(10*KONDO_T_OFFSET+i));
+  }
+
+  // h
   char_array_snprintf(str_identities, "[f%d][f%d]=(1)",10*(KONDO_DIM-1+10*KONDO_H_OFFSET),10*(KONDO_DIM-1+10*KONDO_H_OFFSET));
   for(i=0;i<KONDO_DIM-1;i++){
     char_array_snprintf(str_identities, "+(-1)[f%d][f%d]",10*(i+10*KONDO_H_OFFSET),10*(i+10*KONDO_H_OFFSET));
@@ -809,11 +796,10 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_
 	  if(tmp_poly.length>0){
 	    for(i=0;i<tmp_poly.length;i++){
 	      if(mode==PP_FIELD_SCALAR_MODE){
-		if(offset1!=KONDO_H_OFFSET || offset2!=KONDO_H_OFFSET){
-		  int_array_append(1000*(10*offset1+offset2)+index, tmp_poly.monomials+i);
-		}
+		int_array_append(1000*(10*offset1+offset2)+index, tmp_poly.monomials+i);
 	      }
 	      else{
+		// vector product
 		int_array_append(100*(100*KONDO_A_OFFSET+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+index, tmp_poly.monomials+i);
 	      }
 	    }
@@ -822,28 +808,15 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_
 	  else{
 	    init_Int_Array(&tmp_monomial, MONOMIAL_SIZE);
 	    if(mode==PP_FIELD_SCALAR_MODE){
-	      if(offset1!=KONDO_H_OFFSET || offset2!=KONDO_H_OFFSET){
-		int_array_append(1000*(10*offset1+offset2)+index, &tmp_monomial);
-	      }
+	      int_array_append(1000*(10*offset1+offset2)+index, &tmp_monomial);
 	    }
 	    else{
+	      // vector product
 	      int_array_append(100*(100*KONDO_A_OFFSET+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+index, &tmp_monomial);
 	    }
 	    init_Int_Array(&dummy_factor, 1);
 	    polynomial_append_noinit(tmp_monomial, dummy_factor, number_one(), &tmp_poly);
 	  }
-	  /* // older method in which a scalar product was expanded in A, B and h
-	  // resolve scalar product
-	  kondo_resolve_scalar_prod_symbols(buffer, &scalar_prod_poly);
-	  // add to tmp_poly
-	  if(tmp_poly.length==0){
-	    polynomial_concat(scalar_prod_poly,&tmp_poly);
-	  }
-	  else{
-	    polynomial_prod_chain(scalar_prod_poly,&tmp_poly,fields);
-	  }
-	  free_Polynomial(scalar_prod_poly);
-	  */
 	}
 	// switch back to null mode
 	mode=PP_NULL_MODE;
@@ -941,10 +914,10 @@ int parse_kondo_polynomial(Char_Array kondo_polynomial_str, Polynomial* polynomi
 }
 
 
-// read Aij, Bij, hi where i is a space dimension and j is a box index
-int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){
+// read Aij, Bij, hi, ti where i is a space dimension and j is a box index
+int kondo_resolve_ABht(char* str, Polynomial* output, Fields_Table fields){
   char* ptr;
-  // offset (A,B or H)
+  // offset (A,B, H or T)
   int offset=-1;
   // dimension
   int dim=-1;
@@ -978,6 +951,9 @@ int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){
     case 'h':
       offset=KONDO_H_OFFSET;
       break;
+    case 't':
+      offset=KONDO_T_OFFSET;
+      break;
     default:
       // set index if dim was already set
       if(dim>=0){
@@ -1001,8 +977,8 @@ int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){
     }
   }
 
-  // h's
-  if(offset==KONDO_H_OFFSET){
+  // h's and t's
+  if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){
     // external field
     init_Int_Array(&monomial,1);
     init_Int_Array(&factor,1);
@@ -1062,7 +1038,7 @@ int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){
 // read a Kondo scalar product (generalized to vector products as well)
 int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields){
   char* ptr;
-  // offset of each term (A,B or H)
+  // offset of each term (A,B,H or T)
   int offset=-1;
   // index of each term (0,...,box_count)
   int index=0;
@@ -1091,6 +1067,9 @@ int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields
     case 'h':
       offset=KONDO_H_OFFSET;
       break;
+    case 't':
+      offset=KONDO_T_OFFSET;
+      break;
 
     // scalar product
     case '.':
@@ -1195,8 +1174,8 @@ int kondo_polynomial_vector(int offset, int index, Polynomial (*polys)[3], Field
     init_Polynomial((*polys)+i,POLY_SIZE);
   }
 
-  // h's
-  if(offset==KONDO_H_OFFSET){
+  // h's and t's
+  if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){
     // construct every component field
     for(i=0;i<KONDO_DIM;i++){
       // external field
@@ -1261,7 +1240,7 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output){
   char* ptr;
   // first or second term
   int term=0;
-  // offset of each term (A,B or H)
+  // offset of each term (A,B,H or T)
   int offset[2];
   // index of each term (0,...,box_count)
   int index[2]={0,0};
@@ -1289,6 +1268,9 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output){
     case 'h':
       offset[term]=KONDO_H_OFFSET;
       break;
+    case 't':
+      offset[term]=KONDO_T_OFFSET;
+      break;
     // switch term
     case '.':
       term=1-term;
@@ -1304,13 +1286,13 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output){
     init_Int_Array(&monomial,2);
     init_Int_Array(&factor, 1);
 
-    if(offset[0]==KONDO_H_OFFSET){
+    if(offset[0]==KONDO_H_OFFSET || offset[0]==KONDO_T_OFFSET){
       int_array_append(10*(10*offset[0]+i)+index[0], &monomial);
     }
     else{
       int_array_append(100*(10*offset[0]+i)+index[0], &monomial);
     }
-    if(offset[1]==KONDO_H_OFFSET){
+    if(offset[1]==KONDO_H_OFFSET || offset[1]==KONDO_T_OFFSET){
       int_array_append(10*(10*offset[1]+i)+index[1], &monomial);
     }
     else{
@@ -1344,6 +1326,9 @@ int get_offset_index(char* str, int* offset, int* index){
     case 'h':
       *offset=KONDO_H_OFFSET;
       break;
+    case 't':
+      *offset=KONDO_T_OFFSET;
+      break;
     default:
       // char to int
       *index=*ptr-'0';
@@ -1378,6 +1363,9 @@ int get_offsets_index(char* str, int* offset1, int* offset2, int* index){
     case 'h':
       offset[term]=KONDO_H_OFFSET;
       break;
+    case 't':
+      offset[term]=KONDO_T_OFFSET;
+      break;
     // switch term
     case '.':
       term=1-term;
@@ -1391,6 +1379,11 @@ int get_offsets_index(char* str, int* offset1, int* offset2, int* index){
   *offset1=offset[0];
   *offset2=offset[1];
 
+  // if no A's or B's, then index=0
+  if((offset[0]==KONDO_H_OFFSET || offset[0]==KONDO_T_OFFSET) && (offset[1]==KONDO_H_OFFSET || offset[1]==KONDO_T_OFFSET)){
+    *index=0;
+  }
+
   return(0);
 }
 
@@ -1429,6 +1422,9 @@ int get_symbol_index(char* str){
     case 'h':
       offset=KONDO_H_OFFSET;
       break;
+    case 't':
+      offset=KONDO_T_OFFSET;
+      break;
     default:
       // set index if dim was already set
       if(dim>=0){
@@ -1443,8 +1439,8 @@ int get_symbol_index(char* str){
   if(offset==-1){
     return(-1);
   }
-  // no symbol for h
-  if(offset==KONDO_H_OFFSET){
+  // no symbol for h or t
+  if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){
     return(10*(10*offset+dim));
   }
   else{
diff --git a/src/kondo.h b/src/kondo.h
index 23756ad..c534145 100644
--- a/src/kondo.h
+++ b/src/kondo.h
@@ -55,7 +55,7 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_
 int parse_kondo_polynomial(Char_Array kondo_polynomial_str, Polynomial* polynomial, Fields_Table fields);
 
 // read Aij, Bij, hi where i is a space dimension and j is a box index
-int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields);
+int kondo_resolve_ABht(char* str, Polynomial* output, Fields_Table fields);
 // read a Kondo scalar product
 int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields);
 // compute a scalar product of polynomial vectors
diff --git a/src/mean.c b/src/mean.c
index aad7a5a..39ece56 100644
--- a/src/mean.c
+++ b/src/mean.c
@@ -42,7 +42,7 @@ int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Ma
   *out=polynomial_one();
 
   // sort first
-  monomial_sort(monomial, 0, monomial.length-1, fields, &sign);
+  monomial_sort(monomial, fields, &sign);
   polynomial_multiply_Qscalar(*out, quot(sign,1));
   // get internals
   // (*out).monomials is the first element of out but it only has 1 element
@@ -417,7 +417,7 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po
     if(check_monomial_match(tmp_monomial, fields)==1){
       // sort monomial
       sign=1;
-      monomial_sort(tmp_monomial, 0, tmp_monomial.length-1, fields, &sign);
+      monomial_sort(tmp_monomial, fields, &sign);
       number_Qprod_chain(quot(sign,1), &tmp_num);
 
       // mean
@@ -628,7 +628,7 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol
   init_Polynomial(output, MONOMIAL_SIZE);
 
   // check whether there are symbols
-  // requires the symbols to be at the end of the monomial
+  // IMPORTANT: the symbols must be at the end of the monomial
   if(monomial.length==0 || field_type(monomial.values[monomial.length-1], fields)!=FIELD_SYMBOL){
     // mean
     mean(monomial, &num_mean, fields, propagator);
@@ -639,7 +639,7 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol
     // sort into groups
     if(groups.length>0){
       sign=1;
-      monomial_sort_groups(monomial, 0, monomial.length-1, fields, groups, &sign);
+      monomial_sort_groups(monomial, fields, groups, &sign);
     }
     // construct groups and take mean
     init_Int_Array(&tmp_monomial, MONOMIAL_SIZE);
diff --git a/src/parse_file.c b/src/parse_file.c
index 6054372..19b91c6 100644
--- a/src/parse_file.c
+++ b/src/parse_file.c
@@ -44,16 +44,17 @@ limitations under the License.
 #define PP_EXTERNAL_MODE 8
 #define PP_INTERNAL_MODE 9
 #define PP_FERMIONS_MODE 10
+#define PP_NONCOMMUTING_MODE 11
 // indices
-#define PP_INDEX_MODE 11
+#define PP_INDEX_MODE 12
 // factors or monomials
-#define PP_BRACKET_MODE 12
+#define PP_BRACKET_MODE 13
 // labels
-#define PP_LABEL_MODE 13
+#define PP_LABEL_MODE 14
 // polynomial
-#define PP_POLYNOMIAL_MODE 14
+#define PP_POLYNOMIAL_MODE 15
 // group
-#define PP_GROUP_MODE 15
+#define PP_GROUP_MODE 16
 
 
 // parse fields list
@@ -101,6 +102,11 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
 	  mode=PP_FERMIONS_MODE;
 	}
 	break;
+      case 'a':
+	if(mode==PP_NULL_MODE){
+	  mode=PP_NONCOMMUTING_MODE;
+	}
+	break;
 
       // reset buffer
       case ':':
@@ -123,6 +129,9 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
 	else if(mode==PP_FERMIONS_MODE){
 	  int_array_append(i,&((*fields).fermions));
 	}
+	else if(mode==PP_NONCOMMUTING_MODE){
+	  int_array_append(i,&((*fields).noncommuting));
+	}
 	buffer_ptr=buffer;
 	*buffer_ptr='\0';
 	break;
@@ -142,6 +151,9 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
 	else if(mode==PP_FERMIONS_MODE){
 	  int_array_append(i,&((*fields).fermions));
 	}
+	else if(mode==PP_NONCOMMUTING_MODE){
+	  int_array_append(i,&((*fields).noncommuting));
+	}
 	mode=PP_NULL_MODE;
 	break;
 
@@ -417,7 +429,7 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){
   (*fields).ids.length=0;
   for(i=0;i<tmp;i++){
     sign=1;
-    monomial_sort((*fields).ids.lhs[i], 0, (*fields).ids.lhs[i].length-1, *fields, &sign);
+    monomial_sort((*fields).ids.lhs[i], *fields, &sign);
     polynomial_simplify((*fields).ids.rhs+i, *fields);
     polynomial_multiply_Qscalar((*fields).ids.rhs[i],quot(sign,1));
   }
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);
 }
diff --git a/src/polynomial.h b/src/polynomial.h
index ec50227..6aff76d 100644
--- a/src/polynomial.h
+++ b/src/polynomial.h
@@ -99,14 +99,18 @@ int polynomial_sort(Polynomial* polynomial, int begin, int end);
 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);
+int monomial_sort(Int_Array monomial, Fields_Table fields, int* sign);
+// without noncommuting terms
+int monomial_sort_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign);
 // order fields: parameter, external, internal
 int compare_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields);
 // exchange two fields (with sign)
 int exchange_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields, int* sign);
 
 // 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);
+int monomial_sort_groups(Int_Array monomial, Fields_Table fields, Groups groups, int* sign);
+// without noncommuting terms
+int monomial_sort_groups_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign);
 // order fields: group, then parameter, external, internal
 int compare_monomial_terms_groups(Int_Array monomial, int pos1, int pos2, Fields_Table fields, Groups groups);
 
diff --git a/src/types.h b/src/types.h
index 0452ebc..6b029d5 100644
--- a/src/types.h
+++ b/src/types.h
@@ -154,6 +154,8 @@ typedef struct Fields_Table{
   Symbols symbols;
   // list of anti-commuting variables (fields or symbols)
   Int_Array fermions;
+  // list of non-commuting variables (fields or symbols)
+  Int_Array noncommuting;
 } Fields_Table;
 
 // index labels