/*
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 "parse_file.h"

#include <stdio.h>
#include <stdlib.h>
#include <mpfr.h>
#include "array.h"
#include "fields.h"
#include "rational.h"
#include "number.h"
#include "polynomial.h"
#include "rcc.h"
#include "rcc_mpfr.h"
#include "definitions.cpp"
#include "istring.h"
#include "tools.h"
#include "idtable.h"
#include "tree.h"
#include "symbolic_parser.h"


// parsing modes
#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
#define PP_FERMIONS_MODE 10
#define PP_NONCOMMUTING_MODE 11
// indices
#define PP_INDEX_MODE 12
// factors or monomials
#define PP_BRACKET_MODE 13
// labels
#define PP_LABEL_MODE 14
// polynomial
#define PP_POLYNOMIAL_MODE 15
// group
#define PP_GROUP_MODE 16


// read a positive integer from a string
int read_positive_int(char* str, int* out){
  char* ptr;
  int res;
  for(ptr=str;*ptr!='\0';ptr++){
    if(*ptr-'0'>=10 || *ptr-'0'<0){
      return(-1);
    }
  }
  res=sscanf(str,"%d",out);
  if(res!=1){
    return(-2);
  }
  return(0);
}
// read an integer from a string
int read_int(char* str, int* out){
  char* ptr;
  int res;
  for(ptr=str;*ptr!='\0';ptr++){
    if((*ptr-'0'>=10 || *ptr-'0'<0) && (*ptr!='-')){
      return(-1);
    }
  }
  res=sscanf(str,"%d",out);
  if(res!=1){
    return(-2);
  }
  return(0);
}
// read an long int from a string
int read_long_int(char* str, long int* out){
  char* ptr;
  int res;
  for(ptr=str;*ptr!='\0';ptr++){
    if((*ptr-'0'>=10 || *ptr-'0'<0) && (*ptr!='-')){
      return(-1);
    }
  }
  res=sscanf(str,"%ld",out);
  if(res!=1){
    return(-2);
  }
  return(0);
}

// read a long double
int read_long_double(char* str, long double* out){
  char* ptr;
  int res;
  for(ptr=str;*ptr!='\0';ptr++){
    if((*ptr-'0'>=10 || *ptr-'0'<0) && (*ptr!='-') && (*ptr!='e') && (*ptr!='E') && (*ptr!='.')){
      return(-1);
    }
  }
  res=sscanf(str,"%Lf",out);
  if(res!=1){
    return(-2);
  }
  return(0);
}



// parse fields list
int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
  // buffer
  char* buffer=calloc(str_fields.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  int i,j;
  int mode;
  int comment=0;
  int res;
  int expect_colon=0;

  // allocate memory
  init_Fields_Table(fields);

  // loop over input
  mode=PP_NULL_MODE;
  for(j=0;j<str_fields.length;j++){
    if(comment==1){
      if(str_fields.str[j]=='\n'){
	comment=0;
      }
    }
    else{
      switch(str_fields.str[j]){
      // parameters
      case 'h':
	if(mode==PP_NULL_MODE){
	  mode=PP_PARAMETER_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: fields entry (%d): 'h' should be at the beginning of the line\n",j);
	  exit(-1);
	}
	// flag: colon is expected immediately after letter
	expect_colon=1;
	break;
      // external fields
      case 'x':
	if(mode==PP_NULL_MODE){
	  mode=PP_EXTERNAL_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: fields entry (%d): 'x' should be at the beginning of the line\n",j);
	  exit(-1);
	}
	// flag: colon is expected immediately after letter
	expect_colon=1;
	break;
      // internal fields
      case 'i':
	if(mode==PP_NULL_MODE){
	  mode=PP_INTERNAL_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: fields entry (%d): 'i' should be at the beginning of the line\n",j);
	  exit(-1);
	}
	// flag: colon is expected immediately after letter
	expect_colon=1;
	break;
      case 'f':
	if(mode==PP_NULL_MODE){
	  mode=PP_FERMIONS_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: fields entry (%d): 'f' should be at the beginning of the line\n",j);
	  exit(-1);
	}
	// flag: colon is expected immediately after letter
	expect_colon=1;
	break;
      case 'a':
	if(mode==PP_NULL_MODE){
	  mode=PP_NONCOMMUTING_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: fields entry (%d): 'a' should be at the beginning of the line\n",j);
	  exit(-1);
	}
	// flag: colon is expected immediately after letter
	expect_colon=1;
	break;

      // reset buffer
      case ':':
	// check mode
	if(mode==PP_NULL_MODE || expect_colon==0){
	  fprintf(stderr,"syntax error: fields entry (%d): ':' must be preceded by 'h', 'x', 'i', 'f', or 'a' on the same line\n",j);
	  exit(-1);
	}
	expect_colon=0;
	buffer_ptr=buffer;
	*buffer_ptr='\0';
	break;

      // write to fields
      case ',':
	res=read_positive_int(buffer,&i);

	// check i
	if(res<0){
	  fprintf(stderr,"syntax error: fields entry (%d): fields must be non-negative integers, got '%s'\n",j,buffer);
	  exit(-1);
	}
	if(mode==PP_PARAMETER_MODE || mode==PP_EXTERNAL_MODE || mode==PP_INTERNAL_MODE){
	  if(int_array_find(i,fields->parameter)>=0 || int_array_find(i,fields->external)>=0 || int_array_find(i,fields->internal)>=0){
	    fprintf(stderr,"error: fields entry (%d): field %d is already present in the fields table\n",j,i);
	    exit(-1);
	  }
	}
	else if(mode==PP_FERMIONS_MODE){
	  if(int_array_find(i,fields->fermions)>=0){
	    fprintf(stderr,"error: fields entry (%d): field %d is already present in the Fermions table\n",j,i);
	    exit(-1);
	  }
	}    
	else if(mode==PP_NONCOMMUTING_MODE){
	  if(int_array_find(i,fields->noncommuting)>=0){
	    fprintf(stderr,"error: fields entry (%d): field %d is already present in the noncommuting table\n",j,i);
	    exit(-1);
	  }
	}    

	if(mode==PP_PARAMETER_MODE){
	  int_array_append(i,&((*fields).parameter));
	}
	else if(mode==PP_EXTERNAL_MODE){
	  int_array_append(i,&((*fields).external));
	}
	else if(mode==PP_INTERNAL_MODE){
	  int_array_append(i,&((*fields).internal));
	}
	else if(mode==PP_FERMIONS_MODE){
	  int_array_append(i,&((*fields).fermions));
	}
	else if(mode==PP_NONCOMMUTING_MODE){
	  int_array_append(i,&((*fields).noncommuting));
	}
	else{
	  fprintf(stderr,"syntax error: fields entry (%d): misplaced ','\n",j);
	  exit(-1);
	}
	buffer_ptr=buffer;
	*buffer_ptr='\0';
	break;

      // back to null mode
      case '\n':
	// ignore if in NULL_MODE
	if(mode!=PP_NULL_MODE){
	  res=read_positive_int(buffer,&i);
	  // check i
	  if(res<0){
	    fprintf(stderr,"syntax error: fields entry (%d): fields must be non-negative integers, got '%s'\n",j,buffer);
	    exit(-1);
	  }
	  if(mode==PP_PARAMETER_MODE || mode==PP_EXTERNAL_MODE || mode==PP_INTERNAL_MODE){
	    if(int_array_find(i,fields->parameter)>=0 || int_array_find(i,fields->external)>=0 || int_array_find(i,fields->internal)>=0){
	      fprintf(stderr,"error: fields entry (%d): field %d is already present in the fields table\n",j,i);
	      exit(-1);
	    }
	  }
	  else if(mode==PP_FERMIONS_MODE){
	    if(int_array_find(i,fields->fermions)>=0){
	      fprintf(stderr,"error: fields entry (%d): field %d is already present in the Fermions table\n",j,i);
	      exit(-1);
	    }
	  }    
	  else if(mode==PP_NONCOMMUTING_MODE){
	    if(int_array_find(i,fields->noncommuting)>=0){
	      fprintf(stderr,"error: fields entry (%d): field %d is already present in the noncommuting table\n",j,i);
	      exit(-1);
	    }
	  }    

	  if(mode==PP_PARAMETER_MODE){
	    int_array_append(i,&((*fields).parameter));
	  }
	  else if(mode==PP_EXTERNAL_MODE){
	    int_array_append(i,&((*fields).external));
	  }
	  else if(mode==PP_INTERNAL_MODE){
	    int_array_append(i,&((*fields).internal));
	  }
	  else if(mode==PP_FERMIONS_MODE){
	    int_array_append(i,&((*fields).fermions));
	  }
	  else if(mode==PP_NONCOMMUTING_MODE){
	    int_array_append(i,&((*fields).noncommuting));
	  }
	  mode=PP_NULL_MODE;
	}
	break;

  
      // ignore spaces
      case ' ':break;
      // comment
      case '#':
	comment=1;
	break;

      default:
	if(mode!=PP_NULL_MODE){
	  buffer_ptr=str_addchar(buffer_ptr,str_fields.str[j]);
	}
	else{
	  fprintf(stderr,"syntax error: fields entry (%d): unrecognized character '%c'\n",j,str_fields.str[j]);
	  exit(-1);
	}
	break;

	break;
      }
    }
  }
  free(buffer);

  // check whether there are non-Fermionic internal fields
  for(i=0;i<(*fields).internal.length;i++){
    if(is_fermion((*fields).internal.values[i], *fields)==0){
      fprintf(stderr,"warning: some internal fields are not Fermions: means may be computed using sums over permutations rather than the more efficient LU decomposition\n");
      break;
    }
  }

  return(0);
}


// parse virtual_fields list
// write result to fields
int parse_input_virtual_fields(Char_Array str_virtual_fields, Fields_Table* fields, Variables variables){
  // buffer
  char* buffer=calloc(str_virtual_fields.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  Polynomial polynomial;
  int index;
  int i,j;
  int mode;
  int comment=0;
  int res;

  // loop over input
  mode=PP_INDEX_MODE;
  for(j=0;j<str_virtual_fields.length;j++){
    if(comment==1){
      if(str_virtual_fields.str[j]=='\n'){
	comment=0;
      }
    }
    // stay in polynomial mode until ','
    else if(mode==PP_POLYNOMIAL_MODE){
      if(str_virtual_fields.str[j]==','){
	// parse polynomial
	parse_symbolic_expression_str(buffer, *fields, variables, &polynomial);
	// write index and polynomial
	virtual_fields_append_noinit(index, polynomial, &((*fields).virtual_fields));
	mode=PP_INDEX_MODE;
	buffer_ptr=buffer;
	*buffer_ptr='\0';
      }
      // check there are no '=' signs in polynomial (forgotten ',')
      else if(str_virtual_fields.str[j]=='='){
	fprintf(stderr,"syntax error: virtual fields entry (%d): '=' in polynomial '%s'\n",j, buffer);
	exit(-1);
      }
      else if(str_virtual_fields.str[j]=='#'){
	comment=1;
      }
      else{
	buffer_ptr=str_addchar(buffer_ptr,str_virtual_fields.str[j]);
      }
    }
    else{
      switch(str_virtual_fields.str[j]){
      // polynomial mode
      case '=':
	if(mode==PP_INDEX_MODE){
	  // read index
	  res=read_positive_int(buffer, &index);
	  if(res<0){
	    fprintf(stderr,"syntax error: virtual fields entry (%d): virtual field index must be a non-negative integer, got '%s'\n",j,buffer);
	    exit(-1);
	  }
	  // reset buffer
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	  mode=PP_POLYNOMIAL_MODE;
	}
	break;

      // misplaced ','
      case ',':
	fprintf(stderr,"syntax error: virtual fields entry (%d): misplaced ',' or forgotten virtual field index\n",j);
	exit(-1);
	break;

      // ignore spaces and newlines
      case ' ':break;
      case '\n':break;

      // comment
      case '#':
	comment=1;
	break;

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

  // last step
  parse_symbolic_expression_str(buffer, *fields, variables, &polynomial);
  virtual_fields_append_noinit(index, polynomial, &((*fields).virtual_fields));

  // simplify
  for(i=0;i<(*fields).virtual_fields.length;i++){
    polynomial_simplify((*fields).virtual_fields.expr+i, *fields);
  }

  free(buffer);
  return(0);
}


// parse groups of independent fields
int parse_input_groups(Char_Array str_groups, Groups* groups, Polynomial_Matrix propagator, Fields_Table fields){
  // buffer
  char* buffer=calloc(str_groups.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  int index;
  int i,j;
  Int_Array group;
  int mode;
  int comment=0;
  int res;

  // alloc
  init_Groups(groups, GROUP_SIZE);

  // loop over input
  mode=PP_NULL_MODE;
  for(j=0;j<str_groups.length;j++){
    if(comment==1){
      if(str_groups.str[j]=='\n'){
	comment=0;
      }
    }
    else{
      switch(str_groups.str[j]){
      // group mode
      case '(':
	if(mode==PP_NULL_MODE){
	  // reset buffer
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	  // init
	  init_Int_Array(&group, GROUP_SIZE);
	  mode=PP_GROUP_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: groups entry (%d): nested parentheses are not allowed in this entry\n",j);
	  exit(-1);
	}
	break;
      case')':
	if(mode==PP_GROUP_MODE){
	  res=read_positive_int(buffer,&index);
	  if(res<0){
	    fprintf(stderr,"syntax error: groups entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer);
	    exit(-1);
	  }
	  for(i=0;i<groups->length;i++){
	    if(int_array_find(index,groups->indices[i])>=0){
	      fprintf(stderr,"syntax error: groups entry (%d): index %d is already present in group %d\n",j,index,i);
	      exit(-1);
	    }
	  }

	  int_array_append(index, &group);
	  groups_append_noinit(group, groups);
	  mode=PP_NULL_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: groups entry (%d): extra ')'\n",j);
	  exit(-1);
	}
	break;

      // read index
      case',':
	if(mode==PP_GROUP_MODE){
	  res=read_positive_int(buffer,&index);
	  if(res<0){
	    fprintf(stderr,"syntax errorL groups entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer);
	    exit(-1);
	  }
	  for(i=0;i<groups->length;i++){
	    if(int_array_find(index,groups->indices[i])>=0){
	      fprintf(stderr,"syntax error: groups entry (%d): index %d is already present in group %d\n",j,index,i);
	      exit(-1);
	    }
	  }

	  int_array_append(index, &group);
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	}
	break;

      // skip spaces and newlines
      case ' ':break;
      case '\n':break;

      // comment
      case '#':
	comment=1;
	break;

      default:
	if(mode!=PP_NULL_MODE){
	  buffer_ptr=str_addchar(buffer_ptr,str_groups.str[j]);
	}
	else{
	  fprintf(stderr,"syntax error: groups entry (%d): unrecognized character '%c'\n",j,str_groups.str[j]);
	  exit(-1);
	}
	break;
      }
    }
  }

  // check fields in groups
  check_groups(*groups, propagator, fields);

  free(buffer);
  return(0);
}

// check that the members of groups are independent (assuming the virtual fields and propagator were already parsed)
int check_groups(Groups groups, Polynomial_Matrix propagator, Fields_Table fields){
  Int_Array* group_fields=calloc(groups.length,sizeof(Int_Array));
  int i,j,k,l,a,b;

  // get the lists of fields involved in each group
  for(i=0;i<groups.length;i++){
    // expand virtual_fields
    fields_in_virtual_field_list(groups.indices[i], fields, group_fields+i);
  }

  // check whether pairs of fields are present in the propagator
  for(i=0;i<groups.length;i++){
    for(j=0;j<group_fields[i].length;j++){
      if(field_type(group_fields[i].values[j], fields)==FIELD_INTERNAL){
	for(k=0;k<groups.length;k++){
	  if(k!=i){
	    for(l=0;l<group_fields[k].length;l++){
	      if(field_type(group_fields[k].values[l], fields)==FIELD_INTERNAL){
		if(group_fields[i].values[j]==group_fields[k].values[l]){
		  fprintf(stderr,"error: groups %d and %d both contain internal field %d\n",i,k,group_fields[i].values[j]);
		  exit(-1);
		}
		a=intlist_find(propagator.indices, propagator.length, group_fields[i].values[j]);
		b=intlist_find(propagator.indices, propagator.length, group_fields[k].values[l]);
		if(a>=0 && b>=0 && polynomial_is_zero(propagator.matrix[a][b])==0){
		  fprintf(stderr,"error: groups %d and %d are not independent: they contain %d and %d which have a non-vanishing propagator\n",i,k,group_fields[i].values[j], group_fields[k].values[l]);
		  exit(-1);
		}
	      }
	    }
	  }
	}
      }
    }
  }

  for(i=0;i<groups.length;i++){
    free_Int_Array(group_fields[i]);
  }
  free(group_fields);

  return(0);
}

// list of fields involved in a list of virtual_fields
int fields_in_virtual_field_list(Int_Array indices, Fields_Table fields, Int_Array* output){
  int i,j,k;
  Int_Array tmp_array;

  init_Int_Array(output,FIELDS_SIZE);
  for(k=0;k<indices.length;k++){
    if(field_type(indices.values[k],fields)!=FIELD_VIRTUAL){
      int_array_append(abs(indices.values[k]), output);
    }
    else{
      i=intlist_find_err(fields.virtual_fields.indices, fields.virtual_fields.length, indices.values[k]);
      for(j=0;j<fields.virtual_fields.expr[i].length;j++){
	fields_in_virtual_field_list(fields.virtual_fields.expr[i].monomials[j], fields, &tmp_array);
	int_array_concat_unique(tmp_array, output);
	free_Int_Array(tmp_array);
      }
    }
  }

  return(0);
}



// parse variables list
int parse_input_variables(Char_Array str_variables, Variables* variables){
  // buffer
  char* buffer=calloc(str_variables.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  Tree symbol_tree;
  Char_Array varname;
  int j;
  int mode;
  int comment=0;
  char* ptr;

  init_Variables(variables,VARIABLES_SIZE);

  // loop over input
  mode=PP_INDEX_MODE;
  for(j=0;j<str_variables.length;j++){
    if(comment==1){
      if(str_variables.str[j]=='\n'){
	comment=0;
      }
    }
    // stay in polynomial mode until ','
    else if(mode==PP_POLYNOMIAL_MODE){
      if(str_variables.str[j]==','){
	// parse symbol tree
	str_to_symbol_tree(buffer, &symbol_tree);
	variables_append_noinit(varname , symbol_tree, variables);
	mode=PP_INDEX_MODE;
	buffer_ptr=buffer;
	*buffer_ptr='\0';
      }
      // check there are no '=' signs in polynomial (forgotten ',')
      else if(str_variables.str[j]=='='){
	fprintf(stderr,"syntax error: preprocessor_variables entry (%d): '=' in polynomial '%s'\n",j, buffer);
	exit(-1);
      }
      else if(str_variables.str[j]=='#'){
	comment=1;
      }
      else{
	buffer_ptr=str_addchar(buffer_ptr,str_variables.str[j]);
      }
    }
    else{
      switch(str_variables.str[j]){
      // polynomial mode
      case '=':
	if(mode==PP_INDEX_MODE){
	  if(*buffer=='\0'){
	    fprintf(stderr,"syntax error: preprocessor_variables entry (%d): empty variable name\n",j);
	    exit(-1);
	  }
	  for(ptr=buffer;*ptr!='\0';ptr++){
	    if(*ptr=='<' || *ptr=='>'){
	      fprintf(stderr,"syntax error: preprocessor_variables entry (%d): forbidden '%c' character in variable name '%s'\n",j,*ptr,buffer);
	      exit(-1);
	    }
	  }
	  str_to_char_array(buffer,&varname);
	  // check that the variable name is not 'OUT' or 'FLOW' or 'RCC'
	  if(char_array_cmp_str(varname, "OUT")==1){
	    fprintf(stderr,"error: preprocessor_variables entry (%d): variable name cannot be 'OUT' (this name is reserved)\n",j);
	    exit(-1);
	  }
	  if(char_array_cmp_str(varname, "FLOW")==1){
	    fprintf(stderr,"error: preprocessor_variables entry (%d): variable name cannot be 'FLOW' (this name is reserved)\n",j);
	    exit(-1);
	  }
	  if(char_array_cmp_str(varname, "RCC")==1){
	    fprintf(stderr,"error: preprocessor_variables entry (%d): variable name cannot be 'RCC' (this name is reserved)\n",j);
	    exit(-1);
	  }
	  // reset buffer
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	  mode=PP_POLYNOMIAL_MODE;
	}
	break;

      // misplaced ','
      case ',':
	fprintf(stderr,"syntax error: preprocessor_variables entry (%d): misplaced ',' or forgotten variable name\n",j);
	exit(-1);
	break;

      // ignore characters
      case ' ':break;
      case '\n': break;

      // comment
      case '#':
	comment=1;
	break;

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

  // last step
  if(mode==PP_POLYNOMIAL_MODE){
    str_to_symbol_tree(buffer, &symbol_tree);
    variables_append_noinit(varname , symbol_tree, variables);
  }
  else if(*buffer!='\0'){
    fprintf(stderr,"syntax error: preprocessor_variables entry (%d): mismatched variable/polynomial pair, got '%s'\n",j,buffer);
    exit(-1);
  }

  free(buffer);
  return(0);
}


// parse identities between fields
// write result to fields
int parse_input_identities(Char_Array str_identities, Fields_Table* fields, Variables variables){
  // buffer
  char* buffer=calloc(str_identities.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  Int_Array monomial;
  Polynomial polynomial;
  int i,j;
  int sign;
  int tmp;
  int mode;
  int comment=0;
  int ret;

  init_Int_Array(&monomial, MONOMIAL_SIZE);

  // loop over input
  mode=PP_NULL_MODE;
  for(j=0;j<str_identities.length;j++){
    if(comment==1){
      if(str_identities.str[j]=='\n'){
	comment=0;
      }
    }
    // stay in polynomial mode until ','
    else if(mode==PP_POLYNOMIAL_MODE){
      if(str_identities.str[j]==','){
	// parse polynomial
	parse_symbolic_expression_str(buffer, *fields, variables, &polynomial);
	// write monomial and polynomial
	identities_append_noinit(monomial, polynomial, &((*fields).ids));
	// realloc
	init_Int_Array(&monomial, MONOMIAL_SIZE);
	mode=PP_NULL_MODE;
      }
      // check there are no '=' signs in polynomial (forgotten ',')
      else if(str_identities.str[j]=='='){
	fprintf(stderr,"syntax error: identities entry (%d): '=' in polynomial '%s'\n",j, buffer);
	exit(-1);
      }
      else if(str_identities.str[j]=='#'){
	comment=1;
      }
      else{
	buffer_ptr=str_addchar(buffer_ptr,str_identities.str[j]);
      }
    }
    else{
      switch(str_identities.str[j]){
      // monomial
      case '[':
	if(mode==PP_NULL_MODE){
	  mode=PP_INDEX_MODE;
	}
	if(mode==PP_INDEX_MODE){
	  mode=PP_BRACKET_MODE;
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	}
	else{
	  fprintf(stderr,"syntax error: identities entry (%d): misplaced '['\n",j);
	  exit(-1);
	}
	break;
      // for notational homogeneity
      case 'f':
	if(mode==PP_BRACKET_MODE){
	  mode=PP_FIELD_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: identities entry (%d): misplaced 'f'\n",j);
	  exit(-1);
	}
	break;
      // write monomial term
      case ']':
	if(mode==PP_FIELD_MODE){
	  ret=read_int(buffer,&i);
	  if(ret<0){
	    fprintf(stderr,"syntax error: identities entry (%d): field indices must be integers, got '%s'\n",j,buffer);
	    exit(-1);
	  }
	  int_array_append(i,&monomial);
	  mode=PP_INDEX_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: identities entry (%d): mismatched ']'\n",j);
	  exit(-1);
	}
	break;

      // polynomial mode
      case '=':
	if(mode==PP_INDEX_MODE){
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	  mode=PP_POLYNOMIAL_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: identities entry (%d): misplaced '='\n",j);
	  exit(-1);
	}
	break;

      // misplaced ','
      case ',':
	fprintf(stderr,"syntax error: identities entry (%d): misplaced ',' or forgotten left hand side\n",j);
	exit(-1);
	break;

      // ignore spaces and newlines
      case ' ':break;
      case '\n':break;

      // comment
      case '#':
	comment=1;
	break;

      default:
	if(mode!=PP_NULL_MODE){
	  buffer_ptr=str_addchar(buffer_ptr,str_identities.str[j]);
	}
	else{
	  fprintf(stderr,"syntax error: identities entry (%d): unrecognized character '%c'\n",j,str_identities.str[j]);
	  exit(-1);
	}
	break;
      }
    }
  }

  // last step
  if(mode==PP_POLYNOMIAL_MODE){
    parse_symbolic_expression_str(buffer, *fields, variables, &polynomial);
    identities_append_noinit(monomial, polynomial, &((*fields).ids));
  }
  else if(*buffer!='\0'){
    fprintf(stderr,"syntax error: identities entry (%d): mismatched left/right hand side, got '%s'\n",j,buffer);
    exit(-1);
  }
  else{
    free_Int_Array(monomial);
  }

  // sort
  // don't use the identities to simplify
  tmp=(*fields).ids.length;
  (*fields).ids.length=0;
  for(i=0;i<tmp;i++){
    sign=1;
    monomial_sort((*fields).ids.lhs[i], *fields, &sign);
    polynomial_simplify((*fields).ids.rhs+i, *fields);
    polynomial_multiply_Qscalar((*fields).ids.rhs[i],quot(sign,1));
  }
  (*fields).ids.length=tmp;

  free(buffer);
  return(0);
}

// parse propagator
int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propagator, Fields_Table fields){
  // buffer
  char* buffer=calloc(str_propagator.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  int i,j;
  int index1=-1;
  int index2=-1;
  int mode;
  int comment=0;
  int ret;

  // allocate memory
  init_Polynomial_Matrix(propagator, fields.internal.length);

  // copy indices
  for(i=0;i<fields.internal.length;i++){
    (*propagator).indices[i]=fields.internal.values[i];
  }

  // loop over input
  mode=PP_INDEX_MODE;
  for(j=0;j<str_propagator.length;j++){
    if(comment==1){
      if(str_propagator.str[j]=='\n'){
	comment=0;
      }
    }
    else{
      switch(str_propagator.str[j]){
      // indices
      case ';':
	if(mode==PP_INDEX_MODE){
	  ret=read_positive_int(buffer,&i);
	  if(ret<0){
	    fprintf(stderr,"syntax error: propagator entry (%d): indices must be non-negative integers, got '%s'\n",j,buffer);
	    exit(-1);
	  }
	  index1=intlist_find_err((*propagator).indices, (*propagator).length, i);
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	}
	else{
	  fprintf(stderr,"syntax error: propagator entry (%d): misplaced ';'\n",j);
	  exit(-1);
	}
	break;
      case ':':
	if(mode==PP_INDEX_MODE){
	  ret=read_positive_int(buffer,&i);
	  if(ret<0){
	    fprintf(stderr,"syntax error: propagator entry (%d): indices must be non-negative integers, got '%s'\n",j,buffer);
	    exit(-1);
	  }
	  index2=intlist_find_err((*propagator).indices, (*propagator).length, i);
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	  mode=PP_POLYNOMIAL_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: propagator entry (%d): misplaced ':'\n",j);
	  exit(-1);
	}
	break;

      // num
      case ',':
	if(mode==PP_POLYNOMIAL_MODE && index1>=0 && index2>=0){
	  free_Polynomial((*propagator).matrix[index1][index2]);
	  str_to_Polynomial(buffer,(*propagator).matrix[index1]+index2);
	  index1=-1;
	  index2=-1;
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	  mode=PP_INDEX_MODE;
	}
	else{
	  fprintf(stderr,"syntax error: propagator entry (%d): mismatched index/polynomial pair\n",j);
	  exit(-1);
	}
	break;

      // ignore ' 'and newline
      case ' ':break;
      case '\n':break;

      // comment
      case '#':
	comment=1;
	break;

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

  // last step
  if(mode==PP_POLYNOMIAL_MODE){
    free_Polynomial((*propagator).matrix[index1][index2]);
    str_to_Polynomial(buffer,(*propagator).matrix[index1]+index2);
  }
  else if (*buffer!='\0'){
    fprintf(stderr,"syntax error: propagator entry (%d): mismatched index/polynomial pair, got '%s'\n",j,buffer);
    exit(-1);
  }

  // check whether the polynomial is numeric
  if(polynomial_matrix_is_numeric(*propagator)==0){
    fprintf(stderr,"warning: the propagator contains polynomial entries: means may be computed using sums over permutations rather than the more efficient LU decomposition\n");
  }

  free(buffer);
  return(0);
}


// parse input polynomial
int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields, Variables variables){
  parse_symbolic_expression(str_polynomial, fields, variables, output);
  polynomial_simplify(output, fields);
  return(0);
}



// parse id table
// fields argument for sorting
int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields, Variables variables){
  // buffer
  char* buffer=calloc(str_idtable.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  int index;
  Polynomial polynomial;
  int j;
  int mode;
  int comment=0;
  int ret;

  // allocate memory
  init_Id_Table(idtable,EQUATION_SIZE);

  // loop over input
  mode=PP_INDEX_MODE;
  for(j=0;j<str_idtable.length;j++){
    if(comment==1){
      if(str_idtable.str[j]=='\n'){
	comment=0;
      }
    }
    else{
      switch(str_idtable.str[j]){
      // end polynomial mode
      case ',':
	// write polynomial
	if(mode==PP_POLYNOMIAL_MODE){
	  parse_symbolic_expression_str(buffer, fields, variables, &polynomial);
	  // add to idtable
	  idtable_append_noinit(index,polynomial,idtable);
	  mode=PP_INDEX_MODE;
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	  }
	else{
	  fprintf(stderr,"syntax error: id_table entry (%d): misplaced ','\n",j);
	  exit(-1);
	}
	break;

      case ':':
	if(mode==PP_INDEX_MODE){
	  ret=read_positive_int(buffer,&index);
	  if(ret<0){
	    fprintf(stderr,"syntax error: id_table entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer);
	    exit(-1);
	  }
	  mode=PP_POLYNOMIAL_MODE;
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	}
	else{
	  fprintf(stderr,"syntax error: id_table entry (%d): misplaced ':'\n",j);
	  exit(-1);
	}
	break;

      // ignore ' ' and '\n'
      case ' ':break;
      case '\n':break;

      // comment
      case '#':
	comment=1;
	break;

      default:
	if(mode!=PP_NULL_MODE){
	  buffer_ptr=str_addchar(buffer_ptr,str_idtable.str[j]);
	}
	else{
	  fprintf(stderr,"syntax error: id_table entry (%d): unrecognized character '%c'\n",j, str_idtable.str[j]);
	  exit(-1);
	}
	break;
      }
    }
  }

  //last step
  if(mode==PP_POLYNOMIAL_MODE){
    parse_symbolic_expression_str(buffer, fields, variables, &polynomial);
    idtable_append_noinit(index,polynomial,idtable);
  }
  else if(*buffer!='\0'){
    fprintf(stderr,"syntax error: id_table entry (%d): mismatched index/polynomial pair, got '%s'\n",j,buffer);
    exit(-1);
  }

  // sort
  for(j=0;j<(*idtable).length;j++){
    polynomial_simplify((*idtable).polynomials+j, fields);
  }

  free(buffer);
  return(0);
}

// parse a list of labels
int parse_labels(Char_Array str_labels, Labels* labels){
  // buffer
  char* buffer=calloc(str_labels.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  Char_Array label;
  int index;
  int j;
  int mode;
  int comment=0;
  int ret;

  // allocate memory
  init_Labels(labels,EQUATION_SIZE);

  // loop over input
  mode=PP_INDEX_MODE;
  for(j=0;j<str_labels.length;j++){
    if(comment==1){
      if(str_labels.str[j]=='\n'){
	comment=0;
      }
    }
    else{
      switch(str_labels.str[j]){
      case '"':
	if(mode==PP_INDEX_MODE){
	  mode=PP_LABEL_MODE;
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	}
	// write
	else if(mode==PP_LABEL_MODE){
	  str_to_char_array(buffer,&label);
	  labels_append_noinit(label,index,labels);
	  mode=PP_INDEX_MODE;
	  buffer_ptr=buffer;
	  *buffer_ptr='\0';
	}
	else{
	  fprintf(stderr,"syntax error: labels entry (%d): mismatched '\"'\n",j);
	  exit(-1);
	}
	break;

      case ':':
	// write
	if(mode==PP_INDEX_MODE){
	  sscanf(buffer,"%d",&index);
	  ret=read_positive_int(buffer,&index);
	  if(ret<0){
	    fprintf(stderr,"syntax error: labels entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer);
	    exit(-1);
	  }
	}
	else{
	  fprintf(stderr,"syntax error: labels entry (%d): misplaced ':'\n",j);
	  exit(-1);
	}
	break;

      // characters to ignore
      case ' ':break;
      case '\n':break;
      case ',':break;

      // comment
      case '#':
	comment=1;
	break;

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

  free(buffer);
  return(0);
}


// read initial condition for numerical computation (using either RCC or RCC_mpfr, as specified by mpfr_flag)
int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_flag){
  char* buffer=calloc(init_cd.length+1,sizeof(char));
  char* buffer_ptr=buffer;
  int index=0;
  int i,j;
  int comment_mode=0;
  int dcount=0;
  int ret;

  *buffer_ptr='\0';
  // loop over the input
  for(j=0;j<init_cd.length;j++){
    if(comment_mode==1){
      if(init_cd.str[j]=='\n'){
	comment_mode=0;
      }
    }
    else{
      switch(init_cd.str[j]){
      // new term
      case ',':
	// write init
	if(mpfr_flag==0){
	  ret=read_long_double(buffer,(*init).values+intlist_find_err((*init).indices,(*init).length,index));
	  if(ret<0){
	    fprintf(stderr,"syntax error: initial_condition entry (%d): failed to extract long double from '%s'\n",j,buffer);
	    exit(-1);
	  }
	}
	else{
	  mpfr_strtofr((*init_mpfr).values[intlist_find_err((*init_mpfr).indices,(*init_mpfr).length,index)], buffer, &buffer_ptr, 10, MPFR_RNDN);
	}
	// reset buffer
	buffer_ptr=buffer;
	*buffer_ptr='\0';
	// reset derivatives counter
	dcount=0;
	break;
      // separator
      case ':':
	// write index
	ret=read_int(buffer,&i);
	if(ret<0){
	  fprintf(stderr,"syntax error: initial_condition entry (%d): index must be an integer, got '%s'\n",j,buffer);
	  exit(-1);
	}
	if(i<0){
	  index=i-dcount*DOFFSET;
	}
	else{
	  index=i+dcount*DOFFSET;
	}
	buffer_ptr=buffer;
	*buffer_ptr='\0';
	break;
      // derivatives
      case 'd':
	dcount++;
	break;

      // characters to ignore
      case ' ':break;
      case '\n':break;

      // comments
      case '#':
	comment_mode=1;
	break;

      default:
	// write to buffer
	buffer_ptr=str_addchar(buffer_ptr,init_cd.str[j]);
	break;
      }
    }
  }

  // write init
  if(mpfr_flag==0){
    ret=read_long_double(buffer,(*init).values+intlist_find_err((*init).indices,(*init).length,index));
    if(ret<0){
      fprintf(stderr,"syntax error: initial_condition entry (%d): failed to extract long double from '%s'\n",j,buffer);
      exit(-1);
    }
  }
  else{
    mpfr_strtofr((*init_mpfr).values[intlist_find_err((*init_mpfr).indices,(*init_mpfr).length,index)], buffer, &buffer_ptr, 10, MPFR_RNDN);
  }

  free(buffer);
  return(0);
}


// set indices and length of init
int prepare_init(int* indices, int length, RCC* init){
  int i;
  init_RCC(init, length);
  for(i=0;i<length;i++){
    (*init).indices[i]=indices[i];
    // set constants to 1
    if(indices[i]<0 && indices[i]>-DOFFSET){
      (*init).values[i]=1.;
    }
    else{
      (*init).values[i]=0.;
    }

  }
  return(0);
}
// set indices and length of init for RCC_mpfr
int prepare_init_mpfr(int* indices, int length, RCC_mpfr* init){
  int i;
  init_RCC_mpfr(init, length);
  for(i=0;i<length;i++){
    (*init).indices[i]=indices[i];
    // set constants to 1
    if(indices[i]<0 && indices[i]>-DOFFSET){
      mpfr_set_ui((*init).values[i],1,MPFR_RNDN);
    }
    else{
      mpfr_set_zero((*init).values[i],1);
    }

  }
  return(0);
}