path: root/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);
  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);
  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);
  }
}


// ------------------ 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
int resolve_ids(Polynomial* polynomial, Fields_Table fields){
  int i,j,k,l;
  int sign;
  int fermion_count;
  int at_least_one;
  int security;
  Int_Array 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
	if(int_array_is_subarray_ordered(fields.ids.lhs[j],(*polynomial).monomials[i])==1){
	  init_Int_Array(&monomial, (*polynomial).monomials[i].length);
	  
	  // remove lhs from monomial
	  // sign from moving the fields out of the monomial
	  sign=1;
	  // number of Fermions to remove from the monomial
	  fermion_count=0;
	  for(k=0,l=0;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]){
	      int_array_append((*polynomial).monomials[i].values[k],&monomial);
	      // sign correction
	      if(fermion_count % 2 ==1 && is_fermion((*polynomial).monomials[i].values[k], fields)){
		sign*=-1;
	      }
	    }
	    else{
	      // increment fermion_count
	      if(is_fermion(fields.ids.lhs[j].values[l],fields)){
		fermion_count++;
	      }
	      // 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(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);
	    // re-sort monomial
	    sign=1;
	    monomial_sort((*polynomial).monomials[(*polynomial).length-1],0,(*polynomial).monomials[(*polynomial).length-1].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]=monomial;
	  int_array_concat(fields.ids.rhs[j].monomials[0],(*polynomial).monomials+i);
	  // re-sort monomial
	  sign=1;
	  monomial_sort((*polynomial).monomials[i],0,(*polynomial).monomials[i].length-1,fields,&sign);
	  number_Qprod_chain(quot(sign,1),(*polynomial).nums+i);

	  // free num
	  free_Number(num);

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

  return(0);
}


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