diff options
Diffstat (limited to 'src/fields.c')
| -rw-r--r-- | src/fields.c | 489 |
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); +} |