meankondo/src/fields.c

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

// 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_Symbols(&((*fields).symbols), 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_Symbols(fields.symbols);
  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.symbols.indices, fields.symbols.length, index)>=0){
    return(FIELD_SYMBOL);
  }

  fprintf(stderr,"error: index %d is neither a parameter nor an external or an internal field, nor a symbol\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;

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


// ------------------ Symbols --------------------

// allocate memory
int init_Symbols(Symbols* symbols,int size){
  (*symbols).indices=calloc(size,sizeof(int));
  (*symbols).expr=calloc(size,sizeof(Polynomial));
  (*symbols).length=0;
  (*symbols).memory=size;
  return(0);
}

// free memory
int free_Symbols(Symbols symbols){
  int i;
  for(i=0;i<symbols.length;i++){
    free_Polynomial(symbols.expr[i]);
  }
  free(symbols.indices);
  free(symbols.expr);
  return(0);
}

// resize
int resize_symbols(Symbols* symbols,int new_size){
  Symbols new_symbols;
  int i;

  init_Symbols(&new_symbols,new_size);
  for(i=0;i<(*symbols).length;i++){
    new_symbols.indices[i]=(*symbols).indices[i];
    new_symbols.expr[i]=(*symbols).expr[i];
  }
  new_symbols.length=(*symbols).length;

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

  *symbols=new_symbols;
  return(0);
}

// copy
int symbols_cpy(Symbols input, Symbols* output){
  init_Symbols(output,input.length);
  symbols_cpy_noinit(input,output);
  return(0);
}
int symbols_cpy_noinit(Symbols input, Symbols* output){
  int i;
  if((*output).memory<input.length){
    fprintf(stderr,"error: trying to copy a symbols 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 symbols
int symbols_append(int index, Polynomial expr, Symbols* output){
  int offset=(*output).length;

  if((*output).length>=(*output).memory){
    resize_symbols(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 symbols without allocating memory
int symbols_append_noinit(int index, Polynomial expr, Symbols* output){
  int offset=(*output).length;

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

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

// concatenate two symbolss
int symbols_concat(Symbols input, Symbols* output){
  int i;
  for(i=0;i<input.length;i++){
    symbols_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);
}