Update to v1.5.HEADmaster

  The update to version 1.5 is rather substantial, and introduces some minor
  backward-incompatibilities:
    * The header "#!symbols" has been replaced by "#!virtual_fields"
    * Multiplying polynomials using the '*' symbol is no longer supported (or,
      rather, the symbolic capabilities of meankondo were enhanced, and the
      syntax has been changed).
    * 'meantools exp' has been removed (its functionality is now handled by
      other means)
    * 'meantoolds derive' has been replaced by 'meantools differentiate'

  * The symbolic capabilities were enhanced: polynomials can now be
    multiplied, added, exponentiated, and their logarithms can be taken
    directly in the configuration file.

  * The flow equation can now be processed after being computed using the
    various "#!postprocess_*" entries.

  * Deprecated kondo_preprocess.

  * Compute the mean using an LU decomposition if possible.

  * More detailed checks for syntax errors in configuration file.

  * Check that different '#!group' entries are indeed uncorrelated.

  * New flags in meankondo: '-p' and '-A'.

  * New tool: meantools expand.

  * Improve conversion to LaTeX using meantools-convert

  * Assign terms randomly to different threads.

  * Created vim files to implement syntax highlighting for configuration
    files.

  * Multiple bug fixes

60 files changed, 3070 insertions, 682 deletions

diff --git a/src/array.c b/src/array.c
index 11d14f7..7fb6bd8 100644
--- a/src/array.c
+++ b/src/array.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -72,6 +72,14 @@ int int_array_append(int val, Int_Array* output){
   return(0);
 }
 
+// add a value only if it is not already present
+int int_array_append_unique(int val, Int_Array* output){
+  if(int_array_find(val,*output)<0){
+    int_array_append(val,output);
+  }
+  return(0);
+}
+
 // concatenate
 int int_array_concat(Int_Array input, Int_Array* output){
   int i;
@@ -87,6 +95,15 @@ int int_array_concat(Int_Array input, Int_Array* output){
   return(0);
 }
 
+// concat but only add values that are not already present in the array
+int int_array_concat_unique(Int_Array input, Int_Array* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    int_array_append_unique(input.values[i],output);
+  }
+  return(0);
+}
+
 // find (does not assume the array is sorted)
 int int_array_find(int val, Int_Array array){
   int i;
@@ -199,7 +216,12 @@ int int_array_print(Int_Array array){
   for(i=0;i<array.length-1;i++){
     printf("%d,",array.values[i]);
   }
-  printf("%d)",array.values[array.length-1]);
+  if(array.length>0){
+    printf("%d)",array.values[array.length-1]);
+  }
+  else{
+    printf(")");
+  }
   return(0);
 }
 
@@ -334,6 +356,19 @@ int char_array_concat(Char_Array input, Char_Array* output){
   return(0);
 }
 
+// substring
+int char_array_substring(Char_Array str, int begin, int end, Char_Array* substr){
+  int i;
+  if(begin>end || begin<0 || end>=str.length){
+    fprintf(stderr,"error: cannot extract a substring [%d,%d] from a string of length %d\n", begin, end, str.length);
+    exit(-1);
+  }
+  init_Char_Array(substr,end-begin);
+  for(i=begin;i<=end;i++){
+    char_array_append(str.str[i],substr);
+  }
+  return(0);
+}
 
 
 // convert to char*
@@ -343,7 +378,7 @@ int char_array_to_str(Char_Array input, char** output){
   for(i=0;i<input.length;i++){
     (*output)[i]=input.str[i];
   }
-  if((*output)[input.length-1]!='\0'){
+  if(input.length==0 ||  (*output)[input.length-1]!='\0'){
     (*output)[input.length]='\0';
   }
   return(0);
@@ -371,6 +406,34 @@ int str_to_char_array(char* str, Char_Array* output){
   return(0);
 }
 
+// compare char_array's
+int char_array_cmp(Char_Array char_array1, Char_Array char_array2){
+  int j;
+  if(char_array1.length!=char_array2.length){
+    return(0);
+  }
+  for(j=0;j<char_array1.length && j<char_array2.length;j++){
+    if(char_array1.str[j]!=char_array2.str[j]){
+      return(0);
+    }
+  }
+  return(1);
+}
+
+// compare a char_array and a char*
+int char_array_cmp_str(Char_Array char_array, char* str){
+  int j;
+  for(j=0;j<char_array.length && str[j]!='\0';j++){
+    if(char_array.str[j]!=str[j]){
+      return(0);
+    }
+  }
+  if(j==char_array.length && str[j]=='\0'){
+    return(1);
+  }
+  return(0);
+}
+
 
 // format strings
 int char_array_snprintf(Char_Array* output, char* fmt, ...){
diff --git a/src/array.h b/src/array.h
index fb74e67..480b589 100644
--- a/src/array.h
+++ b/src/array.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -34,8 +34,12 @@ int int_array_resize(Int_Array* array, int newsize);
 
 // add a value
 int int_array_append(int val, Int_Array* output);
+// add a value only if it is not already present
+int int_array_append_unique(int val, Int_Array* output);
 // concatenate
 int int_array_concat(Int_Array input, Int_Array* output);
+// concat but only add values that are not already present in the array
+int int_array_concat_unique(Int_Array input, Int_Array* output);
 
 // find (does not assume the array is sorted)
 int int_array_find(int val, Int_Array array);
@@ -75,6 +79,9 @@ int char_array_append_str(char* str, Char_Array* output);
 // concatenate
 int char_array_concat(Char_Array input, Char_Array* output);
 
+// substring
+int char_array_substring(Char_Array str, int begin, int end, Char_Array* substr);
+
 // convert to char*
 int char_array_to_str(Char_Array input, char** output);
 // noinit (changes the size of input if needed)
@@ -82,6 +89,11 @@ char* char_array_to_str_noinit(Char_Array* input);
 // convert from char*
 int str_to_char_array(char* str, Char_Array* output);
 
+// compare char_array's
+int char_array_cmp(Char_Array char_array1, Char_Array char_array2);
+// compare a char_array and a char*
+int char_array_cmp_str(Char_Array char_array, char* str);
+
 // format strings
 int char_array_snprintf(Char_Array* output, char* fmt, ...);
 
diff --git a/src/cli_parser.c b/src/cli_parser.c
index bdca700..6012e8a 100644
--- a/src/cli_parser.c
+++ b/src/cli_parser.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -121,3 +121,23 @@ int find_str_arg(char* title, Str_Array str_args){
   }
 }
 
+// find a string argument with the specified title
+// namespace support
+int find_str_arg_ns(char* title, Char_Array namespace, Str_Array str_args){
+  Char_Array buffer;
+  int ret;
+
+  // append namespace to title
+  char_array_cpy(namespace,&buffer);
+  char_array_append(':',&buffer);
+  char_array_append_str(title,&buffer);
+
+  // check whether the namespace entry exists
+  ret=find_str_arg(char_array_to_str_noinit(&buffer), str_args);
+  free_Char_Array(buffer);
+  // if not, use global entry
+  if(ret==-1){
+    return(find_str_arg(title, str_args));
+  }
+  return(ret);
+}
diff --git a/src/cli_parser.h b/src/cli_parser.h
index b050e37..b66608c 100644
--- a/src/cli_parser.h
+++ b/src/cli_parser.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -30,6 +30,8 @@ int read_config_file(Str_Array* str_args, const char* file, int read_from_stdin)
 int get_str_arg_title(Char_Array str_arg, Char_Array* out);
 // find a string argument with the specified title
 int find_str_arg(char* title, Str_Array str_args);
+// namespace support
+int find_str_arg_ns(char* title, Char_Array namespace, Str_Array str_args);
 
 
 #endif
diff --git a/src/coefficient.c b/src/coefficient.c
index d4cb6eb..1efe8fd 100644
--- a/src/coefficient.c
+++ b/src/coefficient.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -202,6 +202,235 @@ int coefficient_simplify(Coefficient* coefficient){
   return(0);
 }
 
+// put all terms under a common denominator and simplify the resulting fraction
+int coefficient_simplify_rational(Coefficient constant, Coefficient* coefficient){
+  int ret;
+  Coefficient remainder;
+  Coefficient quotient;
+  Coefficient quotient_prev;
+  Coefficient out;
+  int power;
+  int max_power;
+
+  // common denominator
+  coefficient_common_denominator(constant, coefficient);
+
+  // init
+  init_Coefficient(&out, COEF_SIZE);
+
+  // simplify, one power at a time
+  // largest power (larger powers are at the end)
+  max_power=(*coefficient).denoms[(*coefficient).length-1].power;
+
+  quotient_prev=*coefficient;
+  for(power=max_power;power>=1;power--){
+    ret=coefficient_simplify_fraction(constant, quotient_prev, &remainder, &quotient);
+
+    // if fail to simplify, stop
+    if(ret<0){
+      if(power<max_power){
+	coefficient_concat_noinit(quotient_prev, &out);
+      }
+      else{
+	coefficient_concat(quotient_prev, &out);
+      }
+      break;
+    }
+
+    // add to output
+    coefficient_concat_noinit(remainder, &out);
+  }
+  // if the factorization always succeeded
+  if(max_power>=1 && power==0){
+    coefficient_concat_noinit(quotient, &out);
+  }
+
+  coefficient_simplify(&out);
+
+  // set coefficient to out
+  free_Coefficient(*coefficient);
+  *coefficient=out;
+
+  return 0;
+}
+
+// put all terms under a common denominator
+// only supports coefficients with only one constant
+int coefficient_common_denominator(Coefficient constant, Coefficient* coefficient){
+  int max_power;
+  int i,j;
+  Coefficient tmp;
+  Coefficient out;
+  Coefficient* C_n;
+
+  init_Coefficient(&out, COEF_SIZE);
+
+  // largest power (larger powers are at the end)
+  max_power=(*coefficient).denoms[(*coefficient).length-1].power;
+
+  // store powers of the constant
+  C_n=calloc(sizeof(Coefficient), max_power-1);
+  for(i=0;i<max_power-1;i++){
+    // start from previous product
+    if(i==0){
+      coefficient_cpy(constant, C_n+i);
+    }
+    else{
+      coefficient_cpy(C_n[i-1], C_n+i);
+    }
+    // multiply by constant
+    coefficient_prod_chain(constant, C_n+i);
+  }
+
+  // multiply each term
+  for (i=0;i<(*coefficient).length;i++){
+    init_Coefficient(&tmp, COEF_SIZE);
+    // start with numerator
+    coefficient_append_noinit((*coefficient).factors[i], (*coefficient).nums[i], (*coefficient).denoms[i], &tmp);
+    // multiply
+    if((*coefficient).denoms[i].power<max_power){
+      if((*coefficient).denoms[i].power==max_power-1){
+	coefficient_prod_chain(constant, &tmp);
+      }
+      else{
+	coefficient_prod_chain(C_n[max_power-(*coefficient).denoms[i].power-2], &tmp);
+      }
+    }
+
+    // set denom
+    for(j=0;j<tmp.length;j++){
+      tmp.denoms[j].power=max_power;
+    }
+
+    // add to out
+    coefficient_concat_noinit(tmp, &out);
+  }
+
+  // free C_n
+  for(i=0;i<max_power-1;i++){
+    free_Coefficient(C_n[i]);
+  }
+  free(C_n);
+
+  // free coefficient vectors
+  free((*coefficient).factors);
+  free((*coefficient).nums);
+  free((*coefficient).denoms);
+
+  // set output
+  *coefficient=out;
+
+  coefficient_simplify(coefficient);
+
+  return(0);
+}
+
+// simplify coefficient / constant
+// returns both the remainder and the quotient
+// assumes both coefficient and constant are ordered with the highest order terms last
+int coefficient_simplify_fraction(Coefficient constant, Coefficient coefficient, Coefficient* remainder, Coefficient* out){
+  Coefficient tmp;
+  int step_counter=0;
+  int max_order;
+  int i,j,k;
+  Int_Array rfactors;
+
+  if(constant.length==0){
+    // nothing to do
+    return 0;
+  }
+
+  coefficient_cpy(coefficient, remainder);
+  init_Coefficient(out, COEF_SIZE);
+
+  // continue until (*remainder) is of lower order than constant
+  while((*remainder).length>0 && (*remainder).factors[(*remainder).length-1].length>=constant.factors[constant.length-1].length){
+    step_counter++;
+
+    // interrupt if too long
+    if(step_counter>=coefficient.length*100){
+      free_Coefficient(*remainder);
+      free_Coefficient(*out);
+      return -1;
+    }
+
+    // try to find a term in the constant that divides the last term of the (*remainder)
+    rfactors=(*remainder).factors[(*remainder).length-1];
+
+    // highest order in constant
+    max_order=constant.factors[constant.length-1].length;
+
+    // start from one of the highest order term and look for a common factor
+    for(i=constant.length-1; i>=0; i--){
+      // fail: no highest order terms have been matched
+      if(constant.factors[i].length<max_order){
+	free_Coefficient(*remainder);
+	free_Coefficient(*out);
+	return -2;
+      }
+
+      // check whether the term can be a factor of the last term of the (*remainder)
+      if(int_array_is_subarray_ordered(constant.factors[i], rfactors)==1){
+	// extract the factors that are not in constant
+	init_Coefficient(&tmp, constant.length);
+
+	// init with one term
+	tmp.length=1;
+	init_Int_Array(tmp.factors,MONOMIAL_SIZE);
+
+	for(j=0,k=0;j<rfactors.length;j++){
+	  // check that index is not in constant
+	  if(k<constant.factors[i].length){
+	    if(rfactors.values[j]!=constant.factors[i].values[k]){
+	      int_array_append(rfactors.values[j],tmp.factors);
+	    }
+	    else{
+	      // move to next term in constant
+	      k++;
+	    }
+	  }
+	}
+
+	// numerical prefactor: term in the (*remainder) / term in the constant
+	number_quot((*remainder).nums[(*remainder).length-1], constant.nums[i], tmp.nums);
+	// denominator (dummy)
+	tmp.denoms[0]=(*remainder).denoms[(*remainder).length-1];
+
+	// add to out
+	coefficient_concat(tmp, out);
+
+	// multiply by -1
+	Q minus_1;
+	minus_1.numerator=-1;
+	minus_1.denominator=1;
+	number_Qprod_chain(minus_1, tmp.nums);
+
+	// multiply by constant
+	coefficient_prod_chain(constant, &tmp);
+
+	// add to remainder
+	coefficient_concat(tmp, remainder);
+
+	// free memory
+	free_Coefficient(tmp);
+
+	// simplify
+	coefficient_simplify(remainder);
+
+	break;
+      }
+    }
+  }
+
+  // success!
+  // decrease power of constant
+  for(i=0;i<(*out).length;i++){
+    (*out).denoms[i].power=(*out).denoms[i].power-1;
+  }
+
+  return(0);
+}
+
 // sort the terms in an equation (quicksort algorithm)
 int sort_coefficient(Coefficient* coefficient, int begin, int end){
   int i;
@@ -251,7 +480,7 @@ int exchange_coefficient_terms(int i, int j, Coefficient* coefficient){
   return(0);
 }
 
-// derive a coefficient with respect to an index
+// differentiate a coefficient with respect to an index
 int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){
   int i,j;
   // temp list of indices
@@ -325,7 +554,7 @@ int coefficient_deriv(Coefficient input, int index, Coefficient* output){
 }
 
 /*
-// derive a coefficient with respect to an index (as a polynomial) (does not derive the 1/(1+C)^p )
+// differentiate a coefficient with respect to an index (as a polynomial) (does not differentiate the 1/(1+C)^p )
 int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){
   int i;
   // temp list of indices
@@ -365,7 +594,7 @@ int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){
   return(0);
 }
 
-// derive a monomial with respect to an index
+// differentiate a monomial with respect to an index
 int monomial_deriv(Int_Array factor, int index, Int_Array* out_factor, int* match_count){
   int j;
 
@@ -414,14 +643,14 @@ int coefficient_prod(Coefficient coef1, Coefficient coef2, Coefficient* output){
       int_array_concat(coef2.factors[j],&factor);
       
       // don't throw an error if the power is 0
-      if(coef2.denoms[i].power==0){
-	coef2.denoms[i].index=coef1.denoms[i].index;
+      if(coef2.denoms[j].power==0){
+	coef2.denoms[j].index=coef1.denoms[i].index;
       }
       else if(coef1.denoms[i].power==0){
-	coef1.denoms[i].index=coef2.denoms[i].index;
+	coef1.denoms[i].index=coef2.denoms[j].index;
       }
       if(coef1.denoms[i].index!=coef2.denoms[j].index){
-	fprintf(stderr,"error: cannot multiply flow equations with different constants\n");
+	fprintf(stderr,"error: cannot multiply flow equations with different constants: got %d and %d\n", coef1.denoms[i].index, coef2.denoms[j].index);
 	exit(-1);
       }
       denom=coef1.denoms[i];
diff --git a/src/coefficient.h b/src/coefficient.h
index 43dd313..f8a813a 100644
--- a/src/coefficient.h
+++ b/src/coefficient.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -48,12 +48,20 @@ int coefficient_concat_noinit(Coefficient input, Coefficient* output);
 
 // simplify a Coefficient
 int coefficient_simplify(Coefficient* coefficient);
+
+// put all terms under a common denominator and simplify the resulting fraction
+int coefficient_simplify_rational(Coefficient constant, Coefficient* coefficient);
+// put all terms under a common denominator
+int coefficient_common_denominator(Coefficient constant, Coefficient* coefficient);
+// simplify coefficient / constant
+int coefficient_simplify_fraction(Coefficient constant, Coefficient coefficient, Coefficient* remainder, Coefficient* out);
+
 // sort the terms in an equation (quicksort algorithm)
 int sort_coefficient(Coefficient* coefficient, int begin, int end);
 // exchange two terms (for the sorting algorithm)
 int exchange_coefficient_terms(int i, int j, Coefficient* coefficient);
 
-// derive a coefficient with respect to an index (as a polynomial) (does not derive the 1/(1+C)^p )
+// differentiate a coefficient with respect to an index (as a polynomial) (does not differentiate the 1/(1+C)^p )
 int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output);
 int coefficient_deriv(Coefficient input, int index, Coefficient* output);
 
diff --git a/src/definitions.cpp b/src/definitions.cpp
index 30f021f..db09e57 100644
--- a/src/definitions.cpp
+++ b/src/definitions.cpp
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -33,10 +33,16 @@ limitations under the License.
 #define EQUATION_SIZE 20
 // number of fields
 #define FIELDS_SIZE 50
+// number of variables
+#define VARIABLES_SIZE 10
 // number of elements in numbers
 #define NUMBER_SIZE 5
 // number of elements in a group
 #define GROUP_SIZE 5
+// number of children per node in a symbol_tree
+#define SYMBOL_TREE_SIZE 2
+// size of character string in symbol tree
+#define SYMBOL_TREE_LABEL_SIZE 10
 
 
 // display options
@@ -55,6 +61,6 @@ limitations under the License.
 #define FIELD_PARAMETER 1
 #define FIELD_EXTERNAL 2
 #define FIELD_INTERNAL 3
-#define FIELD_SYMBOL 4
+#define FIELD_VIRTUAL 4
 
 #endif
diff --git a/src/determinant.c b/src/determinant.c
new file mode 100644
index 0000000..906c75f
--- /dev/null
+++ b/src/determinant.c
@@ -0,0 +1,93 @@
+#include "determinant.h"
+
+#include "number.h"
+#include "rational.h"
+#include "definitions.cpp"
+
+// determinant of a matrix
+// replaces the matrix by its LU decomposition
+int determinant_inplace(Number_Matrix M, Number* out){
+  int i;
+  int sign_correction;
+
+  LU_dcmp_inplace(M, &sign_correction);
+
+  if(sign_correction==0){
+    *out=number_zero();
+    return(0);
+  }
+
+  *out=number_one();
+  if(sign_correction==-1){
+    number_Qprod_chain(quot(-1,1), out);
+  }
+
+  for(i=0;i<M.length;i++){
+    number_prod_chain(M.matrix[i][i], out);
+  }
+
+  return(0);
+}
+
+// LU decomposition
+// uses pivoting to avoid dividing by 0
+// the sign_correction should be multiplied to the determinant to obtain the right value
+// if dividing by 0 is unavoidable, then the determinant is 0, and sign_correction is set to 0
+int LU_dcmp_inplace(Number_Matrix M, int* sign_correction){
+  int i,j,k,pivot;
+  Number tmp;
+
+  *sign_correction=1;
+
+  for(j=0;j<M.length;j++){
+    for(i=0;i<=j;i++){
+      for(k=0;k<i;k++){
+	// -M[i][k]*M[k][j]
+	number_prod(M.matrix[i][k], M.matrix[k][j], &tmp);
+	number_Qprod_chain(quot(-1,1), &tmp);
+	number_add_chain(tmp, M.matrix[i]+j);
+	free_Number(tmp);
+      }
+    }
+    for(i=j+1;i<M.length;i++){
+      for(k=0;k<j;k++){
+	// -M[i][k]*M[k][j]
+	number_prod(M.matrix[i][k], M.matrix[k][j], &tmp);
+	number_Qprod_chain(quot(-1,1), &tmp);
+	number_add_chain(tmp, M.matrix[i]+j);
+	free_Number(tmp);
+      }
+    }
+
+    // pivot if M[j][j]==0
+    // find first M[j][j] that is not 0
+    for(pivot=j;pivot<M.length && number_is_zero(M.matrix[pivot][j])==1;pivot++){}
+
+    // no non-zero M[j][j] left: return
+    if(pivot>=M.length){
+      *sign_correction=0;
+      return(0);
+    }
+    // pivot if needed
+    if(pivot!=j){
+      for(k=0;k<M.length;k++){
+	tmp=M.matrix[j][k];
+	M.matrix[j][k]=M.matrix[pivot][k];
+	M.matrix[pivot][k]=tmp;
+      }
+      *sign_correction*=-1;
+
+    }
+
+    for(i=j+1;i<M.length;i++){
+      // do not use the inplace algorithm if M[j][j] has more than one terms, since it would be modified by the inplace function
+      if(M.matrix[j][j].length<=1){
+	number_quot_inplace(M.matrix[i]+j, M.matrix[j]+j);
+      }
+      else{
+	number_quot_chain(M.matrix[i]+j, M.matrix[j][j]);
+      }
+    }
+  }
+  return(0);
+}
diff --git a/src/determinant.h b/src/determinant.h
new file mode 100644
index 0000000..613be87
--- /dev/null
+++ b/src/determinant.h
@@ -0,0 +1,17 @@
+/*
+  Compute the determinant of a number matrix
+*/
+
+#ifndef DETERMINANT_H
+#define DETERMINANT_H
+
+#include "types.h"
+
+// determinant of a matrix
+int determinant_inplace(Number_Matrix M, Number* out);
+
+// LU decomposition
+int LU_dcmp_inplace(Number_Matrix M, int* sign_correction);
+
+
+#endif
diff --git a/src/expansions.c b/src/expansions.c
deleted file mode 100644
index c889829..0000000
--- a/src/expansions.c
+++ /dev/null
@@ -1,28 +0,0 @@
-/*
-Copyright 2015 Ian Jauslin
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-
-    http://www.apache.org/licenses/LICENSE-2.0
-
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
-
-#include "expansions.h"
-
-#include <stdio.h>
-#include <stdlib.h>
-#include "definitions.cpp"
-#include "tools.h"
-#include "array.h"
-#include "polynomial.h"
-#include "number.h"
-#include "rational.h"
-
-
diff --git a/src/expansions.h b/src/expansions.h
deleted file mode 100644
index 85d6c2b..0000000
--- a/src/expansions.h
+++ /dev/null
@@ -1,34 +0,0 @@
-/*
-Copyright 2015 Ian Jauslin
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-
-    http://www.apache.org/licenses/LICENSE-2.0
-
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
-
-/*
-Compute exp(V) and log(1+W)
-*/
-
-#ifndef EXPANSIONS_H
-#define EXPANSIONS_H
-
-#include "polynomial.h"
-#include "fields.h"
-
-// exp(V)
-int expand_exponential(Polynomial input_polynomial,Polynomial* output, Fields_Table fields);
-
-// log(1+W)
-int expand_logarithm(Polynomial input_polynomial, Polynomial* output, Fields_Table fields);
-
-
-#endif
diff --git a/src/fields.c b/src/fields.c
index bfd1f39..7be84df 100644
--- a/src/fields.c
+++ b/src/fields.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -23,6 +23,13 @@ limitations under the License.
 #include "polynomial.h"
 #include "array.h"
 #include "rational.h"
+#include "tree.h"
+
+//---------------------------------------------------------------------
+//
+//                              Fields_Table
+//
+//---------------------------------------------------------------------
 
 // init and free for Fields_Table
 int init_Fields_Table(Fields_Table* fields){
@@ -30,7 +37,7 @@ int init_Fields_Table(Fields_Table* fields){
   init_Int_Array(&((*fields).external),FIELDS_SIZE);
   init_Int_Array(&((*fields).internal),FIELDS_SIZE);
   init_Identities(&((*fields).ids), FIELDS_SIZE);
-  init_Symbols(&((*fields).symbols), FIELDS_SIZE);
+  init_Virtual_fields(&((*fields).virtual_fields), FIELDS_SIZE);
   init_Int_Array(&((*fields).fermions),FIELDS_SIZE);
   init_Int_Array(&((*fields).noncommuting),FIELDS_SIZE);
   return(0);
@@ -40,7 +47,7 @@ int free_Fields_Table(Fields_Table fields){
   free_Int_Array(fields.external);
   free_Int_Array(fields.internal);
   free_Identities(fields.ids);
-  free_Symbols(fields.symbols);
+  free_Virtual_fields(fields.virtual_fields);
   free_Int_Array(fields.fermions);
   free_Int_Array(fields.noncommuting);
   return(0);
@@ -57,11 +64,11 @@ int field_type(int index, Fields_Table fields){
   else if(int_array_find(abs(index), fields.internal)>=0){
     return(FIELD_INTERNAL);
   }
-  else if(intlist_find(fields.symbols.indices, fields.symbols.length, index)>=0){
-    return(FIELD_SYMBOL);
+  else if(intlist_find(fields.virtual_fields.indices, fields.virtual_fields.length, index)>=0){
+    return(FIELD_VIRTUAL);
   }
 
-  fprintf(stderr,"error: index %d is neither a parameter nor an external or an internal field, nor a symbol\n",index);
+  fprintf(stderr,"error: index %d is neither a parameter nor an external or an internal field, nor a virtual field\n",index);
   exit(-1);
 }
 
@@ -86,7 +93,11 @@ int is_noncommuting(int index, Fields_Table fields){
 }
 
 
-// ------------------ Identities --------------------
+//---------------------------------------------------------------------
+//
+//                               Identities
+//
+//---------------------------------------------------------------------
 
 // allocate memory
 int init_Identities(Identities* identities,int size){
@@ -314,6 +325,11 @@ int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fi
   int match_nc;
   int first=-1;
 
+  // cannot fit
+  if(test_array.length<input.length){
+    return(-1);
+  }
+
   // bound noncommuting elements
   while(is_noncommuting(input.values[post_nc], fields)==1){
     post_nc++;
@@ -324,7 +340,7 @@ int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fi
       match_nc=1;
     }
     for(j=1;j<post_nc;j++){
-      if(test_array.values[i+j]!=input.values[j]){
+      if(i+j>=test_array.length || test_array.values[i+j]!=input.values[j]){
 	match_nc=0;
       }
     }
@@ -359,57 +375,61 @@ int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fi
 }
 
 
-// ------------------ Symbols --------------------
+//---------------------------------------------------------------------
+//
+//                                Virtual_fields
+//
+//---------------------------------------------------------------------
 
 // allocate memory
-int init_Symbols(Symbols* symbols,int size){
-  (*symbols).indices=calloc(size,sizeof(int));
-  (*symbols).expr=calloc(size,sizeof(Polynomial));
-  (*symbols).length=0;
-  (*symbols).memory=size;
+int init_Virtual_fields(Virtual_fields* virtual_fields,int size){
+  (*virtual_fields).indices=calloc(size,sizeof(int));
+  (*virtual_fields).expr=calloc(size,sizeof(Polynomial));
+  (*virtual_fields).length=0;
+  (*virtual_fields).memory=size;
   return(0);
 }
 
 // free memory
-int free_Symbols(Symbols symbols){
+int free_Virtual_fields(Virtual_fields virtual_fields){
   int i;
-  for(i=0;i<symbols.length;i++){
-    free_Polynomial(symbols.expr[i]);
+  for(i=0;i<virtual_fields.length;i++){
+    free_Polynomial(virtual_fields.expr[i]);
   }
-  free(symbols.indices);
-  free(symbols.expr);
+  free(virtual_fields.indices);
+  free(virtual_fields.expr);
   return(0);
 }
 
 // resize
-int resize_symbols(Symbols* symbols,int new_size){
-  Symbols new_symbols;
+int resize_virtual_fields(Virtual_fields* virtual_fields,int new_size){
+  Virtual_fields new_virtual_fields;
   int i;
 
-  init_Symbols(&new_symbols,new_size);
-  for(i=0;i<(*symbols).length;i++){
-    new_symbols.indices[i]=(*symbols).indices[i];
-    new_symbols.expr[i]=(*symbols).expr[i];
+  init_Virtual_fields(&new_virtual_fields,new_size);
+  for(i=0;i<(*virtual_fields).length;i++){
+    new_virtual_fields.indices[i]=(*virtual_fields).indices[i];
+    new_virtual_fields.expr[i]=(*virtual_fields).expr[i];
   }
-  new_symbols.length=(*symbols).length;
+  new_virtual_fields.length=(*virtual_fields).length;
 
-  free((*symbols).indices);
-  free((*symbols).expr);
+  free((*virtual_fields).indices);
+  free((*virtual_fields).expr);
 
-  *symbols=new_symbols;
+  *virtual_fields=new_virtual_fields;
   return(0);
 }
 
 // copy
-int symbols_cpy(Symbols input, Symbols* output){
-  init_Symbols(output,input.length);
-  symbols_cpy_noinit(input,output);
+int virtual_fields_cpy(Virtual_fields input, Virtual_fields* output){
+  init_Virtual_fields(output,input.length);
+  virtual_fields_cpy_noinit(input,output);
   return(0);
 }
-int symbols_cpy_noinit(Symbols input, Symbols* output){
+int virtual_fields_cpy_noinit(Virtual_fields input, Virtual_fields* output){
   int i;
   if((*output).memory<input.length){
-    fprintf(stderr,"error: trying to copy a symbols collection of length %d to another with memory %d\n",input.length,(*output).memory);
+    fprintf(stderr,"error: trying to copy a virtual fields collection of length %d to another with memory %d\n",input.length,(*output).memory);
     exit(-1);
   }
   for(i=0;i<input.length;i++){
@@ -421,12 +441,12 @@ int symbols_cpy_noinit(Symbols input, Symbols* output){
   return(0);
 }
 
-// append an element to a symbols
-int symbols_append(int index, Polynomial expr, Symbols* output){
+// append an element to a virtual_fields
+int virtual_fields_append(int index, Polynomial expr, Virtual_fields* output){
   int offset=(*output).length;
 
   if((*output).length>=(*output).memory){
-    resize_symbols(output,2*(*output).memory+1);
+    resize_virtual_fields(output,2*(*output).memory+1);
   }
 
   // copy and allocate
@@ -436,12 +456,12 @@ int symbols_append(int index, Polynomial expr, Symbols* output){
   (*output).length++;
   return(0);
 }
-// append an element to a symbols without allocating memory
-int symbols_append_noinit(int index, Polynomial expr, Symbols* output){
+// append an element to a virtual_fields without allocating memory
+int virtual_fields_append_noinit(int index, Polynomial expr, Virtual_fields* output){
   int offset=(*output).length;
 
   if((*output).length>=(*output).memory){
-    resize_symbols(output,2*(*output).memory+1);
+    resize_virtual_fields(output,2*(*output).memory+1);
   }
 
   // copy without allocating
@@ -452,18 +472,22 @@ int symbols_append_noinit(int index, Polynomial expr, Symbols* output){
   return(0);
 }
 
-// concatenate two symbolss
-int symbols_concat(Symbols input, Symbols* output){
+// concatenate two virtual_fields
+int virtual_fields_concat(Virtual_fields input, Virtual_fields* output){
   int i;
   for(i=0;i<input.length;i++){
-    symbols_append(input.indices[i],input.expr[i],output);
+    virtual_fields_append(input.indices[i],input.expr[i],output);
   }
   return(0);
 }
 
 
 
-// ------------------ Groups --------------------
+//---------------------------------------------------------------------
+//
+//                                 Groups
+//
+//---------------------------------------------------------------------
 
 // allocate memory
 int init_Groups(Groups* groups,int size){
@@ -570,3 +594,164 @@ int find_group(int index, Groups groups){
   }
   return(-1);
 }
+
+
+//---------------------------------------------------------------------
+//
+//                               Variables
+//
+//---------------------------------------------------------------------
+
+// allocate memory
+int init_Variables(Variables* variables,int size){
+  (*variables).var_names=calloc(size,sizeof(Char_Array));
+  (*variables).symbol_trees=calloc(size,sizeof(Tree));
+  (*variables).length=0;
+  (*variables).memory=size;
+  return(0);
+}
+
+// free memory
+int free_Variables(Variables variables){
+  int i;
+  for(i=0;i<variables.length;i++){
+    free_Char_Array(variables.var_names[i]);
+    free_Tree(variables.symbol_trees[i]);
+  }
+  free(variables.var_names);
+  free(variables.symbol_trees);
+  return(0);
+}
+
+// resize
+int resize_variables(Variables* variables,int new_size){
+  Variables new_variables;
+  int i;
+
+  init_Variables(&new_variables,new_size);
+  for(i=0;i<(*variables).length;i++){
+    new_variables.var_names[i]=(*variables).var_names[i];
+    new_variables.symbol_trees[i]=(*variables).symbol_trees[i];
+  }
+  new_variables.length=(*variables).length;
+
+  free((*variables).var_names);
+  free((*variables).symbol_trees);
+
+  *variables=new_variables;
+  return(0);
+}
+
+// copy
+int variables_cpy(Variables input, Variables* output){
+  init_Variables(output,input.length);
+  variables_cpy_noinit(input,output);
+  return(0);
+}
+int variables_cpy_noinit(Variables input, Variables* output){
+  int i;
+  if((*output).memory<input.length){
+    fprintf(stderr,"error: trying to copy a variables collection of length %d to another with memory %d\n",input.length,(*output).memory);
+    exit(-1);
+  }
+  for(i=0;i<input.length;i++){
+    char_array_cpy(input.var_names[i], (*output).var_names+i);
+    tree_cpy(input.symbol_trees[i],(*output).symbol_trees+i);
+  }
+  (*output).length=input.length;
+  
+  return(0);
+}
+
+// append an element to a variables collection
+int variables_append(Char_Array var_name, Tree symbol_tree, Variables* output){
+  int offset=(*output).length;
+
+  if((*output).length>=(*output).memory){
+    resize_variables(output,2*(*output).memory+1);
+  }
+
+  // copy and allocate
+  char_array_cpy(var_name,(*output).var_names+offset);
+  tree_cpy(symbol_tree,(*output).symbol_trees+offset);
+  // increment length
+  (*output).length++;
+  return(0);
+}
+// append an element to a variables collection without allocating memory
+int variables_append_noinit(Char_Array var_name, Tree symbol_tree, Variables* output){
+  int offset=(*output).length;
+
+  if((*output).length>=(*output).memory){
+    resize_variables(output,2*(*output).memory+1);
+  }
+
+  // copy without allocating
+  (*output).var_names[offset]=var_name;
+  (*output).symbol_trees[offset]=symbol_tree;
+  // increment length
+  (*output).length++;
+  return(0);
+}
+
+// concatenate two variables collections
+int variables_concat(Variables input, Variables* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    variables_append(input.var_names[i], input.symbol_trees[i], output);
+  }
+  return(0);
+}
+
+
+// find a variable matching a var_name
+int variables_find_var(Char_Array name, Variables variables, Tree* output){
+  Char_Array varname;
+  int i;
+
+  // drop inital '$'
+  char_array_substring(name, 1, name.length-1, &varname);
+  for(i=0;i<variables.length;i++){
+    if(char_array_cmp(varname, variables.var_names[i])==1){
+      tree_cpy(variables.symbol_trees[i], output);
+      break;
+    }
+  }
+
+  // error if no variable was found
+  if(i==variables.length){
+    fprintf(stderr, "error: variable '$%s' not found\n",char_array_to_str_noinit(&varname));
+    exit(-1);
+  }
+  free_Char_Array(varname);
+
+  return(0);
+}
+
+
+// add a polynomials as a new named variable
+int add_polynomial_to_variables(char* name, Polynomial polynomial, Variables* variables){
+  // save polynomial to string (to convert it to a variable, it must first be a string)
+  Char_Array poly_str;
+  Char_Array out_name;
+  Tree out_tree;
+
+  init_Char_Array(&poly_str, STR_SIZE);
+  polynomial_sprint(polynomial, &poly_str);
+
+  // convert name to Char_Array
+  init_Char_Array(&out_name,STR_SIZE);
+  char_array_append_str(name, &out_name);
+
+  // trivial tree containing the polynomial
+  init_Tree(&out_tree,0,poly_str.length);
+  tree_set_label(poly_str, &out_tree);
+  free_Char_Array(poly_str);
+
+  // add variable
+  variables_append(out_name, out_tree, variables);
+  free_Tree(out_tree);
+  free_Char_Array(out_name);
+
+  return 0;
+}
diff --git a/src/fields.h b/src/fields.h
index b6f05c2..7b7a07a 100644
--- a/src/fields.h
+++ b/src/fields.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -21,6 +21,9 @@ limitations under the License.
 
 #include "types.h"
 
+
+// Fields_Table
+
 // init
 int init_Fields_Table(Fields_Table* fields);
 int free_Fields_Table(Fields_Table fields);
@@ -33,6 +36,8 @@ int is_fermion(int index, Fields_Table fields);
 int is_noncommuting(int index, Fields_Table fields);
 
 
+// Identities
+
 // init
 int init_Identities(Identities* identities,int size);
 int free_Identities(Identities identities);
@@ -57,25 +62,29 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields);
 int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fields_Table fields);
 
 
+// Virtual_fields
+
 // init
-int init_Symbols(Symbols* symbols,int size);
-int free_Symbols(Symbols symbols);
+int init_Virtual_fields(Virtual_fields* virtual_fields,int size);
+int free_Virtual_fields(Virtual_fields virtual_fields);
 
 // resize
-int resize_symbols(Symbols* symbols,int new_size);
+int resize_virtual_fields(Virtual_fields* virtual_fields,int new_size);
 
 // copy
-int symbols_cpy(Symbols input, Symbols* output);
-int symbols_cpy_noinit(Symbols input, Symbols* output);
+int virtual_fields_cpy(Virtual_fields input, Virtual_fields* output);
+int virtual_fields_cpy_noinit(Virtual_fields input, Virtual_fields* output);
 
-// append an element to a symbols
-int symbols_append(int index, Polynomial expr, Symbols* output);
-int symbols_append_noinit(int index, Polynomial expr, Symbols* output);
+// append an element to a virtual_fields
+int virtual_fields_append(int index, Polynomial expr, Virtual_fields* output);
+int virtual_fields_append_noinit(int index, Polynomial expr, Virtual_fields* output);
 
-// concatenate two symbolss
-int symbols_concat(Symbols input, Symbols* output);
+// concatenate two virtual_fields
+int virtual_fields_concat(Virtual_fields input, Virtual_fields* output);
 
 
+// Groups
+
 // init
 int init_Groups(Groups* groups,int size);
 int free_Groups(Groups groups);
@@ -98,5 +107,32 @@ int groups_concat(Groups input, Groups* output);
 int find_group(int index, Groups groups);
 
 
+// Variables
+
+// allocate memory
+int init_Variables(Variables* variables,int size);
+// free memory
+int free_Variables(Variables variables);
+
+// resize
+int resize_variables(Variables* variables,int new_size);
+
+// copy
+int variables_cpy(Variables input, Variables* output);
+int variables_cpy_noinit(Variables input, Variables* output);
+
+// append an element to a variables collection
+int variables_append(Char_Array var_name, Tree symbol_tree, Variables* output);
+int variables_append_noinit(Char_Array var_name, Tree symbol_tree, Variables* output);
+
+// concatenate two variables collections
+int variables_concat(Variables input, Variables* output);
+
+// find a variable matching a var_name
+int variables_find_var(Char_Array name, Variables variables, Tree* output);
+
+// add a polynomials as a new named variable
+int add_polynomial_to_variables(char* name, Polynomial polynomial, Variables* variables);
+
 #define FIELDS_H_DONE
 #endif
diff --git a/src/flow.c b/src/flow.c
index b294bb0..10c5ae0 100644
--- a/src/flow.c
+++ b/src/flow.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -25,14 +25,20 @@ limitations under the License.
 #include "array.h"
 #include "coefficient.h"
 #include "rcc.h"
+#include "grouped_polynomial.h"
 
 
 // compute flow numerically, no exponentials
-int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, int niter, int display_mode){
+int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Grouped_Polynomial postprocess_flow_equation, Labels labels, int niter, int display_mode){
   // running coupling contants
   RCC rccs=init;
+  // for printing
+  RCC rcc_print;
   int i,j;
 
+  // init printing rcc
+  init_RCC(&rcc_print, rccs.length);
+
   if(display_mode==DISPLAY_NUMERICAL){
     // print labels
     printf("%5s  ","n");
@@ -42,9 +48,25 @@ int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, in
     printf("\n\n");
 
     // print initial values
+    RCC_cpy_noinit(rccs,&rcc_print);
+    if(postprocess_flow_equation.length>0){
+      // ignore constants
+      for(j=0;j<rcc_print.length;j++){
+	if(rcc_print.indices[j]<0){
+	  rcc_print.values[j]=1.;
+	}
+      }
+      evaleq(rcc_print, rcc_print, postprocess_flow_equation);
+    }
     printf("%5d  ",0);
-    for(j=0;j<rccs.length;j++){
-      printf("% 14.7Le  ",rccs.values[j]);
+    for(j=0;j<rcc_print.length;j++){
+      // use constants from rcc
+      if(rcc_print.indices[j]<0){
+	printf("% 14.7Le  ",rccs.values[j]);
+      }
+      else{
+	printf("% 14.7Le  ",rcc_print.values[j]);
+      }
     }
     printf("\n");
   }
@@ -52,12 +74,32 @@ int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, in
   for(i=0;i<niter;i++){
     // compute a single step
     step_flow(&rccs, flow_equation);
-    // convert ls to alphas
+
+    // print
+    if(postprocess_flow_equation.length>0){
+      RCC_cpy_noinit(rccs,&rcc_print);
+      // ignore constants
+      for(j=0;j<rcc_print.length;j++){
+	if(rcc_print.indices[j]<0){
+	  rcc_print.values[j]=1.;
+	}
+      }
+      evaleq(rcc_print, rcc_print, postprocess_flow_equation);
+    }
+    else{
+      RCC_cpy_noinit(rccs,&rcc_print);
+    }
     if(display_mode==DISPLAY_NUMERICAL){
       // print the result
       printf("%5d  ",i+1);
-      for(j=0;j<rccs.length;j++){
-	printf("% 14.7Le  ",rccs.values[j]);
+      for(j=0;j<rcc_print.length;j++){
+	// use constants from rcc
+	if(rcc_print.indices[j]<0){
+	  printf("% 14.7Le  ",rccs.values[j]);
+	}
+	else{
+	  printf("% 14.7Le  ",rcc_print.values[j]);
+	}
       }
       printf("\n");
     }
@@ -74,9 +116,24 @@ int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, in
   }
 
   if(display_mode==DISPLAY_FINAL){
-    RCC_print(rccs);
+    if(postprocess_flow_equation.length>0){
+      RCC_cpy_noinit(rccs,&rcc_print);
+      // ignore constants
+      for(j=0;j<rcc_print.length;j++){
+	if(rcc_print.indices[j]<0){
+	  rcc_print.values[j]=1.;
+	}
+      }
+      evaleq(rcc_print, rcc_print, postprocess_flow_equation);
+    }
+    else{
+      RCC_cpy_noinit(rccs,&rcc_print);
+    }
+    RCC_print(rcc_print);
   }
 
+  free_RCC(rcc_print);
+
   return(0);
 }
 
diff --git a/src/flow.h b/src/flow.h
index 9491323..ad43b72 100644
--- a/src/flow.h
+++ b/src/flow.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -24,7 +24,7 @@ Compute flow numerically
 #include "types.h"
 
 // compute flow
-int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, int niter, int display_mode);
+int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Grouped_Polynomial postprocess_flow_equation, Labels labels, int niter, int display_mode);
 // single step
 int step_flow(RCC* rccs, Grouped_Polynomial flow_equation);
 
diff --git a/src/flow_mpfr.c b/src/flow_mpfr.c
index e03c130..1701b15 100644
--- a/src/flow_mpfr.c
+++ b/src/flow_mpfr.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -32,14 +32,21 @@ limitations under the License.
 #include "coefficient.h"
 #include "flow.h"
 #include "rcc_mpfr.h"
+#include "grouped_polynomial.h"
 
 
 
 // compute flow numerically
-int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels labels, int niter, int display_mode){
+int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Grouped_Polynomial postprocess_flow_equation, Labels labels, int niter, int display_mode){
   // running coupling contants
   RCC_mpfr rccs=init;
   int i,j;
+  // for printing
+  RCC_mpfr rcc_print;
+
+
+  // init printing rcc
+  init_RCC_mpfr(&rcc_print, rccs.length);
 
   if(display_mode==DISPLAY_NUMERICAL){
     // print labels
@@ -50,9 +57,25 @@ int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels
     printf("\n\n");
 
     // print initial values
+    RCC_mpfr_cpy_noinit(rccs,&rcc_print);
+    if(postprocess_flow_equation.length>0){
+      // ignore constants
+      for(j=0;j<rcc_print.length;j++){
+	if(rcc_print.indices[j]<0){
+	  mpfr_set_ui(rcc_print.values[j], 1, MPFR_RNDN);
+	}
+      }
+      evaleq_mpfr(rcc_print, rccs, postprocess_flow_equation);
+    }
     printf("%5d  ",0);
-    for(j=0;j<rccs.length;j++){
-      mpfr_printf("% 14.7Re  ",rccs.values[j]);
+    for(j=0;j<rcc_print.length;j++){
+      // use constants from rcc
+      if(rcc_print.indices[j]<0){
+	mpfr_printf("% 14.7Re  ",rccs.values[j]);
+      }
+      else{
+	mpfr_printf("% 14.7Re  ",rcc_print.values[j]);
+      }
     }
     printf("\n");
   }
@@ -60,12 +83,29 @@ int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels
   for(i=0;i<niter;i++){
     // compute a single step
     step_flow_mpfr(&rccs, flow_equation);
-    // convert ls to alphas
+
+    // print
+    RCC_mpfr_cpy_noinit(rccs,&rcc_print);
+    if(postprocess_flow_equation.length>0){
+      // ignore constants
+      for(j=0;j<rcc_print.length;j++){
+	if(rcc_print.indices[j]<0){
+	  mpfr_set_ui(rcc_print.values[j], 1, MPFR_RNDN);
+	}
+      }
+      evaleq_mpfr(rcc_print, rccs, postprocess_flow_equation);
+    }
     if(display_mode==DISPLAY_NUMERICAL){
       // print the result
       printf("%5d  ",i+1);
-      for(j=0;j<rccs.length;j++){
-	mpfr_printf("% 14.7Re  ",rccs.values[j]);
+      for(j=0;j<rcc_print.length;j++){
+	// use constants from rcc
+	if(rcc_print.indices[j]<0){
+	  mpfr_printf("% 14.7Re  ",rccs.values[j]);
+	}
+	else{
+	  mpfr_printf("% 14.7Re  ",rcc_print.values[j]);
+	}
       }
       printf("\n");
     }
@@ -82,9 +122,16 @@ int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels
   }
 
   if(display_mode==DISPLAY_FINAL){
-    RCC_mpfr_print(rccs);
+    if(postprocess_flow_equation.length>0){
+      evaleq_mpfr(rcc_print, rccs, postprocess_flow_equation);
+    }
+    else{
+      rcc_print=rccs;
+    }
+    RCC_mpfr_print(rcc_print);
   }
 
+  free_RCC_mpfr(rcc_print);
   return(0);
 }
 
diff --git a/src/flow_mpfr.h b/src/flow_mpfr.h
index e06cbd8..203770d 100644
--- a/src/flow_mpfr.h
+++ b/src/flow_mpfr.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -25,7 +25,7 @@ Compute flow numerically
 #include "types.h"
 
 // compute flow
-int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels labels, int niter, int display_mode);
+int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Grouped_Polynomial postprocess_flow_equation, Labels labels, int niter, int display_mode);
 // single step
 int step_flow_mpfr(RCC_mpfr* rccs, Grouped_Polynomial flow_equation);
 
diff --git a/src/grouped_polynomial.c b/src/grouped_polynomial.c
index 0d4d75d..9e245a8 100644
--- a/src/grouped_polynomial.c
+++ b/src/grouped_polynomial.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -205,9 +205,9 @@ int group_polynomial(Polynomial polynomial, Grouped_Polynomial* grouped_polynomi
 
     if(index==-2){
       fprintf(stderr,"error: monomial (");
-      for(j=0;j<polynomial.monomials[i].length;j++){
-	fprintf(stderr,"%d", polynomial.monomials[i].values[j]);
-	if(j<polynomial.monomials[i].length-1){
+      for(j=0;j<remainder.monomials[i].length;j++){
+	fprintf(stderr,"%d", remainder.monomials[i].values[j]);
+	if(j<remainder.monomials[i].length-1){
 	  fprintf(stderr,",");
 	}
       }
@@ -436,7 +436,7 @@ int simplify_grouped_polynomial(Grouped_Polynomial* polynomial){
 }
 
 
-// derive a flow equation with respect to an unknown variable
+// differentiate a flow equation with respect to an unknown variable
 // equivalent to DB.dl where dl are symbols for the derivatives of the indices in the flow equation with respect to the unknown variable
 // indices specifies the list of indices that depend on the variable
 int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Grouped_Polynomial* dflow){
@@ -487,7 +487,7 @@ int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Gr
 
 
 /*
-// derive a flow equation with respect to an index
+// differentiate a flow equation with respect to an index
 int flow_equation_deriv(Grouped_Polynomial flow_equation, int index, Grouped_Polynomial* output){
   int i,k;
   // temp list of indices
@@ -603,7 +603,7 @@ int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre
     init_Char_Array(&buffer, STR_SIZE);
     coefficient_sprint(grouped_polynomial.coefs[i],&buffer,9,rhs_pre);
     if(buffer.length>0){
-      printf("%s",buffer.str);
+      printf("%s",char_array_to_str_noinit(&buffer));
     }
     free_Char_Array(buffer);
 
@@ -732,29 +732,33 @@ int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output)
 }
 
 
-// eValuate an equation on a vector
-int evaleq(RCC* rccs, Grouped_Polynomial poly){
+// evaluate an equation on a vector
+int evaleq(RCC out, RCC in, Grouped_Polynomial poly){
   int i;
-  long double* res=calloc((*rccs).length,sizeof(long double));
+  long double* res=calloc(out.length,sizeof(long double));
 
-  if((*rccs).length!=poly.length){
-    fprintf(stderr, "error: trying to evaluate an flow equation with %d components on an rcc with %d\n",poly.length,(*rccs).length);
+  if(in.length!=poly.length){
+    fprintf(stderr, "error: trying to evaluate a flow equation with %d components on an rcc with %d\n",poly.length,in.length);
+    exit(-1);
+  }
+  if(out.length!=poly.length){
+    fprintf(stderr, "error: trying to write the output of a flow equation with %d components on an rcc with %d\n",poly.length,out.length);
     exit(-1);
   }
 
-  // initialize vectors to 0
-  for(i=0;i<(*rccs).length;i++){
+  // initialize vectors to 0 in an auxiliary vector (to allow for out=in without interference)
+  for(i=0;i<in.length;i++){
     res[i]=0.;
   }
 
   // for each equation
   for(i=0;i<poly.length;i++){
-    evalcoef(*rccs, poly.coefs[i], res+i);
+    evalcoef(in, poly.coefs[i], res+i);
   }
 
   // copy res to rccs
-  for(i=0;i<(*rccs).length;i++){
-    (*rccs).values[i]=res[i];
+  for(i=0;i<out.length;i++){
+    out.values[i]=res[i];
   }
 
   // free memory
@@ -763,25 +767,29 @@ int evaleq(RCC* rccs, Grouped_Polynomial poly){
 
 }
 // evaluate an equation on a vector (using mpfr floats)
-int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly){
+int evaleq_mpfr(RCC_mpfr out, RCC_mpfr in, Grouped_Polynomial poly){
   int i;
   mpfr_t* res;
 
-  if((*rccs).length!=poly.length){
-    fprintf(stderr, "error: trying to evaluate an flow equation with %d components on an rcc with %d\n",poly.length,(*rccs).length);
+  if(in.length!=poly.length){
+    fprintf(stderr, "error: trying to evaluate a flow equation with %d components on an rcc with %d\n",poly.length,in.length);
+    exit(-1);
+  }
+  if(out.length!=poly.length){
+    fprintf(stderr, "error: trying to write the output of a flow equation with %d components on an rcc with %d\n",poly.length,out.length);
     exit(-1);
   }
 
-  res=calloc((*rccs).length,sizeof(mpfr_t));
+  res=calloc(out.length,sizeof(mpfr_t));
 
   // for each equation
   for(i=0;i<poly.length;i++){
-    evalcoef_mpfr(*rccs, poly.coefs[i], res[i]);
+    evalcoef_mpfr(in, poly.coefs[i], res[i]);
   }
 
   // copy res to rccs
-  for(i=0;i<(*rccs).length;i++){
-    mpfr_set((*rccs).values[i], res[i], MPFR_RNDN);
+  for(i=0;i<out.length;i++){
+    mpfr_set(out.values[i], res[i], MPFR_RNDN);
     mpfr_clear(res[i]);
   }
 
@@ -791,3 +799,85 @@ int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly){
 }
 
 
+// compose two flow equations (replace the rcc's of flow1 by the right hand side of flow2)
+int compose_flow_equations(Grouped_Polynomial flow1, Grouped_Polynomial flow2, Grouped_Polynomial* out){
+  if(flow1.length!=flow2.length){
+    fprintf(stderr, "error: trying to compose two flow equations of different size\n");
+    exit(-1);
+  }
+  int i,j,k;
+  Coefficient constant;
+
+  // init
+  init_Grouped_Polynomial(out, flow1.length);
+  (*out).length=flow1.length;
+
+  // init constant (so we can tell when the constant was not found)
+  constant.length=0;
+
+  // loop over rcc's
+  for(i=0;i<flow1.length;i++){
+    // set indices
+    (*out).indices[i]=flow1.indices[i];
+
+    // passthrough constant terms
+    if((*out).indices[i]<0){
+      int index=intlist_find_err(flow2.indices,flow2.length,(*out).indices[i]);
+      coefficient_cpy(flow2.coefs[index], (*out).coefs+i);
+      constant=flow2.coefs[index];
+      continue;
+    }
+
+    // init
+    init_Coefficient((*out).coefs+i, COEF_SIZE);
+
+    // loop over terms
+    for(j=0;j<flow1.coefs[i].length;j++){
+      Coefficient tmp_coef;
+
+      // init
+      init_Coefficient(&tmp_coef, COEF_SIZE);
+      // init factor
+      Int_Array tmp_factor;
+      init_Int_Array(&tmp_factor, MONOMIAL_SIZE);
+      // init denom
+      coef_denom denom;
+      // index should be that appearing in flow2
+      if(flow2.coefs[i].length<1){
+	fprintf(stderr,"error: composing two flow equations: the %d-th term in the flow equation is empty\n",flow1.indices[i]);
+	exit(-1);
+      }
+      denom.index=flow2.coefs[i].denoms[0].index;
+      denom.power=0;
+      // init num
+      Number tmp_num;
+      number_cpy(flow1.coefs[i].nums[j], &tmp_num);
+
+      // init coefficient with numerical prefactor
+      coefficient_append_noinit(tmp_factor, tmp_num, denom, &tmp_coef);
+
+      // loop over factors
+      for(k=0;k<flow1.coefs[i].factors[j].length;k++){
+	// multiply factors together
+	coefficient_prod_chain(flow2.coefs[intlist_find_err(flow2.indices,flow2.length,flow1.coefs[i].factors[j].values[k])], &tmp_coef);
+      }
+
+      // add to out
+      coefficient_concat_noinit(tmp_coef, (*out).coefs+i);
+    }
+
+  }
+
+  // simplify fractions
+  if(constant.length!=0){
+    for(i=0;i<(*out).length;i++){
+      if((*out).indices[i]>=0){
+	// reduce them to a common denominator (not much is gained from trying to simplify them)
+	coefficient_common_denominator(constant, (*out).coefs+i);
+	//coefficient_simplify_rational(constant, (*out).coefs+i);
+      }
+    }
+  }
+
+  return(0);
+}
diff --git a/src/grouped_polynomial.h b/src/grouped_polynomial.h
index e369721..7d95a2e 100644
--- a/src/grouped_polynomial.h
+++ b/src/grouped_polynomial.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -59,7 +59,7 @@ int find_id(Int_Array monomial, Id_Table idtable, int start);
 // simplify grouped polynomial
 int simplify_grouped_polynomial(Grouped_Polynomial* polynomial);
 
-// derive a flow equation with respect to an unknown variable
+// differentiate a flow equation with respect to an unknown variable
 int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Grouped_Polynomial* dflow);
 
 // print a grouped polynomial
@@ -69,8 +69,10 @@ int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre
 int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output);
 
 // evaluate an equation on an RCC
-int evaleq(RCC* rccs, Grouped_Polynomial poly);
+int evaleq(RCC out, RCC in, Grouped_Polynomial poly);
 // evaluate an equation on a vector (using mpfr floats)
-int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly);
+int evaleq_mpfr(RCC_mpfr out, RCC_mpfr in, Grouped_Polynomial poly);
 
+// compose two flow equations (replace the rcc's of flow1 by the right hand side of flow2)
+int compose_flow_equations(Grouped_Polynomial flow1, Grouped_Polynomial flow2, Grouped_Polynomial* out);
 #endif
diff --git a/src/idtable.c b/src/idtable.c
index fe409c2..fc72f56 100644
--- a/src/idtable.c
+++ b/src/idtable.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/src/idtable.h b/src/idtable.h
index 38ab249..aa99dca 100644
--- a/src/idtable.h
+++ b/src/idtable.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/src/istring.c b/src/istring.c
index 663c06a..d78f67b 100644
--- a/src/istring.c
+++ b/src/istring.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/src/istring.h b/src/istring.h
index 5b47b9b..5027981 100644
--- a/src/istring.h
+++ b/src/istring.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/src/kondo.c b/src/kondo.c
index 79b7c56..93318d5 100644
--- a/src/kondo.c
+++ b/src/kondo.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -35,7 +35,7 @@ limitations under the License.
 #define KONDO_SPIN 2
 
 // offsets for indices of A, B, h and t
-// order matters for symbols table
+// order matters for virtual_fields table
 #define KONDO_A_OFFSET 1
 #define KONDO_B_OFFSET 2
 #define KONDO_H_OFFSET 3
@@ -72,6 +72,7 @@ limitations under the License.
 int kondo_generate_conf(Str_Array* str_args, int box_count){
   Str_Array new_args;
   Fields_Table fields;
+  Variables variables;
   Char_Array tmp_str;
   int arg_index;
   int i;
@@ -83,16 +84,19 @@ int kondo_generate_conf(Str_Array* str_args, int box_count){
   kondo_fields_table(box_count, &tmp_str, &fields);
   str_array_append_noinit(tmp_str, &new_args);
 
-  // symbols
-  kondo_symbols(&tmp_str, box_count, &fields);
-  arg_index=find_str_arg("symbols", *str_args);
+  // dummy variables
+  init_Variables(&variables,1);
+
+  // virtual fields
+  kondo_virtual_fields(&tmp_str, box_count, &fields);
+  arg_index=find_str_arg("virtual_fields", *str_args);
   if(arg_index>=0){
     if(tmp_str.length>0){
       char_array_snprintf(&tmp_str,",\n");
     }
     char_array_concat((*str_args).strs[arg_index], &tmp_str);
   }
-  parse_input_symbols(tmp_str, &fields);
+  parse_input_virtual_fields(tmp_str, &fields, variables);
   str_array_append_noinit(tmp_str, &new_args);
 
   // identities
@@ -104,7 +108,7 @@ int kondo_generate_conf(Str_Array* str_args, int box_count){
     }
     char_array_concat((*str_args).strs[arg_index], &tmp_str);
   }
-  parse_input_identities(tmp_str, &fields);
+  parse_input_identities(tmp_str, &fields, variables);
   str_array_append_noinit(tmp_str, &new_args);
 
   // groups
@@ -136,7 +140,7 @@ int kondo_generate_conf(Str_Array* str_args, int box_count){
   // copy remaining entries
   for(i=0;i<(*str_args).length;i++){
     get_str_arg_title((*str_args).strs[i], &title);
-    if(str_cmp(title.str, "symbols")==0 &&\
+    if(str_cmp(title.str, "virtual_fields")==0 &&\
        str_cmp(title.str, "identities")==0 &&\
        str_cmp(title.str, "propagator")==0 &&\
        str_cmp(title.str, "input_polynomial")==0 &&\
@@ -149,6 +153,7 @@ int kondo_generate_conf(Str_Array* str_args, int box_count){
   }
 
   free_Fields_Table(fields);
+  free_Variables(variables);
   free_Str_Array(*str_args);
   *str_args=new_args;
 
@@ -250,15 +255,15 @@ int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fiel
 }
 
 
-// generate Kondo symbols
-int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){
+// generate Kondo virtual_fields
+int kondo_virtual_fields(Char_Array* str_virtual_fields, int box_count, Fields_Table* fields){
   int i,j,k,l;
   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");
+  init_Char_Array(str_virtual_fields, STR_SIZE);
+  char_array_snprintf(str_virtual_fields, "#!virtual_fields\n");
 
   // loop over box index
   for(i=1;i<=box_count;i++){
@@ -267,7 +272,7 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){
       // loop over space dimension
       for(k=0;k<KONDO_DIM;k++){
 	// write index
-	char_array_snprintf(str_symbols, "%d=", 100*(10*(KONDO_A_OFFSET+j)+k)+i);
+	char_array_snprintf(str_virtual_fields, "%d=", 100*(10*(KONDO_A_OFFSET+j)+k)+i);
 	// write the name of the scalar product
 	init_Char_Array(&tmp_str, 6);
 	char_array_snprintf(&tmp_str, "%c%d%d", letters[j], k, i);
@@ -275,10 +280,10 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){
 	kondo_resolve_ABht(tmp_str.str, &poly, *fields);
 	free_Char_Array(tmp_str);
 	// write to output
-	polynomial_sprint(poly, str_symbols);
+	polynomial_sprint(poly, str_virtual_fields);
 	free_Polynomial(poly);
 	// add ,
-	char_array_snprintf(str_symbols,",\n");
+	char_array_snprintf(str_virtual_fields,",\n");
       }
     }
   }
@@ -290,46 +295,46 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){
     for(j=0;j<3;j++){
       for(k=0;k<3;k++){
 	// write index
-	char_array_snprintf(str_symbols, "%d=", 1000*(10*(KONDO_A_OFFSET+j)+KONDO_A_OFFSET+k)+i);
+	char_array_snprintf(str_virtual_fields, "%d=", 1000*(10*(KONDO_A_OFFSET+j)+KONDO_A_OFFSET+k)+i);
 	for(l=0;l<KONDO_DIM;l++){
-	  char_array_snprintf(str_symbols, "(1)");
+	  char_array_snprintf(str_virtual_fields, "(1)");
 	  if(j<2){
-	    char_array_snprintf(str_symbols,"[f%d]", 100*(10*(KONDO_A_OFFSET+j)+l)+i);
+	    char_array_snprintf(str_virtual_fields,"[f%d]", 100*(10*(KONDO_A_OFFSET+j)+l)+i);
 	  }
 	  else{
-	    char_array_snprintf(str_symbols,"[f%d]", 10*(10*(KONDO_A_OFFSET+j)+l));
+	    char_array_snprintf(str_virtual_fields,"[f%d]", 10*(10*(KONDO_A_OFFSET+j)+l));
 	  }
 	  if(k<2){
-	    char_array_snprintf(str_symbols,"[f%d]", 100*(10*(KONDO_A_OFFSET+k)+l)+i);
+	    char_array_snprintf(str_virtual_fields,"[f%d]", 100*(10*(KONDO_A_OFFSET+k)+l)+i);
 	  }
 	  else{
-	    char_array_snprintf(str_symbols,"[f%d]", 10*(10*(KONDO_A_OFFSET+k)+l));
+	    char_array_snprintf(str_virtual_fields,"[f%d]", 10*(10*(KONDO_A_OFFSET+k)+l));
 	  }
 	    
 	  if(l<KONDO_DIM-1){
-	    char_array_append('+',str_symbols);
+	    char_array_append('+',str_virtual_fields);
 	  }
 	}
 	// add ,
-	char_array_snprintf(str_symbols,",\n");
+	char_array_snprintf(str_virtual_fields,",\n");
       }
     }
   }
 
   // vector products
   for(i=1;i<=box_count;i++){
-    char_array_snprintf(str_symbols, "%d=", 100*(100*(KONDO_A_OFFSET)+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+i);
+    char_array_snprintf(str_virtual_fields, "%d=", 100*(100*(KONDO_A_OFFSET)+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+i);
     for(l=0;l<KONDO_DIM;l++){
       // remember (-1 %3 = -1)
-      char_array_snprintf(str_symbols, "(1)[f%d][f%d][f%d]+(-1)[f%d][f%d][f%d]", 100*(10*KONDO_A_OFFSET+((l+1)%KONDO_DIM))+i, 100*(10*KONDO_B_OFFSET+((l+2)%KONDO_DIM))+i, 10*(10*KONDO_H_OFFSET+l), 100*(10*KONDO_A_OFFSET+((l+2)%KONDO_DIM))+i, 100*(10*KONDO_B_OFFSET+((l+1)%KONDO_DIM))+i, 10*(10*KONDO_H_OFFSET+l));
+      char_array_snprintf(str_virtual_fields, "(1)[f%d][f%d][f%d]+(-1)[f%d][f%d][f%d]", 100*(10*KONDO_A_OFFSET+((l+1)%KONDO_DIM))+i, 100*(10*KONDO_B_OFFSET+((l+2)%KONDO_DIM))+i, 10*(10*KONDO_H_OFFSET+l), 100*(10*KONDO_A_OFFSET+((l+2)%KONDO_DIM))+i, 100*(10*KONDO_B_OFFSET+((l+1)%KONDO_DIM))+i, 10*(10*KONDO_H_OFFSET+l));
       if(l<KONDO_DIM-1){
-	char_array_append('+',str_symbols);
+	char_array_append('+',str_virtual_fields);
       }
     }
 
     // add ,
     if(i<box_count){
-      char_array_snprintf(str_symbols,",\n");
+      char_array_snprintf(str_virtual_fields,",\n");
     }
   }
 
@@ -778,13 +783,13 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_
 	  // if polynomial exists, add to each monomial
 	  if(tmp_poly.length>0){
 	    for(i=0;i<tmp_poly.length;i++){
-	      int_array_append(get_symbol_index(buffer), tmp_poly.monomials+i);
+	      int_array_append(get_virtual_field_index(buffer), tmp_poly.monomials+i);
 	    }
 	  }
 	  // if not, create a new term in the polynomial
 	  else{
 	    init_Int_Array(&tmp_monomial, MONOMIAL_SIZE);
-	    int_array_append(get_symbol_index(buffer), &tmp_monomial);
+	    int_array_append(get_virtual_field_index(buffer), &tmp_monomial);
 	    init_Int_Array(&dummy_factor, 1);
 	    polynomial_append_noinit(tmp_monomial, dummy_factor, number_one(), &tmp_poly);
 	  }
@@ -1387,8 +1392,8 @@ int get_offsets_index(char* str, int* offset1, int* offset2, int* index){
   return(0);
 }
 
-// get the index of the symbol corresponding to a given string
-int get_symbol_index(char* str){
+// get the index of the virtual_field corresponding to a given string
+int get_virtual_field_index(char* str){
   char* ptr;
   int offset=-1;
   int index=0;
@@ -1439,7 +1444,7 @@ int get_symbol_index(char* str){
   if(offset==-1){
     return(-1);
   }
-  // no symbol for h or t
+  // no virtual field for h or t
   if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){
     return(10*(10*offset+dim));
   }
diff --git a/src/kondo.h b/src/kondo.h
index c534145..18488ae 100644
--- a/src/kondo.h
+++ b/src/kondo.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -27,10 +27,10 @@ int kondo_generate_conf(Str_Array* str_args, int box_count);
 // generate the Kondo fields table
 int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fields);
 
-// generate Kondo symbols
-int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields);
-// 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);
+// generate Kondo virtual_fields
+int kondo_virtual_fields(Char_Array* str_virtual_fields, int box_count, Fields_Table* fields);
+// generate Kondo virtual_fields (older method: one virtual_field for each scalar product)
+int kondo_virtual_fields_scalarprod(Char_Array* str_virtual_fields, int box_count, Fields_Table* fields);
 
 // generate Kondo groups (groups of independent variables)
 int kondo_groups(Char_Array* str_groups, int box_count);
@@ -71,6 +71,6 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output);
 int get_offset_index(char* str, int* offset, int* index);
 // get the offsets and index of a scalar product
 int get_offsets_index(char* str, int* offset1, int* offset2, int* index);
-// get the index of the symbol corresponding to a given string
-int get_symbol_index(char* str);
+// get the index of the virtual_field corresponding to a given string
+int get_virtual_field_index(char* str);
 #endif
diff --git a/src/kondo_preprocess.c b/src/kondo_preprocess.c
index 3371014..c472e9d 100644
--- a/src/kondo_preprocess.c
+++ b/src/kondo_preprocess.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -104,6 +104,11 @@ for(i=1;i<argc;i++){
     // number of dimensions
     else if (flag==CP_FLAG_DIMENSION){
       sscanf(argv[i],"%d",&((*opts).dimension));
+      // check value of the dimension
+      if((*opts).dimension<=0 || (*opts).dimension>=4){
+	fprintf(stderr,"error: kondo_preprocess only supports dimensions 1, 2 and 3 (got %d)\n",(*opts).dimension);
+	exit(-1);
+      }
       flag=0;
     }
     // read file name from command-line
diff --git a/src/mean.c b/src/mean.c
index 39ece56..0c5707e 100644
--- a/src/mean.c
+++ b/src/mean.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -30,6 +30,7 @@ As of version 1.0, the mean of a monomial is computed directly
 #include "array.h"
 #include "fields.h"
 #include "number.h"
+#include "determinant.h"
 
 // mean of a monomial
 int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Matrix propagator){
@@ -61,10 +62,79 @@ int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Ma
 
 // compute the mean of a monomial of internal fields (with split + and -)
 int mean_internal(Int_Array internal_plus, Int_Array internal_minus, Polynomial* out, Polynomial_Matrix propagator, Fields_Table fields){
+  int ret;
+  Number num;
+
   if(internal_plus.length!=internal_minus.length){
     fprintf(stderr,"error: monomial contains unmatched +/- fields\n");
     exit(-1);
   }
+
+  ret=mean_determinant(internal_plus, internal_minus, &num, propagator, fields);
+  // cannot compute the mean as a determinant, use permutations
+  // can be because some fields are not Fermions
+  // can be because the propagator has non-numeric values (inverting polynomials is not implemented, and would be required for the computation of the determinant)
+  if(ret==-1){
+    mean_permutations(internal_plus, internal_minus, out, propagator, fields);
+  }
+  else{
+    polynomial_multiply_scalar(*out, num);
+    free_Number(num);
+  }
+
+  return(0);
+}
+
+// compute the mean of a monomial by computing a determinant
+// can only be used if all of the propagators are numbers
+int mean_determinant(Int_Array internal_plus, Int_Array internal_minus, Number* out, Polynomial_Matrix propagator, Fields_Table fields){
+  Number_Matrix M;
+  int n=internal_minus.length;
+  int i,j;
+  int a,b;
+  int sign;
+
+  init_Number_Matrix(&M,n);
+  
+  // extra sign: the monomial is sorted in such a way that minus fields are on the left of plus fields, but the determinant formula requires the fields to be alternated +-
+  if((n+1)/2%2==1){
+    sign=-1;
+  }
+  else{
+    sign=1;
+  }
+
+  // construct matrix
+  for(i=0;i<n;i++){
+    a=intlist_find_err(propagator.indices, propagator.length, internal_plus.values[i]);
+    for(j=0;j<n;j++){
+      b=intlist_find_err(propagator.indices, propagator.length, -internal_minus.values[j]);
+      // ignore 0
+      if(propagator.matrix[a][b].length!=0){
+	// check whether the fields are Fermions, and whether the entry is a number
+	if(is_fermion(internal_plus.values[i], fields)==0 || is_fermion(internal_minus.values[j], fields)==0 || polynomial_is_number(propagator.matrix[a][b])==0){
+	  free_Number_Matrix(M);
+	  return(-1);
+	}
+
+	number_add_chain(propagator.matrix[a][b].nums[0], M.matrix[i]+j);
+      }
+    }
+  }
+
+  // compute determinant
+  determinant_inplace(M, out);
+
+  number_Qprod_chain(quot(sign,1), out);
+
+  free_Number_Matrix(M);
+
+  return(0);
+}
+
+
+// compute the mean of a monomial by summing over permutations
+int mean_permutations(Int_Array internal_plus, Int_Array internal_minus, Polynomial* out, Polynomial_Matrix propagator, Fields_Table fields){
   int n=internal_minus.length;
   // pairing as an array of positions
   int* pairing=calloc(n,sizeof(int));
@@ -118,7 +188,7 @@ int mean_internal(Int_Array internal_plus, Int_Array internal_minus, Polynomial*
     pairing_sign=permutation_signature(pairing,n);
   }
 
-  // only simplify in mean_symbols
+  // only simplify in mean_virtual_fields
   polynomial_prod_chain_nosimplify(num_summed,out,fields);
   free_Polynomial(num_summed);
   free(pairing);
@@ -346,10 +416,10 @@ int get_internals(Int_Array monomial, Int_Array* internal_plus, Int_Array* inter
 }
 
 
-// compute the mean of a monomial containing symbolic expressions
+// compute the mean of a monomial containing virtual fields
 // keep track of which means were already computed
-int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed){
-  Int_Array symbol_list;
+int mean_virtual_fields(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed){
+  Int_Array virtual_field_list;
   int i;
   int power;
   int* current_term;
@@ -375,43 +445,43 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po
     }
   }
 
-  init_Int_Array(&symbol_list, monomial.length);
+  init_Int_Array(&virtual_field_list, monomial.length);
   init_Int_Array(&base_monomial, monomial.length);
 
-  // generate symbols list
+  // generate virtual_fields list
   for(i=0;i<monomial.length;i++){
-    if(field_type(monomial.values[i], fields)==FIELD_SYMBOL){
-      int_array_append(intlist_find_err(fields.symbols.indices, fields.symbols.length, monomial.values[i]), &symbol_list);
+    if(field_type(monomial.values[i], fields)==FIELD_VIRTUAL){
+      int_array_append(intlist_find_err(fields.virtual_fields.indices, fields.virtual_fields.length, monomial.values[i]), &virtual_field_list);
     }
     else{
       int_array_append(monomial.values[i], &base_monomial);
     }
   }
-  power=symbol_list.length;
+  power=virtual_field_list.length;
 
   // trivial case
   if(power==0){
     mean(monomial, &mean_num, fields, propagator);
     polynomial_concat_noinit(mean_num, output);
     
-    free_Int_Array(symbol_list);
+    free_Int_Array(virtual_field_list);
     free_Int_Array(base_monomial);
     return(0);
   }
   else{
     // initialize current term to a position that has no repetitions
     current_term=calloc(power,sizeof(int));
-    exists_next=init_prod(current_term, symbol_list, fields, power, base_monomial)+1;
+    exists_next=init_prod(current_term, virtual_field_list, fields, power, base_monomial)+1;
   }
 
-  // loop over terms; the loop stops when all the pointers are at the end of the first symbol
+  // loop over terms; the loop stops when all the pointers are at the end of the first virtual field
   while(exists_next==1){
     // construct monomial
     int_array_cpy(base_monomial, &tmp_monomial);
     tmp_num=number_one();
     for(i=0;i<power;i++){
-      int_array_concat(fields.symbols.expr[symbol_list.values[i]].monomials[current_term[i]], &tmp_monomial);
-      number_prod_chain(fields.symbols.expr[symbol_list.values[i]].nums[current_term[i]], &tmp_num);
+      int_array_concat(fields.virtual_fields.expr[virtual_field_list.values[i]].monomials[current_term[i]], &tmp_monomial);
+      number_prod_chain(fields.virtual_fields.expr[virtual_field_list.values[i]].nums[current_term[i]], &tmp_num);
     }
     // check whether the monomial vanishes
     if(check_monomial_match(tmp_monomial, fields)==1){
@@ -432,7 +502,7 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po
     free_Int_Array(tmp_monomial);
 
     // next term
-    exists_next=next_prod(current_term, symbol_list, fields, power, base_monomial)+1;
+    exists_next=next_prod(current_term, virtual_field_list, fields, power, base_monomial)+1;
 
     
     // simplfiy every 25 steps (improves both memory usage and performance)
@@ -451,13 +521,13 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po
 
   // free memory
   free(current_term);
-  free_Int_Array(symbol_list);
+  free_Int_Array(virtual_field_list);
   free_Int_Array(base_monomial);
   return(0);
 }
 
 // first term in product with no repetitions
-int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial){
+int init_prod(int* current_term, Int_Array virtual_field_list, Fields_Table fields, int power, Int_Array base_monomial){
   // index we want to increment
   int move=0;
   // tmp monomial
@@ -474,7 +544,7 @@ int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int
   // loop until move is out of range
   while(move>=0 && move<power){
     // move
-    current_term[move]=next_term_norepeat(current_term[move], fields.symbols.expr[symbol_list.values[move]], &monomial, fields);
+    current_term[move]=next_term_norepeat(current_term[move], fields.virtual_fields.expr[virtual_field_list.values[move]], &monomial, fields);
     // if the next term does not exist, then move previous index
     if(current_term[move]==-1){
       move--;
@@ -495,7 +565,7 @@ int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int
 }
 
 // next term in product with no repetitions
-int next_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial){
+int next_prod(int* current_term, Int_Array virtual_field_list, Fields_Table fields, int power, Int_Array base_monomial){
   // index we want to increment
   int move=power-1;
   // tmp monomial
@@ -506,13 +576,13 @@ int next_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int
   init_Int_Array(&monomial, base_monomial.length+5*power);
   int_array_cpy_noinit(base_monomial, &monomial);
   for(i=0;i<=move;i++){
-    int_array_concat(fields.symbols.expr[symbol_list.values[i]].monomials[current_term[i]],&monomial);
+    int_array_concat(fields.virtual_fields.expr[virtual_field_list.values[i]].monomials[current_term[i]],&monomial);
   }
 
   // loop until move is out of range
   while(move>=0 && move<power){
     // move
-    current_term[move]=next_term_norepeat(current_term[move], fields.symbols.expr[symbol_list.values[move]], &monomial, fields);
+    current_term[move]=next_term_norepeat(current_term[move], fields.virtual_fields.expr[virtual_field_list.values[move]], &monomial, fields);
     // if the next term does not exist, then move previous index
     if(current_term[move]==-1){
       move--;
@@ -623,13 +693,13 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol
   int group=-2;
   int next_group=-2;
   Polynomial tmp_poly;
-  int sign;
+  int sign=1;
 
   init_Polynomial(output, MONOMIAL_SIZE);
 
-  // check whether there are symbols
-  // 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){
+  // check whether there are virtual fields
+  // IMPORTANT: the virtual fields must be at the end of the monomial
+  if(monomial.length==0 || field_type(monomial.values[monomial.length-1], fields)!=FIELD_VIRTUAL){
     // mean
     mean(monomial, &num_mean, fields, propagator);
     // add to output
@@ -650,7 +720,7 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol
       }
       // if group changes, take mean
       if((i>0 && next_group!=group) || i==monomial.length){
-	mean_symbols(tmp_monomial, &tmp_poly, fields, propagator, groups, computed);
+	mean_virtual_fields(tmp_monomial, &tmp_poly, fields, propagator, groups, computed);
 	// if zero
 	if(polynomial_is_zero(tmp_poly)==1){
 	  // set output to 0
@@ -702,10 +772,11 @@ struct mean_args{
   Fields_Table fields;
   Polynomial_Matrix propagator;
   Groups groups;
+  int print_progress;
 };
 
 // multithreaded
-int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads){
+int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, int print_progress){
   int i;
   Polynomial thread_polys[threads];
   pthread_t thread_ids[threads];
@@ -720,11 +791,15 @@ int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Pol
     args[i].fields=fields;
     args[i].propagator=propagator;
     args[i].groups=groups;
+    args[i].print_progress=print_progress;
   }
 
   // split polynomial
+  // randomly choose the thread
+  // see random number generator
+  srand(time(NULL));
   for(i=0;i<len;i++){
-    polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], thread_polys+(i % threads));
+    polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], thread_polys+(rand() % threads));
   }
 
   // start threads
@@ -750,12 +825,12 @@ int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Pol
 // mean for one of the threads
 void* polynomial_mean_thread(void* mean_args){
   struct mean_args *args=mean_args;
-  polynomial_mean((*args).polynomial,(*args).fields,(*args).propagator,(*args).groups);
+  polynomial_mean((*args).polynomial,(*args).fields,(*args).propagator,(*args).groups, (*args).print_progress);
   return(NULL);
 }
 
 // single threaded version
-int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups){
+int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int print_progress){
   int i,j;
   Polynomial output;
   Polynomial tmp_poly;
@@ -769,7 +844,9 @@ int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matr
 
   // mean of each monomial
   for(i=0;i<(*polynomial).length;i++){
-    fprintf(stderr,"computing %d of %d means\n",i,(*polynomial).length-1);
+    if(print_progress==1){
+      fprintf(stderr,"computing %d of %d means\n",i,(*polynomial).length-1);
+    }
     mean_groups((*polynomial).monomials[i], &tmp_poly, fields, propagator, groups, &computed);
 
     // write factors
diff --git a/src/mean.h b/src/mean.h
index 34ec1d3..2726992 100644
--- a/src/mean.h
+++ b/src/mean.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -28,6 +28,12 @@ int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Ma
 
 // compute the mean of a monomial of internal fields (with split + and -)
 int mean_internal(Int_Array internal_plus, Int_Array internal_minus, Polynomial* out, Polynomial_Matrix propagator, Fields_Table fields);
+
+// compute the mean of a monomial by computing a determinant
+int mean_determinant(Int_Array internal_plus, Int_Array internal_minus, Number* out, Polynomial_Matrix propagator,Fields_Table fields);
+
+// compute the mean of a monomial by summing over permutations
+int mean_permutations(Int_Array internal_plus, Int_Array internal_minus, Polynomial* out, Polynomial_Matrix propagator, Fields_Table fields);
 // first pairing with a non-vanishing propagator
 int init_pairing(int* pairing, int* mask, int n, Polynomial_Matrix propagator, int* indices_plus, int* indices_minus);
 // next pairing with a non-vanishing propagator
@@ -43,12 +49,12 @@ int mean_internal_slow(Int_Array internal_plus, Int_Array internal_minus, Number
 // requires the monomial to be sorted (for the sign to be correct)
 int get_internals(Int_Array monomial, Int_Array* internal_plus, Int_Array* internal_minus, Int_Array* others, Fields_Table fields);
 
-// compute the mean of a monomial containing symbolic expressions
-int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed);
+// compute the mean of a monomial containing virtual fields
+int mean_virtual_fields(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed);
 // first term in product with no repetitions
-int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial);
+int init_prod(int* current_term, Int_Array virtual_field_list, Fields_Table fields, int power, Int_Array base_monomial);
 // next term in product with no repetitions
-int next_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial);
+int next_prod(int* current_term, Int_Array virtual_field_list, Fields_Table fields, int power, Int_Array base_monomial);
 // find the next term in a polynomial that can be multiplied to a given monomial and add it to the monomial
 int next_term_norepeat(int start, Polynomial polynomial, Int_Array* monomial, Fields_Table fields);
 
@@ -62,9 +68,9 @@ int sort_fermions(int* array, int begin, int end, int* sign);
 int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed);
 
 // compute the mean of a polynomial
-int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups);
+int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int print_progress);
 // multithreaded
-int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads);
+int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, int print_progress);
 // single thread mean
 void* polynomial_mean_thread(void* mean_args);
 #endif
diff --git a/src/meankondo.c b/src/meankondo.c
index 1c976d4..8b5c197 100644
--- a/src/meankondo.c
+++ b/src/meankondo.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -17,7 +17,7 @@ limitations under the License.
 /*
 meankondo
 
-A simple tool to compute the renormalization group flow for Fermionic hierarchical models
+A tool to compute the renormalization group flow for Fermionic hierarchical models
 
 */
 
@@ -49,15 +49,19 @@ A simple tool to compute the renormalization group flow for Fermionic hierarchic
 #include "mean.h"
 // various string operations
 #include "istring.h"
+// symbolic trees
+# include "tree.h"
 
 // read cli arguments
 int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meankondo_Options* opts);
 // print usage message
 int print_usage_meankondo();
+// check consistency of options
+int check_meankondo_opts(Meankondo_Options opts);
 // compute flow
 int compute_flow(Str_Array str_args, Meankondo_Options opts);
-// compute the flow equation
-int compute_flow_equation(Polynomial init_poly, Id_Table idtable, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, Grouped_Polynomial* flow_equation);
+// compute average
+int compute_average(Polynomial init_poly, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, int print_progress, Polynomial* exp_poly);
 
 
 int main (int argc, const char* argv[]){
@@ -69,6 +73,9 @@ int main (int argc, const char* argv[]){
   // read command-line arguments
   read_args_meankondo(argc,argv,&str_args,&opts);
 
+  // check command-line arguments
+  check_meankondo_opts(opts);
+
   // warning message if representing rational numbers as floats
 #ifdef RATIONAL_AS_FLOAT
   fprintf(stderr,"info: representing rational numbers using floats\n");
@@ -101,6 +108,10 @@ int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meanko
   (*opts).threads=1;
   // do not chain
   (*opts).chain=0;
+  // do not print progress
+  (*opts).print_progress=0;
+  // print the flow equation
+  (*opts).group_poly=1;
 
   // loop over arguments
   for(i=1;i<argc;i++){
@@ -116,6 +127,13 @@ int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meanko
 	case 'C':
 	  (*opts).chain=1;
 	  break;
+	// print progress
+	case 'p':
+	  (*opts).print_progress=1;
+	  break;
+	case 'A':
+	  (*opts).group_poly=0;
+	  break;
 	// print version
 	case 'v':
 	  printf("meankondo " VERSION "\n");
@@ -147,7 +165,16 @@ int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meanko
 
 // print usage message
 int print_usage_meankondo(){
-  printf("\nusage:\n   meankondo [-t threads] [-C] <filename>\n\n");
+  printf("\nusage:\n   meankondo [-t threads] [-C] [-p] [-A] <filename>\n\n");
+  return(0);
+}
+
+// check consistency of options
+int check_meankondo_opts(Meankondo_Options opts){
+  if(opts.chain==1 && opts.group_poly==0){
+    fprintf(stderr,"aborting: the '-C' and '-A' options are incompatible\n");
+    exit(-1);
+  }
   return(0);
 }
 
@@ -163,6 +190,8 @@ int compute_flow(Str_Array str_args, Meankondo_Options opts){
   Fields_Table fields;
   // their propagator
   Polynomial_Matrix propagator;
+  // preprocessor variables
+  Variables variables;
   // initial polynomial
   Polynomial init_poly;
   // list of rccs
@@ -171,6 +200,8 @@ int compute_flow(Str_Array str_args, Meankondo_Options opts){
   Groups groups;
   // flow equation
   Grouped_Polynomial flow_equation;
+  // polynomial produced by the averaging operation
+  Polynomial exp_poly;
 
 
   // parse fields
@@ -183,29 +214,38 @@ int compute_flow(Str_Array str_args, Meankondo_Options opts){
     parse_input_fields(str_args.strs[arg_index],&fields);
   }
 
-  // parse id table
-  arg_index=find_str_arg("id_table", str_args);
-  if(arg_index<0){
-    fprintf(stderr,"error: no id table entry in the configuration file\n");
-    exit(-1);
+  // parse variables
+  // must precede id_table, virtual_fields, identities and input_polynomial
+  arg_index=find_str_arg("preprocessor_variables", str_args);
+  if(arg_index>=0){
+    parse_input_variables(str_args.strs[arg_index],&variables);
   }
   else{
-    parse_input_id_table(str_args.strs[arg_index],&idtable, fields);
+    init_Variables(&variables,1);
   }
 
-  // parse symbols
-  arg_index=find_str_arg("symbols", str_args);
-  if(arg_index>=0){
-    parse_input_symbols(str_args.strs[arg_index],&fields);
+  // parse id table
+  if(opts.group_poly==1){
+    arg_index=find_str_arg("id_table", str_args);
+    if(arg_index<0){
+      fprintf(stderr,"error: no id table entry in the configuration file\n");
+      exit(-1);
+    }
+    else{
+      parse_input_id_table(str_args.strs[arg_index],&idtable, fields, variables);
+    }
   }
-  else{
-    init_Symbols(&(fields.symbols),1);
+
+  // parse virtual_fields
+  arg_index=find_str_arg("virtual_fields", str_args);
+  if(arg_index>=0){
+    parse_input_virtual_fields(str_args.strs[arg_index], &fields, variables);
   }
 
   // parse input polynomial
   arg_index=find_str_arg("input_polynomial", str_args);
   if(arg_index>=0){
-    parse_input_polynomial(str_args.strs[arg_index],&init_poly, fields);
+    parse_input_polynomial(str_args.strs[arg_index],&init_poly, fields, variables);
   }
   else{
     fprintf(stderr,"error: no input polynomial entry in the configuration file\n");
@@ -225,41 +265,90 @@ int compute_flow(Str_Array str_args, Meankondo_Options opts){
   // parse identities
   arg_index=find_str_arg("identities", str_args);
   if(arg_index>=0){
-    parse_input_identities(str_args.strs[arg_index],&fields);
-  }
-  else{
-    init_Identities(&(fields.ids),1);
+    parse_input_identities(str_args.strs[arg_index],&fields, variables);
   }
 
-  // parse groups
+  // parse groups (must come after virtual_fields and propagator)
   arg_index=find_str_arg("groups", str_args);
   if(arg_index>=0){
-    parse_input_groups(str_args.strs[arg_index],&groups);
+    parse_input_groups(str_args.strs[arg_index],&groups, propagator, fields);
   }
   else{
     init_Groups(&groups, 1);
   }
 
-  // flow equation
-  compute_flow_equation(init_poly, idtable, fields, propagator, groups, opts.threads, &flow_equation);
+  // compute the average
+  compute_average(init_poly, fields, propagator, groups, opts.threads, opts.print_progress, &exp_poly);
+
   free_Polynomial(init_poly);
   free_Polynomial_Matrix(propagator);
-  free_Fields_Table(fields);
   free_Groups(groups);
 
+  // parse postprocessing entry
+  arg_index=find_str_arg("postprocess_operation", str_args);
+  if(arg_index>=0){
+    add_polynomial_to_variables("OUT", exp_poly, &variables);
+    // parse postprocess entry
+    Polynomial postprocess_operation;
+    parse_input_polynomial(str_args.strs[arg_index], &postprocess_operation, fields, variables);
+    // replace exp_poly
+    free_Polynomial(exp_poly);
+    exp_poly=postprocess_operation;
+  }
+
+  if(opts.group_poly==1){
+    // flow equation
+    group_polynomial(exp_poly, &flow_equation, idtable, fields);
+  }
+
+  // postprocess flow equation
+  arg_index=find_str_arg("postprocess_flow_equation", str_args);
+  if(arg_index>=0){
+    Polynomial flow_polynomial;
+    // polynomial made of the rcc's multiplied by the corresponding fields (parsed from idtable)
+    idtable_to_polynomial(idtable, &flow_polynomial);
+
+    // add to variables
+    add_polynomial_to_variables("FLOW", flow_polynomial, &variables);
+    free_Polynomial(flow_polynomial);
+
+    // parse postprocess entry
+    Polynomial postprocess_polynomial;
+    parse_input_polynomial(str_args.strs[arg_index], &postprocess_polynomial, fields, variables);
+
+    // convert to flow equation
+    Grouped_Polynomial postprocess_flow_equation;
+    group_polynomial(postprocess_polynomial, &postprocess_flow_equation, idtable, fields);
+    free_Polynomial(postprocess_polynomial);
+
+    // apply postprocessing to flow equation
+    Grouped_Polynomial new_flow;
+    compose_flow_equations(postprocess_flow_equation, flow_equation, &new_flow);
+    free_Grouped_Polynomial(postprocess_flow_equation);
+
+    // replace flow_equation
+    free_Grouped_Polynomial(flow_equation);
+    flow_equation=new_flow;
+  }
+
+
   // if chain then print config file
   if(opts.chain==1){
     for(i=0;i<str_args.length;i++){
       // check whether to print the str_arg
       get_str_arg_title(str_args.strs[i], &arg_header);
       if (\
-	str_cmp(arg_header.str, "symbols")==0 &&\
+	str_cmp(arg_header.str, "variables")==0 &&\
+	str_cmp(arg_header.str, "virtual_fields")==0 &&\
 	str_cmp(arg_header.str, "groups")==0 &&\
 	str_cmp(arg_header.str, "fields")==0 &&\
 	str_cmp(arg_header.str, "identities")==0 &&\
 	str_cmp(arg_header.str, "propagator")==0 &&\
 	str_cmp(arg_header.str, "input_polynomial")==0 &&\
-	str_cmp(arg_header.str, "id_table")==0 ){
+	str_cmp(arg_header.str, "id_table")==0 &&\
+	str_cmp(arg_header.str, "postprocess_operation")==0 &&\
+	str_cmp(arg_header.str, "numerical_postprocess_operation")==0
+	){
 	  printf("%s\n&\n",str_args.strs[i].str);
       }
       free_Char_Array(arg_header);
@@ -268,34 +357,66 @@ int compute_flow(Str_Array str_args, Meankondo_Options opts){
     printf("#!flow_equation\n");
   }
 
-  // print flow equation
-  grouped_polynomial_print(flow_equation,'%','%');
+  // print result
+  if(opts.group_poly==1){
+    grouped_polynomial_print(flow_equation,'%','%');
+
+    // free memory
+    free_Grouped_Polynomial(flow_equation);
+  }
+  else{
+    polynomial_print(exp_poly);
+  }
 
   // free memory
-  free_Id_Table(idtable);
-  free_Grouped_Polynomial(flow_equation);
+  free_Polynomial(exp_poly);
+
+  // parse numerical_postprocessing entry
+  arg_index=find_str_arg("numerical_postprocess_operation", str_args);
+  if(arg_index>=0){
+    Polynomial rcc_polynomial;
+    // polynomial made of the rcc's multiplied by the corresponding fields (parsed from idtable)
+    idtable_to_polynomial(idtable, &rcc_polynomial);
+
+    // add to variables
+    add_polynomial_to_variables("RCC", rcc_polynomial, &variables);
+    free_Polynomial(rcc_polynomial);
+
+    // parse postprocess entry
+    Polynomial numerical_postprocess_operation;
+    parse_input_polynomial(str_args.strs[arg_index], &numerical_postprocess_operation, fields, variables);
+
+    // convert to flow equation
+    Grouped_Polynomial numerical_postprocess_flow_equation;
+    group_polynomial(numerical_postprocess_operation, &numerical_postprocess_flow_equation, idtable, fields);
+    free_Polynomial(numerical_postprocess_operation);
+
+    // print postprocessing flow equation
+    printf("\n&\n#!postprocess_operation\n");
+    grouped_polynomial_print(numerical_postprocess_flow_equation,'%','%');
+    free_Grouped_Polynomial(numerical_postprocess_flow_equation);
+  }
+
+  if(opts.group_poly==1){
+    free_Id_Table(idtable);
+  }
+  free_Fields_Table(fields);
+  free_Variables(variables);
+
   return(0);
 }
 
 
-// compute the flow equation
-int compute_flow_equation(Polynomial init_poly, Id_Table idtable, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, Grouped_Polynomial* flow_equation){
-  // expectation
-  Polynomial exp_poly;
-
-  polynomial_cpy(init_poly,&exp_poly);
+// compute average
+int compute_average(Polynomial init_poly, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, int print_progress, Polynomial* exp_poly){
+  polynomial_cpy(init_poly,exp_poly);
 
   // average
   if(threads>1){
-    polynomial_mean_multithread(&exp_poly, fields, propagator, groups, threads);
+    polynomial_mean_multithread(exp_poly, fields, propagator, groups, threads, print_progress);
   }
   else{
-    polynomial_mean(&exp_poly, fields, propagator, groups);
+    polynomial_mean(exp_poly, fields, propagator, groups, print_progress);
   }
-    
-  // grouped representation of expanded_poly
-  group_polynomial(exp_poly,flow_equation,idtable, fields);
-  free_Polynomial(exp_poly);
-
   return(0);
 }
diff --git a/src/meantools.c b/src/meantools.c
index e3e7247..64cf3f1 100644
--- a/src/meantools.c
+++ b/src/meantools.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -39,13 +39,13 @@ Utility to perform various operations on flow equations
 // string functions
 #include "istring.h"
 // tools
-#include "meantools_exp.h"
 #include "meantools_deriv.h"
 #include "meantools_eval.h"
+#include "meantools_expand.h"
 
-#define EXP_COMMAND 1
-#define DERIV_COMMAND 2
-#define EVAL_COMMAND 3
+#define DERIV_COMMAND 1
+#define EVAL_COMMAND 2
+#define EXPAND_COMMAND 3
 
 // read cli arguments
 int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meantools_Options* opts);
@@ -63,12 +63,12 @@ int main (int argc, const char* argv[]){
   read_args_meantools(argc,argv,&str_args, &opts);
 
   switch(opts.command){
-  case EXP_COMMAND: tool_exp(str_args);
-    break;
   case DERIV_COMMAND: tool_deriv(str_args,opts);
     break;
   case EVAL_COMMAND: tool_eval(str_args,opts);
     break;
+  case EXPAND_COMMAND: tool_expand(str_args,opts);
+    break;
   }
 
   //free memory
@@ -86,11 +86,7 @@ int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meanto
     exit(-1);
   }
 
-  if(str_cmp((char*)argv[1],"exp")==1){
-    (*opts).command=EXP_COMMAND;
-    tool_exp_read_args(argc, argv, str_args);
-  }
-  else if(str_cmp((char*)argv[1],"derive")==1){
+  if(str_cmp((char*)argv[1],"differentiate")==1){
     (*opts).command=DERIV_COMMAND;
     tool_deriv_read_args(argc, argv, str_args, opts);
   }
@@ -98,6 +94,10 @@ int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meanto
     (*opts).command=EVAL_COMMAND;
     tool_eval_read_args(argc, argv, str_args, opts);
   }
+  else if(str_cmp((char*)argv[1],"expand")==1){
+    (*opts).command=EXPAND_COMMAND;
+    tool_expand_read_args(argc, argv, str_args, opts);
+  }
   else{
     print_usage_meantools();
     exit(-1);
@@ -108,9 +108,6 @@ int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meanto
 
 // print usage message
 int print_usage_meantools(){
-  printf("\nusage:\n   meantools exp [config_file]\n   meantools derive [-d derivatives] [-V variables] [-C] [config_file]\n   meantools eval [-R values] [-P precision] [-E max_exponent] [config_file]\n\n");
+  printf("\nusage:\n   meantools differentiate [-d derivatives] [-V variables] [-C] [config_file]\n   meantools eval [-R values] [-P precision] [-E max_exponent] [config_file]\n   meantools expand [-N namespace] [config_file]\n\n");
   return(0);
 }
-
-
-
diff --git a/src/meantools_deriv.c b/src/meantools_deriv.c
index b096380..3530808 100644
--- a/src/meantools_deriv.c
+++ b/src/meantools_deriv.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -42,9 +42,9 @@ int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Mean
   char* ptr;
 
   // defaults
-  // derive once
+  // differentiate once
   (*opts).deriv_derivs=1;
-  // derive with respect to all variables
+  // differentiate with respect to all variables
   (*opts).deriv_vars.length=-1;
   // do not chain
   (*opts).chain=0;
@@ -77,7 +77,7 @@ int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Mean
     }
     // variables
     else if(flag==CP_FLAG_VARS){
-      // if the argument is "all" then derive wrt all variables
+      // if the argument is "all" then differentiate wrt all variables
       if(str_cmp((char*)argv[i],"all")){
 	(*opts).deriv_vars.length=-2;
       }
@@ -101,7 +101,7 @@ int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Mean
 }
 
 
-// derive a flow equation
+// differentiate a flow equation
 int tool_deriv(Str_Array str_args, Meantools_Options opts){
   // index of the entry in the input file
   int arg_index;
@@ -158,7 +158,7 @@ int tool_deriv(Str_Array str_args, Meantools_Options opts){
       get_str_arg_title(str_args.strs[i], &arg_header);
       if (\
 	str_cmp(arg_header.str, "flow_equation")==0 &&\
-	str_cmp(arg_header.str, "symbols")==0 &&\
+	str_cmp(arg_header.str, "virtual_fields")==0 &&\
 	str_cmp(arg_header.str, "groups")==0 &&\
 	str_cmp(arg_header.str, "fields")==0 &&\
 	str_cmp(arg_header.str, "identities")==0 &&\
@@ -208,7 +208,7 @@ int flow_equation_derivative(int n, Int_Array variables, Grouped_Polynomial flow
     // free tmpflow
     free_Grouped_Polynomial(tmpflow);
 
-    // add the derived indices as variables for the next derivative
+    // add the differentiated indices as variables for the next derivative
     for(i=0;i<variables.length;i++){
       if(variables.values[i]>=0){
 	int_array_append((j+1)*DOFFSET+variables.values[i], &indices);
diff --git a/src/meantools_deriv.h b/src/meantools_deriv.h
index 4ea36a9..7425fd1 100644
--- a/src/meantools_deriv.h
+++ b/src/meantools_deriv.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -21,7 +21,7 @@ limitations under the License.
 
 // read arguments
 int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts);
-// derive a flow equation
+// differentiate a flow equation
 int tool_deriv(Str_Array str_args, Meantools_Options opts);
 // n first derivatives of a flow equation wrt to variables
 int flow_equation_derivative(int n, Int_Array variables, Grouped_Polynomial flow_equation, Grouped_Polynomial* flow_equation_derivs);
diff --git a/src/meantools_eval.c b/src/meantools_eval.c
index 4f166cd..027d18d 100644
--- a/src/meantools_eval.c
+++ b/src/meantools_eval.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -163,12 +163,12 @@ int tool_eval(Str_Array str_args, Meantools_Options opts){
 
   // evaluate
   if(mpfr_flag==0){
-    evaleq(&rccs, flow_equation);
+    evaleq(rccs, rccs, flow_equation);
     RCC_print(rccs);
     free_RCC(rccs);
   }
   else{
-    evaleq_mpfr(&rccs_mpfr, flow_equation);
+    evaleq_mpfr(rccs_mpfr, rccs_mpfr, flow_equation);
     RCC_mpfr_print(rccs_mpfr);
     free_RCC_mpfr(rccs_mpfr);
   }
diff --git a/src/meantools_eval.h b/src/meantools_eval.h
index 518049d..54f6e51 100644
--- a/src/meantools_eval.h
+++ b/src/meantools_eval.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/src/meantools_exp.c b/src/meantools_exp.c
deleted file mode 100644
index 1aca928..0000000
--- a/src/meantools_exp.c
+++ /dev/null
@@ -1,130 +0,0 @@
-/*
-Copyright 2015 Ian Jauslin
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-
-    http://www.apache.org/licenses/LICENSE-2.0
-
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
-
-#include "meantools_exp.h"
-
-#include <stdio.h>
-#include <stdlib.h>
-#include "parse_file.h"
-#include "cli_parser.h"
-#include "polynomial.h"
-#include "fields.h"
-#include "grouped_polynomial.h"
-#include "idtable.h"
-
-// read command line arguments
-int tool_exp_read_args(int argc, const char* argv[], Str_Array* str_args){
-  // file to read the polynomial from in flow mode
-  const char* file="";
-  // whether a file was specified on the command-line
-  int exists_file=0;
-
-  if(argc>=3){
-    file=argv[2];
-    exists_file=1;
-  }
-  read_config_file(str_args, file, 1-exists_file);
-
-  return(0);
-}
-
-
-// compute the exponential of the input polynomial
-int tool_exp(Str_Array str_args){
-  // index of the entry in the input file
-  int arg_index;
-  // list of fields
-  Fields_Table fields;
-  // input polynomial
-  Polynomial poly;
-  // exp as a polynomial
-  Polynomial exp_poly;
-  // list of rccs
-  Id_Table idtable;
-  // exp
-  Grouped_Polynomial exp;
-  int i,j;
-
-  // parse fields
-  arg_index=find_str_arg("fields", str_args);
-  if(arg_index<0){
-    fprintf(stderr,"error: no fields entry in the configuration file\n");
-    exit(-1);
-  }
-  else{
-    parse_input_fields(str_args.strs[arg_index],&fields);
-  }
-
-  // parse id table
-  arg_index=find_str_arg("id_table", str_args);
-  if(arg_index<0){
-    fprintf(stderr,"error: no id table entry in the configuration file\n");
-    exit(-1);
-  }
-  else{
-    parse_input_id_table(str_args.strs[arg_index],&idtable, fields);
-  }
-
-  // parse input polynomial
-  arg_index=find_str_arg("input_polynomial", str_args);
-  if(arg_index>=0){
-    parse_input_polynomial(str_args.strs[arg_index],&poly, fields);
-  }
-  else{
-    fprintf(stderr,"error: no input polynomial entry in the configuration file\n");
-    exit(-1);
-  }
-
-  // parse symbols
-  arg_index=find_str_arg("symbols", str_args);
-  if(arg_index>=0){
-    parse_input_symbols(str_args.strs[arg_index],&fields);
-  }
-  else{
-    init_Symbols(&(fields.symbols),1);
-  }
-
-  // parse identities
-  arg_index=find_str_arg("identities", str_args);
-  if(arg_index>=0){
-    parse_input_identities(str_args.strs[arg_index],&fields);
-  }
-  else{
-    init_Identities(&(fields.ids),1);
-  }
-
-  // exp(V)
-  polynomial_exponential(poly,&exp_poly, fields);
-  // grouped representation
-  group_polynomial(exp_poly, &exp, idtable, fields);
-  free_Polynomial(exp_poly);
-  free_Polynomial(poly);
-
-  // no denominators
-  for(i=0;i<exp.length;i++){
-    for(j=0;j<exp.coefs[i].length;j++){
-      exp.coefs[i].denoms[j].power=0;
-    }
-  }
-
-  grouped_polynomial_print(exp,'%','%');
-
-  // free memory
-  free_Fields_Table(fields);
-  free_Id_Table(idtable);
-  free_Grouped_Polynomial(exp);
-  return(0);
-}
diff --git a/src/meantools_expand.c b/src/meantools_expand.c
new file mode 100644
index 0000000..2586ff5
--- /dev/null
+++ b/src/meantools_expand.c
@@ -0,0 +1,166 @@
+/*
+Copyright 2015-2022 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "meantools_expand.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "cli_parser.h"
+#include "parse_file.h"
+#include "polynomial.h"
+#include "array.h"
+#include "fields.h"
+
+
+#define CP_FLAG_NAMESPACE 1
+// read command line arguments
+int tool_expand_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts){
+  // file to read the polynomial from in flow mode
+  const char* file="";
+  // whether a file was specified on the command-line
+  int exists_file=0;
+  // flag
+  int flag=0;
+  int i;
+  char* ptr;
+
+  // defaults
+  // no namespace
+  (*opts).namespace.length=-1;
+
+
+  // loop over arguments
+  for(i=2;i<argc;i++){
+    // flag
+    if(argv[i][0]=='-'){
+      for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
+	switch(*ptr){
+	// number of derivatives
+	case 'N':
+	  flag=CP_FLAG_NAMESPACE;
+	  break;
+	}
+      }
+    }
+    // namespace
+    else if(flag==CP_FLAG_NAMESPACE){
+      str_to_char_array((char*)argv[i], &(opts->namespace));
+      flag=0;
+    }
+    // read file name from command-line
+    else{
+      file=argv[i];
+      exists_file=1;
+    }
+  }
+
+  read_config_file(str_args, file, 1-exists_file);
+
+  return(0);
+}
+
+
+// expand a polynomial
+int tool_expand(Str_Array str_args, Meantools_Options opts){
+  // index of the entry in the input file
+  int arg_index;
+  // input polynomial
+  Polynomial polynomial;
+  // fields table
+  Fields_Table fields;
+  // preprocessor variables
+  Variables variables;
+
+  // parse fields
+  if(opts.namespace.length>=0){
+    arg_index=find_str_arg_ns("fields", opts.namespace, str_args);
+  }
+  else{
+    arg_index=find_str_arg("fields", str_args);
+  }
+  if(arg_index<0){
+    fprintf(stderr,"error: no fields entry in the configuration file\n");
+    exit(-1);
+  }
+  else{
+    parse_input_fields(str_args.strs[arg_index],&fields);
+  }
+
+  // parse variables
+  // must precede id_table, virtual_fields, identities and input_polynomial
+  if(opts.namespace.length>=0){
+    arg_index=find_str_arg_ns("preprocessor_variables", opts.namespace, str_args);
+  }
+  else{
+    arg_index=find_str_arg("preprocessor_variables", str_args);
+  }
+  if(arg_index>=0){
+    parse_input_variables(str_args.strs[arg_index],&variables);
+  }
+  else{
+    init_Variables(&variables,1);
+  }
+
+  // parse virtual_fields
+  if(opts.namespace.length>=0){
+    arg_index=find_str_arg_ns("virtual_fields", opts.namespace, str_args);
+  }
+  else{
+    arg_index=find_str_arg("virtual_fields", str_args);
+  }
+  if(arg_index>=0){
+    parse_input_virtual_fields(str_args.strs[arg_index], &fields, variables);
+  }
+
+  // parse identities
+  if(opts.namespace.length>=0){
+    arg_index=find_str_arg_ns("identities", opts.namespace, str_args);
+  }
+  else{
+    arg_index=find_str_arg("identities", str_args);
+  }
+  if(arg_index>=0){
+    parse_input_identities(str_args.strs[arg_index],&fields, variables);
+  }
+
+  // parse input polynomial
+  if(opts.namespace.length>=0){
+    arg_index=find_str_arg_ns("input_polynomial", opts.namespace, str_args);
+  }
+  else{
+    arg_index=find_str_arg("input_polynomial", str_args);
+  }
+  if(arg_index>=0){
+    parse_input_polynomial(str_args.strs[arg_index], &polynomial, fields, variables);
+  }
+  else{
+    fprintf(stderr,"error: no input polynomial entry in the configuration file\n");
+    exit(-1);
+  }
+
+  // print polynomial
+  polynomial_print(polynomial);
+
+  // free memory
+  free_Variables(variables);
+  free_Fields_Table(fields);
+  free_Polynomial(polynomial);
+
+  if(opts.namespace.length>=0){
+    free_Char_Array(opts.namespace);
+  }
+  return(0);
+}
diff --git a/src/meantools_exp.h b/src/meantools_expand.h
index f3f6666..9d988b1 100644
--- a/src/meantools_exp.h
+++ b/src/meantools_expand.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -14,14 +14,16 @@ See the License for the specific language governing permissions and
 limitations under the License.
 */
 
-#ifndef MEANTOOLS_EXP_H
-#define MEANTOOLS_EXP_H
+#ifndef MEANTOOLS_EXPAND_H
+#define MEANTOOLS_EXPAND_H
 
 #include "types.h"
 
 // read arguments
-int tool_exp_read_args(int argc, const char* argv[], Str_Array* str_args);
-// compute the exponential of the input polynomial
-int tool_exp(Str_Array str_args);
+int tool_expand_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts);
+// expand a flow equation
+int tool_expand(Str_Array str_args, Meantools_Options opts);
 
 #endif
+
+
diff --git a/src/number.c b/src/number.c
index e63427f..5bc87de 100644
--- a/src/number.c
+++ b/src/number.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -27,6 +27,7 @@ limitations under the License.
 #include "tools.h"
 #include "rational.h"
 #include "array.h"
+#include "parse_file.h"
 
 // init
 int init_Number(Number* number, int memory){
@@ -271,50 +272,107 @@ Number number_Qprod_ret(Q q, Number x){
   return(ret);
 }
 
-// inverse
-int number_inverse_inplace(Number* inout){
-  int i;
-  for(i=0;i<(*inout).length;i++){
-    if((*inout).base[i]>0){
-      (*inout).scalars[i]=Q_inverse((*inout).scalars[i]);
-      (*inout).scalars[i].denominator*=(*inout).base[i];
+// quotient of two numbers
+// recursively simplify numerator/denominator until denominator only has one term
+// the output is set to 'numerator'
+// both numerator and denominator may be changed by this function
+// this algorithm is not optimal in cases where denominator has several terms, in particular if their bases are large
+// it is optimal if denominator only has one term
+int number_quot_inplace(Number* numerator, Number* denominator){
+  Number tmp;
+  int i,factor;
+
+  switch((*denominator).length){
+  // error
+  case 0:
+    fprintf(stderr,"error: attempting to invert 0\n");
+    exit(-1);
+    break;
+
+  // trivial case
+  case 1:
+    // trivial base
+    if((*denominator).base[0]==1){
+      number_Qprod_chain(Q_inverse((*denominator).scalars[0]), numerator);
     }
-    else if((*inout).base[i]<0){
-      (*inout).scalars[i]=Q_inverse((*inout).scalars[i]);
-      (*inout).scalars[i].denominator*=-(*inout).base[i];
-      (*inout).scalars[i].numerator*=-1;
+    // non-trivial base
+    else{
+      // set tmp=1/denominator
+      tmp=number_one();
+      tmp.base[0]=(*denominator).base[0];
+      tmp.scalars[0].denominator=(*denominator).scalars[0].numerator*abs((*denominator).base[0]);
+      if((*denominator).base[0]>0){
+	tmp.scalars[0].numerator=(*denominator).scalars[0].denominator;
+      }
+      else if((*denominator).base[0]<0){
+	tmp.scalars[0].numerator=-(*denominator).scalars[0].denominator;
+      }
+      else{
+	fprintf(stderr,"error: attempting to invert 0\n");
+	exit(-1);
+      }
+      number_prod_chain(tmp, numerator);
+    }
+    break;
+
+  default:
+    // find non-trivial basis
+    for(i=0;(*denominator).base[i]==1;i++){}
+
+    // smallest prime factor of the base
+    if((*denominator).base[i]<0){
+      factor=-1;
+    }
+    else if((*denominator).base[i]==2){
+      factor=2;
     }
     else{
-      fprintf(stderr,"error: attempting to invert 0\n");
-      exit(-1);
+      // COMMENT: this is not optimal, but, provided the basis is not too large, this should not be a problem
+      for(factor=3;is_factor(factor, (*denominator).base[i])==0;factor=factor+2){}
+    }
+
+    // tmp is set to ((terms that k does not divide)-(terms that k divides))
+    // tmp will be multiplied to numerator and denominator
+    init_Number(&tmp,(*denominator).length);
+
+    // find all terms whose base is a multiple of k
+    for(i=0;i<(*denominator).length;i++){
+      if(is_factor(factor, (*denominator).base[i])){
+	// add to tmp with a - sign
+	number_add_elem(quot(-(*denominator).scalars[i].numerator,(*denominator).scalars[i].denominator), (*denominator).base[i], &tmp);
+      }
+      else{
+	// add to tmp
+	number_add_elem((*denominator).scalars[i], (*denominator).base[i], &tmp);
+      }
     }
+
+    number_prod_chain(tmp, numerator);
+    number_prod_chain(tmp, denominator);
+    free_Number(tmp);
+
+    // recurse
+    number_quot_inplace(numerator, denominator);
   }
+
   return(0);
 }
-// write to output
-int number_inverse(Number input, Number* output){
-  number_cpy(input,output);
-  number_inverse_inplace(output);
-  return(0);
-}
-// return result
-Number number_inverse_ret(Number x){
-  Number ret;
-  number_inverse(x,&ret);
-  return(ret);
-}
 
-// quotient
+// not inplace
 int number_quot(Number x1, Number x2, Number* output){
-  Number inv;
-  number_inverse(x2, &inv);
-  number_prod(x1, inv, output);
-  free_Number(inv);
+  Number numerator, denominator;
+  number_cpy(x1, &numerator);
+  number_cpy(x2, &denominator);
+  number_quot_inplace(&numerator, &denominator);
+  *output=numerator;
+  free_Number(denominator);
   return(0);
 }
-int number_quot_chain(Number x1, Number* inout){
-  number_inverse_inplace(inout);
-  number_prod_chain(x1, inout);
+int number_quot_chain(Number* inout, Number x2){
+  Number tmp;
+  number_quot(*inout,x2,&tmp);
+  free_Number(*inout);
+  *inout=tmp;
   return(0);
 }
 Number number_quot_ret(Number x1, Number x2){
@@ -419,7 +477,7 @@ int number_print(Number number){
   Char_Array buffer;
   init_Char_Array(&buffer,5*number.length);
   number_sprint(number, &buffer);
-  printf("%s",buffer.str);
+  printf("%s",char_array_to_str_noinit(&buffer));
   return(0);
 }
 
@@ -434,18 +492,54 @@ int str_to_Number(char* str, Number* number){
   char* buffer_ptr=buffer;
   Q num;
   int base;
-  // whether there are parentheses in the string
-  int exist_parenthesis=0;
+  int ret;
+  int j;
+  char* aux_str;
+  int aux_free=0;
 
   init_Number(number, NUMBER_SIZE);
 
+  // check whether the string is blank (return 0 in that case)
+  for(ptr=str;*ptr!='\0';ptr++){
+    if(*ptr!=' ' && *ptr!='\n'){
+      break;
+    }
+  }
+  // blank string
+  if(*ptr=='\0'){
+    number_add_elem(quot(0,1), 1, number);
+    free(buffer);
+    return(0);
+  }
+
   // init num and base
-  // init to 0 so that if str is empty, then the number is set to 0
-  num=quot(0,1);
+  num=quot(1,1);
   base=1;
 
+  // check whether the str only contains a rational number, and add parentheses
+  // keep rtack of the length of str
+  for(j=0,ptr=str;*ptr!='\0';j++,ptr++){
+    if((*ptr<'0' || *ptr>'9') && *ptr!='-' && *ptr!='/'){
+      break;
+    }
+  }
+  // only rational
+  if(*ptr=='\0'){
+    aux_str=calloc(j+3,sizeof(char));
+    aux_str[0]='(';
+    for(j=0,ptr=str;*ptr!='\0';ptr++,j++){
+      aux_str[j+1]=*ptr;
+    }
+    aux_str[j+1]=')';
+    aux_str[j+2]='\0';
+    aux_free=1;
+  }
+  else{
+    aux_str=str;
+  }
+
   mode=PP_NULL_MODE;
-  for(ptr=str;*ptr!='\0';ptr++){
+  for(ptr=aux_str;*ptr!='\0';ptr++){
     switch(*ptr){
     // read number
     case '(':
@@ -453,7 +547,10 @@ int str_to_Number(char* str, Number* number){
 	// init base
 	base=1;
 	mode=PP_NUM_MODE;
-	exist_parenthesis=1;
+      }
+      else{
+	fprintf(stderr,"syntax error: misplaced '(' in number '%s'\n",str);
+	exit(-1);
       }
       break;
     case ')':
@@ -463,6 +560,10 @@ int str_to_Number(char* str, Number* number){
 	*buffer_ptr='\0';
 	mode=PP_NULL_MODE;
       }
+      else{
+	fprintf(stderr,"syntax error: mismatched ')' in number '%s'\n",str);
+	exit(-1);
+      }
       break;
 
     // read sqrt
@@ -474,16 +575,26 @@ int str_to_Number(char* str, Number* number){
       if(mode==PP_NULL_MODE){
 	mode=PP_SQRT_MODE;
       }
-      // if there is a square root, then do not read a fraction (see end of loop)
-      exist_parenthesis=1;
+      else{
+	fprintf(stderr,"syntax error: misplaced '{' in number '%s'\n",str);
+	exit(-1);
+      }
       break;
     case '}':
       if(mode==PP_SQRT_MODE){
-	sscanf(buffer,"%d",&base);
+	ret=read_int(buffer,&base);
+	if(ret<0){
+	  fprintf(stderr,"syntax error: number base should be an integer, got '%s' in '%s'",buffer,str);
+	  exit(-1);
+	}
 	buffer_ptr=buffer;
 	*buffer_ptr='\0';
 	mode=PP_NULL_MODE;
       }
+      else{
+	fprintf(stderr,"syntax error: mismatched '}' in number '%s'\n",str);
+	exit(-1);
+      }
       break;
 
     // write num
@@ -494,24 +605,40 @@ int str_to_Number(char* str, Number* number){
 	num=quot(0,1);
 	base=1;
 	}
+      else{
+	fprintf(stderr,"syntax error: misplaced '+' in number '%s'\n",str);
+	exit(-1);
+      }
       break;
 
+    // ignore 's', ' ' and '\n'
+    case 's':break;
+    case ' ':break;
+    case '\n':break;
     default:
       if(mode!=PP_NULL_MODE){
 	buffer_ptr=str_addchar(buffer_ptr,*ptr);
       }
+      else{
+	fprintf(stderr,"syntax error: unrecognized character '%c' in number '%s'\n",*ptr,str);
+	exit(-1);
+      }
       break;
     }
   }
 
-  // last step
   if(mode==PP_NULL_MODE){
-    if(exist_parenthesis==0){
-      str_to_Q(str, &num);
-    }
     number_add_elem(num, base, number);
   }
+  else{
+    fprintf(stderr,"syntax error: mismatched '(' in number '%s'\n",str);
+    exit(-1);
+  }
 
+
+  if(aux_free==1){
+    free(aux_str);
+  }
   free(buffer);
   return(0);
 }
diff --git a/src/number.h b/src/number.h
index a33d425..5eacdaa 100644
--- a/src/number.h
+++ b/src/number.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -78,16 +78,11 @@ int number_Qprod(Q q, Number x, Number* inout);
 // return result
 Number number_Qprod_ret(Q q, Number x);
 
-// inverse
-int number_inverse_inplace(Number* inout);
-// write to output
-int number_inverse(Number input, Number* output);
-// return result
-Number number_inverse_ret(Number x);
-
-// quotient
+// quotient of two numbers
+int number_quot_inplace(Number* numerator, Number* denominator);
+// not inplace
 int number_quot(Number x1, Number x2, Number* output);
-int number_quot_chain(Number x1, Number* inout);
+int number_quot_chain(Number* inout, Number x2);
 Number number_quot_ret(Number x1, Number x2);
 
 // remove 0's
diff --git a/src/numkondo.c b/src/numkondo.c
index 0568861..91e6be2 100644
--- a/src/numkondo.c
+++ b/src/numkondo.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -86,12 +86,6 @@ int read_args_numkondo(int argc,const char* argv[], Str_Array* str_args, Numkond
   // whether a file was specified on the command-line
   int exists_file=0;
 
-  // if there are no arguments
-  if(argc==1){
-    print_usage_numkondo();
-    exit(-1);
-  }
-
   // defaults
   // display entire flow
   (*opts).display_mode=DISPLAY_NUMERICAL;
@@ -193,6 +187,8 @@ int numflow(Str_Array str_args, Numkondo_Options opts){
   Grouped_Polynomial flow_equation;
   // whether or not to use mpfr floats
   int mpfr_flag=0;
+  // postprocess flow equation
+  Grouped_Polynomial postprocess_flow_equation;
 
   // set mpfr defaults
   if(opts.mpfr_prec!=0){
@@ -205,7 +201,7 @@ int numflow(Str_Array str_args, Numkondo_Options opts){
   }
 
 
-  // parse id table
+  // parse labels
   arg_index=find_str_arg("labels", str_args);
   if(arg_index<0){
     fprintf(stderr,"error: no labels entry in the configuration file\n");
@@ -225,6 +221,16 @@ int numflow(Str_Array str_args, Numkondo_Options opts){
     char_array_to_Grouped_Polynomial(str_args.strs[arg_index], &flow_equation);
   }
 
+  // parse postprocess operation
+  arg_index=find_str_arg("postprocess_operation", str_args);
+  if(arg_index>=0){
+    char_array_to_Grouped_Polynomial(str_args.strs[arg_index], &postprocess_flow_equation);
+  }
+  else{
+    init_Grouped_Polynomial(&postprocess_flow_equation,1);
+  }
+
+
   // initial conditions
   // check they were not specified on the command line
   if(opts.eval_rccstring.length==-1){
@@ -252,17 +258,18 @@ int numflow(Str_Array str_args, Numkondo_Options opts){
   }
 
   if(mpfr_flag==0){
-    numerical_flow(flow_equation, init_cd, labels, opts.niter, opts.display_mode);
+    numerical_flow(flow_equation, init_cd, postprocess_flow_equation, labels, opts.niter, opts.display_mode);
     free_RCC(init_cd);
   }
   else{
-    numerical_flow_mpfr(flow_equation, init_cd_mpfr, labels, opts.niter, opts.display_mode);
+    numerical_flow_mpfr(flow_equation, init_cd_mpfr, postprocess_flow_equation, labels, opts.niter, opts.display_mode);
     free_RCC_mpfr(init_cd_mpfr);
   }
 
 
   // free memory
   free_Labels(labels);
+  free_Grouped_Polynomial(postprocess_flow_equation);
   free_Grouped_Polynomial(flow_equation);
   return(0);
 }
diff --git a/src/parse_file.c b/src/parse_file.c
index ae3a9af..e998803 100644
--- a/src/parse_file.c
+++ b/src/parse_file.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -30,6 +30,8 @@ limitations under the License.
 #include "istring.h"
 #include "tools.h"
 #include "idtable.h"
+#include "tree.h"
+#include "symbolic_parser.h"
 
 
 // parsing modes
@@ -59,6 +61,70 @@ limitations under the License.
 #define PP_GROUP_MODE 16
 
 
+// read a positive integer from a string
+int read_positive_int(char* str, int* out){
+  char* ptr;
+  int res;
+  for(ptr=str;*ptr!='\0';ptr++){
+    if(*ptr-'0'>=10 || *ptr-'0'<0){
+      return(-1);
+    }
+  }
+  res=sscanf(str,"%d",out);
+  if(res!=1){
+    return(-2);
+  }
+  return(0);
+}
+// read an integer from a string
+int read_int(char* str, int* out){
+  char* ptr;
+  int res;
+  for(ptr=str;*ptr!='\0';ptr++){
+    if((*ptr-'0'>=10 || *ptr-'0'<0) && (*ptr!='-')){
+      return(-1);
+    }
+  }
+  res=sscanf(str,"%d",out);
+  if(res!=1){
+    return(-2);
+  }
+  return(0);
+}
+// read an long int from a string
+int read_long_int(char* str, long int* out){
+  char* ptr;
+  int res;
+  for(ptr=str;*ptr!='\0';ptr++){
+    if((*ptr-'0'>=10 || *ptr-'0'<0) && (*ptr!='-')){
+      return(-1);
+    }
+  }
+  res=sscanf(str,"%ld",out);
+  if(res!=1){
+    return(-2);
+  }
+  return(0);
+}
+
+// read a long double
+int read_long_double(char* str, long double* out){
+  char* ptr;
+  int res;
+  for(ptr=str;*ptr!='\0';ptr++){
+    if((*ptr-'0'>=10 || *ptr-'0'<0) && (*ptr!='-') && (*ptr!='e') && (*ptr!='E') && (*ptr!='.')){
+      return(-1);
+    }
+  }
+  res=sscanf(str,"%Lf",out);
+  if(res!=1){
+    return(-2);
+  }
+  return(0);
+}
+
+
+
 // parse fields list
 int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
   // buffer
@@ -67,6 +133,8 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
   int i,j;
   int mode;
   int comment=0;
+  int res;
+  int expect_colon=0;
 
   // allocate memory
   init_Fields_Table(fields);
@@ -86,39 +154,100 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
 	if(mode==PP_NULL_MODE){
 	  mode=PP_PARAMETER_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: fields entry (%d): 'h' should be at the beginning of the line\n",j);
+	  exit(-1);
+	}
+	// flag: colon is expected immediately after letter
+	expect_colon=1;
 	break;
       // external fields
       case 'x':
 	if(mode==PP_NULL_MODE){
 	  mode=PP_EXTERNAL_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: fields entry (%d): 'x' should be at the beginning of the line\n",j);
+	  exit(-1);
+	}
+	// flag: colon is expected immediately after letter
+	expect_colon=1;
 	break;
       // internal fields
       case 'i':
 	if(mode==PP_NULL_MODE){
 	  mode=PP_INTERNAL_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: fields entry (%d): 'i' should be at the beginning of the line\n",j);
+	  exit(-1);
+	}
+	// flag: colon is expected immediately after letter
+	expect_colon=1;
 	break;
       case 'f':
 	if(mode==PP_NULL_MODE){
 	  mode=PP_FERMIONS_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: fields entry (%d): 'f' should be at the beginning of the line\n",j);
+	  exit(-1);
+	}
+	// flag: colon is expected immediately after letter
+	expect_colon=1;
 	break;
       case 'a':
 	if(mode==PP_NULL_MODE){
 	  mode=PP_NONCOMMUTING_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: fields entry (%d): 'a' should be at the beginning of the line\n",j);
+	  exit(-1);
+	}
+	// flag: colon is expected immediately after letter
+	expect_colon=1;
 	break;
 
       // reset buffer
       case ':':
+	// check mode
+	if(mode==PP_NULL_MODE || expect_colon==0){
+	  fprintf(stderr,"syntax error: fields entry (%d): ':' must be preceded by 'h', 'x', 'i', 'f', or 'a' on the same line\n",j);
+	  exit(-1);
+	}
+	expect_colon=0;
 	buffer_ptr=buffer;
 	*buffer_ptr='\0';
 	break;
 
       // write to fields
       case ',':
-	sscanf(buffer,"%d",&i);
+	res=read_positive_int(buffer,&i);
+
+	// check i
+	if(res<0){
+	  fprintf(stderr,"syntax error: fields entry (%d): fields must be non-negative integers, got '%s'\n",j,buffer);
+	  exit(-1);
+	}
+	if(mode==PP_PARAMETER_MODE || mode==PP_EXTERNAL_MODE || mode==PP_INTERNAL_MODE){
+	  if(int_array_find(i,fields->parameter)>=0 || int_array_find(i,fields->external)>=0 || int_array_find(i,fields->internal)>=0){
+	    fprintf(stderr,"error: fields entry (%d): field %d is already present in the fields table\n",j,i);
+	    exit(-1);
+	  }
+	}
+	else if(mode==PP_FERMIONS_MODE){
+	  if(int_array_find(i,fields->fermions)>=0){
+	    fprintf(stderr,"error: fields entry (%d): field %d is already present in the Fermions table\n",j,i);
+	    exit(-1);
+	  }
+	}    
+	else if(mode==PP_NONCOMMUTING_MODE){
+	  if(int_array_find(i,fields->noncommuting)>=0){
+	    fprintf(stderr,"error: fields entry (%d): field %d is already present in the noncommuting table\n",j,i);
+	    exit(-1);
+	  }
+	}    
+
 	if(mode==PP_PARAMETER_MODE){
 	  int_array_append(i,&((*fields).parameter));
 	}
@@ -134,31 +263,65 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
 	else if(mode==PP_NONCOMMUTING_MODE){
 	  int_array_append(i,&((*fields).noncommuting));
 	}
+	else{
+	  fprintf(stderr,"syntax error: fields entry (%d): misplaced ','\n",j);
+	  exit(-1);
+	}
 	buffer_ptr=buffer;
 	*buffer_ptr='\0';
 	break;
 
       // back to null mode
       case '\n':
-	sscanf(buffer,"%d",&i);
-	if(mode==PP_PARAMETER_MODE){
-	  int_array_append(i,&((*fields).parameter));
-	}
-	else if(mode==PP_EXTERNAL_MODE){
-	  int_array_append(i,&((*fields).external));
-	}
-	else if(mode==PP_INTERNAL_MODE){
-	  int_array_append(i,&((*fields).internal));
-	}
-	else if(mode==PP_FERMIONS_MODE){
-	  int_array_append(i,&((*fields).fermions));
-	}
-	else if(mode==PP_NONCOMMUTING_MODE){
-	  int_array_append(i,&((*fields).noncommuting));
+	// ignore if in NULL_MODE
+	if(mode!=PP_NULL_MODE){
+	  res=read_positive_int(buffer,&i);
+	  // check i
+	  if(res<0){
+	    fprintf(stderr,"syntax error: fields entry (%d): fields must be non-negative integers, got '%s'\n",j,buffer);
+	    exit(-1);
+	  }
+	  if(mode==PP_PARAMETER_MODE || mode==PP_EXTERNAL_MODE || mode==PP_INTERNAL_MODE){
+	    if(int_array_find(i,fields->parameter)>=0 || int_array_find(i,fields->external)>=0 || int_array_find(i,fields->internal)>=0){
+	      fprintf(stderr,"error: fields entry (%d): field %d is already present in the fields table\n",j,i);
+	      exit(-1);
+	    }
+	  }
+	  else if(mode==PP_FERMIONS_MODE){
+	    if(int_array_find(i,fields->fermions)>=0){
+	      fprintf(stderr,"error: fields entry (%d): field %d is already present in the Fermions table\n",j,i);
+	      exit(-1);
+	    }
+	  }    
+	  else if(mode==PP_NONCOMMUTING_MODE){
+	    if(int_array_find(i,fields->noncommuting)>=0){
+	      fprintf(stderr,"error: fields entry (%d): field %d is already present in the noncommuting table\n",j,i);
+	      exit(-1);
+	    }
+	  }    
+
+	  if(mode==PP_PARAMETER_MODE){
+	    int_array_append(i,&((*fields).parameter));
+	  }
+	  else if(mode==PP_EXTERNAL_MODE){
+	    int_array_append(i,&((*fields).external));
+	  }
+	  else if(mode==PP_INTERNAL_MODE){
+	    int_array_append(i,&((*fields).internal));
+	  }
+	  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;
 	}
-	mode=PP_NULL_MODE;
 	break;
 
+  
+      // ignore spaces
+      case ' ':break;
       // comment
       case '#':
 	comment=1;
@@ -168,57 +331,85 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
 	if(mode!=PP_NULL_MODE){
 	  buffer_ptr=str_addchar(buffer_ptr,str_fields.str[j]);
 	}
+	else{
+	  fprintf(stderr,"syntax error: fields entry (%d): unrecognized character '%c'\n",j,str_fields.str[j]);
+	  exit(-1);
+	}
+	break;
+
 	break;
       }
     }
   }
   free(buffer);
+
+  // check whether there are non-Fermionic internal fields
+  for(i=0;i<(*fields).internal.length;i++){
+    if(is_fermion((*fields).internal.values[i], *fields)==0){
+      fprintf(stderr,"warning: some internal fields are not Fermions: means may be computed using sums over permutations rather than the more efficient LU decomposition\n");
+      break;
+    }
+  }
+
   return(0);
 }
 
 
-// parse symbols list
+// parse virtual_fields list
 // write result to fields
-int parse_input_symbols(Char_Array str_symbols, Fields_Table* fields){
+int parse_input_virtual_fields(Char_Array str_virtual_fields, Fields_Table* fields, Variables variables){
   // buffer
-  char* buffer=calloc(str_symbols.length+1,sizeof(char));
+  char* buffer=calloc(str_virtual_fields.length+1,sizeof(char));
   char* buffer_ptr=buffer;
   Polynomial polynomial;
   int index;
   int i,j;
   int mode;
   int comment=0;
+  int res;
 
   // loop over input
   mode=PP_INDEX_MODE;
-  for(j=0;j<str_symbols.length;j++){
+  for(j=0;j<str_virtual_fields.length;j++){
     if(comment==1){
-      if(str_symbols.str[j]=='\n'){
+      if(str_virtual_fields.str[j]=='\n'){
 	comment=0;
       }
     }
     // stay in polynomial mode until ','
     else if(mode==PP_POLYNOMIAL_MODE){
-      if(str_symbols.str[j]==','){
+      if(str_virtual_fields.str[j]==','){
 	// parse polynomial
-	str_to_Polynomial(buffer, &polynomial);
+	parse_symbolic_expression_str(buffer, *fields, variables, &polynomial);
 	// write index and polynomial
-	symbols_append_noinit(index, polynomial, &((*fields).symbols));
+	virtual_fields_append_noinit(index, polynomial, &((*fields).virtual_fields));
 	mode=PP_INDEX_MODE;
 	buffer_ptr=buffer;
 	*buffer_ptr='\0';
       }
+      // check there are no '=' signs in polynomial (forgotten ',')
+      else if(str_virtual_fields.str[j]=='='){
+	fprintf(stderr,"syntax error: virtual fields entry (%d): '=' in polynomial '%s'\n",j, buffer);
+	exit(-1);
+      }
+      else if(str_virtual_fields.str[j]=='#'){
+	comment=1;
+      }
       else{
-	buffer_ptr=str_addchar(buffer_ptr,str_symbols.str[j]);
+	buffer_ptr=str_addchar(buffer_ptr,str_virtual_fields.str[j]);
       }
     }
     else{
-      switch(str_symbols.str[j]){
+      switch(str_virtual_fields.str[j]){
       // polynomial mode
       case '=':
 	if(mode==PP_INDEX_MODE){
 	  // read index
-	  sscanf(buffer,"%d",&index);
+	  res=read_positive_int(buffer, &index);
+	  if(res<0){
+	    fprintf(stderr,"syntax error: virtual fields entry (%d): virtual field index must be a non-negative integer, got '%s'\n",j,buffer);
+	    exit(-1);
+	  }
 	  // reset buffer
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
@@ -226,43 +417,53 @@ int parse_input_symbols(Char_Array str_symbols, Fields_Table* fields){
 	}
 	break;
 
+      // misplaced ','
+      case ',':
+	fprintf(stderr,"syntax error: virtual fields entry (%d): misplaced ',' or forgotten virtual field index\n",j);
+	exit(-1);
+	break;
+
+      // ignore spaces and newlines
+      case ' ':break;
+      case '\n':break;
+
       // comment
       case '#':
 	comment=1;
 	break;
 
       default:
-	buffer_ptr=str_addchar(buffer_ptr,str_symbols.str[j]);
+	buffer_ptr=str_addchar(buffer_ptr,str_virtual_fields.str[j]);
 	break;
       }
     }
   }
 
   // last step
-  if(polynomial.length>0){
-    str_to_Polynomial(buffer, &polynomial);
-    symbols_append_noinit(index, polynomial, &((*fields).symbols));
-  }
+  parse_symbolic_expression_str(buffer, *fields, variables, &polynomial);
+  virtual_fields_append_noinit(index, polynomial, &((*fields).virtual_fields));
 
   // simplify
-  for(i=0;i<(*fields).symbols.length;i++){
-    polynomial_simplify((*fields).symbols.expr+i, *fields);
+  for(i=0;i<(*fields).virtual_fields.length;i++){
+    polynomial_simplify((*fields).virtual_fields.expr+i, *fields);
   }
 
   free(buffer);
   return(0);
 }
 
+
 // parse groups of independent fields
-int parse_input_groups(Char_Array str_groups, Groups* groups){
+int parse_input_groups(Char_Array str_groups, Groups* groups, Polynomial_Matrix propagator, Fields_Table fields){
   // buffer
   char* buffer=calloc(str_groups.length+1,sizeof(char));
   char* buffer_ptr=buffer;
   int index;
-  int j;
+  int i,j;
   Int_Array group;
   int mode;
   int comment=0;
+  int res;
 
   // alloc
   init_Groups(groups, GROUP_SIZE);
@@ -287,26 +488,60 @@ int parse_input_groups(Char_Array str_groups, Groups* groups){
 	  init_Int_Array(&group, GROUP_SIZE);
 	  mode=PP_GROUP_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: groups entry (%d): nested parentheses are not allowed in this entry\n",j);
+	  exit(-1);
+	}
 	break;
       case')':
 	if(mode==PP_GROUP_MODE){
-	  sscanf(buffer,"%d",&index);
+	  res=read_positive_int(buffer,&index);
+	  if(res<0){
+	    fprintf(stderr,"syntax error: groups entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer);
+	    exit(-1);
+	  }
+	  for(i=0;i<groups->length;i++){
+	    if(int_array_find(index,groups->indices[i])>=0){
+	      fprintf(stderr,"syntax error: groups entry (%d): index %d is already present in group %d\n",j,index,i);
+	      exit(-1);
+	    }
+	  }
+
 	  int_array_append(index, &group);
 	  groups_append_noinit(group, groups);
 	  mode=PP_NULL_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: groups entry (%d): extra ')'\n",j);
+	  exit(-1);
+	}
 	break;
 
       // read index
       case',':
 	if(mode==PP_GROUP_MODE){
-	  sscanf(buffer,"%d",&index);
+	  res=read_positive_int(buffer,&index);
+	  if(res<0){
+	    fprintf(stderr,"syntax errorL groups entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer);
+	    exit(-1);
+	  }
+	  for(i=0;i<groups->length;i++){
+	    if(int_array_find(index,groups->indices[i])>=0){
+	      fprintf(stderr,"syntax error: groups entry (%d): index %d is already present in group %d\n",j,index,i);
+	      exit(-1);
+	    }
+	  }
+
 	  int_array_append(index, &group);
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
 	}
 	break;
 
+      // skip spaces and newlines
+      case ' ':break;
+      case '\n':break;
+
       // comment
       case '#':
 	comment=1;
@@ -316,11 +551,204 @@ int parse_input_groups(Char_Array str_groups, Groups* groups){
 	if(mode!=PP_NULL_MODE){
 	  buffer_ptr=str_addchar(buffer_ptr,str_groups.str[j]);
 	}
+	else{
+	  fprintf(stderr,"syntax error: groups entry (%d): unrecognized character '%c'\n",j,str_groups.str[j]);
+	  exit(-1);
+	}
+	break;
+      }
+    }
+  }
+
+  // check fields in groups
+  check_groups(*groups, propagator, fields);
+
+  free(buffer);
+  return(0);
+}
+
+// check that the members of groups are independent (assuming the virtual fields and propagator were already parsed)
+int check_groups(Groups groups, Polynomial_Matrix propagator, Fields_Table fields){
+  Int_Array* group_fields=calloc(groups.length,sizeof(Int_Array));
+  int i,j,k,l,a,b;
+
+  // get the lists of fields involved in each group
+  for(i=0;i<groups.length;i++){
+    // expand virtual_fields
+    fields_in_virtual_field_list(groups.indices[i], fields, group_fields+i);
+  }
+
+  // check whether pairs of fields are present in the propagator
+  for(i=0;i<groups.length;i++){
+    for(j=0;j<group_fields[i].length;j++){
+      if(field_type(group_fields[i].values[j], fields)==FIELD_INTERNAL){
+	for(k=0;k<groups.length;k++){
+	  if(k!=i){
+	    for(l=0;l<group_fields[k].length;l++){
+	      if(field_type(group_fields[k].values[l], fields)==FIELD_INTERNAL){
+		if(group_fields[i].values[j]==group_fields[k].values[l]){
+		  fprintf(stderr,"error: groups %d and %d both contain internal field %d\n",i,k,group_fields[i].values[j]);
+		  exit(-1);
+		}
+		a=intlist_find(propagator.indices, propagator.length, group_fields[i].values[j]);
+		b=intlist_find(propagator.indices, propagator.length, group_fields[k].values[l]);
+		if(a>=0 && b>=0 && polynomial_is_zero(propagator.matrix[a][b])==0){
+		  fprintf(stderr,"error: groups %d and %d are not independent: they contain %d and %d which have a non-vanishing propagator\n",i,k,group_fields[i].values[j], group_fields[k].values[l]);
+		  exit(-1);
+		}
+	      }
+	    }
+	  }
+	}
+      }
+    }
+  }
+
+  for(i=0;i<groups.length;i++){
+    free_Int_Array(group_fields[i]);
+  }
+  free(group_fields);
+
+  return(0);
+}
+
+// list of fields involved in a list of virtual_fields
+int fields_in_virtual_field_list(Int_Array indices, Fields_Table fields, Int_Array* output){
+  int i,j,k;
+  Int_Array tmp_array;
+
+  init_Int_Array(output,FIELDS_SIZE);
+  for(k=0;k<indices.length;k++){
+    if(field_type(indices.values[k],fields)!=FIELD_VIRTUAL){
+      int_array_append(abs(indices.values[k]), output);
+    }
+    else{
+      i=intlist_find_err(fields.virtual_fields.indices, fields.virtual_fields.length, indices.values[k]);
+      for(j=0;j<fields.virtual_fields.expr[i].length;j++){
+	fields_in_virtual_field_list(fields.virtual_fields.expr[i].monomials[j], fields, &tmp_array);
+	int_array_concat_unique(tmp_array, output);
+	free_Int_Array(tmp_array);
+      }
+    }
+  }
+
+  return(0);
+}
+
+
+
+// parse variables list
+int parse_input_variables(Char_Array str_variables, Variables* variables){
+  // buffer
+  char* buffer=calloc(str_variables.length+1,sizeof(char));
+  char* buffer_ptr=buffer;
+  Tree symbol_tree;
+  Char_Array varname;
+  int j;
+  int mode;
+  int comment=0;
+  char* ptr;
+
+  init_Variables(variables,VARIABLES_SIZE);
+
+  // loop over input
+  mode=PP_INDEX_MODE;
+  for(j=0;j<str_variables.length;j++){
+    if(comment==1){
+      if(str_variables.str[j]=='\n'){
+	comment=0;
+      }
+    }
+    // stay in polynomial mode until ','
+    else if(mode==PP_POLYNOMIAL_MODE){
+      if(str_variables.str[j]==','){
+	// parse symbol tree
+	str_to_symbol_tree(buffer, &symbol_tree);
+	variables_append_noinit(varname , symbol_tree, variables);
+	mode=PP_INDEX_MODE;
+	buffer_ptr=buffer;
+	*buffer_ptr='\0';
+      }
+      // check there are no '=' signs in polynomial (forgotten ',')
+      else if(str_variables.str[j]=='='){
+	fprintf(stderr,"syntax error: preprocessor_variables entry (%d): '=' in polynomial '%s'\n",j, buffer);
+	exit(-1);
+      }
+      else if(str_variables.str[j]=='#'){
+	comment=1;
+      }
+      else{
+	buffer_ptr=str_addchar(buffer_ptr,str_variables.str[j]);
+      }
+    }
+    else{
+      switch(str_variables.str[j]){
+      // polynomial mode
+      case '=':
+	if(mode==PP_INDEX_MODE){
+	  if(*buffer=='\0'){
+	    fprintf(stderr,"syntax error: preprocessor_variables entry (%d): empty variable name\n",j);
+	    exit(-1);
+	  }
+	  for(ptr=buffer;*ptr!='\0';ptr++){
+	    if(*ptr=='<' || *ptr=='>'){
+	      fprintf(stderr,"syntax error: preprocessor_variables entry (%d): forbidden '%c' character in variable name '%s'\n",j,*ptr,buffer);
+	      exit(-1);
+	    }
+	  }
+	  str_to_char_array(buffer,&varname);
+	  // check that the variable name is not 'OUT' or 'FLOW' or 'RCC'
+	  if(char_array_cmp_str(varname, "OUT")==1){
+	    fprintf(stderr,"error: preprocessor_variables entry (%d): variable name cannot be 'OUT' (this name is reserved)\n",j);
+	    exit(-1);
+	  }
+	  if(char_array_cmp_str(varname, "FLOW")==1){
+	    fprintf(stderr,"error: preprocessor_variables entry (%d): variable name cannot be 'FLOW' (this name is reserved)\n",j);
+	    exit(-1);
+	  }
+	  if(char_array_cmp_str(varname, "RCC")==1){
+	    fprintf(stderr,"error: preprocessor_variables entry (%d): variable name cannot be 'RCC' (this name is reserved)\n",j);
+	    exit(-1);
+	  }
+	  // reset buffer
+	  buffer_ptr=buffer;
+	  *buffer_ptr='\0';
+	  mode=PP_POLYNOMIAL_MODE;
+	}
+	break;
+
+      // misplaced ','
+      case ',':
+	fprintf(stderr,"syntax error: preprocessor_variables entry (%d): misplaced ',' or forgotten variable name\n",j);
+	exit(-1);
+	break;
+
+      // ignore characters
+      case ' ':break;
+      case '\n': break;
+
+      // comment
+      case '#':
+	comment=1;
+	break;
+
+      default:
+	buffer_ptr=str_addchar(buffer_ptr,str_variables.str[j]);
 	break;
       }
     }
   }
 
+  // last step
+  if(mode==PP_POLYNOMIAL_MODE){
+    str_to_symbol_tree(buffer, &symbol_tree);
+    variables_append_noinit(varname , symbol_tree, variables);
+  }
+  else if(*buffer!='\0'){
+    fprintf(stderr,"syntax error: preprocessor_variables entry (%d): mismatched variable/polynomial pair, got '%s'\n",j,buffer);
+    exit(-1);
+  }
+
   free(buffer);
   return(0);
 }
@@ -328,7 +756,7 @@ int parse_input_groups(Char_Array str_groups, Groups* groups){
 
 // parse identities between fields
 // write result to fields
-int parse_input_identities(Char_Array str_identities, Fields_Table* fields){
+int parse_input_identities(Char_Array str_identities, Fields_Table* fields, Variables variables){
   // buffer
   char* buffer=calloc(str_identities.length+1,sizeof(char));
   char* buffer_ptr=buffer;
@@ -339,6 +767,7 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){
   int tmp;
   int mode;
   int comment=0;
+  int ret;
 
   init_Int_Array(&monomial, MONOMIAL_SIZE);
 
@@ -354,13 +783,21 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){
     else if(mode==PP_POLYNOMIAL_MODE){
       if(str_identities.str[j]==','){
 	// parse polynomial
-	str_to_Polynomial(buffer, &polynomial);
+	parse_symbolic_expression_str(buffer, *fields, variables, &polynomial);
 	// write monomial and polynomial
 	identities_append_noinit(monomial, polynomial, &((*fields).ids));
 	// realloc
 	init_Int_Array(&monomial, MONOMIAL_SIZE);
 	mode=PP_NULL_MODE;
       }
+      // check there are no '=' signs in polynomial (forgotten ',')
+      else if(str_identities.str[j]=='='){
+	fprintf(stderr,"syntax error: identities entry (%d): '=' in polynomial '%s'\n",j, buffer);
+	exit(-1);
+      }
+      else if(str_identities.str[j]=='#'){
+	comment=1;
+      }
       else{
 	buffer_ptr=str_addchar(buffer_ptr,str_identities.str[j]);
       }
@@ -377,20 +814,36 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
 	}
+	else{
+	  fprintf(stderr,"syntax error: identities entry (%d): misplaced '['\n",j);
+	  exit(-1);
+	}
 	break;
       // for notational homogeneity
       case 'f':
 	if(mode==PP_BRACKET_MODE){
 	  mode=PP_FIELD_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: identities entry (%d): misplaced 'f'\n",j);
+	  exit(-1);
+	}
 	break;
       // write monomial term
       case ']':
 	if(mode==PP_FIELD_MODE){
-	  sscanf(buffer,"%d",&i);
+	  ret=read_int(buffer,&i);
+	  if(ret<0){
+	    fprintf(stderr,"syntax error: identities entry (%d): field indices must be integers, got '%s'\n",j,buffer);
+	    exit(-1);
+	  }
 	  int_array_append(i,&monomial);
 	  mode=PP_INDEX_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: identities entry (%d): mismatched ']'\n",j);
+	  exit(-1);
+	}
 	break;
 
       // polynomial mode
@@ -400,8 +853,22 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){
 	  *buffer_ptr='\0';
 	  mode=PP_POLYNOMIAL_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: identities entry (%d): misplaced '='\n",j);
+	  exit(-1);
+	}
 	break;
 
+      // misplaced ','
+      case ',':
+	fprintf(stderr,"syntax error: identities entry (%d): misplaced ',' or forgotten left hand side\n",j);
+	exit(-1);
+	break;
+
+      // ignore spaces and newlines
+      case ' ':break;
+      case '\n':break;
+
       // comment
       case '#':
 	comment=1;
@@ -411,6 +878,10 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){
 	if(mode!=PP_NULL_MODE){
 	  buffer_ptr=str_addchar(buffer_ptr,str_identities.str[j]);
 	}
+	else{
+	  fprintf(stderr,"syntax error: identities entry (%d): unrecognized character '%c'\n",j,str_identities.str[j]);
+	  exit(-1);
+	}
 	break;
       }
     }
@@ -418,9 +889,13 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){
 
   // last step
   if(mode==PP_POLYNOMIAL_MODE){
-    str_to_Polynomial(buffer, &polynomial);
+    parse_symbolic_expression_str(buffer, *fields, variables, &polynomial);
     identities_append_noinit(monomial, polynomial, &((*fields).ids));
   }
+  else if(*buffer!='\0'){
+    fprintf(stderr,"syntax error: identities entry (%d): mismatched left/right hand side, got '%s'\n",j,buffer);
+    exit(-1);
+  }
   else{
     free_Int_Array(monomial);
   }
@@ -451,6 +926,7 @@ int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propaga
   int index2=-1;
   int mode;
   int comment=0;
+  int ret;
 
   // allocate memory
   init_Polynomial_Matrix(propagator, fields.internal.length);
@@ -473,20 +949,36 @@ int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propaga
       // indices
       case ';':
 	if(mode==PP_INDEX_MODE){
-	  sscanf(buffer,"%d",&i);
+	  ret=read_positive_int(buffer,&i);
+	  if(ret<0){
+	    fprintf(stderr,"syntax error: propagator entry (%d): indices must be non-negative integers, got '%s'\n",j,buffer);
+	    exit(-1);
+	  }
 	  index1=intlist_find_err((*propagator).indices, (*propagator).length, i);
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
 	}
+	else{
+	  fprintf(stderr,"syntax error: propagator entry (%d): misplaced ';'\n",j);
+	  exit(-1);
+	}
 	break;
       case ':':
 	if(mode==PP_INDEX_MODE){
-	  sscanf(buffer,"%d",&i);
+	  ret=read_positive_int(buffer,&i);
+	  if(ret<0){
+	    fprintf(stderr,"syntax error: propagator entry (%d): indices must be non-negative integers, got '%s'\n",j,buffer);
+	    exit(-1);
+	  }
 	  index2=intlist_find_err((*propagator).indices, (*propagator).length, i);
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
 	  mode=PP_POLYNOMIAL_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: propagator entry (%d): misplaced ':'\n",j);
+	  exit(-1);
+	}
 	break;
 
       // num
@@ -494,12 +986,22 @@ int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propaga
 	if(mode==PP_POLYNOMIAL_MODE && index1>=0 && index2>=0){
 	  free_Polynomial((*propagator).matrix[index1][index2]);
 	  str_to_Polynomial(buffer,(*propagator).matrix[index1]+index2);
+	  index1=-1;
+	  index2=-1;
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
 	  mode=PP_INDEX_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: propagator entry (%d): mismatched index/polynomial pair\n",j);
+	  exit(-1);
+	}
 	break;
 
+      // ignore ' 'and newline
+      case ' ':break;
+      case '\n':break;
+
       // comment
       case '#':
 	comment=1;
@@ -517,6 +1019,15 @@ int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propaga
     free_Polynomial((*propagator).matrix[index1][index2]);
     str_to_Polynomial(buffer,(*propagator).matrix[index1]+index2);
   }
+  else if (*buffer!='\0'){
+    fprintf(stderr,"syntax error: propagator entry (%d): mismatched index/polynomial pair, got '%s'\n",j,buffer);
+    exit(-1);
+  }
+
+  // check whether the polynomial is numeric
+  if(polynomial_matrix_is_numeric(*propagator)==0){
+    fprintf(stderr,"warning: the propagator contains polynomial entries: means may be computed using sums over permutations rather than the more efficient LU decomposition\n");
+  }
 
   free(buffer);
   return(0);
@@ -524,49 +1035,9 @@ int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propaga
 
 
 // parse input polynomial
-int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields){
-  int j;
-  // buffer
-  char* buffer=calloc(str_polynomial.length+1,sizeof(char));
-  char* buffer_ptr=buffer;
-  Polynomial tmp_poly;
-
-  // allocate memory
-  init_Polynomial(output,POLY_SIZE);
-
-  for(j=0;j<str_polynomial.length;j++){
-    switch(str_polynomial.str[j]){
-    case '*':
-      str_to_Polynomial(buffer, &tmp_poly);
-      if((*output).length==0){
-	polynomial_concat(tmp_poly, output);
-      }
-      else{
-	polynomial_prod_chain(tmp_poly, output, fields);
-      }
-      free_Polynomial(tmp_poly);
-
-      buffer_ptr=buffer;
-      *buffer_ptr='\0';
-      break;
-
-    default:
-      buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]);
-      break;
-    }
-  }
-
-  //last step
-  str_to_Polynomial(buffer, &tmp_poly);
-  if((*output).length==0){
-    polynomial_concat(tmp_poly, output);
-  }
-  else{
-    polynomial_prod_chain(tmp_poly, output, fields);
-  }
-  free_Polynomial(tmp_poly);
-
-  free(buffer);
+int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields, Variables variables){
+  parse_symbolic_expression(str_polynomial, fields, variables, output);
+  polynomial_simplify(output, fields);
   return(0);
 }
 
@@ -574,7 +1045,7 @@ int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields
 
 // parse id table
 // fields argument for sorting
-int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields){
+int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields, Variables variables){
   // buffer
   char* buffer=calloc(str_idtable.length+1,sizeof(char));
   char* buffer_ptr=buffer;
@@ -583,6 +1054,7 @@ int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table
   int j;
   int mode;
   int comment=0;
+  int ret;
 
   // allocate memory
   init_Id_Table(idtable,EQUATION_SIZE);
@@ -601,24 +1073,40 @@ int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table
       case ',':
 	// write polynomial
 	if(mode==PP_POLYNOMIAL_MODE){
-	  str_to_Polynomial(buffer,&polynomial);
+	  parse_symbolic_expression_str(buffer, fields, variables, &polynomial);
 	  // add to idtable
 	  idtable_append_noinit(index,polynomial,idtable);
 	  mode=PP_INDEX_MODE;
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
 	  }
+	else{
+	  fprintf(stderr,"syntax error: id_table entry (%d): misplaced ','\n",j);
+	  exit(-1);
+	}
 	break;
 
       case ':':
 	if(mode==PP_INDEX_MODE){
-	  sscanf(buffer,"%d",&index);
+	  ret=read_positive_int(buffer,&index);
+	  if(ret<0){
+	    fprintf(stderr,"syntax error: id_table entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer);
+	    exit(-1);
+	  }
 	  mode=PP_POLYNOMIAL_MODE;
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
 	}
+	else{
+	  fprintf(stderr,"syntax error: id_table entry (%d): misplaced ':'\n",j);
+	  exit(-1);
+	}
 	break;
 
+      // ignore ' ' and '\n'
+      case ' ':break;
+      case '\n':break;
+
       // comment
       case '#':
 	comment=1;
@@ -628,6 +1116,10 @@ int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table
 	if(mode!=PP_NULL_MODE){
 	  buffer_ptr=str_addchar(buffer_ptr,str_idtable.str[j]);
 	}
+	else{
+	  fprintf(stderr,"syntax error: id_table entry (%d): unrecognized character '%c'\n",j, str_idtable.str[j]);
+	  exit(-1);
+	}
 	break;
       }
     }
@@ -635,9 +1127,13 @@ int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table
 
   //last step
   if(mode==PP_POLYNOMIAL_MODE){
-    str_to_Polynomial(buffer,&polynomial);
+    parse_symbolic_expression_str(buffer, fields, variables, &polynomial);
     idtable_append_noinit(index,polynomial,idtable);
   }
+  else if(*buffer!='\0'){
+    fprintf(stderr,"syntax error: id_table entry (%d): mismatched index/polynomial pair, got '%s'\n",j,buffer);
+    exit(-1);
+  }
 
   // sort
   for(j=0;j<(*idtable).length;j++){
@@ -658,6 +1154,7 @@ int parse_labels(Char_Array str_labels, Labels* labels){
   int j;
   int mode;
   int comment=0;
+  int ret;
 
   // allocate memory
   init_Labels(labels,EQUATION_SIZE);
@@ -686,18 +1183,30 @@ int parse_labels(Char_Array str_labels, Labels* labels){
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
 	}
+	else{
+	  fprintf(stderr,"syntax error: labels entry (%d): mismatched '\"'\n",j);
+	  exit(-1);
+	}
 	break;
 
       case ':':
 	// write
 	if(mode==PP_INDEX_MODE){
 	  sscanf(buffer,"%d",&index);
+	  ret=read_positive_int(buffer,&index);
+	  if(ret<0){
+	    fprintf(stderr,"syntax error: labels entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer);
+	    exit(-1);
+	  }
+	}
+	else{
+	  fprintf(stderr,"syntax error: labels entry (%d): misplaced ':'\n",j);
+	  exit(-1);
 	}
 	break;
 
       // characters to ignore
       case ' ':break;
-      case '&':break;
       case '\n':break;
       case ',':break;
 
@@ -726,6 +1235,7 @@ int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_f
   int i,j;
   int comment_mode=0;
   int dcount=0;
+  int ret;
 
   *buffer_ptr='\0';
   // loop over the input
@@ -741,7 +1251,11 @@ int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_f
       case ',':
 	// write init
 	if(mpfr_flag==0){
-	  sscanf(buffer,"%Lf",(*init).values+intlist_find_err((*init).indices,(*init).length,index));
+	  ret=read_long_double(buffer,(*init).values+intlist_find_err((*init).indices,(*init).length,index));
+	  if(ret<0){
+	    fprintf(stderr,"syntax error: initial_condition entry (%d): failed to extract long double from '%s'\n",j,buffer);
+	    exit(-1);
+	  }
 	}
 	else{
 	  mpfr_strtofr((*init_mpfr).values[intlist_find_err((*init_mpfr).indices,(*init_mpfr).length,index)], buffer, &buffer_ptr, 10, MPFR_RNDN);
@@ -755,7 +1269,11 @@ int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_f
       // separator
       case ':':
 	// write index
-	sscanf(buffer,"%d",&i);
+	ret=read_int(buffer,&i);
+	if(ret<0){
+	  fprintf(stderr,"syntax error: initial_condition entry (%d): index must be an integer, got '%s'\n",j,buffer);
+	  exit(-1);
+	}
 	if(i<0){
 	  index=i-dcount*DOFFSET;
 	}
@@ -772,7 +1290,6 @@ int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_f
 
       // characters to ignore
       case ' ':break;
-      case '&':break;
       case '\n':break;
 
       // comments
@@ -790,7 +1307,11 @@ int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_f
 
   // write init
   if(mpfr_flag==0){
-    sscanf(buffer,"%Lf",(*init).values+intlist_find_err((*init).indices,(*init).length,index));
+    ret=read_long_double(buffer,(*init).values+intlist_find_err((*init).indices,(*init).length,index));
+    if(ret<0){
+      fprintf(stderr,"syntax error: initial_condition entry (%d): failed to extract long double from '%s'\n",j,buffer);
+      exit(-1);
+    }
   }
   else{
     mpfr_strtofr((*init_mpfr).values[intlist_find_err((*init_mpfr).indices,(*init_mpfr).length,index)], buffer, &buffer_ptr, 10, MPFR_RNDN);
diff --git a/src/parse_file.h b/src/parse_file.h
index 0430763..e22393e 100644
--- a/src/parse_file.h
+++ b/src/parse_file.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -23,26 +23,43 @@ Parse the input file
 
 #include "types.h"
 
+// read a positive integer from a string
+int read_positive_int(char* str, int* out);
+// read an integer from a string
+int read_int(char* str, int* out);
+// read an long int from a string
+int read_long_int(char* str, long int* out);
+// read a long double
+int read_long_double(char* str, long double* out);
+
+
 // parse fields list
 int parse_input_fields(Char_Array str_fields, Fields_Table* fields);
 
-// parse symbols list
-int parse_input_symbols(Char_Array str_symbols, Fields_Table* fields);
+// parse virtual_fields list
+int parse_input_virtual_fields(Char_Array str_virtual_fields, Fields_Table* fields, Variables variables);
 
 // parse groups of independent fields
-int parse_input_groups(Char_Array str_groups, Groups* groups);
+int parse_input_groups(Char_Array str_groups, Groups* groups, Polynomial_Matrix propagator, Fields_Table fields);
+// check that the members of groups are independent (assuming the virtual_fields and propagator were already parsed)
+int check_groups(Groups groups, Polynomial_Matrix propagator, Fields_Table fields);
+// list of fields involved in a list of virtual_fields
+int fields_in_virtual_field_list(Int_Array indices, Fields_Table fields, Int_Array* output);
+
+// parse variables list
+int parse_input_variables(Char_Array str_variables, Variables* variables);
 
 // parse identities between fields
-int parse_input_identities(Char_Array str_identities, Fields_Table* fields);
+int parse_input_identities(Char_Array str_identities, Fields_Table* fields, Variables variables);
 
 // parse propagator
 int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propagator, Fields_Table fields);
 
 // parse input polynomial
-int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields);
+int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields, Variables variables);
 
 // parse id table
-int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields);
+int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields, Variables variables);
 
 // parse a list of labels
 int parse_labels(Char_Array str_labels, Labels* labels);
diff --git a/src/polynomial.c b/src/polynomial.c
index d082b60..fae0c08 100644
--- a/src/polynomial.c
+++ b/src/polynomial.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -27,6 +27,7 @@ limitations under the License.
 #include "array.h"
 #include "number.h"
 #include "fields.h"
+#include "parse_file.h"
 
 
 // allocate memory
@@ -368,7 +369,7 @@ int polynomial_prod_chain_nosimplify(Polynomial input, Polynomial* inout, Fields
   Int_Array out_monomial;
   Int_Array out_factor;
   Number out_num;
-  // save length of inout (which changes during the loop
+  // save length of inout (which changes during the loop)
   int inout_length=(*inout).length;
   // first position in input which can multiply a term of inout without vanishing
   int firstpos;
@@ -630,8 +631,8 @@ int remove_unmatched_plusminus(Polynomial* polynomial, Fields_Table fields){
 	  match_internals--;
 	}
       }
-      // don't remove a term containing symbols
-      else if(type==FIELD_SYMBOL){
+      // don't remove a term containing virtual_field
+      else if(type==FIELD_VIRTUAL){
 	match_internals=0;
 	break;
       }
@@ -1093,7 +1094,7 @@ int polynomial_print(Polynomial polynomial){
   Char_Array buffer;
   init_Char_Array(&buffer, STR_SIZE);
   polynomial_sprint(polynomial, &buffer);
-  printf("%s",buffer.str);
+  printf("%s",char_array_to_str_noinit(&buffer));
   free_Char_Array(buffer);
   return(0);
 }
@@ -1104,7 +1105,7 @@ int polynomial_print(Polynomial polynomial){
 #define PP_MONOMIAL_MODE 2
 #define PP_FACTOR_MODE 3
 #define PP_NUMBER_MODE 4
-int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){
+int Char_Array_to_Polynomial(Char_Array str_polynomial, Polynomial* output){
   // buffer
   char* buffer=calloc(str_polynomial.length+1,sizeof(char));
   char* buffer_ptr=buffer;
@@ -1115,6 +1116,7 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){
   int comment=0;
   int i,j;
   int parenthesis_count=0;
+  int ret;
 
   // allocate memory
   init_Polynomial(output,POLY_SIZE);
@@ -1174,6 +1176,10 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){
 	  init_Int_Array(&factor, MONOMIAL_SIZE);
 	  num=number_one();
 	}
+	else{
+	  fprintf(stderr,"syntax error: misplaced '+' in polynomial\n");
+	  exit(-1);
+	}
 	break;
 
       // enter monomial or factor mode
@@ -1181,6 +1187,10 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){
 	if(mode==PP_NULL_MODE){
 	  mode=PP_BRACKET_MODE;
 	}
+	else{
+	  fprintf(stderr,"syntax error: misplaced '[' in polynomial\n");
+	  exit(-1);
+	}
 	break;
       // factor mode
       case 'l':
@@ -1189,6 +1199,10 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
 	}
+	else{
+	  fprintf(stderr,"syntax error: misplaced 'l' in polynomial\n");
+	  exit(-1);
+	}
 	break;
       // monomial mode
       case 'f':
@@ -1197,16 +1211,30 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){
 	  buffer_ptr=buffer;
 	  *buffer_ptr='\0';
 	}
+	else{
+	  fprintf(stderr,"syntax error: misplaced 'j' in polynomial\n");
+	  exit(-1);
+	}
 	break;
       // read monomial or factor
       case ']':
 	sscanf(buffer,"%d",&i);
+	ret=read_int(buffer,&i);
+	if(ret<0){
+	  fprintf(stderr,"syntax error: in polynomial, expected integer field or factor index, got '%s'\n",buffer);
+	  exit(-1);
+	}
+	  
 	if(mode==PP_FACTOR_MODE){
 	  int_array_append(i,&factor);
 	}
 	else if(mode==PP_MONOMIAL_MODE){
 	  int_array_append(i,&monomial);
 	}
+	else{
+	  fprintf(stderr,"syntax error: mismatched ']' in polynomial\n");
+	  exit(-1);
+	}
 	// switch back to null mode
 	mode=PP_NULL_MODE;
 	break;
@@ -1224,6 +1252,10 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){
 	  parenthesis_count++;
 	  buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]);
 	}
+	else{
+	  fprintf(stderr,"syntax error: misplaced '(' in polynomial\n");
+	  exit(-1);
+	}
 	break;
       case ')':
 	if(mode==PP_NUMBER_MODE){
@@ -1240,11 +1272,14 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){
 	    buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]);
 	  }
 	}
+	else{
+	  fprintf(stderr,"syntax error: mismatched ')' in polynomial\n");
+	  exit(-1);
+	}
 	break;
 
       // characters to ignore
       case ' ':break;
-      case '&':break;
       case '\n':break;
       
       // comments
@@ -1257,6 +1292,10 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){
 	  // write to buffer
 	  buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]);
 	}
+	else{
+	  fprintf(stderr,"syntax error: in polynomial, unrecognized character '%c'\n",str_polynomial.str[j]);
+	  exit(-1);
+	}
 	break;
       }
     }
@@ -1270,7 +1309,7 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){
 }
 
 // with str input
-int str_to_Polynomial(char* str_polynomial,Polynomial* output){
+int str_to_Polynomial(char* str_polynomial, Polynomial* output){
   Char_Array buffer;
   str_to_char_array(str_polynomial, &buffer);
   Char_Array_to_Polynomial(buffer, output);
@@ -1278,6 +1317,15 @@ int str_to_Polynomial(char* str_polynomial,Polynomial* output){
   return(0);
 }
 
+// check whether the polynomial is a constant
+int polynomial_is_number(Polynomial poly){
+  if(poly.length==0 || (poly.length==1 && poly.monomials[0].length==0 && poly.factors[0].length==0)){
+    return(1);
+  }
+  else{
+    return(0);
+  }
+}
 
 // -------------------- Polynomial_Matrix ---------------------
 
@@ -1307,3 +1355,16 @@ int free_Polynomial_Matrix(Polynomial_Matrix matrix){
   free(matrix.indices);
   return(0);
 }
+
+// check whether the entries are numbers
+int polynomial_matrix_is_numeric(Polynomial_Matrix matrix){
+  int i,j;
+  for(i=0;i<matrix.length;i++){
+    for(j=0;j<matrix.length;j++){
+      if(polynomial_is_number(matrix.matrix[i][j])==0){
+	return(0);
+      }
+    }
+  }
+  return(1);
+}
diff --git a/src/polynomial.h b/src/polynomial.h
index 6aff76d..f3488f5 100644
--- a/src/polynomial.h
+++ b/src/polynomial.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -124,12 +124,18 @@ int replace_factors(Grouped_Polynomial equations, Polynomial* polynomial);
 int polynomial_sprint(Polynomial polynomial, Char_Array* output);
 int polynomial_print(Polynomial polynomial);
 // read a polynomial
-int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output);
-int str_to_Polynomial(char* str_polynomial,Polynomial* output);
+int Char_Array_to_Polynomial(Char_Array str_polynomial, Polynomial* output);
+int str_to_Polynomial(char* str_polynomial, Polynomial* output);
+
+// check whether the polynomial is a constant
+int polynomial_is_number(Polynomial poly);
 
 //------------------------ Polynomial_Matrix --------------------------
 // init
 int init_Polynomial_Matrix(Polynomial_Matrix* matrix, int length);
 int free_Polynomial_Matrix(Polynomial_Matrix matrix);
 
+// check whether the entries are numbers
+int polynomial_matrix_is_numeric(Polynomial_Matrix matrix);
+
 #endif
diff --git a/src/rational.h b/src/rational.h
index 5beab8d..aea0625 100644
--- a/src/rational.h
+++ b/src/rational.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/src/rational_float.c b/src/rational_float.c
index d6a6a4c..2538d41 100644
--- a/src/rational_float.c
+++ b/src/rational_float.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/src/rational_float.h b/src/rational_float.h
index 02c4a22..ef5541b 100644
--- a/src/rational_float.h
+++ b/src/rational_float.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/src/rational_int.c b/src/rational_int.c
index 7aab561..f3e5475 100644
--- a/src/rational_int.c
+++ b/src/rational_int.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -25,6 +25,7 @@ limitations under the License.
 #include <mpfr.h>
 #include "istring.h"
 #include "array.h"
+#include "parse_file.h"
 
 Q quot(long int p, long int q){
   Q ret;
@@ -171,13 +172,18 @@ int str_to_Q(char* str, Q* num){
   int mode;
   char* buffer=calloc(str_len(str)+1,sizeof(char));
   char* buffer_ptr=buffer;
+  int ret;
 
   *num=quot(0,1);
 
   mode=PP_NUMERATOR_MODE;
   for(ptr=str;*ptr!='\0';ptr++){
     if(*ptr=='/'){
-      sscanf(buffer,"%ld",&((*num).numerator));
+      ret=read_long_int(buffer,&((*num).numerator));
+      if(ret<0){
+	fprintf(stderr,"syntax error: numerator should be an integer, got '%s' in '%s'\n",buffer,str);
+	exit(-1);
+      }
       buffer_ptr=buffer;
       *buffer_ptr='\0';
       mode=PP_DENOMINATOR_MODE;
@@ -189,10 +195,18 @@ int str_to_Q(char* str, Q* num){
 
   // last step
   if(mode==PP_NUMERATOR_MODE){
-    sscanf(buffer,"%ld",&((*num).numerator));
+    ret=read_long_int(buffer,&((*num).numerator));
+    if(ret<0){
+      fprintf(stderr,"syntax error: numerator should be an integer, got '%s' in '%s'\n",buffer,str);
+      exit(-1);
+    }
   }
   else if(mode==PP_DENOMINATOR_MODE){
-    sscanf(buffer,"%ld",&((*num).denominator));
+    ret=read_long_int(buffer,&((*num).denominator));
+    if(ret<0){
+      fprintf(stderr,"syntax error: numerator should be an integer, got '%s' in '%s'\n",buffer,str);
+      exit(-1);
+    }
   }
 
   free(buffer);
diff --git a/src/rational_int.h b/src/rational_int.h
index a4dd000..1aba15d 100644
--- a/src/rational_int.h
+++ b/src/rational_int.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/src/rcc.c b/src/rcc.c
index 484e20a..89d2c96 100644
--- a/src/rcc.c
+++ b/src/rcc.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -49,6 +49,13 @@ int RCC_cpy(RCC input,RCC* output){
   }
   return(0);
 }
+int RCC_cpy_noinit(RCC input,RCC* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    RCC_set_elem(input.values[i], input.indices[i], output, i);
+  }
+  return(0);
+}
 
 // concatenate rccs
 int RCC_concat(RCC rccs1, RCC rccs2, RCC* output){
diff --git a/src/rcc.h b/src/rcc.h
index 770309c..6beeee2 100644
--- a/src/rcc.h
+++ b/src/rcc.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -28,6 +28,7 @@ int free_RCC(RCC rccs);
 int RCC_set_elem(long double value, int index, RCC* rcc, int pos);
 // copy
 int RCC_cpy(RCC input,RCC* output);
+int RCC_cpy_noinit(RCC input,RCC* output);
 // concatenate 2 rccs
 int RCC_concat(RCC rccs1, RCC rccs2, RCC* output);
 // append an rcc to another
diff --git a/src/rcc_mpfr.c b/src/rcc_mpfr.c
index 8a362e3..7703b1e 100644
--- a/src/rcc_mpfr.c
+++ b/src/rcc_mpfr.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -65,6 +65,13 @@ int RCC_mpfr_cpy(RCC_mpfr input,RCC_mpfr* output){
   }
   return(0);
 }
+int RCC_mpfr_cpy_noinit(RCC_mpfr input,RCC_mpfr* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    RCC_mpfr_set_elem(input.values[i], input.indices[i], output, i);
+  }
+  return(0);
+}
 
 // concatenate rcc_mpfr
 int RCC_mpfr_concat(RCC_mpfr rcc_mpfr1, RCC_mpfr rcc_mpfr2, RCC_mpfr* output){
diff --git a/src/rcc_mpfr.h b/src/rcc_mpfr.h
index 18ebed3..940fbb1 100644
--- a/src/rcc_mpfr.h
+++ b/src/rcc_mpfr.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -32,6 +32,7 @@ int free_RCC_mpfr(RCC_mpfr rcc_mpfr);
 int RCC_mpfr_set_elem(mpfr_t value, int index, RCC_mpfr* rcc_mpfr, int pos);
 // copy
 int RCC_mpfr_cpy(RCC_mpfr input,RCC_mpfr* output);
+int RCC_mpfr_cpy_noinit(RCC_mpfr input,RCC_mpfr* output);
 // concatenate 2 rcc_mpfr_mpfr
 int RCC_mpfr_concat(RCC_mpfr rcc_mpfr_mpfr1, RCC_mpfr rcc_mpfr_mpfr2, RCC_mpfr* output);
 // append an rcc_mpfr to another
diff --git a/src/symbolic_parser.c b/src/symbolic_parser.c
new file mode 100644
index 0000000..08cae88
--- /dev/null
+++ b/src/symbolic_parser.c
@@ -0,0 +1,330 @@
+/*
+Copyright 2015-2022 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "symbolic_parser.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include "tree.h"
+#include "definitions.cpp"
+#include "array.h"
+#include "istring.h"
+#include "fields.h"
+#include "polynomial.h"
+
+#define SP_NULL_MODE 0
+#define SP_FUNCTION_MODE 1
+
+// parse a symbolic expression from a char_array
+int parse_symbolic_expression(Char_Array str, Fields_Table fields, Variables variables, Polynomial* polynomial){
+  Tree symbol_tree;
+  char_array_to_symbol_tree(str, &symbol_tree);
+  resolve_symbol_tree(symbol_tree, fields, variables, polynomial);
+  free_Tree(symbol_tree);
+  return(0);
+}
+// from char*
+int parse_symbolic_expression_str(char* str, Fields_Table fields, Variables variables, Polynomial* polynomial){
+  Char_Array char_array;
+  str_to_char_array(str,&char_array);
+  parse_symbolic_expression(char_array, fields, variables, polynomial);
+  free_Char_Array(char_array);
+  return(0);
+}
+
+// compute the symbol tree from a string
+int char_array_to_symbol_tree(Char_Array str, Tree* symbol_tree){
+  // buffer
+  char* buffer=calloc(str.length+1,sizeof(char));
+  char* buffer_ptr=buffer;
+  Tree child;
+  int match;
+  Char_Array nodestr;
+  Char_Array label;
+  Char_Array str_clean;
+  int mode;
+  int comment=0;
+  int j;
+  int gotanode=0;
+
+  // allocate memory
+  init_Tree(symbol_tree,SYMBOL_TREE_SIZE, SYMBOL_TREE_LABEL_SIZE);
+
+  // remove comments, ' ' and '\n'
+  init_Char_Array(&str_clean,str.length);
+  for(j=0;j<str.length;j++){
+    if(comment==1){
+      if(str.str[j]=='\n'){
+	comment=0;
+      }
+    }
+    else{
+      switch(str.str[j]){
+      case ' ':break;
+      case '\n':break;
+      // comments
+      case '#':
+	comment=1;
+	break;
+      default:
+	char_array_append(str.str[j],&str_clean);
+	break;
+      }
+    }
+  }
+
+  // if the string contains no '<', then trivial tree
+  for(j=0;j<str_clean.length;j++){
+    if(str_clean.str[j]=='<'){
+      break;
+    }
+  }
+  // no '<': trivial tree
+  if(j==str_clean.length){
+    tree_set_label(str_clean, symbol_tree);
+    free(buffer);
+    free_Char_Array(str_clean);
+    return(0);
+  }
+
+  *buffer_ptr='\0';
+  // loop over the input string
+  // start in null mode
+  mode=SP_NULL_MODE;
+  for(j=0;j<str_clean.length;j++){
+    switch(str_clean.str[j]){
+    // new node
+    case '<':
+      // find matching bracket
+      match=matching_bracket(str_clean,j);
+      // check whether it exists
+      if(match<0){
+	fprintf(stderr,"syntax error: unmatched brackets in %s\n",char_array_to_str_noinit(&str));
+	exit(-1);
+      }
+      // extract substring until bracket
+      char_array_substring(str_clean,j+1,match-1,&nodestr);
+
+      // check whether node is trivial
+      if(j==0 && match==str_clean.length-1){
+	free_Tree(*symbol_tree);
+	char_array_to_symbol_tree(nodestr, symbol_tree);
+	free_Char_Array(nodestr);
+	j=match;
+	break;
+      }
+
+      // parse subexpression
+      char_array_to_symbol_tree(nodestr, &child);
+      free_Char_Array(nodestr);
+      // add child to tree
+      tree_append_child_noinit(child, symbol_tree);
+      // boolean indicating a node has been found
+      gotanode=1;
+      // set next position after the node
+      j=match;
+
+      // if function mode, then check that the match is at the end of the node
+      if(mode==SP_FUNCTION_MODE){
+	if(match<str_clean.length-1){
+	  fprintf(stderr,"syntax error: functions must occupy an entire node (e.g. <%%exp<...>>), got %s\n",char_array_to_str_noinit(&str));
+	  exit(-1);
+	}
+	else{
+	  // set label
+	  str_to_char_array(buffer,&label);
+	  tree_set_label(label,symbol_tree);
+	  free_Char_Array(label);
+	}
+      }
+      break;
+
+    // function
+    case '%':
+      if(j>0){
+	fprintf(stderr,"syntax error: functions must occupy an entire node (e.g. <%%exp<...>>), got %s\n",char_array_to_str_noinit(&str));
+	exit(-1);
+      }
+      mode=SP_FUNCTION_MODE;
+      break;
+
+    // product
+    case '*':
+      if(gotanode==0){
+	fprintf(stderr,"syntax error: '*' is not preceded by a node in %s\n",char_array_to_str_noinit(&str));
+	exit(-1);
+      }
+      if(j>=str_clean.length-1){
+	fprintf(stderr,"syntax error: '*' cannot be at the end of an expression, got %s\n",char_array_to_str_noinit(&str));
+	exit(-1);
+      }
+      // set label
+      init_Char_Array(&label,1);
+      char_array_append('*',&label);
+      tree_set_label(label,symbol_tree);
+      free_Char_Array(label);
+      // next child
+      char_array_substring(str_clean,j+1,str_clean.length-1,&nodestr);
+      // parse subexpression
+      char_array_to_symbol_tree(nodestr, &child);
+      free_Char_Array(nodestr);
+      // append next child
+      tree_append_child_noinit(child, symbol_tree);
+
+      // make it stop
+      j=str_clean.length-1;
+      break;
+
+    // sum
+    case '+':
+      if(gotanode==0){
+	fprintf(stderr,"syntax error: '+' is not preceded by a node in %s\n",char_array_to_str_noinit(&str));
+	exit(-1);
+      }
+      if(j>=str_clean.length-1){
+	fprintf(stderr,"syntax error: '+' cannot be at the end of an expression, got %s\n",char_array_to_str_noinit(&str));
+	exit(-1);
+      }
+      // set label
+      init_Char_Array(&label,1);
+      char_array_append('+',&label);
+      tree_set_label(label,symbol_tree);
+      free_Char_Array(label);
+      // next child
+      char_array_substring(str_clean,j+1,str_clean.length-1,&nodestr);
+      // parse subexpression
+      char_array_to_symbol_tree(nodestr, &child);
+      free_Char_Array(nodestr);
+      // append next child
+      tree_append_child_noinit(child, symbol_tree);
+
+      // make it stop
+      j=str_clean.length-1;
+      break;
+
+    default:
+      if(mode!=SP_NULL_MODE){
+	// write to buffer
+	buffer_ptr=str_addchar(buffer_ptr,str_clean.str[j]);
+      }
+      break;
+    }
+  }
+
+  free_Char_Array(str_clean);
+  free(buffer);
+  return(0);
+}
+// from char*
+int str_to_symbol_tree(char* str, Tree* symbol_tree){
+  Char_Array char_array;
+  str_to_char_array(str,&char_array);
+  char_array_to_symbol_tree(char_array,symbol_tree);
+  free_Char_Array(char_array);
+  return(0);
+}
+
+
+// find matching '<' and '>'
+int matching_bracket(Char_Array str, int start){
+  int bracket_count=0;
+  int i;
+  for(i=start;i<str.length;i++){
+    if(str.str[i]=='<'){
+      bracket_count++;
+    }
+    else if(str.str[i]=='>'){
+      bracket_count--;
+      if(bracket_count==0){
+	return(i);
+      }
+    }
+  }
+  // if the function has not returned, then no matching bracket
+  return(-1);
+}
+
+
+// resolve a symbol tree to its corresponding polynomial
+int resolve_symbol_tree(Tree symbol_tree, Fields_Table fields, Variables variables, Polynomial* output){
+  Polynomial poly;
+  Tree variable_tree;
+
+  // trivial tree
+  if(symbol_tree.length==0){
+    // variable
+    if(symbol_tree.root_label.length>0 && symbol_tree.root_label.str[0]=='$'){
+      variables_find_var(symbol_tree.root_label, variables, &variable_tree);
+      resolve_symbol_tree(variable_tree, fields, variables, output);
+      free_Tree(variable_tree);
+    }
+    //polynomial
+    else{
+      Char_Array_to_Polynomial(symbol_tree.root_label, output);
+    }
+  }
+
+  // exp
+  else if (char_array_cmp_str(symbol_tree.root_label,"exp")==1){
+    if(symbol_tree.length!=1){
+      fprintf(stderr,"syntax error: exp must have 1 argument\n");
+      exit(-1);
+    }
+    resolve_symbol_tree(symbol_tree.children[0], fields, variables, &poly);
+    polynomial_exponential(poly, output, fields);
+    free_Polynomial(poly);
+  }
+
+  // log
+  else if (char_array_cmp_str(symbol_tree.root_label,"log_1")==1){
+    if(symbol_tree.length!=1){
+      fprintf(stderr,"syntax error: log_1 must have 1 argument\n");
+      exit(-1);
+    }
+    resolve_symbol_tree(symbol_tree.children[0], fields, variables, &poly);
+    polynomial_logarithm(poly, output, fields);
+    free_Polynomial(poly);
+  }
+
+  // product
+  else if (char_array_cmp_str(symbol_tree.root_label,"*")==1){
+    if(symbol_tree.length!=2){
+      fprintf(stderr,"syntax error: '*' must have 2 arguments\n");
+      exit(-1);
+    }
+    resolve_symbol_tree(symbol_tree.children[0], fields, variables, output);
+    resolve_symbol_tree(symbol_tree.children[1], fields, variables, &poly);
+    polynomial_prod_chain(poly, output, fields);
+    free_Polynomial(poly);
+  }
+
+  // sum
+  else if (char_array_cmp_str(symbol_tree.root_label,"+")==1){
+    if(symbol_tree.length!=2){
+      fprintf(stderr,"syntax error: '+' must have 2 arguments\n");
+      exit(-1);
+    }
+    resolve_symbol_tree(symbol_tree.children[0], fields, variables, output);
+    resolve_symbol_tree(symbol_tree.children[1], fields, variables, &poly);
+    polynomial_add_chain_noinit(poly, output, fields);
+  }
+
+  else{
+    fprintf(stderr,"syntax error: unrecognized operation '%s'\n",char_array_to_str_noinit(&(symbol_tree.root_label)));
+    exit(-1);
+  }
+
+  return(0);
+}
diff --git a/src/symbolic_parser.h b/src/symbolic_parser.h
new file mode 100644
index 0000000..e91c902
--- /dev/null
+++ b/src/symbolic_parser.h
@@ -0,0 +1,42 @@
+/*
+Copyright 2015-2022 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+/*
+  parse symbolic expressions
+*/
+
+#ifndef SYMBOLIC_PARSER_H
+#define SYMBOLIC_PARSER_H
+
+#include "types.h"
+
+// parse a symbolic expression from a char_array
+int parse_symbolic_expression(Char_Array str, Fields_Table fields, Variables variables, Polynomial* polynomial);
+// from char*
+int parse_symbolic_expression_str(char* str, Fields_Table fields, Variables variables, Polynomial* polynomial);
+
+// compute the symbol tree from a string
+int char_array_to_symbol_tree(Char_Array str, Tree* symbol_tree);
+// from char*
+int str_to_symbol_tree(char* str, Tree* symbol_tree);
+
+// find matching '<' and '>'
+int matching_bracket(Char_Array str, int start);
+
+// resolve a symbol tree to its corresponding polynomial
+int resolve_symbol_tree(Tree symbol_tree, Fields_Table fields, Variables variables, Polynomial* output);
+
+#endif
diff --git a/src/tools.c b/src/tools.c
index 2b7e85d..5912819 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -140,3 +140,13 @@ int min(int x1, int x2){
     return(x2);
   }
 }
+
+// check whether a divides b
+int is_factor(int a, int b){
+  if(b-a*(b/a)==0){
+    return(1);
+  }
+  else{
+    return(0);
+  }
+}
diff --git a/src/tools.h b/src/tools.h
index e5f319f..8f29cf8 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -46,4 +46,7 @@ int intlist_find_err(int* list, int size, int x);
 int max(int x1, int x2);
 int min(int x1, int x2);
 
+// check whether a divides b
+int is_factor(int a, int b);
+
 #endif
diff --git a/src/tree.c b/src/tree.c
new file mode 100644
index 0000000..0edf558
--- /dev/null
+++ b/src/tree.c
@@ -0,0 +1,117 @@
+/*
+Copyright 2015-2022 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "tree.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include "array.h"
+
+// init
+int init_Tree(Tree* tree, int memory_children, int memory_label){
+  init_Char_Array(&(tree->root_label),memory_label);
+  (*tree).children=calloc(memory_children,sizeof(Tree));
+  (*tree).memory=memory_children;
+  (*tree).length=0;
+  return(0);
+}
+int free_Tree(Tree tree){
+  int i;
+  free_Char_Array(tree.root_label);
+  for(i=0;i<tree.length;i++){
+    free_Tree(tree.children[i]);
+  }
+  free(tree.children);
+  return(0);
+}
+
+// copy
+int tree_cpy(Tree input, Tree* output){
+  init_Tree(output,input.length, input.root_label.length);
+  tree_cpy_noinit(input,output);
+  return(0);
+}
+int tree_cpy_noinit(Tree input, Tree* output){
+  int i;
+  if((*output).memory<input.length){
+    fprintf(stderr,"error: trying to copy a tree of length %d to another with memory %d\n",input.length, (*output).memory);
+    exit(-1);
+  }
+  char_array_cpy_noinit(input.root_label,&(output->root_label));
+  for(i=0;i<input.length;i++){
+    tree_cpy(input.children[i],(*output).children+i);
+  }
+
+  (*output).length=input.length;
+  return(0);
+}
+
+// resize memory
+int tree_resize(Tree* tree, int newsize){
+  Tree new_tree;
+  init_Tree(&new_tree,newsize,tree->root_label.memory);
+  tree_cpy_noinit(*tree,&new_tree);
+  free_Tree(*tree);
+  *tree=new_tree;
+  return(0);
+}
+
+// set label
+int tree_set_label(Char_Array label, Tree* tree){
+  if(label.length > tree->root_label.memory){
+    char_array_resize(&(tree->root_label),label.length);
+  }
+  char_array_cpy_noinit(label,&(tree->root_label));
+  return(0);
+}
+
+// add a child to a tree
+int tree_append_child(Tree child, Tree* output){
+  if((*output).length>=(*output).memory){
+    tree_resize(output,2*(*output).memory+1);
+  }
+  tree_cpy(child,(*output).children+(*output).length);
+  (*output).length++;
+  return(0);
+}
+// add a child to a tree without allocating memory for the new child
+int tree_append_child_noinit(Tree child, Tree* output){
+  if((*output).length>=(*output).memory){
+    tree_resize(output,2*(*output).memory+1);
+  }
+  (*output).children[(*output).length]=child;
+  (*output).length++;
+  return(0);
+}
+
+// concatenate the children of two trees
+int tree_concat_children(Tree input, Tree* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    tree_append_child(input.children[i],output);
+  }
+  return(0);
+}
+// noinit
+int tree_concat_children_noinit(Tree input, Tree* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    tree_append_child_noinit(input.children[i],output);
+  }
+  // free input array
+  free(input.children);
+  return(0);
+}
+
diff --git a/src/tree.h b/src/tree.h
new file mode 100644
index 0000000..60b222c
--- /dev/null
+++ b/src/tree.h
@@ -0,0 +1,48 @@
+/*
+Copyright 2015-2022 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+/* Trees */
+
+#ifndef TREE_H
+#define TREE_H
+
+#include "types.h"
+
+int init_Tree(Tree* tree, int memory_children, int memory_label);
+int free_Tree(Tree tree);
+
+// copy
+int tree_cpy(Tree input, Tree* output);
+int tree_cpy_noinit(Tree input, Tree* output);
+
+// resize memory
+int tree_resize(Tree* tree, int newsize);
+
+// set label
+int tree_set_label(Char_Array label, Tree* tree);
+
+// add a child to a tree
+int tree_append_child(Tree tree, Tree* output);
+// add a child to a tree without allocating memory for the new child
+int tree_append_child_noinit(Tree tree, Tree* output);
+
+// concatenate the children of two trees
+int tree_concat_children(Tree input, Tree* output);
+// noinit
+int tree_concat_children_noinit(Tree input, Tree* output);
+
+#endif
+
diff --git a/src/types.h b/src/types.h
index d30840e..00a34ac 100644
--- a/src/types.h
+++ b/src/types.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -71,6 +71,15 @@ typedef struct Str_Array{
   int memory;
 } Str_Array;
 
+// tree
+typedef struct Tree Tree;
+struct Tree{
+  Char_Array root_label;
+  Tree* children;
+  int length;
+  int memory;
+};
+
 // polynomial
 typedef struct Polynomial{
   Int_Array* monomials;
@@ -133,13 +142,21 @@ typedef struct Identities{
   int memory;
 } Identities;
 
-// symbolic expressions
-typedef struct Symbols{
+// virtual_fields
+typedef struct Virtual_fields{
   int* indices;
   Polynomial* expr;
   int length;
   int memory;
-} Symbols;
+} Virtual_fields;
+
+// variables used in symbolic expressions
+typedef struct Variables{
+  Char_Array* var_names;
+  Tree* symbol_trees;
+  int length;
+  int memory;
+} Variables;
 
 // groups of independent fields
 typedef struct Groups{
@@ -159,10 +176,10 @@ typedef struct Fields_Table{
   // identities between fields
   Identities ids;
   // symbolic expressions (commuting)
-  Symbols symbols;
-  // list of anti-commuting variables (fields or symbols)
+  Virtual_fields virtual_fields;
+  // list of anti-commuting variables
   Int_Array fermions;
-  // list of non-commuting variables (fields or symbols)
+  // list of non-commuting variables
   Int_Array noncommuting;
 } Fields_Table;
 
@@ -187,6 +204,8 @@ typedef struct Id_Table{
 typedef struct Meankondo_Options{
   int threads;
   int chain;
+  int print_progress;
+  int group_poly;
 } Meankondo_Options;
 
 typedef struct Numkondo_Options{
@@ -203,6 +222,7 @@ typedef struct Meantools_Options{
   Int_Array deriv_vars;
   Char_Array eval_rccstring;
   int chain;
+  Char_Array namespace;
   mpfr_prec_t mpfr_prec;
   mpfr_exp_t mpfr_emax;
 } Meantools_Options;