Initial commitv1.2

1 files changed, 489 insertions, 0 deletions

diff --git a/src/fields.c b/src/fields.c
new file mode 100644
index 0000000..1b221f2
--- /dev/null
+++ b/src/fields.c
@@ -0,0 +1,489 @@
+/*
+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);
+}