path: root/src/fields.c

/*
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 "fields.h"
#include "definitions.cpp"
#include <stdio.h>
#include <stdlib.h>
#include "number.h"
#include "tools.h"
#include "polynomial.h"
#include "array.h"
#include "rational.h"
#include "tree.h"

//---------------------------------------------------------------------
//
//                              Fields_Table
//
//---------------------------------------------------------------------

// init and free for Fields_Table
int init_Fields_Table(Fields_Table* fields){
  init_Int_Array(&((*fields).parameter),FIELDS_SIZE);
  init_Int_Array(&((*fields).external),FIELDS_SIZE);
  init_Int_Array(&((*fields).internal),FIELDS_SIZE);
  init_Identities(&((*fields).ids), FIELDS_SIZE);
  init_Virtual_fields(&((*fields).virtual_fields), FIELDS_SIZE);
  init_Int_Array(&((*fields).fermions),FIELDS_SIZE);
  init_Int_Array(&((*fields).noncommuting),FIELDS_SIZE);
  return(0);
}
int free_Fields_Table(Fields_Table fields){
  free_Int_Array(fields.parameter);
  free_Int_Array(fields.external);
  free_Int_Array(fields.internal);
  free_Identities(fields.ids);
  free_Virtual_fields(fields.virtual_fields);
  free_Int_Array(fields.fermions);
  free_Int_Array(fields.noncommuting);
  return(0);
}

// determine field type
int field_type(int index, Fields_Table fields){
  if(int_array_find(abs(index), fields.parameter)>=0){
    return(FIELD_PARAMETER);
  }
  else if(int_array_find(abs(index), fields.external)>=0){
    return(FIELD_EXTERNAL);
  }
  else if(int_array_find(abs(index), fields.internal)>=0){
    return(FIELD_INTERNAL);
  }
  else if(intlist_find(fields.virtual_fields.indices, fields.virtual_fields.length, index)>=0){
    return(FIELD_VIRTUAL);
  }

  fprintf(stderr,"error: index %d is neither a parameter nor an external or an internal field, nor a virtual field\n",index);
  exit(-1);
}

// check whether a field anticommutes
int is_fermion(int index, Fields_Table fields){
  if(int_array_find(abs(index), fields.fermions)>=0){
    return(1);
  }
  else{
    return(0);
  }
}

// check whether a field is non-commuting
int is_noncommuting(int index, Fields_Table fields){
  if(int_array_find(abs(index), fields.noncommuting)>=0){
    return(1);
  }
  else{
    return(0);
  }
}


//---------------------------------------------------------------------
//
//                               Identities
//
//---------------------------------------------------------------------

// allocate memory
int init_Identities(Identities* identities,int size){
  (*identities).lhs=calloc(size,sizeof(Int_Array));
  (*identities).rhs=calloc(size,sizeof(Polynomial));
  (*identities).length=0;
  (*identities).memory=size;
  return(0);
}

// free memory
int free_Identities(Identities identities){
  int i;
  for(i=0;i<identities.length;i++){
    free_Int_Array(identities.lhs[i]);
    free_Polynomial(identities.rhs[i]);
  }
  free(identities.lhs);
  free(identities.rhs);
  return(0);
}

// resize
int resize_identities(Identities* identities,int new_size){
  Identities new_identities;
  int i;

  init_Identities(&new_identities,new_size);
  for(i=0;i<(*identities).length;i++){
    new_identities.lhs[i]=(*identities).lhs[i];
    new_identities.rhs[i]=(*identities).rhs[i];
  }
  new_identities.length=(*identities).length;

  free((*identities).lhs);
  free((*identities).rhs);

  *identities=new_identities;
  return(0);
}

// copy
int identities_cpy(Identities input, Identities* output){
  init_Identities(output,input.length);
  identities_cpy_noinit(input,output);
  return(0);
}
int identities_cpy_noinit(Identities input, Identities* output){
  int i;
  if((*output).memory<input.length){
    fprintf(stderr,"error: trying to copy an identities collection of length %d to another with memory %d\n",input.length,(*output).memory);
    exit(-1);
  }
  for(i=0;i<input.length;i++){
    int_array_cpy(input.lhs[i],(*output).lhs+i);
    polynomial_cpy(input.rhs[i],(*output).rhs+i);
  }
  (*output).length=input.length;
  
  return(0);
}

// append an element to a identities
int identities_append(Int_Array lhs, Polynomial rhs, Identities* output){
  int offset=(*output).length;

  if((*output).length>=(*output).memory){
    resize_identities(output,2*(*output).memory+1);
  }

  // copy and allocate
  int_array_cpy(lhs,(*output).lhs+offset);
  polynomial_cpy(rhs,(*output).rhs+offset);
  // increment length
  (*output).length++;
  return(0);
}
// append an element to a identities without allocating memory
int identities_append_noinit(Int_Array lhs, Polynomial rhs, Identities* output){
  int offset=(*output).length;

  if((*output).length>=(*output).memory){
    resize_identities(output,2*(*output).memory+1);
  }

  // copy without allocating
  (*output).lhs[offset]=lhs;
  (*output).rhs[offset]=rhs;
  // increment length
  (*output).length++;
  return(0);
}

// concatenate two identitiess
int identities_concat(Identities input, Identities* output){
  int i;
  for(i=0;i<input.length;i++){
    identities_append(input.lhs[i],input.rhs[i],output);
  }
  return(0);
}


// resolve the identities
// requires both the monomials in polynomial and the ids in fields to be sorted
// IMPORTANT: the sorting must be such that noncommuting fields must come before the other fields
int resolve_ids(Polynomial* polynomial, Fields_Table fields){
  int i,j,k,l;
  int sign;
  int fermion_count;
  int first_field;
  int at_least_one;
  int security;
  Int_Array pre_monomial;
  Int_Array post_monomial;
  Number num;
  Number tmp_num;

  // loop over monomials
  for(i=0;i<(*polynomial).length;i++){
    at_least_one=1;
    security=0;
    // repeat the simplification until the monomial is fully simplified
    while(at_least_one>0){
      at_least_one=0;

      // prevent infinite loops
      security++;
      if(security>1000000){
	fprintf(stderr,"error: 1000000 iterations reached when trying to simplify a monomial\n");
	exit(-1);
      }

      // loop over ids
      for(j=0;j<fields.ids.length;j++){

	// check whether the monomial matches the id
	first_field=int_array_is_subarray_noncommuting(fields.ids.lhs[j],(*polynomial).monomials[i],fields);
	if(first_field>=0){
	  init_Int_Array(&pre_monomial, (*polynomial).monomials[i].length);
	  init_Int_Array(&post_monomial, (*polynomial).monomials[i].length);

	  // add whatever is before the first field to pre
	  for(k=0;k<first_field;k++){
	    int_array_append((*polynomial).monomials[i].values[k],&pre_monomial);
	  }

	  // find others and move them together
	  // sign from moving the fields
	  sign=1;
	  // number of Fermions to jump over
	  fermion_count=0;
	  for(l=1,k=first_field+1;k<(*polynomial).monomials[i].length;k++){
	    // check whether the field is identical to the "current" one in the id
	    // if l is too large, then keep the field
	    if(l>=fields.ids.lhs[j].length || (*polynomial).monomials[i].values[k]!=fields.ids.lhs[j].values[l]){
	      // add to post
	      int_array_append((*polynomial).monomials[i].values[k],&post_monomial);
	      // count Fermions to jump
	      if(is_fermion((*polynomial).monomials[i].values[k],fields)){
		fermion_count++;
	      }
	    }
	    else{
	      // sign correction
	      if(is_fermion(fields.ids.lhs[j].values[l],fields) && fermion_count % 2 == 1){
		sign*=-1;
	      }
	      // increment "current" field in the id
	      l++;
	    }
	  }

	  num=number_Qprod_ret(quot(sign,1),(*polynomial).nums[i]);
	  // add extra monomials (if there are more than 1)
	  for(k=1;k<fields.ids.rhs[j].length;k++){
	    number_prod(num, fields.ids.rhs[j].nums[k], &tmp_num);
	    polynomial_append(pre_monomial, (*polynomial).factors[i], tmp_num, polynomial);
	    free_Number(tmp_num);
	    int_array_concat(fields.ids.rhs[j].monomials[k],(*polynomial).monomials+(*polynomial).length-1);
	    int_array_concat(post_monomial,(*polynomial).monomials+(*polynomial).length-1);
	    // re-sort monomial
	    sign=1;
	    monomial_sort((*polynomial).monomials[(*polynomial).length-1],fields,&sign);
	    number_Qprod_chain(quot(sign,1),(*polynomial).nums+(*polynomial).length-1);
	  } 
	  // correct i-th monomial
	  free_Number((*polynomial).nums[i]);
	  (*polynomial).nums[i]=number_prod_ret(num,fields.ids.rhs[j].nums[0]);
	  free_Int_Array((*polynomial).monomials[i]);
	  (*polynomial).monomials[i]=pre_monomial;
	  int_array_concat(fields.ids.rhs[j].monomials[0],(*polynomial).monomials+i);
	  int_array_concat(post_monomial,(*polynomial).monomials+i);
	  free_Int_Array(post_monomial);
	  // re-sort monomial
	  sign=1;
	  monomial_sort((*polynomial).monomials[i],fields,&sign);
	  number_Qprod_chain(quot(sign,1),(*polynomial).nums+i);

	  // free num
	  free_Number(num);

	  // repeat the replacement (in order to perform all of the replacements if several are necessary)
	  j=0;
	  if(at_least_one==0){
	    at_least_one=1;
	  }
	}
      }
    }
  }

  return(0);
}

// check whether an array is a sub-array of another
// requires noncommuting elements to be next to each other
// other elements may be separated, but the order must be respected
// returns the first index of the sub-array
// IMPORTANT: the noncommuting elements must precede all others in input and in test_array
int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fields_Table fields){
  int i,j;
  int matches=0;
  int post_nc=0;
  int match_nc;
  int first=-1;

  // cannot fit
  if(test_array.length<input.length){
    return(-1);
  }

  // bound noncommuting elements
  while(is_noncommuting(input.values[post_nc], fields)==1){
    post_nc++;
  }

  for(i=0,match_nc=0;i<test_array.length;i++){
    if(test_array.values[i]==input.values[0]){
      match_nc=1;
    }
    for(j=1;j<post_nc;j++){
      if(i+j>=test_array.length || test_array.values[i+j]!=input.values[j]){
	match_nc=0;
      }
    }
    if(match_nc==1){
      first=i;
      break;
    }
  }

  if(first<0){
    return(-1);
  }

  if(post_nc>0){
    matches=post_nc;
  }
  else{
    matches=1;
  }

  for(i=first+1;i<test_array.length && matches<input.length;i++){
    if(input.values[matches]==test_array.values[i]){
      matches++;
    }
  }
  if(matches==input.length){
    return(first);
  }
  else{
    return(-1);
  }
}


//---------------------------------------------------------------------
//
//                                Virtual_fields
//
//---------------------------------------------------------------------

// allocate memory
int init_Virtual_fields(Virtual_fields* virtual_fields,int size){
  (*virtual_fields).indices=calloc(size,sizeof(int));
  (*virtual_fields).expr=calloc(size,sizeof(Polynomial));
  (*virtual_fields).length=0;
  (*virtual_fields).memory=size;
  return(0);
}

// free memory
int free_Virtual_fields(Virtual_fields virtual_fields){
  int i;
  for(i=0;i<virtual_fields.length;i++){
    free_Polynomial(virtual_fields.expr[i]);
  }
  free(virtual_fields.indices);
  free(virtual_fields.expr);
  return(0);
}

// resize
int resize_virtual_fields(Virtual_fields* virtual_fields,int new_size){
  Virtual_fields new_virtual_fields;
  int i;

  init_Virtual_fields(&new_virtual_fields,new_size);
  for(i=0;i<(*virtual_fields).length;i++){
    new_virtual_fields.indices[i]=(*virtual_fields).indices[i];
    new_virtual_fields.expr[i]=(*virtual_fields).expr[i];
  }
  new_virtual_fields.length=(*virtual_fields).length;

  free((*virtual_fields).indices);
  free((*virtual_fields).expr);

  *virtual_fields=new_virtual_fields;
  return(0);
}

// copy
int virtual_fields_cpy(Virtual_fields input, Virtual_fields* output){
  init_Virtual_fields(output,input.length);
  virtual_fields_cpy_noinit(input,output);
  return(0);
}
int virtual_fields_cpy_noinit(Virtual_fields input, Virtual_fields* output){
  int i;
  if((*output).memory<input.length){
    fprintf(stderr,"error: trying to copy a virtual fields collection of length %d to another with memory %d\n",input.length,(*output).memory);
    exit(-1);
  }
  for(i=0;i<input.length;i++){
    (*output).indices[i]=input.indices[i];
    polynomial_cpy(input.expr[i],(*output).expr+i);
  }
  (*output).length=input.length;
  
  return(0);
}

// append an element to a virtual_fields
int virtual_fields_append(int index, Polynomial expr, Virtual_fields* output){
  int offset=(*output).length;

  if((*output).length>=(*output).memory){
    resize_virtual_fields(output,2*(*output).memory+1);
  }

  // copy and allocate
  (*output).indices[offset]=index;
  polynomial_cpy(expr,(*output).expr+offset);
  // increment length
  (*output).length++;
  return(0);
}
// append an element to a virtual_fields without allocating memory
int virtual_fields_append_noinit(int index, Polynomial expr, Virtual_fields* output){
  int offset=(*output).length;

  if((*output).length>=(*output).memory){
    resize_virtual_fields(output,2*(*output).memory+1);
  }

  // copy without allocating
  (*output).indices[offset]=index;
  (*output).expr[offset]=expr;
  // increment length
  (*output).length++;
  return(0);
}

// concatenate two virtual_fields
int virtual_fields_concat(Virtual_fields input, Virtual_fields* output){
  int i;
  for(i=0;i<input.length;i++){
    virtual_fields_append(input.indices[i],input.expr[i],output);
  }
  return(0);
}



//---------------------------------------------------------------------
//
//                                 Groups
//
//---------------------------------------------------------------------

// allocate memory
int init_Groups(Groups* groups,int size){
  (*groups).indices=calloc(size,sizeof(Int_Array));
  (*groups).length=0;
  (*groups).memory=size;
  return(0);
}

// free memory
int free_Groups(Groups groups){
  int i;
  for(i=0;i<groups.length;i++){
    free_Int_Array(groups.indices[i]);
  }
  free(groups.indices);
  return(0);
}

// resize
int resize_groups(Groups* groups,int new_size){
  Groups new_groups;
  int i;

  init_Groups(&new_groups,new_size);
  for(i=0;i<(*groups).length;i++){
    new_groups.indices[i]=(*groups).indices[i];
  }
  new_groups.length=(*groups).length;

  free((*groups).indices);

  *groups=new_groups;
  return(0);
}

// copy
int groups_cpy(Groups input, Groups* output){
  init_Groups(output,input.length);
  groups_cpy_noinit(input,output);
  return(0);
}
int groups_cpy_noinit(Groups input, Groups* output){
  int i;
  if((*output).memory<input.length){
    fprintf(stderr,"error: trying to copy a groups collection of length %d to another with memory %d\n",input.length,(*output).memory);
    exit(-1);
  }
  for(i=0;i<input.length;i++){
    int_array_cpy(input.indices[i],(*output).indices+i);
  }
  (*output).length=input.length;
  
  return(0);
}

// append an element to a groups
int groups_append(Int_Array indices, Groups* output){
  int offset=(*output).length;

  if((*output).length>=(*output).memory){
    resize_groups(output,2*(*output).memory+1);
  }

  // copy and allocate
  int_array_cpy(indices,(*output).indices+offset);
  // increment length
  (*output).length++;
  return(0);
}
// append an element to a groups without allocating memory
int groups_append_noinit(Int_Array indices, Groups* output){
  int offset=(*output).length;

  if((*output).length>=(*output).memory){
    resize_groups(output,2*(*output).memory+1);
  }

  // copy without allocating
  (*output).indices[offset]=indices;
  // increment length
  (*output).length++;
  return(0);
}

// concatenate two groupss
int groups_concat(Groups input, Groups* output){
  int i;
  for(i=0;i<input.length;i++){
    groups_append(input.indices[i],output);
  }
  return(0);
}

// find which group an index belongs to
int find_group(int index, Groups groups){
  int i,j;
  for(i=0;i<groups.length;i++){
    for(j=0;j<groups.indices[i].length;j++){
      if(groups.indices[i].values[j]==index){
	return(i);
      }
    }
  }
  return(-1);
}


//---------------------------------------------------------------------
//
//                               Variables
//
//---------------------------------------------------------------------

// allocate memory
int init_Variables(Variables* variables,int size){
  (*variables).var_names=calloc(size,sizeof(Char_Array));
  (*variables).symbol_trees=calloc(size,sizeof(Tree));
  (*variables).length=0;
  (*variables).memory=size;
  return(0);
}

// free memory
int free_Variables(Variables variables){
  int i;
  for(i=0;i<variables.length;i++){
    free_Char_Array(variables.var_names[i]);
    free_Tree(variables.symbol_trees[i]);
  }
  free(variables.var_names);
  free(variables.symbol_trees);
  return(0);
}

// resize
int resize_variables(Variables* variables,int new_size){
  Variables new_variables;
  int i;

  init_Variables(&new_variables,new_size);
  for(i=0;i<(*variables).length;i++){
    new_variables.var_names[i]=(*variables).var_names[i];
    new_variables.symbol_trees[i]=(*variables).symbol_trees[i];
  }
  new_variables.length=(*variables).length;

  free((*variables).var_names);
  free((*variables).symbol_trees);

  *variables=new_variables;
  return(0);
}

// copy
int variables_cpy(Variables input, Variables* output){
  init_Variables(output,input.length);
  variables_cpy_noinit(input,output);
  return(0);
}
int variables_cpy_noinit(Variables input, Variables* output){
  int i;
  if((*output).memory<input.length){
    fprintf(stderr,"error: trying to copy a variables collection of length %d to another with memory %d\n",input.length,(*output).memory);
    exit(-1);
  }
  for(i=0;i<input.length;i++){
    char_array_cpy(input.var_names[i], (*output).var_names+i);
    tree_cpy(input.symbol_trees[i],(*output).symbol_trees+i);
  }
  (*output).length=input.length;
  
  return(0);
}

// append an element to a variables collection
int variables_append(Char_Array var_name, Tree symbol_tree, Variables* output){
  int offset=(*output).length;

  if((*output).length>=(*output).memory){
    resize_variables(output,2*(*output).memory+1);
  }

  // copy and allocate
  char_array_cpy(var_name,(*output).var_names+offset);
  tree_cpy(symbol_tree,(*output).symbol_trees+offset);
  // increment length
  (*output).length++;
  return(0);
}
// append an element to a variables collection without allocating memory
int variables_append_noinit(Char_Array var_name, Tree symbol_tree, Variables* output){
  int offset=(*output).length;

  if((*output).length>=(*output).memory){
    resize_variables(output,2*(*output).memory+1);
  }

  // copy without allocating
  (*output).var_names[offset]=var_name;
  (*output).symbol_trees[offset]=symbol_tree;
  // increment length
  (*output).length++;
  return(0);
}

// concatenate two variables collections
int variables_concat(Variables input, Variables* output){
  int i;
  for(i=0;i<input.length;i++){
    variables_append(input.var_names[i], input.symbol_trees[i], output);
  }
  return(0);
}


// find a variable matching a var_name
int variables_find_var(Char_Array name, Variables variables, Tree* output){
  Char_Array varname;
  int i;

  // drop inital '$'
  char_array_substring(name, 1, name.length-1, &varname);
  for(i=0;i<variables.length;i++){
    if(char_array_cmp(varname, variables.var_names[i])==1){
      tree_cpy(variables.symbol_trees[i], output);
      break;
    }
  }

  // error if no variable was found
  if(i==variables.length){
    fprintf(stderr, "error: variable '$%s' not found\n",char_array_to_str_noinit(&varname));
    exit(-1);
  }
  free_Char_Array(varname);

  return(0);
}


// add a polynomials as a new named variable
int add_polynomial_to_variables(char* name, Polynomial polynomial, Variables* variables){
  // save polynomial to string (to convert it to a variable, it must first be a string)
  Char_Array poly_str;
  Char_Array out_name;
  Tree out_tree;

  init_Char_Array(&poly_str, STR_SIZE);
  polynomial_sprint(polynomial, &poly_str);

  // convert name to Char_Array
  init_Char_Array(&out_name,STR_SIZE);
  char_array_append_str(name, &out_name);

  // trivial tree containing the polynomial
  init_Tree(&out_tree,0,poly_str.length);
  tree_set_label(poly_str, &out_tree);
  free_Char_Array(poly_str);

  // add variable
  variables_append(out_name, out_tree, variables);
  free_Tree(out_tree);
  free_Char_Array(out_name);

  return 0;
}