Initial commitv1.2

1 files changed, 1263 insertions, 0 deletions

diff --git a/src/polynomial.c b/src/polynomial.c
new file mode 100644
index 0000000..639728a
--- /dev/null
+++ b/src/polynomial.c
@@ -0,0 +1,1263 @@
+/*
+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 "polynomial.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "definitions.cpp"
+#include "rational.h"
+#include "tools.h"
+#include "mean.h"
+#include "coefficient.h"
+#include "istring.h"
+#include "array.h"
+#include "number.h"
+#include "fields.h"
+
+
+// allocate memory
+int init_Polynomial(Polynomial* polynomial,int size){
+  (*polynomial).monomials=calloc(size,sizeof(Int_Array));
+  (*polynomial).factors=calloc(size,sizeof(Int_Array));
+  (*polynomial).nums=calloc(size,sizeof(Number));
+  (*polynomial).length=0;
+  (*polynomial).memory=size;
+  return(0);
+}
+
+// free memory
+int free_Polynomial(Polynomial polynomial){
+  int i;
+  for(i=0;i<polynomial.length;i++){
+    free_Int_Array(polynomial.monomials[i]);
+    free_Int_Array(polynomial.factors[i]);
+    free_Number(polynomial.nums[i]);
+  }
+  free(polynomial.monomials);
+  free(polynomial.factors);
+  free(polynomial.nums);
+
+  return(0);
+}
+
+// resize the memory allocated to a polynomial
+int resize_polynomial(Polynomial* polynomial,int new_size){
+  Polynomial new_poly;
+  int i;
+
+  init_Polynomial(&new_poly,new_size);
+  for(i=0;i<(*polynomial).length;i++){
+    new_poly.monomials[i]=(*polynomial).monomials[i];
+    new_poly.factors[i]=(*polynomial).factors[i];
+    new_poly.nums[i]=(*polynomial).nums[i];
+  }
+  new_poly.length=(*polynomial).length;
+
+  free((*polynomial).monomials);
+  free((*polynomial).factors);
+  free((*polynomial).nums);
+
+  *polynomial=new_poly;
+  return(0);
+}
+
+// copy a polynomial
+int polynomial_cpy(Polynomial input, Polynomial* output){
+  init_Polynomial(output,input.length);
+  polynomial_cpy_noinit(input,output);
+  return(0);
+}
+int polynomial_cpy_noinit(Polynomial input, Polynomial* output){
+  int i;
+  if((*output).memory<input.length){
+    fprintf(stderr,"error: trying to copy a polynomial of length %d to another with memory %d\n",input.length,(*output).memory);
+    exit(-1);
+  }
+  for(i=0;i<input.length;i++){
+    int_array_cpy(input.monomials[i],(*output).monomials+i);
+    int_array_cpy(input.factors[i],(*output).factors+i);
+    number_cpy(input.nums[i],(*output).nums+i);
+  }
+  (*output).length=input.length;
+  
+  return(0);
+}
+
+// append an element to a polynomial
+int polynomial_append(Int_Array monomial, Int_Array factor, Number num, Polynomial* output){
+  int offset=(*output).length;
+
+  if((*output).length>=(*output).memory){
+    resize_polynomial(output,2*(*output).memory+1);
+  }
+
+  // copy and allocate
+  int_array_cpy(monomial,(*output).monomials+offset);
+  int_array_cpy(factor,(*output).factors+offset);
+  number_cpy(num,(*output).nums+offset);
+  //increment length
+  (*output).length++;
+
+  return(0);
+}
+// append an element to a polynomial without allocating memory
+int polynomial_append_noinit(Int_Array monomial, Int_Array factor, Number num, Polynomial* output){
+  int offset=(*output).length;
+
+  if((*output).length>=(*output).memory){
+    resize_polynomial(output,2*(*output).memory+1);
+  }
+
+  // copy without allocating
+  (*output).monomials[offset]=monomial;
+  (*output).factors[offset]=factor;
+  (*output).nums[offset]=num;
+  // increment length
+  (*output).length++;
+  return(0);
+}
+// noinit, and if there already exists an element with the same monomial and factor, then just add numbers
+int polynomial_append_noinit_inplace(Int_Array monomial, Int_Array factor, Number num, Polynomial* output){
+  int i;
+  int foundit=0;
+  for(i=0;i<(*output).length;i++){
+    if(int_array_cmp(monomial, (*output).monomials[i])==0 && int_array_cmp(factor, (*output).factors[i])==0){
+      number_add_chain(num,(*output).nums+i);
+      foundit=1;
+      free_Number(num);
+      free_Int_Array(monomial);
+      free_Int_Array(factor);
+      break;
+    }
+  }
+  if(foundit==0){
+    polynomial_append_noinit(monomial, factor, num, output);
+  }
+  return(0);
+}
+
+// concatenate two polynomials
+int polynomial_concat(Polynomial input, Polynomial* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    polynomial_append(input.monomials[i],input.factors[i],input.nums[i],output);
+  }
+  return(0);
+}
+// noinit
+int polynomial_concat_noinit(Polynomial input, Polynomial* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    polynomial_append_noinit(input.monomials[i],input.factors[i],input.nums[i],output);
+  }
+
+  // free input arrays
+  free(input.monomials);
+  free(input.factors);
+  free(input.nums);
+  return(0);
+}
+// noinit, inplace
+int polynomial_concat_noinit_inplace(Polynomial input, Polynomial* output){
+  int i;
+  for(i=0;i<input.length;i++){
+    polynomial_append_noinit_inplace(input.monomials[i],input.factors[i],input.nums[i],output);
+  }
+
+  // free input arrays
+  free(input.monomials);
+  free(input.factors);
+  free(input.nums);
+  return(0);
+}
+
+// add polynomials
+int polynomial_add_chain(Polynomial input, Polynomial* inout, Fields_Table fields){
+  polynomial_concat(input,inout);
+  polynomial_simplify(inout, fields);
+  return(0);
+}
+// add polynomials (noinit)
+int polynomial_add_chain_noinit(Polynomial input, Polynomial* inout, Fields_Table fields){
+  polynomial_concat_noinit(input,inout);
+  polynomial_simplify(inout, fields);
+  return(0);
+}
+
+// multiply a polynomial by a scalar
+int polynomial_multiply_scalar(Polynomial polynomial, Number num){
+  int i;
+  for(i=0;i<polynomial.length;i++){
+    number_prod_chain(num,polynomial.nums+i);
+  }
+  return(0);
+}
+// multiply a polynomial by a rational number
+int polynomial_multiply_Qscalar(Polynomial polynomial, Q q){
+  int i;
+  for(i=0;i<polynomial.length;i++){
+    number_Qprod_chain(q,polynomial.nums+i);
+  }
+  return(0);
+}
+
+// change the sign of the monomials in a polynomial
+int polynomial_conjugate(Polynomial polynomial){
+  int i,j;
+  for(i=0;i<polynomial.length;i++){
+    for(j=0;j<polynomial.monomials[i].length;j++){
+      polynomial.monomials[i].values[j]*=-1;
+    }
+  }
+  return(0);
+}
+
+
+// returns an initialized polynomial, equal to 1
+Polynomial polynomial_one(){
+  Polynomial ret;
+  Int_Array dummy_monomial;
+  Int_Array dummy_factor;
+
+  init_Polynomial(&ret,1);
+  init_Int_Array(&dummy_monomial,1);
+  init_Int_Array(&dummy_factor,1);
+  polynomial_append_noinit(dummy_monomial,dummy_factor,number_one(),&ret);
+
+  return(ret);
+}
+
+// returns an initialized polynomial, equal to 0
+Polynomial polynomial_zero(){
+  Polynomial ret;
+  Int_Array dummy_monomial;
+  Int_Array dummy_factor;
+
+  init_Polynomial(&ret,1);
+  init_Int_Array(&dummy_monomial,1);
+  init_Int_Array(&dummy_factor,1);
+  polynomial_append_noinit(dummy_monomial,dummy_factor,number_zero(),&ret);
+
+  return(ret);
+}
+
+// check whether a polynomial is 0
+int polynomial_is_zero(Polynomial poly){
+  int i;
+  for(i=0;i<poly.length;i++){
+    if(number_is_zero(poly.nums[i])==0){
+      return(0);
+    }
+  }
+  return(1);
+}
+
+// compute V^p
+int polynomial_power(Polynomial input_polynomial, Polynomial* output, int power, Fields_Table fields){
+  int i,j;
+  int* current_term;
+  // out buffers: since we first check whether the monomial vanishes before adding it to the output,
+  // it is more efficient to concatenate the monomials in a separate variable, and then add it to the ouput
+  // instead of incrementally adding terms to the output
+  Int_Array out_monomial;
+  int monomial_size;
+  Int_Array out_factor;
+  int factor_size;
+  Number out_num;
+
+  init_Polynomial(output, POLY_SIZE);
+
+  // trivial case
+  if(power==0){
+    init_Int_Array(&out_monomial, 1);
+    init_Int_Array(&out_factor, 1);
+    polynomial_append_noinit(out_monomial, out_factor, number_one(), output);
+    return(0);
+  }
+
+  // initialize current term
+  current_term=calloc(power,sizeof(int));
+  for(i=0;i<power;i++){
+    current_term[i]=0;
+  }
+
+  // loop over terms; the loop stops when all the pointers in the list
+  // mentioned above are at the end of the input polynomial
+  while(current_term[0]<input_polynomial.length){
+    // compute the amount of memory to allocate
+    for(i=0,monomial_size=0,factor_size=0;i<power;i++){
+      monomial_size+=input_polynomial.monomials[current_term[i]].length;
+      factor_size+=input_polynomial.factors[current_term[i]].length;
+    }
+    // allocate
+    init_Int_Array(&out_monomial, monomial_size);
+    init_Int_Array(&out_factor, factor_size);
+    out_num=number_one();
+
+    // concatenate monomial and factor
+    for(i=0;i<power;i++){
+      int_array_concat(input_polynomial.monomials[current_term[i]],&out_monomial);
+      int_array_concat(input_polynomial.factors[current_term[i]],&out_factor);
+      number_prod_chain(input_polynomial.nums[current_term[i]], &out_num);
+    }
+
+    // check whether the monomial vanishes
+    if(check_monomial(out_monomial, fields)==1){
+      // multinomial combinatorial factor n!/(m1!m2!m3!...)
+      number_Qprod_chain(quot(factorial(power),multinomial(power,current_term)),&out_num);
+      // write monomials and factors
+      polynomial_append_noinit(out_monomial, out_factor, out_num, output);
+    }
+    else{
+      free_Int_Array(out_monomial);
+      free_Int_Array(out_factor);
+      free_Number(out_num);
+    }
+
+    // move to the next term of V^p by advancing the last pointer if possible,
+    // the next to last if not, and so forth, until all the pointers have
+    // reached the end of the input polynomial
+    (current_term[power-1])++;
+    // loop until the last pointer is in an adequate place or the first
+    // pointer has reached the end
+    for(i=power-2;current_term[0]<input_polynomial.length && current_term[power-1]>=input_polynomial.length;i--){
+      // try to advance pointer
+      (current_term[i])++;
+      // if the step fails, keep on looping, if not, set all the following
+      // pointers to the latest position
+      if(current_term[i]<input_polynomial.length){
+	for(j=i+1;j<power;j++){
+	  current_term[j]=current_term[i];
+	}
+      }
+    }      
+  }
+
+  polynomial_simplify(output, fields);
+  
+  // free memory
+  free(current_term);
+  return(0);
+}
+
+
+// compute V*W
+int polynomial_prod(Polynomial input1, Polynomial input2, Polynomial* output, Fields_Table fields){
+  polynomial_cpy(input2,output);
+  polynomial_prod_chain(input1,output, fields);
+  return(0);
+}
+// chain
+int polynomial_prod_chain_nosimplify(Polynomial input, Polynomial* inout, Fields_Table fields){
+  // position in V and W
+  int pos1, pos2;
+  Int_Array out_monomial;
+  Int_Array out_factor;
+  Number out_num;
+  // 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;
+
+  // loop over terms
+  for(pos1=0;pos1<inout_length;pos1++){
+    // multiply first term of input to inout[pos1]
+    // search for the first possible product
+    firstpos=-1;
+    for(pos2=0; pos2<input.length; pos2++){
+      if(check_monomial_willnot_vanish((*inout).monomials[pos1],input.monomials[pos2],fields)==1){
+	firstpos=pos2;
+	pos2++;
+	break;
+      }
+    }
+
+    // if there was no possible product
+    if(firstpos==-1){
+      number_Qprod_chain(quot(0,1),(*inout).nums+pos1);
+    }
+    else{
+      // add other terms at the end of inout
+      for(;pos2<input.length;pos2++){
+	// check whether the term will vanish
+	if(check_monomial_willnot_vanish((*inout).monomials[pos1],input.monomials[pos2],fields)==1){
+	  // allocate
+	  init_Int_Array(&out_monomial, (*inout).monomials[pos1].length+input.monomials[pos2].length);
+	  init_Int_Array(&out_factor, (*inout).factors[pos1].length+input.factors[pos2].length);
+
+	  // concatenate monomial and factor
+	  int_array_concat((*inout).monomials[pos1],&out_monomial);
+	  int_array_concat(input.monomials[pos2],&out_monomial);
+	  int_array_concat((*inout).factors[pos1],&out_factor);
+	  int_array_concat(input.factors[pos2],&out_factor);
+	  number_prod((*inout).nums[pos1],input.nums[pos2],&out_num);
+
+	  // write monomials and factors
+	  polynomial_append_noinit(out_monomial, out_factor, out_num, inout);
+	}
+      }
+      // first term
+      int_array_concat(input.monomials[firstpos],(*inout).monomials+pos1);
+      int_array_concat(input.factors[firstpos],(*inout).factors+pos1);
+      number_prod_chain(input.nums[firstpos],(*inout).nums+pos1);
+    }
+  }
+
+  return(0);
+}
+// simplify
+int polynomial_prod_chain(Polynomial input, Polynomial* inout, Fields_Table fields){
+  polynomial_prod_chain_nosimplify(input, inout, fields);
+  polynomial_simplify(inout, fields);
+  return(0);
+}
+
+// exp(V)
+int polynomial_exponential(Polynomial input_polynomial, Polynomial* output, Fields_Table fields){
+  // a buffer for the result of a given power
+  Polynomial tmp_poly;
+  // a buffer for the previous power
+  Polynomial previous_power;
+  // power
+  int power=1;
+  Int_Array out_monomial;
+  Int_Array out_factor;
+
+  // allocate memory
+  init_Polynomial(output,POLY_SIZE);
+
+  // 1
+  init_Int_Array(&out_monomial, 1);
+  init_Int_Array(&out_factor, 1);
+  polynomial_append_noinit(out_monomial, out_factor, number_one(), output);
+  
+  while(1){
+    if(power>1){
+      if(power>33){
+	fprintf(stderr,"error: trying to take a power of a polynomial that is too high (>33)\n");
+	exit(-1);
+      }
+      // next power
+      polynomial_prod(input_polynomial,previous_power,&tmp_poly, fields);
+
+      // free
+      free_Polynomial(previous_power);
+    }
+    else{
+      polynomial_cpy(input_polynomial,&tmp_poly);
+    }
+
+    // if the power is high enough that V^p=0, then stop
+    if(tmp_poly.length==0){
+      free_Polynomial(tmp_poly);
+      break;
+    }
+
+    // copy for next power
+    polynomial_cpy(tmp_poly,&previous_power);
+
+    // 1/p!
+    polynomial_multiply_Qscalar(tmp_poly,quot(1,factorial(power)));
+    // append tmp to the output
+    polynomial_concat_noinit(tmp_poly,output);
+
+    // increase power
+    power++;
+  }
+
+  return(0);
+}
+
+
+// log(1+W)
+int polynomial_logarithm(Polynomial input_polynomial,Polynomial* output, Fields_Table fields){
+  // a buffer for the result of a given power
+  Polynomial tmp_poly;
+  // a buffer for the previous power
+  Polynomial previous_power;
+  // power
+  int power=1;
+
+  // allocate memory
+  init_Polynomial(output,POLY_SIZE);
+
+  while(1){
+    if(power>1){
+      // next power
+      polynomial_prod(input_polynomial,previous_power,&tmp_poly, fields);
+
+      // free
+      free_Polynomial(previous_power);
+    }
+    else{
+      polynomial_cpy(input_polynomial,&tmp_poly);
+    }
+
+    // if the power is high enough that V^p=0, then stop
+    if(tmp_poly.length==0){
+      free_Polynomial(tmp_poly);
+      break;
+    }
+
+    // copy for next power
+    polynomial_cpy(tmp_poly,&previous_power);
+
+    // (-1)^{p-1}/p
+    polynomial_multiply_Qscalar(tmp_poly,quot(ipower(-1,power-1),power));
+    // append tmp to the output
+    polynomial_concat_noinit(tmp_poly,output);
+
+    // increase power
+    power++;
+  }
+
+  return(0);
+}
+
+// check whether a monomial vanishes
+int check_monomial(Int_Array monomial, Fields_Table fields){
+  int i,j;
+  for(i=0;i<monomial.length;i++){
+    // check for repetitions
+    if(is_fermion(monomial.values[i], fields)==1){
+      for(j=i+1;j<monomial.length;j++){
+	if(monomial.values[j]==monomial.values[i]){
+	  return(0);
+	}
+      }
+    }
+  }
+  return(1);
+}
+// check whether the product of two monomials will vanish
+int check_monomial_willnot_vanish(Int_Array monomial1, Int_Array monomial2, Fields_Table fields){
+  int i,j;
+  for(i=0;i<monomial1.length;i++){
+    // check for repetitions
+    if(is_fermion(monomial1.values[i], fields)==1){
+      for(j=0;j<monomial2.length;j++){
+	if(monomial2.values[j]==monomial1.values[i]){
+	  return(0);
+	}
+      }
+    }
+  }
+  return(1);
+}
+
+// check whether one can add a term to a monomial without creating repetitions
+int check_monomial_addterm(Int_Array monomial, Int_Array term, Fields_Table fields){
+  int i,j;
+  for(i=0;i<term.length;i++){
+    // check for repetitions
+    if(is_fermion(term.values[i], fields)==1){
+      for(j=0;j<monomial.length;j++){
+	if(monomial.values[j]==term.values[i]){
+	  return(0);
+	}
+      }
+    }
+  }
+  return(1);
+}
+
+// check whether a monomial vanishes or has unmatched +/- fields
+int check_monomial_match(Int_Array monomial, Fields_Table fields){
+  int i,j;
+  int match=0;
+  for(i=0;i<monomial.length;i++){
+    // count match
+    if(field_type(monomial.values[i], fields)==FIELD_INTERNAL){
+      if(monomial.values[i]>0){
+	match++;
+      }
+      else if(monomial.values[i]<0){
+	match--;
+      }
+    }
+
+    // check for repetitions
+    if(is_fermion(monomial.values[i], fields)==1){
+      for(j=i+1;j<monomial.length;j++){
+	if(monomial.values[j]==monomial.values[i]){
+	  return(0);
+	}
+      }
+    }
+  }
+  if(match==0){
+    return(1);
+  }
+  else{
+    // different return codes depending on why the monomial was rejected
+    return(-1);
+  }
+}
+
+// remove terms with more plus internal fields than minus ones
+int remove_unmatched_plusminus(Polynomial* polynomial, Fields_Table fields){
+  int i,j;
+  int match_internals;
+  int type;
+  Polynomial output;
+
+  init_Polynomial(&output, (*polynomial).length);
+
+  for(i=0;i<(*polynomial).length;i++){
+    match_internals=0;
+    for(j=0;j<(*polynomial).monomials[i].length;j++){
+      // check for unmatched internal fields
+      type=field_type((*polynomial).monomials[i].values[j],fields);
+      if(type==FIELD_INTERNAL){
+	if((*polynomial).monomials[i].values[j]>0){
+	  match_internals++;
+	}
+	else if((*polynomial).monomials[i].values[j]<0){
+	  match_internals--;
+	}
+      }
+      // don't remove a term containing symbols
+      else if(type==FIELD_SYMBOL){
+	match_internals=0;
+	break;
+      }
+    }
+    if(match_internals==0){
+      polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], &output);
+    }
+  }
+
+  free_Polynomial(*polynomial);
+  *polynomial=output;
+  return(0);
+}
+
+
+// denominator of a multinomal factor: m1!m2!...
+// requires terms to be sorted
+int multinomial(int power,int* terms){
+  int multiple=1;
+  int ret=1;
+  int i;
+  // the number of numbers to be multiplied in the multinomial is
+  // equal to power-1 (the first is 1)
+  for(i=1;i<power;i++){
+    // if there is a degeneracy, then increment the multiple by
+    // which the multinomial is multiplied
+    if(terms[i-1]==terms[i]){
+      multiple++;
+    }
+    // if there is no degeneracy, reset it to 1
+    else{
+      multiple=1;
+    }
+    // multiply the result by the multiple
+    ret*=multiple;
+  }
+  return(ret);
+}
+
+
+// simplify a Polynomial
+// the fields table is there to compute the sign coming from re-arranging monomials
+int polynomial_simplify(Polynomial* polynomial, Fields_Table fields){
+  int i;
+  int monomial_cmp;
+  int factor_cmp;
+  int sign;
+  Polynomial output;
+  init_Polynomial(&output,(*polynomial).length);
+  // the combination of numerical factors
+  Number new_num;
+  init_Number(&new_num,NUMBER_SIZE);
+
+  // sort monomials and factors
+  for(i=0;i<(*polynomial).length;i++){
+    sign=1;
+    monomial_sort((*polynomial).monomials[i],0,(*polynomial).monomials[i].length-1,fields,&sign);
+    number_Qprod_chain(quot(sign,1),(*polynomial).nums+i);
+    int_array_sort((*polynomial).factors[i],0,(*polynomial).factors[i].length-1);
+  }
+
+  // resolve the identities specified in the fields table
+  resolve_ids(polynomial, fields);
+
+  // in order to perform a simplification, the list of terms must be
+  // sorted (so that terms that are proportional are next to each other)
+  polynomial_sort(polynomial,0,(*polynomial).length-1);
+
+  for(i=0;i<(*polynomial).length;i++){
+    // if the term actually exists
+    if(number_is_zero((*polynomial).nums[i])!=1){
+      // combine numerical factors
+      number_add_chain((*polynomial).nums[i], &new_num);
+    }
+    // if the numerator is 0, the previous terms that may have the same factors should still be added, hence the 'if' ends here
+
+    // if either the monomial or the factor is different from the next then add term
+    if(i<(*polynomial).length-1){
+      monomial_cmp=int_array_cmp((*polynomial).monomials[i],(*polynomial).monomials[i+1]);
+      factor_cmp=int_array_cmp((*polynomial).factors[i],(*polynomial).factors[i+1]);
+    }
+    if(i>=(*polynomial).length-1 || monomial_cmp!=0 || factor_cmp!=0 ){
+      // check that the polynomial is not trivial
+      if(number_is_zero(new_num)!=1){
+	polynomial_append((*polynomial).monomials[i],(*polynomial).factors[i],new_num,&output);
+      }
+
+      // reset new numerical factor
+      free_Number(new_num);
+      init_Number(&new_num,NUMBER_SIZE);
+    }
+  }
+
+  free_Number(new_num);
+  free_Polynomial(*polynomial);
+  *polynomial=output;
+  return(0);
+}
+
+// sort a polynomial
+// requires the monomials and factors to be sorted
+int polynomial_sort(Polynomial* polynomial, int begin, int end){
+  int i;
+  int index;
+  int monomial_cmp;
+  int factor_cmp;
+
+  // the pivot: middle of the array
+  int pivot=(begin+end)/2;
+  // if the array is non trivial
+  if(begin<end){
+    // send pivot to the end
+    exchange_polynomial_terms(pivot,end,polynomial);
+    // loop over the others
+    for(i=begin, index=begin;i<end;i++){
+      // compare with pivot
+      monomial_cmp=int_array_cmp((*polynomial).monomials[i],(*polynomial).monomials[end]);
+      factor_cmp=int_array_cmp((*polynomial).factors[i],(*polynomial).factors[end]);
+      if(monomial_cmp<0 ||  (monomial_cmp==0 && factor_cmp<0)){
+	// if smaller, exchange with reference index
+	exchange_polynomial_terms(i,index,polynomial);
+	// move reference index
+	index++;
+      }
+    }
+    // put pivot (which we had sent to the end) in the right place
+    exchange_polynomial_terms(index,end,polynomial);
+    // recurse
+    polynomial_sort(polynomial, begin, index-1);
+    polynomial_sort(polynomial, index+1, end);
+  }
+  return(0);
+}
+
+// exchange two terms (for the sorting algorithm)
+int exchange_polynomial_terms(int i, int j, Polynomial* polynomial){
+  Int_Array ptmp;
+  Number tmp;
+
+  ptmp=(*polynomial).monomials[j];
+  (*polynomial).monomials[j]=(*polynomial).monomials[i];
+  (*polynomial).monomials[i]=ptmp;
+
+  ptmp=(*polynomial).factors[j];
+  (*polynomial).factors[j]=(*polynomial).factors[i];
+  (*polynomial).factors[i]=ptmp;
+
+  tmp=(*polynomial).nums[j];
+  (*polynomial).nums[j]=(*polynomial).nums[i];
+  (*polynomial).nums[i]=tmp;
+
+  return(0);
+}
+
+// sort a monomial (with sign coming from exchanging two Fermions)
+int monomial_sort(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign){
+  int i;
+  int index;
+  // the pivot: middle of the monomial
+  int pivot=(begin+end)/2;
+
+  // if the monomial is non trivial
+  if(begin<end){
+    // send pivot to the end
+    exchange_monomial_terms(monomial, pivot, end, fields, sign);
+
+    // loop over the others
+    for(i=begin, index=begin;i<end;i++){
+      // compare with pivot
+      if(compare_monomial_terms(monomial, i, end, fields)<0){
+	// if smaller, exchange with reference index
+	exchange_monomial_terms(monomial, index, i, fields, sign);
+	// move reference index
+	index++;
+      }
+    }
+    // put pivot (which we had sent to the end) in the right place
+    exchange_monomial_terms(monomial, index, end, fields, sign);
+
+    // recurse
+    monomial_sort(monomial, begin, index-1, fields, sign);
+    monomial_sort(monomial, index+1, end, fields, sign);
+  }
+  return(0);
+}
+
+// order fields: parameter, external, internal
+int compare_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields){
+  int type1=field_type(monomial.values[pos1], fields);
+  int type2=field_type(monomial.values[pos2],fields);
+
+  if(type1 < type2){
+    return(-1);
+  }
+  else if(type1 > type2){
+    return(1);
+  }
+  
+  if(monomial.values[pos1] < monomial.values[pos2]){
+    return(-1);
+  }
+  else if (monomial.values[pos1] > monomial.values[pos2]){
+    return(1);
+  }
+
+  return(0);
+}
+
+// exchange two fields (with sign)
+int exchange_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields, int* sign){
+  int tmp=monomial.values[pos1];
+  int f1,f2;
+  int i;
+
+  monomial.values[pos1]=monomial.values[pos2];
+  monomial.values[pos2]=tmp;
+
+  // sign change
+  // only if the exchange is not trivial
+  if(pos1!=pos2){
+    f1=is_fermion(monomial.values[pos1],fields);
+    f2=is_fermion(monomial.values[pos2],fields);
+    // if both Fermions then sign
+    if(f1==1 && f2==1){
+      *sign*=-1;
+    }
+    // if only one of them is a Fermion, then count the number of Fermions between them
+    else if(f1==1 || f2==1){
+      for(i=min(pos1,pos2)+1;i<max(pos1,pos2);i++){
+	if(is_fermion(monomial.values[i],fields)==1){
+	  *sign*=-1;
+	}
+      }
+    }
+  }
+  return(0);
+}
+
+
+// sort a monomial by putting each group together
+int monomial_sort_groups(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign){
+  int i;
+  int index;
+  // the pivot: middle of the monomial
+  int pivot=(begin+end)/2;
+
+  // if the monomial is non trivial
+  if(begin<end){
+    // send pivot to the end
+    exchange_monomial_terms(monomial, pivot, end, fields, sign);
+
+    // loop over the others
+    for(i=begin, index=begin;i<end;i++){
+      // compare with pivot
+      if(compare_monomial_terms_groups(monomial, i, end, fields, groups)<0){
+	// if smaller, exchange with reference index
+	exchange_monomial_terms(monomial, index, i, fields, sign);
+	// move reference index
+	index++;
+      }
+    }
+    // put pivot (which we had sent to the end) in the right place
+    exchange_monomial_terms(monomial, index, end, fields, sign);
+
+    // recurse
+    monomial_sort(monomial, begin, index-1, fields, sign);
+    monomial_sort(monomial, index+1, end, fields, sign);
+  }
+  return(0);
+}
+
+// order fields: group, then parameter, external, internal
+int compare_monomial_terms_groups(Int_Array monomial, int pos1, int pos2, Fields_Table fields, Groups groups){
+  int group1=find_group(monomial.values[pos1], groups);
+  int group2=find_group(monomial.values[pos2], groups);
+
+  if(group1 < group2){
+    return(-1);
+  }
+  else if(group1 > group2){
+    return(1);
+  }
+  
+  int type1=field_type(monomial.values[pos1], fields);
+  int type2=field_type(monomial.values[pos2],fields);
+
+  if(type1 < type2){
+    return(-1);
+  }
+  else if(type1 > type2){
+    return(1);
+  }
+  
+  if(monomial.values[pos1] < monomial.values[pos2]){
+    return(-1);
+  }
+  else if (monomial.values[pos1] > monomial.values[pos2]){
+    return(1);
+  }
+
+  return(0);
+}
+
+
+// convert and idtable to a polynomial
+int idtable_to_polynomial(Id_Table idtable, Polynomial* polynomial){
+  int i,j;
+  int start=0;
+
+  // allocate memory
+  init_Polynomial(polynomial,POLY_SIZE);
+
+  for(i=0;i<idtable.length;i++){
+    polynomial_concat(idtable.polynomials[i],polynomial);
+    // set factors
+    for(j=start;j<(*polynomial).length;j++){
+      int_array_append(idtable.indices[i],(*polynomial).factors+j);
+    }
+    // shift start point
+    start+=idtable.polynomials[i].length;
+  }
+
+  return(0);
+}
+
+
+// replace the factors in a polynomial as prescribed by an equation in the form of a Grouped_Polynomial
+int replace_factors(Grouped_Polynomial equations, Polynomial* polynomial){
+  int i,j;
+  int index;
+  Polynomial output;
+  Coefficient coef;
+
+  init_Polynomial(&output, POLY_SIZE);
+
+  // loop over monomials
+  for(i=0;i<(*polynomial).length;i++){
+    if((*polynomial).factors[i].length>0){
+      // initialize coef to store the product of factors
+      init_Coefficient(&coef, POLY_SIZE);
+
+      // first term
+      index=intlist_find_err(equations.indices, equations.length, (*polynomial).factors[i].values[0]);
+      if(index>=0){
+	coefficient_cpy_noinit(equations.coefs[index],&coef);
+      }
+
+      // other terms
+      for(j=1;j<(*polynomial).factors[i].length;j++){
+	index=intlist_find_err(equations.indices, equations.length, (*polynomial).factors[i].values[j]);
+	if(index>=0){
+	  coefficient_prod_chain(equations.coefs[index],&coef);
+	}
+      }
+
+      // new polynomial terms
+      for(j=0;j<coef.length;j++){
+	// add to output
+	polynomial_append((*polynomial).monomials[i], coef.factors[j], number_prod_ret((*polynomial).nums[i],coef.nums[j]), &output);
+      }
+
+      // free memory
+      free_Coefficient(coef);
+    }
+    // if no factors
+    else{
+      polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], &output);
+    }
+  }
+  
+  // replace output
+  free_Polynomial(*polynomial);
+  *polynomial=output;
+
+  return(0);
+}
+
+
+// print a polynomial to a string
+int polynomial_sprint(Polynomial polynomial, Char_Array* output){
+  int i,j;
+
+  // for each monomial
+  for(i=0;i<polynomial.length;i++){
+    if(i==0){
+      char_array_snprintf(output, "  ");
+    }
+    else{
+      char_array_snprintf(output, "+ ",i);
+    }
+
+    // print num
+    char_array_append('(',output);
+    number_sprint(polynomial.nums[i],output);
+    char_array_append(')',output);
+
+    // print factors
+    for(j=0;j<polynomial.factors[i].length;j++){
+      char_array_snprintf(output,"[l%d]",polynomial.factors[i].values[j]);
+    }
+
+    // print monomials
+    for(j=0;j<polynomial.monomials[i].length;j++){
+      char_array_snprintf(output,"[f%d]",polynomial.monomials[i].values[j]);
+    }
+
+    char_array_append('\n',output);
+  }
+  return(0);
+}
+// print
+int polynomial_print(Polynomial polynomial){
+  Char_Array buffer;
+  init_Char_Array(&buffer, STR_SIZE);
+  polynomial_sprint(polynomial, &buffer);
+  printf("%s",buffer.str);
+  free_Char_Array(buffer);
+  return(0);
+}
+
+// read a polynomial from a Char_Array
+#define PP_NULL_MODE 0
+#define PP_BRACKET_MODE 1
+#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){
+  // buffer
+  char* buffer=calloc(str_polynomial.length+1,sizeof(char));
+  char* buffer_ptr=buffer;
+  Int_Array monomial;
+  Int_Array factor;
+  Number num, tmp1_num;
+  int mode;
+  int comment=0;
+  int i,j;
+  int parenthesis_count=0;
+
+  // allocate memory
+  init_Polynomial(output,POLY_SIZE);
+
+  // init
+  init_Int_Array(&monomial, MONOMIAL_SIZE);
+  init_Int_Array(&factor, MONOMIAL_SIZE);
+  num=number_one();
+
+  // if the string contains no '[', then read a number
+  for(j=0;j<str_polynomial.length;j++){
+    // ignore comments
+    if(comment==1){
+      if(str_polynomial.str[j]=='\n'){
+	comment=0;
+      }
+    }
+    else{
+      if(str_polynomial.str[j]=='['){
+	break;
+      }
+      if(str_polynomial.str[j]=='#'){
+	comment=1;
+      }
+    }
+  }
+  // no '[': read a number
+  if(j==str_polynomial.length){
+    free_Number(num);
+    char_array_to_Number(str_polynomial,&num);
+    polynomial_append_noinit(monomial, factor, num, output);
+    free(buffer);
+    return(0);
+  }
+
+  // reset comment flag
+  comment=0;
+
+  *buffer_ptr='\0';
+  // loop over the input polynomial
+  // start in null mode
+  mode=PP_NULL_MODE;
+  for(j=0;j<str_polynomial.length;j++){
+    if(comment==1){
+      if(str_polynomial.str[j]=='\n'){
+	comment=0;
+      }
+    }
+    else{
+      switch(str_polynomial.str[j]){
+      // new monomial
+      case '+':
+	if(mode==PP_NULL_MODE){
+	  polynomial_append_noinit(monomial, factor, num, output);
+	  // reset monomial, factor, num
+	  init_Int_Array(&monomial, MONOMIAL_SIZE);
+	  init_Int_Array(&factor, MONOMIAL_SIZE);
+	  num=number_one();
+	}
+	break;
+
+      // enter monomial or factor mode
+      case '[':
+	if(mode==PP_NULL_MODE){
+	  mode=PP_BRACKET_MODE;
+	}
+	break;
+      // factor mode
+      case 'l':
+	if(mode==PP_BRACKET_MODE){
+	  mode=PP_FACTOR_MODE;
+	  buffer_ptr=buffer;
+	  *buffer_ptr='\0';
+	}
+	break;
+      // monomial mode
+      case 'f':
+	if(mode==PP_BRACKET_MODE){
+	  mode=PP_MONOMIAL_MODE;
+	  buffer_ptr=buffer;
+	  *buffer_ptr='\0';
+	}
+	break;
+      // read monomial or factor
+      case ']':
+	sscanf(buffer,"%d",&i);
+	if(mode==PP_FACTOR_MODE){
+	  int_array_append(i,&factor);
+	}
+	else if(mode==PP_MONOMIAL_MODE){
+	  int_array_append(i,&monomial);
+	}
+	// switch back to null mode
+	mode=PP_NULL_MODE;
+	break;
+	
+      // numerical pre-factor
+      case '(':
+	if(mode==PP_NULL_MODE){
+	  mode=PP_NUMBER_MODE;
+	  parenthesis_count=0;
+	  buffer_ptr=buffer;
+	  *buffer_ptr='\0';
+	}
+	else if(mode==PP_NUMBER_MODE){
+	  // match parentheses
+	  parenthesis_count++;
+	  buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]);
+	}
+	break;
+      case ')':
+	if(mode==PP_NUMBER_MODE){
+	  if(parenthesis_count==0){
+	    // write num
+	    str_to_Number(buffer,&tmp1_num);
+	    number_prod_chain(tmp1_num,&num);
+	    free_Number(tmp1_num);
+	    // back to null mode
+	    mode=PP_NULL_MODE;
+	  }
+	  else{
+	    parenthesis_count--;
+	    buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]);
+	  }
+	}
+	break;
+
+      // characters to ignore
+      case ' ':break;
+      case '&':break;
+      case '\n':break;
+      
+      // comments
+      case '#':
+	comment=1;
+	break;
+
+      default:
+	if(mode!=PP_NULL_MODE){
+	  // write to buffer
+	  buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]);
+	}
+	break;
+      }
+    }
+  }
+
+  // last term
+  polynomial_append_noinit(monomial, factor, num, output);
+
+  free(buffer);
+  return(0);
+}
+
+// with str input
+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);
+  free_Char_Array(buffer);
+  return(0);
+}
+
+
+// -------------------- Polynomial_Matrix ---------------------
+
+// init
+int init_Polynomial_Matrix(Polynomial_Matrix* matrix, int length){
+  int i,j;
+  (*matrix).matrix=calloc(length,sizeof(Polynomial*));
+  (*matrix).indices=calloc(length,sizeof(int));
+  for(i=0;i<length;i++){
+    (*matrix).matrix[i]=calloc(length,sizeof(Polynomial));
+    for(j=0;j<length;j++){
+      (*matrix).matrix[i][j]=polynomial_zero();
+    }
+  }
+  (*matrix).length=length;
+  return(0);
+}
+int free_Polynomial_Matrix(Polynomial_Matrix matrix){
+  int i,j;
+  for(i=0;i<matrix.length;i++){
+    for(j=0;j<matrix.length;j++){
+      free_Polynomial(matrix.matrix[i][j]);
+    }
+    free(matrix.matrix[i]);
+  }
+  free(matrix.matrix);
+  free(matrix.indices);
+  return(0);
+}