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

1 files changed, 237 insertions, 8 deletions

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];