From aa0f3ae2988d372b190b9bde2e75a6d17e744e93 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian.jauslin@roma1.infn.it>
Date: Sun, 14 Jun 2015 00:52:45 +0000
Subject: Initial commit

---
 src/kondo.c | 1449 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 1449 insertions(+)
 create mode 100644 src/kondo.c

(limited to 'src/kondo.c')

diff --git a/src/kondo.c b/src/kondo.c
new file mode 100644
index 0000000..c86b4b3
--- /dev/null
+++ b/src/kondo.c
@@ -0,0 +1,1449 @@
+/*
+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 "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 and h (must be <10)
+#define KONDO_DIM 3
+// number of spin components
+#define KONDO_SPIN 2
+
+// offsets for indices of A, B and h
+// order matters for symbols table
+#define KONDO_A_OFFSET 1
+#define KONDO_B_OFFSET 2
+#define KONDO_H_OFFSET 3
+
+// 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;
+  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);
+
+  // symbols
+  kondo_symbols(&tmp_str, box_count, &fields);
+  arg_index=find_str_arg("symbols", *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_symbols(tmp_str, &fields);
+  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);
+  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, "symbols")==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_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);
+  for(i=0;i<KONDO_DIM;i++){
+    char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_H_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);
+
+  // parse fields table
+  parse_input_fields(*str_fields, fields);
+
+  return(0);
+}
+
+
+// generate Kondo symbols
+int kondo_symbols(Char_Array* str_symbols, 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_symbols, STR_SIZE);
+  char_array_snprintf(str_symbols, "#!symbols\n");
+
+  // loop over box index
+  for(i=1;i<=box_count;i++){
+    // loop over letters
+    for(j=0;j<2;j++){
+      // loop over space dimension
+      for(k=0;k<KONDO_DIM;k++){
+	// write index
+	char_array_snprintf(str_symbols, "%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_ABh(tmp_str.str, &poly, *fields);
+	free_Char_Array(tmp_str);
+	// write to output
+	polynomial_sprint(poly, str_symbols);
+	free_Polynomial(poly);
+	// add ,
+	char_array_snprintf(str_symbols,",\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_symbols, "%d=", 1000*(10*(KONDO_A_OFFSET+j)+KONDO_A_OFFSET+k)+i);
+	for(l=0;l<KONDO_DIM;l++){
+	  char_array_snprintf(str_symbols, "(1)");
+	  if(j<2){
+	    char_array_snprintf(str_symbols,"[f%d]", 100*(10*(KONDO_A_OFFSET+j)+l)+i);
+	  }
+	  else{
+	    char_array_snprintf(str_symbols,"[f%d]", 10*(10*(KONDO_A_OFFSET+j)+l));
+	  }
+	  if(k<2){
+	    char_array_snprintf(str_symbols,"[f%d]", 100*(10*(KONDO_A_OFFSET+k)+l)+i);
+	  }
+	  else{
+	    char_array_snprintf(str_symbols,"[f%d]", 10*(10*(KONDO_A_OFFSET+k)+l));
+	  }
+	    
+	  if(l<KONDO_DIM-1){
+	    char_array_append('+',str_symbols);
+	  }
+	}
+	// add ,
+	char_array_snprintf(str_symbols,",\n");
+      }
+    }
+  }
+
+  // vector products
+  for(i=1;i<=box_count;i++){
+    char_array_snprintf(str_symbols, "%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_symbols, "(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_symbols);
+      }
+    }
+
+    // add ,
+    if(i<box_count){
+      char_array_snprintf(str_symbols,",\n");
+    }
+  }
+
+
+  return(0);
+}
+
+// generate Kondo symbols (older method: one symbol for each scalar product)
+int kondo_symbols_scalarprod(Char_Array* str_symbols, int box_count, Fields_Table* fields){
+  int i,j,k;
+  Char_Array tmp_str;
+  Polynomial poly;
+  char letters[3]={'A','B','h'};
+
+  init_Char_Array(str_symbols, STR_SIZE);
+  char_array_snprintf(str_symbols, "#!symbols\n");
+
+  // 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_symbols, "%d=", 100*(10*(KONDO_A_OFFSET+j)+KONDO_A_OFFSET+k)+i);
+	// write the name of the scalar product
+	init_Char_Array(&tmp_str, 6);
+	char_array_snprintf(&tmp_str, "%c%d.%c%d", letters[j], i, letters[k], i);
+	// compute corresponding polynomial
+	kondo_resolve_scalar_prod(tmp_str.str, &poly, *fields);
+	free_Char_Array(tmp_str);
+	// write to output
+	polynomial_sprint(poly, str_symbols);
+	free_Polynomial(poly);
+	// add ,
+	if(i<box_count || j<2 || k<2){
+	  char_array_snprintf(str_symbols,",\n");
+	}
+      }
+    }
+  }
+
+  parse_input_symbols(*str_symbols, fields);
+
+  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
+int kondo_identities(Char_Array* str_identities){
+  int i;
+
+  init_Char_Array(str_identities,STR_SIZE);
+  char_array_snprintf(str_identities, "#!identities\n");
+
+  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_symbol_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_symbol_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){
+		if(offset1!=KONDO_H_OFFSET || offset2!=KONDO_H_OFFSET){
+		  int_array_append(1000*(10*offset1+offset2)+index, tmp_poly.monomials+i);
+		}
+	      }
+	      else{
+		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){
+	      if(offset1!=KONDO_H_OFFSET || offset2!=KONDO_H_OFFSET){
+		int_array_append(1000*(10*offset1+offset2)+index, &tmp_monomial);
+	      }
+	    }
+	    else{
+	      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);
+	  }
+	  /* // older method in which a scalar product was expanded in A, B and h
+	  // resolve scalar product
+	  kondo_resolve_scalar_prod_symbols(buffer, &scalar_prod_poly);
+	  // 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);
+	  */
+	}
+	// 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 where i is a space dimension and j is a box index
+int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){
+  char* ptr;
+  // offset (A,B or H)
+  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;
+    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
+  if(offset==KONDO_H_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);
+	// 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 or H)
+  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;
+
+    // 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
+  if(offset==KONDO_H_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);
+	  // 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 or H)
+  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;
+    // 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){
+      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){
+      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;
+    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;
+    // switch term
+    case '.':
+      term=1-term;
+      break;
+    default:
+      // char to int
+      *index=*ptr-'0';
+    }
+  }
+
+  *offset1=offset[0];
+  *offset2=offset[1];
+
+  return(0);
+}
+
+// get the index of the symbol corresponding to a given string
+int get_symbol_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;
+    default:
+      // set index if dim was already set
+      if(dim>=0){
+	index=*ptr-'0';
+      }
+      else{
+	dim=*ptr-'0';
+      }
+    }
+  }
+
+  if(offset==-1){
+    return(-1);
+  }
+  // no symbol for h
+  if(offset==KONDO_H_OFFSET){
+    return(10*(10*offset+dim));
+  }
+  else{
+    return(100*(10*offset+dim)+index);
+  }
+}
-- 
cgit v1.2.3-70-g09d2