/*
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 "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"
#include "parse_file.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 virtual_field
      else if(type==FIELD_VIRTUAL){
	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],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)
// if the monomial contains noncommuting elements, put them at the beginning of the monomial
int monomial_sort(Int_Array monomial, Fields_Table fields, int* sign){
  int i,j;
  int tmp;
  // first index after noncommuting indices
  int post_nc=0;

  for(i=0;i<monomial.length;i++){
    if(is_noncommuting(monomial.values[i], fields)){
      tmp=monomial.values[i];
      for(j=i;j>post_nc;j--){
	monomial.values[j]=monomial.values[j-1];
      }
      monomial.values[post_nc]=tmp;
      post_nc++;
    }
  }

  monomial_sort_nonc(monomial, post_nc, monomial.length-1, fields, sign);

  return(0);
}
// without noncommuting terms
int monomial_sort_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign){
  int i;
  int index;
  // the pivot: middle of the monomial
  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_nonc(monomial, begin, index-1, fields, sign);
    monomial_sort_nonc(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
// if the monomial contains noncommuting elements, put them at the beginning of the monomial
int monomial_sort_groups(Int_Array monomial, Fields_Table fields, Groups groups, int* sign){
  int i,j;
  int tmp;
  // first index after noncommuting indices
  int post_nc=0;

  for(i=0;i<monomial.length;i++){
    if(is_noncommuting(monomial.values[i], fields)){
      tmp=monomial.values[i];
      for(j=post_nc;j<i;j++){
	monomial.values[j+1]=monomial.values[j];
      }
      monomial.values[post_nc]=tmp;
      post_nc++;
    }
  }

  monomial_sort_groups_nonc(monomial, post_nc, monomial.length-1, fields, groups, sign);

  return(0);
}
// without noncommuting terms
int monomial_sort_groups_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign){
  int i;
  int index;
  // the pivot: middle of the monomial
  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_groups_nonc(monomial, begin, index-1, fields, groups, sign);
    monomial_sort_groups_nonc(monomial, index+1, end, fields, groups, 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",char_array_to_str_noinit(&buffer));
  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;
  int ret;

  // 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();
	}
	else{
	  fprintf(stderr,"syntax error: misplaced '+' in polynomial\n");
	  exit(-1);
	}
	break;

      // enter monomial or factor mode
      case '[':
	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':
	if(mode==PP_BRACKET_MODE){
	  mode=PP_FACTOR_MODE;
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	}
	else{
	  fprintf(stderr,"syntax error: misplaced 'l' in polynomial\n");
	  exit(-1);
	}
	break;
      // monomial mode
      case 'f':
	if(mode==PP_BRACKET_MODE){
	  mode=PP_MONOMIAL_MODE;
	  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;
	
      // 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]);
	}
	else{
	  fprintf(stderr,"syntax error: misplaced '(' in polynomial\n");
	  exit(-1);
	}
	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]);
	  }
	}
	else{
	  fprintf(stderr,"syntax error: mismatched ')' in polynomial\n");
	  exit(-1);
	}
	break;

      // characters to ignore
      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]);
	}
	else{
	  fprintf(stderr,"syntax error: in polynomial, unrecognized character '%c'\n",str_polynomial.str[j]);
	  exit(-1);
	}
	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);
}

// 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 ---------------------

// 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);
}

// 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);
}