/*
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 "kondo.h"
#include <stdlib.h>
#include <stdio.h>

#include "idtable.h"
#include "array.h"
#include "number.h"
#include "istring.h"
#include "cli_parser.h"
#include "fields.h"
#include "parse_file.h"
#include "polynomial.h"
#include "definitions.cpp"
#include "rational.h"

// dimension of A, B, h and t (must be <10)
#define KONDO_DIM 3
// number of spin components
#define KONDO_SPIN 2

// offsets for indices of A, B, h and t
// order matters for virtual_fields table
#define KONDO_A_OFFSET 1
#define KONDO_B_OFFSET 2
#define KONDO_H_OFFSET 3
#define KONDO_T_OFFSET 4

// parsing modes (from parse_file.c)
#define PP_NULL_MODE 0
// when reading a factor
#define PP_FACTOR_MODE 1
// reading a monomial
#define PP_MONOMIAL_MODE 2
// reading a numerator and denominator
#define PP_NUMBER_MODE 3
// types of fields
#define PP_FIELD_MODE 6
#define PP_PARAMETER_MODE 7
#define PP_EXTERNAL_MODE 8
#define PP_INTERNAL_MODE 9
// indices
#define PP_INDEX_MODE 10
// factors or monomials
#define PP_BRACKET_MODE 11
// labels
#define PP_LABEL_MODE 12
// polynomial
#define PP_POLYNOMIAL_MODE 13
// field scalar product
#define PP_FIELD_SCALAR_MODE 14
#define PP_FIELD_VECTOR_PROD_MODE 15



// generate configuration file
int kondo_generate_conf(Str_Array* str_args, int box_count){
  Str_Array new_args;
  Fields_Table fields;
  Variables variables;
  Char_Array tmp_str;
  int arg_index;
  int i;
  Char_Array title;

  init_Str_Array(&new_args,8);

  // fields
  kondo_fields_table(box_count, &tmp_str, &fields);
  str_array_append_noinit(tmp_str, &new_args);

  // dummy variables
  init_Variables(&variables,1);

  // virtual fields
  kondo_virtual_fields(&tmp_str, box_count, &fields);
  arg_index=find_str_arg("virtual_fields", *str_args);
  if(arg_index>=0){
    if(tmp_str.length>0){
      char_array_snprintf(&tmp_str,",\n");
    }
    char_array_concat((*str_args).strs[arg_index], &tmp_str);
  }
  parse_input_virtual_fields(tmp_str, &fields, variables);
  str_array_append_noinit(tmp_str, &new_args);

  // identities
  kondo_identities(&tmp_str);
  arg_index=find_str_arg("identities", *str_args);
  if(arg_index>=0){
    if(tmp_str.length>0){
      char_array_snprintf(&tmp_str,",\n");
    }
    char_array_concat((*str_args).strs[arg_index], &tmp_str);
  }
  parse_input_identities(tmp_str, &fields, variables);
  str_array_append_noinit(tmp_str, &new_args);

  // groups
  kondo_groups(&tmp_str, box_count);
  str_array_append_noinit(tmp_str, &new_args);


  // propagator
  arg_index=find_str_arg("propagator", *str_args);
  if(arg_index>=0){
    kondo_propagator((*str_args).strs[arg_index], &tmp_str);
    str_array_append_noinit(tmp_str, &new_args);
  }

  // input polynomial
  arg_index=find_str_arg("input_polynomial", *str_args);
  if(arg_index>=0){
    kondo_input_polynomial((*str_args).strs[arg_index], &tmp_str, fields, box_count);
    str_array_append_noinit(tmp_str, &new_args);
  }

  // id table
  arg_index=find_str_arg("id_table", *str_args);
  if(arg_index>=0){
    kondo_idtable((*str_args).strs[arg_index], &tmp_str, fields);
    str_array_append_noinit(tmp_str, &new_args);
  }

  // copy remaining entries
  for(i=0;i<(*str_args).length;i++){
    get_str_arg_title((*str_args).strs[i], &title);
    if(str_cmp(title.str, "virtual_fields")==0 &&\
       str_cmp(title.str, "identities")==0 &&\
       str_cmp(title.str, "propagator")==0 &&\
       str_cmp(title.str, "input_polynomial")==0 &&\
       str_cmp(title.str, "id_table")==0 ){
       
      char_array_cpy((*str_args).strs[i], &tmp_str);
      str_array_append_noinit(tmp_str, &new_args);
    }
    free_Char_Array(title);
  }

  free_Fields_Table(fields);
  free_Variables(variables);
  free_Str_Array(*str_args);
  *str_args=new_args;

  return(0);
}


// generate the Kondo fields table
int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fields){
  int i,j;

  init_Char_Array(str_fields,STR_SIZE);
  char_array_snprintf(str_fields, "#!fields\n");

  // external fields
  char_array_append_str("x:",str_fields);
  for(i=0;i<KONDO_SPIN;i++){
    char_array_snprintf(str_fields, "%d,%d", 10*(i+10*KONDO_A_OFFSET), 10*(i+10*KONDO_B_OFFSET));
    if(i<KONDO_SPIN-1){
      char_array_append(',',str_fields);
    }
  }
  char_array_append('\n',str_fields);

  // internal fields: A
  char_array_append_str("i:",str_fields);
  for(i=0;i<KONDO_SPIN;i++){
    for(j=1;j<=box_count;j++){
      char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_A_OFFSET)+j);
      char_array_append(',',str_fields);
    }
  }
  // B
  for(i=0;i<KONDO_SPIN;i++){
    for(j=1;j<=2;j++){
      char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_B_OFFSET)+j);
      if(i<KONDO_SPIN-1 || j<2){
	char_array_append(',',str_fields);
      }
    }
  }
  char_array_append('\n',str_fields);

  //  parameters
  char_array_append_str("h:",str_fields);
  // h
  for(i=0;i<KONDO_DIM;i++){
    char_array_snprintf(str_fields, "%d,", 10*(i+10*KONDO_H_OFFSET));
  }
  // t
  for(i=0;i<KONDO_DIM;i++){
    char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_T_OFFSET));
    if(i<KONDO_DIM-1){
      char_array_append(',', str_fields);
    }
  }
  char_array_append('\n', str_fields);


  // declare Fermions
  char_array_append_str("f:",str_fields);
  // external fields
  for(i=0;i<KONDO_SPIN;i++){
    char_array_snprintf(str_fields, "%d,%d", 10*(i+10*KONDO_A_OFFSET), 10*(i+10*KONDO_B_OFFSET));
    char_array_append(',',str_fields);
  }
  // internal fields: A
  for(i=0;i<KONDO_SPIN;i++){
    for(j=1;j<=box_count;j++){
      char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_A_OFFSET)+j);
      char_array_append(',',str_fields);
    }
  }
  // B
  for(i=0;i<KONDO_SPIN;i++){
    for(j=1;j<=2;j++){
      char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_B_OFFSET)+j);
      if(i<KONDO_SPIN-1 || j<2){
	char_array_append(',',str_fields);
      }
    }
  }
  char_array_append('\n',str_fields);

  // declare noncommuting
  char_array_append_str("a:",str_fields);
  for(i=0;i<KONDO_DIM;i++){
    char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_T_OFFSET));
    if(i<KONDO_DIM-1){
      char_array_append(',',str_fields);
    }
  }
  char_array_append('\n', str_fields);

  // parse fields table
  parse_input_fields(*str_fields, fields);

  return(0);
}


// generate Kondo virtual_fields
int kondo_virtual_fields(Char_Array* str_virtual_fields, int box_count, Fields_Table* fields){
  int i,j,k,l;
  Char_Array tmp_str;
  Polynomial poly;
  char letters[3]={'A','B','h'};

  init_Char_Array(str_virtual_fields, STR_SIZE);
  char_array_snprintf(str_virtual_fields, "#!virtual_fields\n");

  // loop over box index
  for(i=1;i<=box_count;i++){
    // loop over letters (A and B)
    for(j=0;j<2;j++){
      // loop over space dimension
      for(k=0;k<KONDO_DIM;k++){
	// write index
	char_array_snprintf(str_virtual_fields, "%d=", 100*(10*(KONDO_A_OFFSET+j)+k)+i);
	// write the name of the scalar product
	init_Char_Array(&tmp_str, 6);
	char_array_snprintf(&tmp_str, "%c%d%d", letters[j], k, i);
	// compute corresponding polynomial
	kondo_resolve_ABht(tmp_str.str, &poly, *fields);
	free_Char_Array(tmp_str);
	// write to output
	polynomial_sprint(poly, str_virtual_fields);
	free_Polynomial(poly);
	// add ,
	char_array_snprintf(str_virtual_fields,",\n");
      }
    }
  }

  // scalar products
  // loop over box index
  for(i=1;i<=box_count;i++){
    // loop over letters
    for(j=0;j<3;j++){
      for(k=0;k<3;k++){
	// write index
	char_array_snprintf(str_virtual_fields, "%d=", 1000*(10*(KONDO_A_OFFSET+j)+KONDO_A_OFFSET+k)+i);
	for(l=0;l<KONDO_DIM;l++){
	  char_array_snprintf(str_virtual_fields, "(1)");
	  if(j<2){
	    char_array_snprintf(str_virtual_fields,"[f%d]", 100*(10*(KONDO_A_OFFSET+j)+l)+i);
	  }
	  else{
	    char_array_snprintf(str_virtual_fields,"[f%d]", 10*(10*(KONDO_A_OFFSET+j)+l));
	  }
	  if(k<2){
	    char_array_snprintf(str_virtual_fields,"[f%d]", 100*(10*(KONDO_A_OFFSET+k)+l)+i);
	  }
	  else{
	    char_array_snprintf(str_virtual_fields,"[f%d]", 10*(10*(KONDO_A_OFFSET+k)+l));
	  }
	    
	  if(l<KONDO_DIM-1){
	    char_array_append('+',str_virtual_fields);
	  }
	}
	// add ,
	char_array_snprintf(str_virtual_fields,",\n");
      }
    }
  }

  // vector products
  for(i=1;i<=box_count;i++){
    char_array_snprintf(str_virtual_fields, "%d=", 100*(100*(KONDO_A_OFFSET)+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+i);
    for(l=0;l<KONDO_DIM;l++){
      // remember (-1 %3 = -1)
      char_array_snprintf(str_virtual_fields, "(1)[f%d][f%d][f%d]+(-1)[f%d][f%d][f%d]", 100*(10*KONDO_A_OFFSET+((l+1)%KONDO_DIM))+i, 100*(10*KONDO_B_OFFSET+((l+2)%KONDO_DIM))+i, 10*(10*KONDO_H_OFFSET+l), 100*(10*KONDO_A_OFFSET+((l+2)%KONDO_DIM))+i, 100*(10*KONDO_B_OFFSET+((l+1)%KONDO_DIM))+i, 10*(10*KONDO_H_OFFSET+l));
      if(l<KONDO_DIM-1){
	char_array_append('+',str_virtual_fields);
      }
    }

    // add ,
    if(i<box_count){
      char_array_snprintf(str_virtual_fields,",\n");
    }
  }


  return(0);
}


// generate Kondo groups (groups of independent variables)
int kondo_groups(Char_Array* str_groups, int box_count){
  int i,j;

  init_Char_Array(str_groups, STR_SIZE);
  char_array_snprintf(str_groups, "#!groups\n");
  char_array_append('(',str_groups);
  for(i=0;i<KONDO_DIM;i++){
    for(j=1;j<=box_count;j++){
      char_array_snprintf(str_groups, "%d",100*(10*KONDO_A_OFFSET+i)+j);
      if(j<box_count || i<KONDO_DIM-1){
	char_array_append(',',str_groups);
      }
    }
  }
  char_array_append(')',str_groups);
  char_array_append('\n',str_groups);

  char_array_append('(',str_groups);
  for(i=0;i<KONDO_DIM;i++){
    for(j=1;j<=box_count;j++){
      char_array_snprintf(str_groups, "%d",100*(10*KONDO_B_OFFSET+i)+j);
      if(j<box_count || i<KONDO_DIM-1){
	char_array_append(',',str_groups);
      }
    }
  }
  char_array_append(')',str_groups);
  char_array_append('\n',str_groups);
  return(0);
}


// generate Kondo identities
// h_3^2=1-h_1^2-h_2^2
// and Pauli matrices
int kondo_identities(Char_Array* str_identities){
  int i;

  init_Char_Array(str_identities,STR_SIZE);
  char_array_snprintf(str_identities, "#!identities\n");

  // Pauli matrices
  for(i=0;i<KONDO_DIM;i++){
    char_array_snprintf(str_identities,"[f%d][f%d]=(1),\n",10*(10*KONDO_T_OFFSET+i),10*(10*KONDO_T_OFFSET+i));
    char_array_snprintf(str_identities,"[f%d][f%d]=(s{-1})[f%d],\n",10*(10*KONDO_T_OFFSET+i),10*(10*KONDO_T_OFFSET+(i+1)%3),10*(10*KONDO_T_OFFSET+(i+2)%3));
    char_array_snprintf(str_identities,"[f%d][f%d]=((-1)s{-1})[f%d],\n",10*(10*KONDO_T_OFFSET+(i+2)%3),10*(10*KONDO_T_OFFSET+(i+1)%3),10*(10*KONDO_T_OFFSET+i));
  }

  // h
  char_array_snprintf(str_identities, "[f%d][f%d]=(1)",10*(KONDO_DIM-1+10*KONDO_H_OFFSET),10*(KONDO_DIM-1+10*KONDO_H_OFFSET));
  for(i=0;i<KONDO_DIM-1;i++){
    char_array_snprintf(str_identities, "+(-1)[f%d][f%d]",10*(i+10*KONDO_H_OFFSET),10*(i+10*KONDO_H_OFFSET));
  }

  return(0);
}


// convert the Kondo propagator
int kondo_propagator(Char_Array str_kondo_propagator, Char_Array* str_propagator){
  int i,j;
  // buffer
  char* buffer=calloc(str_kondo_propagator.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  // offset and index for each element
  int offset[2]={-1,-1};
  int index[2]={-1,-1};
  int mode;
  int comment=0;

  // allocate memory
  init_Char_Array(str_propagator,STR_SIZE);

  // reproduce the loop from parse_input_propagatore but merely copy values, and replace indices
  mode=PP_INDEX_MODE;
  for(j=0;j<str_kondo_propagator.length;j++){
    if(comment==1){
      // write comments to str
      char_array_append(str_kondo_propagator.str[j],str_propagator);
      if(str_kondo_propagator.str[j]=='\n'){
	comment=0;
      }
    }
    else{
      switch(str_kondo_propagator.str[j]){
      // indices
      case ';':
	if(mode==PP_INDEX_MODE){
	  get_offset_index(buffer,offset,index);
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	}
	break;
      case ':':
	if(mode==PP_INDEX_MODE){
	  get_offset_index(buffer,offset+1,index+1);
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	  mode=PP_NUMBER_MODE;
	}
	break;

      // num
      case ',':
	if(mode==PP_NUMBER_MODE && offset[0]>=0 && offset[1]>=0 && index[0]>=0 && index[1]>=0){
	  // write indices and num
	  for(i=0;i<KONDO_SPIN;i++){
	    char_array_snprintf(str_propagator,"%d;%d:%s,",10*(i+10*offset[0])+index[0], 10*(i+10*offset[1])+index[1], buffer);
	  }
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	  mode=PP_INDEX_MODE;
	}
	break;

      // comment
      case '#':
	comment=1;
	char_array_append(str_kondo_propagator.str[j],str_propagator);
	break;

      // ignore line breaks
      case '\n': break;

      default:
	buffer_ptr=str_addchar(buffer_ptr,str_kondo_propagator.str[j]);
	break;
      }
    }
  }

  // last step
  if(mode==PP_NUMBER_MODE){
    for(i=0;i<KONDO_SPIN;i++){
      char_array_snprintf(str_propagator,"%d;%d:%s",10*(i+10*offset[0])+index[0], 10*(i+10*offset[1])+index[1], buffer);
      if(i<KONDO_SPIN-1){
	char_array_append(',',str_propagator);
      }
    }
  }

  free(buffer);
  return(0);
  
}


// convert Kondo input polynomial
int kondo_input_polynomial(Char_Array str_kondo_polynomial, Char_Array* str_polynomial, Fields_Table fields, int box_count){
  Polynomial tmp_poly;
  Polynomial out_poly;
  Char_Array tmp_str_kondo_polynomial;
  int i;
  // whether there is a '%' in the input polynomial
  int star=0;

  init_Char_Array(str_polynomial, STR_SIZE);
  char_array_snprintf(str_polynomial, "#!input_polynomial\n");

  // check for a '%'
  for(i=0;i<str_kondo_polynomial.length;i++){
    if(str_kondo_polynomial.str[i]=='%'){
      star=1;
      break;
    }
  }

  // if there is a '%', then take a product over boxes
  if(star==1){
    // product over i from 1 to box_count
    for(i=1;i<=box_count;i++){
      // replace '%' with the appropriate index
      replace_star('0'+i,str_kondo_polynomial, &tmp_str_kondo_polynomial);
      // read polynomial
      parse_kondo_polynomial_factors(tmp_str_kondo_polynomial, &tmp_poly, fields);

      // product
      if(i==1){
	polynomial_cpy(tmp_poly,&out_poly);
      }
      else{
	polynomial_prod_chain(tmp_poly,&out_poly, fields);
      }

      free_Polynomial(tmp_poly);
      free_Char_Array(tmp_str_kondo_polynomial);
    }
  }
  // if no '%' then read polynomial as is
  else{
    parse_kondo_polynomial_factors(str_kondo_polynomial, &out_poly, fields);
  }

  // useful simplification
  remove_unmatched_plusminus(&out_poly, fields);
  polynomial_sprint(out_poly, str_polynomial);
  free_Polynomial(out_poly);
  return(0);
}


// convert the Kondo idtable
int kondo_idtable(Char_Array str_kondo_idtable, Char_Array* str_idtable, Fields_Table fields){
  int j;
  // buffer
  char* buffer=calloc(str_kondo_idtable.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  Polynomial tmp_poly;
  int mode;

  // allocate memory
  init_Char_Array(str_idtable,STR_SIZE);

  // reproduce the loop from parse_input_id_table but merely copy labels and indices, and replace Kondo polynomials
  mode=PP_INDEX_MODE;
  for(j=0;j<str_kondo_idtable.length;j++){
    // unless inside a polynomial write to output
    if(mode!=PP_POLYNOMIAL_MODE){
      char_array_append(str_kondo_idtable.str[j],str_idtable);
    }

    switch(str_kondo_idtable.str[j]){
    // end polynomial mode
    case ',':
      if(mode==PP_POLYNOMIAL_MODE){
	mode=PP_INDEX_MODE;
	// write polynomial
	parse_kondo_polynomial_str(buffer, &tmp_poly, fields);
	polynomial_sprint(tmp_poly, str_idtable);
	free_Polynomial(tmp_poly);
	char_array_append(',',str_idtable);
	}
      break;

    case ':':
      if(mode==PP_INDEX_MODE){
	mode=PP_POLYNOMIAL_MODE;
	buffer_ptr=buffer;
	*buffer_ptr='\0';
      }
      break;

    default:
      if(mode==PP_POLYNOMIAL_MODE){
	buffer_ptr=str_addchar(buffer_ptr,str_kondo_idtable.str[j]);
      }
      break;
    }
  }

  //last step
  if(mode==PP_POLYNOMIAL_MODE){
    parse_kondo_polynomial_str(buffer, &tmp_poly, fields);
    polynomial_sprint(tmp_poly, str_idtable);
    free_Polynomial(tmp_poly);
  }

  free(buffer);
  return(0);
}

// read a product of polynomials
int parse_kondo_polynomial_factors(Char_Array str_polynomial, Polynomial* output, Fields_Table fields){
  int j;
  // buffer
  char* buffer=calloc(str_polynomial.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  Polynomial tmp_poly;

  // allocate memory
  init_Polynomial(output,POLY_SIZE);

  for(j=0;j<str_polynomial.length;j++){
    switch(str_polynomial.str[j]){
    case '*':
      parse_kondo_polynomial_str(buffer, &tmp_poly, fields);
      if((*output).length==0){
	polynomial_concat(tmp_poly, output);
      }
      else{
	polynomial_prod_chain(tmp_poly, output, fields);
      }
      free_Polynomial(tmp_poly);

      buffer_ptr=buffer;
      *buffer_ptr='\0';
      break;

    default:
      buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]);
      break;
    }
  }

  //last step
  parse_kondo_polynomial_str(buffer, &tmp_poly, fields);
  if((*output).length==0){
    polynomial_concat(tmp_poly, output);
  }
  else{
    polynomial_prod_chain(tmp_poly, output, fields);
  }
  free_Polynomial(tmp_poly);

  free(buffer);
  return(0);
}


// read a kondo polynomial and convert it to a polynomial expressed in terms of the fields in the fields table
int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_Table fields){
  // input pointer
  char* polynomial_ptr;
  // buffer
  char* buffer=calloc(str_len(str_polynomial),sizeof(char));
  char* buffer_ptr=buffer;
  int mode;
  int comment=0;
  int parenthesis_count=0;
  int i;
  int offset1, offset2;
  int index;
  Polynomial tmp_poly;
  Number tmp_num, tmp1_num;
  Int_Array tmp_factor, tmp_monomial, dummy_factor;
  Polynomial scalar_prod_poly;

  // allocate memory
  init_Polynomial(output,POLY_SIZE);

  init_Polynomial(&tmp_poly,MONOMIAL_SIZE);
  tmp_num=number_one();
  init_Int_Array(&tmp_factor, MONOMIAL_SIZE);

  *buffer_ptr='\0';
  // loop over the input polynomial
  // start in null mode
  mode=PP_NULL_MODE;
  for(polynomial_ptr=str_polynomial;*polynomial_ptr!='\0';polynomial_ptr++){
    if(comment==1){
      if(*polynomial_ptr=='\n'){
	comment=0;
      }
    }
    else{
      switch(*polynomial_ptr){
      // new monomial
      case '+':
	if(mode==PP_NULL_MODE){
	  // if not a constant
	  if(tmp_poly.length>0){
	    // write num
	    polynomial_multiply_scalar(tmp_poly, tmp_num);
	    // replace factor
	    for(i=0;i<tmp_poly.length;i++){
	      free_Int_Array(tmp_poly.factors[i]);
	      int_array_cpy(tmp_factor,tmp_poly.factors+i);
	    }
	  }
	  // if constant
	  else{
	    init_Int_Array(&tmp_monomial,1);
	    polynomial_append(tmp_monomial,tmp_factor,tmp_num,&tmp_poly);
	    free_Int_Array(tmp_monomial);
	  }
	  free_Int_Array(tmp_factor);
	  free_Number(tmp_num);
	  // write polynomial
	  polynomial_concat_noinit(tmp_poly, output);
	  // reset tmp_poly
	  init_Polynomial(&tmp_poly,MONOMIAL_SIZE);
	  tmp_num=number_one();
	  init_Int_Array(&tmp_factor,MONOMIAL_SIZE);
	}
	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++;
	}
	break;
      case ')':
	if(mode==PP_NUMBER_MODE){
	  if(parenthesis_count==0){
	    // write num
	    str_to_Number(buffer,&tmp1_num);
	    number_prod_chain(tmp1_num,&tmp_num);
	    free_Number(tmp1_num);
	    // back to null mode
	    mode=PP_NULL_MODE;
	  }
	  else{
	    parenthesis_count--;
	  }
	}
	break;

      // enter 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;
      // symbol mode
      case 'f':
	if(mode==PP_BRACKET_MODE){
	  mode=PP_FIELD_MODE;
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	}
	break;
      // read factor
      case ']':
	// factor
	if(mode==PP_FACTOR_MODE){
	  sscanf(buffer,"%d",&i);
	  int_array_append(i,&tmp_factor);
	}
	// symbol
	else if(mode==PP_FIELD_MODE){
	  // if polynomial exists, add to each monomial
	  if(tmp_poly.length>0){
	    for(i=0;i<tmp_poly.length;i++){
	      int_array_append(get_virtual_field_index(buffer), tmp_poly.monomials+i);
	    }
	  }
	  // if not, create a new term in the polynomial
	  else{
	    init_Int_Array(&tmp_monomial, MONOMIAL_SIZE);
	    int_array_append(get_virtual_field_index(buffer), &tmp_monomial);
	    init_Int_Array(&dummy_factor, 1);
	    polynomial_append_noinit(tmp_monomial, dummy_factor, number_one(), &tmp_poly);
	  }
	}
	// scalar product of symbols
	else if(mode==PP_FIELD_SCALAR_MODE || mode==PP_FIELD_VECTOR_PROD_MODE){
	  get_offsets_index(buffer, &offset1, &offset2, &index);
	  // if polynomial exists, add to each monomial
	  if(tmp_poly.length>0){
	    for(i=0;i<tmp_poly.length;i++){
	      if(mode==PP_FIELD_SCALAR_MODE){
		int_array_append(1000*(10*offset1+offset2)+index, tmp_poly.monomials+i);
	      }
	      else{
		// vector product
		int_array_append(100*(100*KONDO_A_OFFSET+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+index, tmp_poly.monomials+i);
	      }
	    }
	  }
	  // if not, create a new term in the polynomial
	  else{
	    init_Int_Array(&tmp_monomial, MONOMIAL_SIZE);
	    if(mode==PP_FIELD_SCALAR_MODE){
	      int_array_append(1000*(10*offset1+offset2)+index, &tmp_monomial);
	    }
	    else{
	      // vector product
	      int_array_append(100*(100*KONDO_A_OFFSET+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+index, &tmp_monomial);
	    }
	    init_Int_Array(&dummy_factor, 1);
	    polynomial_append_noinit(tmp_monomial, dummy_factor, number_one(), &tmp_poly);
	  }
	}
	// switch back to null mode
	mode=PP_NULL_MODE;
	break;

      // symbol scalar product
      case '.':
	if(mode==PP_FIELD_MODE){
	  mode=PP_FIELD_SCALAR_MODE;
	}
	buffer_ptr=str_addchar(buffer_ptr,*polynomial_ptr);
	break;
      case 'x':
	if(mode==PP_FIELD_MODE){
	  mode=PP_FIELD_VECTOR_PROD_MODE;
	}
	buffer_ptr=str_addchar(buffer_ptr,*polynomial_ptr);
	break;

      // scalar product
      case '<':
	if(mode==PP_NULL_MODE){
	  mode=PP_MONOMIAL_MODE;
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	}
	break;
      case '>':
	if(mode==PP_MONOMIAL_MODE){
	  // resolve scalar product
	  kondo_resolve_scalar_prod(buffer, &scalar_prod_poly, fields);
	  // add to tmp_poly
	  if(tmp_poly.length==0){
	    polynomial_concat(scalar_prod_poly,&tmp_poly);
	  }
	  else{
	    polynomial_prod_chain(scalar_prod_poly,&tmp_poly,fields);
	  }
	  free_Polynomial(scalar_prod_poly);

	  mode=PP_NULL_MODE;
	}
	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,*polynomial_ptr);
	}
	break;
      }
    }
  }

  // last term
  if(tmp_poly.length>0){
    polynomial_multiply_scalar(tmp_poly,tmp_num);
    for(i=0;i<tmp_poly.length;i++){
      free_Int_Array(tmp_poly.factors[i]);
      int_array_cpy(tmp_factor,tmp_poly.factors+i);
    }
  }
  else{
    init_Int_Array(&tmp_monomial,1);
    polynomial_append(tmp_monomial,tmp_factor,tmp_num,&tmp_poly);
  }
  free_Int_Array(tmp_factor);
  free_Number(tmp_num);
  polynomial_concat_noinit(tmp_poly, output);

  // simplify
  polynomial_simplify(output, fields);

  free(buffer);
  return(0);
}

// as Char_Array
int parse_kondo_polynomial(Char_Array kondo_polynomial_str, Polynomial* polynomial, Fields_Table fields){
  char* str;
  char_array_to_str(kondo_polynomial_str, &str);
  parse_kondo_polynomial_str(str, polynomial, fields);
  free(str);
  return(0);
}


// read Aij, Bij, hi, ti where i is a space dimension and j is a box index
int kondo_resolve_ABht(char* str, Polynomial* output, Fields_Table fields){
  char* ptr;
  // offset (A,B, H or T)
  int offset=-1;
  // dimension
  int dim=-1;
  // box index
  int index=-1;
  // polynomial for each term
  Polynomial psi[KONDO_SPIN];
  Polynomial poly_conjugate;
  Int_Array monomial;
  Int_Array factor;
  Number_Matrix pauli_mat;
  int i,a,b;

  // memory
  init_Polynomial(output, MONOMIAL_SIZE);

  for(ptr=str;*ptr!='\0';ptr++){
    switch(*ptr){
    case 'A':
      offset=KONDO_A_OFFSET;
      break;
    case 'a':
      offset=KONDO_A_OFFSET;
      break;
    case 'B':
      offset=KONDO_B_OFFSET;
      break;
    case 'b':
      offset=KONDO_B_OFFSET;
      break;
    case 'h':
      offset=KONDO_H_OFFSET;
      break;
    case 't':
      offset=KONDO_T_OFFSET;
      break;
    default:
      // set index if dim was already set
      if(dim>=0){
	index=*ptr-'0';
      }
      else{
	dim=*ptr-'0';
      }
    }
  }

  // turn B3 into B2 and B4 into B1
  if(offset==KONDO_B_OFFSET){
    switch(index){
    case 3:
      index=2;
      break;
    case 4:
      index=1;
      break;
    }
  }

  // h's and t's
  if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){
    // external field
    init_Int_Array(&monomial,1);
    init_Int_Array(&factor,1);
    int_array_append(10*(dim+10*offset), &monomial);
    polynomial_append_noinit(monomial, factor, number_one(), output);
  }
  // psi's
  else{
    // construct spin indices
    for(i=0;i<KONDO_SPIN;i++){
      init_Polynomial(psi+i,2);

      // external field
      init_Int_Array(&monomial,1);
      init_Int_Array(&factor,1);
      int_array_append(10*(i+10*offset), &monomial);
      polynomial_append_noinit(monomial, factor, number_one(), psi+i);

      // internal field if applicable
      if(index>0){
	init_Int_Array(&monomial,1);
	init_Int_Array(&factor,1);

	int_array_append(10*(i+10*offset)+index, &monomial);
	polynomial_append_noinit(monomial, factor, number_one(), psi+i);
      }
    }

    // multiply by Pauli matrices
    Pauli_matrix(dim+1,&pauli_mat);
    for(a=0;a<KONDO_SPIN;a++){
      for(b=0;b<KONDO_SPIN;b++){
	polynomial_cpy(psi[b],&poly_conjugate);
	polynomial_conjugate(poly_conjugate);
	polynomial_multiply_scalar(poly_conjugate, pauli_mat.matrix[a][b]);
	polynomial_prod_chain(psi[a],&poly_conjugate,fields);
	// correct sign: psi[a]^+ should be on the left of psi[b]^-
	polynomial_multiply_Qscalar(poly_conjugate,quot(-1,1));
	// add to poly
	polynomial_concat_noinit(poly_conjugate, output);
      }
    }

    free_Number_Matrix(pauli_mat);

    // free spin indices
    for(i=0;i<KONDO_SPIN;i++){
      free_Polynomial(psi[i]);
    }
  }

  return(0);
}

#define K_VECT_PROD 1
#define K_SCALAR_PROD 2
// read a Kondo scalar product (generalized to vector products as well)
int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields){
  char* ptr;
  // offset of each term (A,B,H or T)
  int offset=-1;
  // index of each term (0,...,box_count)
  int index=0;
  int i;
  int operation=0;
  Polynomial poly_vect1[KONDO_DIM];
  Polynomial poly_vect2[KONDO_DIM];

  // memory
  init_Polynomial(output, MONOMIAL_SIZE);

  for(ptr=str;*ptr!='\0';ptr++){
    switch(*ptr){
    case 'A':
      offset=KONDO_A_OFFSET;
      break;
    case 'a':
      offset=KONDO_A_OFFSET;
      break;
    case 'B':
      offset=KONDO_B_OFFSET;
      break;
    case 'b':
      offset=KONDO_B_OFFSET;
      break;
    case 'h':
      offset=KONDO_H_OFFSET;
      break;
    case 't':
      offset=KONDO_T_OFFSET;
      break;

    // scalar product
    case '.':
      // if no previous vector product
      if(operation!=K_VECT_PROD){
	kondo_polynomial_vector(offset, index, &poly_vect1, fields);
      }
      // compute vector product
      else{
	kondo_polynomial_vector(offset, index, &poly_vect2, fields);
	kondo_polynomial_vector_product(&poly_vect1, poly_vect2, fields);
      }
      operation=K_SCALAR_PROD;
      break;

    // vector product
    case 'x':
      if(offset>=0){
	kondo_polynomial_vector(offset, index, &poly_vect1, fields);
	operation=K_VECT_PROD;
      }
      break;

    // index
    default:
      // char to int
      index=*ptr-'0';
    }
  }

  // final scalar product
  if(operation==K_SCALAR_PROD){
    if(offset>=0){
      kondo_polynomial_vector(offset, index, &poly_vect2, fields);
      kondo_polynomial_scalar_product(poly_vect1, poly_vect2, output, fields);
    }
  }

  // free memory
  for(i=0;i<KONDO_DIM;i++){
    free_Polynomial(poly_vect1[i]);
    free_Polynomial(poly_vect2[i]);
  }

  return(0);
}

// compute a scalar product of polynomial vectors
int kondo_polynomial_scalar_product(Polynomial poly_vect1[3], Polynomial poly_vect2[3], Polynomial* output, Fields_Table fields){
  int i;
  Polynomial tmp_poly;

  for(i=0;i<KONDO_DIM;i++){
    polynomial_prod(poly_vect1[i],poly_vect2[i],&tmp_poly,fields);

    // add to output
    polynomial_concat_noinit(tmp_poly, output);
  }

  polynomial_simplify(output, fields);
  
  return(0);
}

// compute a vector product of polynomial vectors
int kondo_polynomial_vector_product(Polynomial (*poly_vect1)[3], Polynomial poly_vect2[3], Fields_Table fields){
  int i;
  Polynomial out[3];
  Polynomial tmp_poly;

  for(i=0;i<3;i++){
    init_Polynomial(out+i, POLY_SIZE);

    polynomial_prod((*poly_vect1)[(i+1)%3],poly_vect2[(i+2)%3], &tmp_poly, fields);
    polynomial_concat_noinit(tmp_poly, out+i);

    polynomial_prod((*poly_vect1)[(i+2)%3],poly_vect2[(i+1)%3], &tmp_poly, fields);
    polynomial_multiply_Qscalar(tmp_poly,quot(-1,1));
    polynomial_concat_noinit(tmp_poly, out+i);
  }

  for(i=0;i<3;i++){
    free_Polynomial((*poly_vect1)[i]);
    (*poly_vect1)[i]=out[i];
  }

  return(0);
}

// compute the 3 components of a kondo vector
int kondo_polynomial_vector(int offset, int index, Polynomial (*polys)[3], Fields_Table fields){
  int i,a,b;
  // polynomial for each term
  Polynomial psi[KONDO_SPIN];
  Polynomial poly_conjugate;
  Int_Array monomial;
  Int_Array factor;
  Number_Matrix pauli_mat;

  for(i=0;i<KONDO_DIM;i++){
    // memory
    init_Polynomial((*polys)+i,POLY_SIZE);
  }

  // h's and t's
  if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){
    // construct every component field
    for(i=0;i<KONDO_DIM;i++){
      // external field
      init_Int_Array(&monomial,1);
      init_Int_Array(&factor,1);
      int_array_append(10*(i+10*offset), &monomial);
      polynomial_append_noinit(monomial, factor, number_one(), (*polys)+i);
    }
  }
  // psi's
  else{
    // construct spin indices
    for(i=0;i<KONDO_SPIN;i++){
      init_Polynomial(psi+i,2);

      // external field
      init_Int_Array(&monomial,1);
      init_Int_Array(&factor,1);
      int_array_append(10*(i+10*offset), &monomial);
      polynomial_append_noinit(monomial, factor, number_one(), psi+i);

      // internal field if applicable
      if(index>0){
	init_Int_Array(&monomial,1);
	init_Int_Array(&factor,1);

	int_array_append(10*(i+10*offset)+index, &monomial);
	polynomial_append_noinit(monomial, factor, number_one(), psi+i);
      }
    }

    // multiply by Pauli matrices
    for(i=0;i<KONDO_DIM;i++){
      Pauli_matrix(i+1,&pauli_mat);
      for(a=0;a<KONDO_SPIN;a++){
	for(b=0;b<KONDO_SPIN;b++){
	  polynomial_cpy(psi[b],&poly_conjugate);
	  polynomial_conjugate(poly_conjugate);
	  polynomial_multiply_scalar(poly_conjugate, pauli_mat.matrix[a][b]);
	  polynomial_prod_chain(psi[a],&poly_conjugate,fields);
	  // correct sign: psi[a]^+ should be on the left of psi[b]^-
	  polynomial_multiply_Qscalar(poly_conjugate,quot(-1,1));
	  // add to polys[j]
	  polynomial_concat_noinit(poly_conjugate, (*polys)+i);
	}
      }

      free_Number_Matrix(pauli_mat);
    }
    
    // free spin indices
    for(i=0;i<KONDO_SPIN;i++){
      free_Polynomial(psi[i]);
    }
  }

  return(0);
}

// read a scalar product of symbols
int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output){
  char* ptr;
  // first or second term
  int term=0;
  // offset of each term (A,B,H or T)
  int offset[2];
  // index of each term (0,...,box_count)
  int index[2]={0,0};
  Int_Array monomial;
  Int_Array factor;
  int i;

  // memory
  init_Polynomial(output, KONDO_DIM);

  for(ptr=str;*ptr!='\0';ptr++){
    switch(*ptr){
    case 'A':
      offset[term]=KONDO_A_OFFSET;
      break;
    case 'a':
      offset[term]=KONDO_A_OFFSET;
      break;
    case 'B':
      offset[term]=KONDO_B_OFFSET;
      break;
    case 'b':
      offset[term]=KONDO_B_OFFSET;
      break;
    case 'h':
      offset[term]=KONDO_H_OFFSET;
      break;
    case 't':
      offset[term]=KONDO_T_OFFSET;
      break;
    // switch term
    case '.':
      term=1-term;
      break;
    default:
      // char to int
      index[term]=*ptr-'0';
    }
  }

  // scalar product
  for(i=0;i<KONDO_DIM;i++){
    init_Int_Array(&monomial,2);
    init_Int_Array(&factor, 1);

    if(offset[0]==KONDO_H_OFFSET || offset[0]==KONDO_T_OFFSET){
      int_array_append(10*(10*offset[0]+i)+index[0], &monomial);
    }
    else{
      int_array_append(100*(10*offset[0]+i)+index[0], &monomial);
    }
    if(offset[1]==KONDO_H_OFFSET || offset[1]==KONDO_T_OFFSET){
      int_array_append(10*(10*offset[1]+i)+index[1], &monomial);
    }
    else{
      int_array_append(100*(10*offset[1]+i)+index[1], &monomial);
    }

    polynomial_append_noinit(monomial, factor, number_one(), output);
  }
  return(0);
}

// get the offset and index of a monomial term
// (e.g. A1 yields KONDO_A_OFFSET and 1)
int get_offset_index(char* str, int* offset, int* index){
  char* ptr;

  for(ptr=str;*ptr!='\0';ptr++){
    switch(*ptr){
    case 'A':
      *offset=KONDO_A_OFFSET;
      break;
    case 'a':
      *offset=KONDO_A_OFFSET;
      break;
    case 'B':
      *offset=KONDO_B_OFFSET;
      break;
    case 'b':
      *offset=KONDO_B_OFFSET;
      break;
    case 'h':
      *offset=KONDO_H_OFFSET;
      break;
    case 't':
      *offset=KONDO_T_OFFSET;
      break;
    default:
      // char to int
      *index=*ptr-'0';
    }
  }

  return(0);
}

// get the offsets and index of a scalar product
int get_offsets_index(char* str, int* offset1, int* offset2, int* index){
  int offset[2]={-1,-1};
  char* ptr;
  int term=0;

  *index=-1;

  for(ptr=str;*ptr!='\0';ptr++){
    switch(*ptr){
    case 'A':
      offset[term]=KONDO_A_OFFSET;
      break;
    case 'a':
      offset[term]=KONDO_A_OFFSET;
      break;
    case 'B':
      offset[term]=KONDO_B_OFFSET;
      break;
    case 'b':
      offset[term]=KONDO_B_OFFSET;
      break;
    case 'h':
      offset[term]=KONDO_H_OFFSET;
      break;
    case 't':
      offset[term]=KONDO_T_OFFSET;
      break;
    // switch term
    case '.':
      term=1-term;
      break;
    default:
      // char to int
      *index=*ptr-'0';
    }
  }

  *offset1=offset[0];
  *offset2=offset[1];

  // if no A's or B's, then index=0
  if((offset[0]==KONDO_H_OFFSET || offset[0]==KONDO_T_OFFSET) && (offset[1]==KONDO_H_OFFSET || offset[1]==KONDO_T_OFFSET)){
    *index=0;
  }

  return(0);
}

// get the index of the virtual_field corresponding to a given string
int get_virtual_field_index(char* str){
  char* ptr;
  int offset=-1;
  int index=0;
  int dim=-1;

  // first check whether the field already is an index
  for(ptr=str;*ptr!='\0';ptr++){
    if((*ptr-'0'>=10 || *ptr-'0'<0) && (*ptr!='-')){
      break;
    }
  }
  if(*ptr=='\0'){
    sscanf(str,"%d",&index);
    return(index);
  }

  for(ptr=str;*ptr!='\0';ptr++){
    switch(*ptr){
    case 'A':
      offset=KONDO_A_OFFSET;
      break;
    case 'a':
      offset=KONDO_A_OFFSET;
      break;
    case 'B':
      offset=KONDO_B_OFFSET;
      break;
    case 'b':
      offset=KONDO_B_OFFSET;
      break;
    case 'h':
      offset=KONDO_H_OFFSET;
      break;
    case 't':
      offset=KONDO_T_OFFSET;
      break;
    default:
      // set index if dim was already set
      if(dim>=0){
	index=*ptr-'0';
      }
      else{
	dim=*ptr-'0';
      }
    }
  }

  if(offset==-1){
    return(-1);
  }
  // no virtual field for h or t
  if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){
    return(10*(10*offset+dim));
  }
  else{
    return(100*(10*offset+dim)+index);
  }
}