Ian Jauslin
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorIan Jauslin <ian.jauslin@roma1.infn.it>2015-06-14 00:52:45 +0000
committerIan Jauslin <ian.jauslin@roma1.infn.it>2015-06-14 00:52:45 +0000
commitaa0f3ae2988d372b190b9bde2e75a6d17e744e93 (patch)
tree14482245c2fca27fcdad3078e97d0871352d52a7 /src/fields.c
Initial commitv1.2
Diffstat (limited to 'src/fields.c')
-rw-r--r--src/fields.c489
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);
+}