Initial commitv1.2

1 files changed, 765 insertions, 0 deletions

diff --git a/src/grouped_polynomial.c b/src/grouped_polynomial.c
new file mode 100644
index 0000000..b7bea42
--- /dev/null
+++ b/src/grouped_polynomial.c
@@ -0,0 +1,765 @@
+/*
+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 "grouped_polynomial.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "definitions.cpp"
+#include "rational.h"
+#include "istring.h"
+#include "coefficient.h"
+#include "polynomial.h"
+#include "array.h"
+#include "number.h"
+#include "tools.h"
+
+
+// allocate memory
+int init_Grouped_Polynomial(Grouped_Polynomial* gpolynomial, int size){
+  (*gpolynomial).coefs=calloc(size,sizeof(Coefficient));
+  (*gpolynomial).indices=calloc(size,sizeof(int));
+  (*gpolynomial).length=0;
+  (*gpolynomial).memory=size;
+
+  return(0);
+}
+
+// free memory
+int free_Grouped_Polynomial(Grouped_Polynomial gpolynomial){
+  int i;
+  for(i=0;i<gpolynomial.length;i++){
+    free_Coefficient(gpolynomial.coefs[i]);
+  }
+  free(gpolynomial.coefs);
+  free(gpolynomial.indices);
+
+  return(0);
+}
+
+// resize the memory allocated to a grouped_polynomial
+int resize_grouped_polynomial(Grouped_Polynomial* grouped_polynomial,int new_size){
+  Grouped_Polynomial new_poly;
+  int i;
+
+  init_Grouped_Polynomial(&new_poly,new_size);
+  for(i=0;i<(*grouped_polynomial).length;i++){
+    new_poly.indices[i]=(*grouped_polynomial).indices[i];
+    new_poly.coefs[i]=(*grouped_polynomial).coefs[i];
+  }
+  new_poly.length=(*grouped_polynomial).length;
+
+  free((*grouped_polynomial).indices);
+  free((*grouped_polynomial).coefs);
+
+  *grouped_polynomial=new_poly;
+  return(0);
+}
+
+// copy a grouped_polynomial
+int grouped_polynomial_cpy(Grouped_Polynomial input, Grouped_Polynomial* output){
+  init_Grouped_Polynomial(output,input.length);
+  grouped_polynomial_cpy_noinit(input,output);
+  return(0);
+}
+int grouped_polynomial_cpy_noinit(Grouped_Polynomial input, Grouped_Polynomial* output){
+  int i;
+  if((*output).memory<input.length){
+    fprintf(stderr,"error: trying to copy a grouped polynomial of length %d to another with memory %d\n",input.length,(*output).memory);
+    exit(-1);
+  }
+  for(i=0;i<input.length;i++){
+    (*output).indices[i]=input.indices[i];
+    coefficient_cpy(input.coefs[i], (*output).coefs+i);
+  }
+  (*output).length=input.length;
+  
+  return(0);
+}
+
+// append an element to a grouped_polynomial
+int grouped_polynomial_append(int index, Coefficient coef, Grouped_Polynomial* output){
+  int offset=(*output).length;
+
+  if((*output).length>=(*output).memory){
+    resize_grouped_polynomial(output,2*(*output).memory+1);
+  }
+
+  // copy and allocate
+  (*output).indices[offset]=index;
+  coefficient_cpy(coef, (*output).coefs+offset);
+  //increment length
+  (*output).length++;
+
+  return(0);
+}
+// append an element to a grouped_polynomial without allocating memory
+int grouped_polynomial_append_noinit(int index, Coefficient coef, Grouped_Polynomial* output){
+  int offset=(*output).length;
+
+  if((*output).length>=(*output).memory){
+    resize_grouped_polynomial(output,2*(*output).memory+1);
+  }
+
+  // copy without allocating
+  (*output).indices[offset]=index;
+  (*output).coefs[offset]=coef;
+  // increment length
+  (*output).length++;
+  return(0);
+}
+
+// concatenate two grouped_polynomials
+int grouped_polynomial_concat(Grouped_Polynomial input, Grouped_Polynomial* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    grouped_polynomial_append(input.indices[i],input.coefs[i],output);
+  }
+  return(0);
+}
+int grouped_polynomial_concat_noinit(Grouped_Polynomial input, Grouped_Polynomial* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    grouped_polynomial_append_noinit(input.indices[i],input.coefs[i],output);
+  }
+
+  // free input arrays
+  free(input.indices);
+  free(input.coefs);
+  return(0);
+}
+
+
+// construct a grouped polynomial from a polynomial, grouping together the terms specified in the id table
+// robust algorithm: allows for a term in the polynomial to contribute to several id table entries
+int group_polynomial(Polynomial polynomial, Grouped_Polynomial* grouped_polynomial, Id_Table idtable, Fields_Table fields){
+  int i,j;
+  // what is left to group
+  Polynomial remainder;
+  int index;
+  Number ratio;
+  Number num;
+  Int_Array factor;
+  int security=0;
+  int pos;
+  coef_denom denom;
+  
+  // init remainder
+  polynomial_cpy(polynomial, &remainder);
+
+  // allocate memory
+  init_Grouped_Polynomial(grouped_polynomial, idtable.length+1);
+
+  // copy indices from idtable and allocate
+  // the constant term
+  (*grouped_polynomial).indices[0]=-1;
+  init_Coefficient((*grouped_polynomial).coefs, COEF_SIZE);
+  for(i=1;i<=idtable.length;i++){
+    (*grouped_polynomial).indices[i]=idtable.indices[i-1];
+    init_Coefficient((*grouped_polynomial).coefs+i, COEF_SIZE);
+  }
+  (*grouped_polynomial).length=idtable.length+1;
+
+  // keep on going as long as there are terms in the remainder
+  while(remainder.length>0){
+    // stop if the number of iterations exceeds 100 times the length of the polynomial
+    if(security >= 100*polynomial.length){
+      fprintf(stderr,"error: polynomial could not be grouped in less than %d groupings\n", 100*polynomial.length);
+      exit(-1);
+    }
+    security++;
+
+    // index of the last element
+    i=remainder.length-1;
+
+    // find entry
+    if(remainder.monomials[i].length==0){
+      // constant
+      index=-1;
+    }
+    else{
+      // loop over entries
+      for(j=0,index=-2;j<idtable.length && index==-2;j++){
+	// loop over terms in the polynomial
+	for(pos=0;pos<idtable.polynomials[j].length;pos++){
+	  if(int_array_cmp(idtable.polynomials[j].monomials[pos],remainder.monomials[i])==0){
+	    index=j;
+	    break;
+	  }
+	}
+      }
+    }
+
+    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){
+	  fprintf(stderr,",");
+	}
+      }
+      fprintf(stderr,") not found in idtable\n");
+      exit(-1);
+    }
+
+    // if not constant
+    if(index>=0){
+      ratio=number_quot_ret(remainder.nums[i],idtable.polynomials[index].nums[pos]);
+      factor=remainder.factors[i];
+      // add to coefficient
+      denom.index=-1;
+      denom.power=1;
+      coefficient_append(factor, ratio, denom, (*grouped_polynomial).coefs+index+1);
+
+      // remove from remainder
+      free_Int_Array(remainder.monomials[i]);
+      // do not free factor yet
+      free_Number(remainder.nums[i]);
+      remainder.length--;
+
+      // add terms from idtable with minus sign
+      for(j=0;j<idtable.polynomials[index].length;j++){
+	if(j!=pos){
+	  num=number_prod_ret(ratio, idtable.polynomials[index].nums[j]);
+	  number_Qprod_chain(quot(-1,1),&num);
+	  polynomial_append(idtable.polynomials[index].monomials[j], factor, num, &remainder);
+	  free_Number(num);
+	}
+      }
+
+      free_Int_Array(factor);
+      free_Number(ratio);
+
+      // simplify remainder
+      polynomial_simplify(&remainder, fields);
+    }
+    // constant
+    else if(index==-1){
+      // add to coefficient
+      denom.index=-1;
+      denom.power=0;
+      coefficient_append(remainder.factors[i], remainder.nums[i], denom, (*grouped_polynomial).coefs);
+      // remove from remainder
+      free_Int_Array(remainder.monomials[i]);
+      free_Int_Array(remainder.factors[i]);
+      free_Number(remainder.nums[i]);
+      remainder.length--;
+    }
+  }
+
+  // simplify the result
+  simplify_grouped_polynomial(grouped_polynomial);
+
+  free_Polynomial(remainder);
+  return(0);
+}
+
+
+// construct a grouped polynomial from a polynomial, grouping together the terms specified in the id table.
+// identifies sub-polynomials in the polynomial corresponding to the entire rhs of an entry in the id table.
+// requires the polynomial and the idtable to be sorted
+// can only treat cases in which monomials in different polynomials of the idtable are distinct
+int group_polynomial_pickandchoose(Polynomial polynomial, Grouped_Polynomial* grouped_polynomial, Id_Table idtable){
+  int i,j,k;
+  // a mask specifying which terms of the polynomial have already been grouped
+  int* mask=calloc(polynomial.length, sizeof(int));
+  int index;
+  Number ratio, ratio_check;
+  // whether ratio was ever allocated
+  int alloc_ratio=0;
+  // whether the correct index was found
+  int found_index;
+  int start_index_search;
+  Int_Array mask_tmp_flips;
+  coef_denom denom;
+  
+  // allocate memory
+  init_Grouped_Polynomial(grouped_polynomial, idtable.length+1);
+
+  // copy indices from idtable and allocate
+  // the constant term
+  (*grouped_polynomial).indices[0]=-1;
+  init_Coefficient((*grouped_polynomial).coefs, COEF_SIZE);
+  for(i=1;i<=idtable.length;i++){
+    (*grouped_polynomial).indices[i]=idtable.indices[i-1];
+    init_Coefficient((*grouped_polynomial).coefs+i, COEF_SIZE);
+  }
+  (*grouped_polynomial).length=idtable.length+1;
+
+  // loop over monomials
+  for(i=0;i<polynomial.length;i++){
+    // check that the term hasn't already been added
+    if(mask[i]==0){
+      // loop until the correct index is found (the polynomial must contain all the terms in the index and the numerical factors must match)
+      found_index=0;
+      start_index_search=0;
+      while(found_index==0){
+	found_index=1;
+	// find entry
+	index=find_id(polynomial.monomials[i], idtable,start_index_search);
+	// easier to debug if the error is here instead of inside find_id
+	if(index==-2){
+	  fprintf(stderr,"error: monomial not found in idtable\n");
+	  exit(-1);
+	}
+
+	// if not constant
+	if(index>=0){
+	  // a vector in which to store the indices that were masked
+	  init_Int_Array(&mask_tmp_flips,idtable.polynomials[index].length);
+
+	  // loop over all monomials in that entry of the idtable
+	  for(j=0;j<idtable.polynomials[index].length && found_index==1;j++){
+	    // find the monomial in the polynomial
+	    for(k=i;k<polynomial.length;k++){
+	      // only check if mask==0
+	      // only check if the factors are correct
+	      if(mask[k]==0 && int_array_cmp(polynomial.factors[i],polynomial.factors[k])==0 && int_array_cmp(idtable.polynomials[index].monomials[j],polynomial.monomials[k])==0){
+		ratio_check=number_quot_ret(polynomial.nums[k],idtable.polynomials[index].nums[j]);
+
+		// if the factors don't factor
+		if(alloc_ratio!=0 && number_compare(ratio,ratio_check)==0){
+		  found_index=0;
+		  break;
+		}
+		// check that ratio was allocated
+		if(alloc_ratio!=0){
+		  free_Number(ratio);
+		}
+		ratio=ratio_check;
+		alloc_ratio=1;
+
+		// added to polynomial
+		mask[k]=1;
+		// keep track of the flips so that they can be undone if the index turns out to be incorrect
+		int_array_append(k,&mask_tmp_flips);
+		break;
+	      }
+	    }
+
+	    // error if the monomial could not be found
+	    if(k==polynomial.length){
+	      found_index=0;
+	    }
+	  }
+
+	  // if the index was incorrect
+	  if(found_index==0){
+	    // reset mask
+	    for(j=0;j<mask_tmp_flips.length;j++){
+	      mask[mask_tmp_flips.values[j]]=0;
+	    }
+	    // start index search at next item
+	    start_index_search=index+1;
+	  }
+	  else{
+	    // add to grouped polynomial
+	    denom.index=-1;
+	    denom.power=1;
+	    coefficient_append(polynomial.factors[i], ratio, denom, (*grouped_polynomial).coefs+index+1);
+	  }
+
+	  if(alloc_ratio==1){
+	    free_Number(ratio);
+	    alloc_ratio=0;
+	  }
+	  free_Int_Array(mask_tmp_flips);
+	}
+	// constant
+	else if(index==-1){
+	  mask[i]=1;
+	  denom.index=-1;
+	  denom.power=0;
+	  coefficient_append(polynomial.factors[i], polynomial.nums[i], denom, (*grouped_polynomial).coefs);
+	}
+      }
+    }
+  }
+
+  // check all the terms were grouped
+  for(i=0;i<polynomial.length;i++){
+    if(mask[i]==0){
+      fprintf(stderr,"error: this polynomial could not be grouped: no matches were found for some of the terms\n");
+      exit(-1);
+    }
+  }
+
+  free(mask);
+  return(0);
+}
+
+
+// find the entry in the idtable containing monomial
+// start search at the specified index
+int find_id(Int_Array monomial, Id_Table idtable, int start){
+  int i,j;
+
+  // constant
+  if(monomial.length==0){
+    return(-1);
+  }
+
+  // loop over entries
+  for(i=start;i<idtable.length;i++){
+    // loop over terms in the polynomial
+    for(j=0;j<idtable.polynomials[i].length;j++){
+      if(int_array_cmp(idtable.polynomials[i].monomials[j],monomial)==0){
+	return(i);
+      }
+    }
+  }
+
+  return(-2);
+}
+
+
+// simplify grouped polynomial
+int simplify_grouped_polynomial(Grouped_Polynomial* polynomial){
+  int i;
+  for(i=0;i<(*polynomial).length;i++){
+    coefficient_simplify((*polynomial).coefs+i);
+  }
+  return(0);
+}
+
+
+// derive 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){
+  int i,j,k;
+  Coefficient tmp_coef;
+
+  // alloc
+  init_Grouped_Polynomial(dflow, flow_equation.length);
+
+  // for each equation
+  for(i=0;i<flow_equation.length;i++){
+    // copy indices
+    if(flow_equation.indices[i]>=0){
+      (*dflow).indices[i]=flow_equation.indices[i]+DOFFSET;
+    }
+    else{
+      (*dflow).indices[i]=flow_equation.indices[i]-DOFFSET;
+    }
+
+    init_Coefficient((*dflow).coefs+i, COEF_SIZE);
+    // for each index
+    for(j=0;j<indices.length;j++){
+      coefficient_deriv(flow_equation.coefs[i], indices.values[j], &tmp_coef);
+      // multiply each coefficient by the appropriate dl[j]
+      for(k=0;k<tmp_coef.length;k++){
+	// only in non-trivial cases
+	if(number_is_zero(tmp_coef.nums[k])==0){
+	  // non-constants
+	  if(indices.values[j]>=0){
+	    int_array_append(DOFFSET + indices.values[j], tmp_coef.factors+k);
+	  }
+	  // constants are offset with -doffset (so that the derivatives of constants also have a negative index)
+	  else{
+	    int_array_append(-DOFFSET + indices.values[j], tmp_coef.factors+k);
+	  }
+	}
+      }
+      
+      // add to output
+      coefficient_concat_noinit(tmp_coef, (*dflow).coefs+i);
+    }
+  }
+
+  (*dflow).length=flow_equation.length;
+
+  return(0);
+}
+
+
+/*
+// derive 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
+  Int_Array factor;
+  // number of times index was found
+  int match_count;
+  coef_denom denom;
+  // store the computation of the derivative of the constant
+  int previous_constant_index=0;
+  Coefficient dC;
+  Coefficient tmp_coef;
+
+  init_Grouped_Polynomial(output, flow_equation.length);
+
+  // loop over equations
+  for(k=0;k<flow_equation.length;k++){
+    init_Coefficient((*output).coefs+k, COEF_SIZE);
+
+    // loop over monomials
+    for(i=0;i<flow_equation.coefs[k].length;i++){
+
+      // derivative of the numerator
+      monomial_deriv(flow_equation.coefs[k].factors[i], index, &factor, &match_count);
+
+      // if the derivative doesn't vanish, add it to the coefficient
+      if(match_count>0){
+	coefficient_append_noinit(factor,number_Qprod_ret(quot(match_count,1),flow_equation.coefs[k].nums[i]), flow_equation.coefs[k].denoms[i], (*output).coefs+k);
+      }
+      else{
+	free_Int_Array(factor);
+      }
+
+      // derivative of the denominator
+      if(flow_equation.coefs[k].denoms[i].power>0){
+	// check whether the derivative was already computed
+	if(flow_equation.coefs[k].denoms[i].index!=previous_constant_index){
+	  // if not first, then free
+	  if(previous_constant_index!=0){
+	    free_Coefficient(dC);
+	    previous_constant_index=0;
+	  }
+	  init_Coefficient(&dC,COEF_SIZE);
+	  coefficient_deriv_noinit(flow_equation.coefs[intlist_find_err(flow_equation.indices, flow_equation.length, flow_equation.coefs[k].denoms[i].index)], index, &dC);
+	  previous_constant_index=flow_equation.coefs[k].denoms[i].index;
+	}
+
+	init_Coefficient(&tmp_coef, dC.length);
+	coefficient_append(flow_equation.coefs[k].factors[i], number_Qprod_ret(quot(-flow_equation.coefs[k].denoms[i].power,1), flow_equation.coefs[k].nums[i]), flow_equation.coefs[k].denoms[i], &tmp_coef);
+	(tmp_coef.denoms[0].power)++;
+
+	coefficient_prod_chain(dC, &tmp_coef);
+
+	coefficient_concat_noinit(tmp_coef, (*output).coefs+k);
+      }
+    }
+
+    // memory safe
+    if((*output).coefs[k].length>0){
+      coefficient_simplify((*output).coefs+k);
+    }
+    else{
+      // add a trivial element to the coefficient
+      init_Int_Array(&factor,0);
+      denom.index=-1;
+      denom.power=0;
+      coefficient_append_noinit(factor,number_zero(),denom,(*output).coefs+k);
+    }
+  }
+
+  free_Coefficient(dC);
+  return(0);
+}
+*/
+
+
+// print a grouped polynomial
+// prepend the indices on the left side with lhs_pre, and those on the right by rhs_pre
+int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre, char rhs_pre){
+  int i,j;
+  Char_Array buffer;
+  int dcount;
+
+  // for each equation
+  for(i=0;i<grouped_polynomial.length;i++){
+    //print lhs
+    // negative indices are constants
+    if(grouped_polynomial.indices[i]<0){
+      // count derivatives
+      dcount=-grouped_polynomial.indices[i]/DOFFSET;
+      for(j=0;j<3-dcount;j++){
+	printf(" ");
+      }
+      printf("[");
+      for(j=0;j<dcount;j++){
+	printf("d");
+      }
+      printf("C%d] =",-grouped_polynomial.indices[i]-dcount*DOFFSET);
+    }
+    else{
+      // count derivatives
+      dcount=grouped_polynomial.indices[i]/DOFFSET;
+      for(j=0;j<2-dcount;j++){
+	printf(" ");
+      }
+      printf("[");
+      for(j=0;j<dcount;j++){
+	printf("d");
+      }
+      printf("%c%2d] =",lhs_pre,grouped_polynomial.indices[i]-dcount*DOFFSET);
+    }
+
+    // rhs
+    init_Char_Array(&buffer, STR_SIZE);
+    coefficient_sprint(grouped_polynomial.coefs[i],&buffer,9,rhs_pre);
+    if(buffer.length>0){
+      printf("%s",buffer.str);
+    }
+    free_Char_Array(buffer);
+
+    if(i<grouped_polynomial.length-1){
+      printf(",");
+    }
+    // extra \n
+    printf("\n");
+  }
+  return(0);
+}
+
+// read from string
+#define PP_NULL_MODE 0
+#define PP_COEF_MODE 1
+#define PP_INDEX_MODE 3
+#define PP_COMMENT_MODE 4
+#define PP_BRACKET_MODE 5
+#define PP_CONSTANT_MODE 6
+int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output){
+  // buffer
+  char* buffer=calloc(str.length+1,sizeof(char));
+  char* buffer_ptr=buffer;
+  int index=-2;
+  Coefficient coef;
+  int i,j;
+  int mode;
+  int dcount=0;
+
+  init_Grouped_Polynomial(output, EQUATION_SIZE);
+
+  // loop over input
+  mode=PP_NULL_MODE;
+  for(j=0;j<str.length;j++){
+    if(mode==PP_COMMENT_MODE){
+      if(str.str[j]=='\n'){
+	mode=PP_NULL_MODE;
+      }
+    }
+    // stay in polynomial mode until ','
+    else if(mode==PP_COEF_MODE){
+      if(str.str[j]==','){
+	// parse polynomial
+	str_to_Coefficient(buffer, &coef);
+	// write index and polynomial
+	grouped_polynomial_append_noinit(index, coef, output);
+	mode=PP_NULL_MODE;
+      }
+      else{
+	buffer_ptr=str_addchar(buffer_ptr,str.str[j]);
+      }
+    }
+    else{
+      switch(str.str[j]){
+      // index
+      case '[':
+	if(mode==PP_NULL_MODE){
+	  mode=PP_BRACKET_MODE;
+	  buffer_ptr=buffer;
+	  *buffer_ptr='\0';
+	  // reset derivatives count
+	  dcount=0;
+	}
+	break;
+      case '%':
+	if(mode==PP_BRACKET_MODE){
+	  mode=PP_INDEX_MODE;
+	}
+	break;
+      // constant term
+      case 'C':
+	if(mode==PP_BRACKET_MODE){
+	  mode=PP_CONSTANT_MODE;
+	}
+	break;
+      // derivatives
+      case 'd':
+	if(mode==PP_BRACKET_MODE || mode==PP_INDEX_MODE || mode==PP_CONSTANT_MODE){
+	  dcount++;
+	}
+	break;
+      // write index
+      case ']':
+	sscanf(buffer,"%d",&i);
+	if(mode==PP_INDEX_MODE){
+	  index=i+dcount*DOFFSET;
+	}
+	else if(mode==PP_CONSTANT_MODE){
+	  index=-i-dcount*DOFFSET;
+	}
+	mode=PP_NULL_MODE;
+	  
+	break;
+
+      // coef mode
+      case '=':
+	if(mode==PP_NULL_MODE){
+	  buffer_ptr=buffer;
+	  *buffer_ptr='\0';
+	  mode=PP_COEF_MODE;
+	}
+	break;
+
+      // comment
+      case '#':
+	mode=PP_COMMENT_MODE;
+	break;
+
+      default:
+	if(mode!=PP_NULL_MODE){
+	  buffer_ptr=str_addchar(buffer_ptr,str.str[j]);
+	}
+	break;
+      }
+    }
+  }
+
+  // last step
+  if(mode==PP_COEF_MODE){
+    str_to_Coefficient(buffer, &coef);
+    grouped_polynomial_append_noinit(index, coef, output);
+  }
+
+  free(buffer);
+  return(0);
+}
+
+
+// evaluate an equation on a vector
+int evaleq(RCC* rccs, Grouped_Polynomial poly){
+  int i;
+  long double* res=calloc((*rccs).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);
+    exit(-1);
+  }
+
+  // initialize vectors to 0
+  for(i=0;i<(*rccs).length;i++){
+    res[i]=0.;
+  }
+
+  // for each equation
+  for(i=0;i<poly.length;i++){
+    evalcoef(*rccs, poly.coefs[i], res+i);
+  }
+
+  // copy res to rccs
+  for(i=0;i<(*rccs).length;i++){
+    (*rccs).values[i]=res[i];
+  }
+
+  // free memory
+  free(res);
+  return(0);
+
+}
+