diff options
Diffstat (limited to 'src/mean.c')
| -rw-r--r-- | src/mean.c | 141 |
1 files changed, 109 insertions, 32 deletions
@@ -1,5 +1,5 @@ /* -Copyright 2015 Ian Jauslin +Copyright 2015-2022 Ian Jauslin Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -30,6 +30,7 @@ As of version 1.0, the mean of a monomial is computed directly #include "array.h" #include "fields.h" #include "number.h" +#include "determinant.h" // mean of a monomial int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Matrix propagator){ @@ -61,10 +62,79 @@ int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Ma // compute the mean of a monomial of internal fields (with split + and -) int mean_internal(Int_Array internal_plus, Int_Array internal_minus, Polynomial* out, Polynomial_Matrix propagator, Fields_Table fields){ + int ret; + Number num; + if(internal_plus.length!=internal_minus.length){ fprintf(stderr,"error: monomial contains unmatched +/- fields\n"); exit(-1); } + + ret=mean_determinant(internal_plus, internal_minus, &num, propagator, fields); + // cannot compute the mean as a determinant, use permutations + // can be because some fields are not Fermions + // can be because the propagator has non-numeric values (inverting polynomials is not implemented, and would be required for the computation of the determinant) + if(ret==-1){ + mean_permutations(internal_plus, internal_minus, out, propagator, fields); + } + else{ + polynomial_multiply_scalar(*out, num); + free_Number(num); + } + + return(0); +} + +// compute the mean of a monomial by computing a determinant +// can only be used if all of the propagators are numbers +int mean_determinant(Int_Array internal_plus, Int_Array internal_minus, Number* out, Polynomial_Matrix propagator, Fields_Table fields){ + Number_Matrix M; + int n=internal_minus.length; + int i,j; + int a,b; + int sign; + + init_Number_Matrix(&M,n); + + // extra sign: the monomial is sorted in such a way that minus fields are on the left of plus fields, but the determinant formula requires the fields to be alternated +- + if((n+1)/2%2==1){ + sign=-1; + } + else{ + sign=1; + } + + // construct matrix + for(i=0;i<n;i++){ + a=intlist_find_err(propagator.indices, propagator.length, internal_plus.values[i]); + for(j=0;j<n;j++){ + b=intlist_find_err(propagator.indices, propagator.length, -internal_minus.values[j]); + // ignore 0 + if(propagator.matrix[a][b].length!=0){ + // check whether the fields are Fermions, and whether the entry is a number + if(is_fermion(internal_plus.values[i], fields)==0 || is_fermion(internal_minus.values[j], fields)==0 || polynomial_is_number(propagator.matrix[a][b])==0){ + free_Number_Matrix(M); + return(-1); + } + + number_add_chain(propagator.matrix[a][b].nums[0], M.matrix[i]+j); + } + } + } + + // compute determinant + determinant_inplace(M, out); + + number_Qprod_chain(quot(sign,1), out); + + free_Number_Matrix(M); + + return(0); +} + + +// compute the mean of a monomial by summing over permutations +int mean_permutations(Int_Array internal_plus, Int_Array internal_minus, Polynomial* out, Polynomial_Matrix propagator, Fields_Table fields){ int n=internal_minus.length; // pairing as an array of positions int* pairing=calloc(n,sizeof(int)); @@ -118,7 +188,7 @@ int mean_internal(Int_Array internal_plus, Int_Array internal_minus, Polynomial* pairing_sign=permutation_signature(pairing,n); } - // only simplify in mean_symbols + // only simplify in mean_virtual_fields polynomial_prod_chain_nosimplify(num_summed,out,fields); free_Polynomial(num_summed); free(pairing); @@ -346,10 +416,10 @@ int get_internals(Int_Array monomial, Int_Array* internal_plus, Int_Array* inter } -// compute the mean of a monomial containing symbolic expressions +// compute the mean of a monomial containing virtual fields // keep track of which means were already computed -int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed){ - Int_Array symbol_list; +int mean_virtual_fields(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed){ + Int_Array virtual_field_list; int i; int power; int* current_term; @@ -375,43 +445,43 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po } } - init_Int_Array(&symbol_list, monomial.length); + init_Int_Array(&virtual_field_list, monomial.length); init_Int_Array(&base_monomial, monomial.length); - // generate symbols list + // generate virtual_fields list for(i=0;i<monomial.length;i++){ - if(field_type(monomial.values[i], fields)==FIELD_SYMBOL){ - int_array_append(intlist_find_err(fields.symbols.indices, fields.symbols.length, monomial.values[i]), &symbol_list); + if(field_type(monomial.values[i], fields)==FIELD_VIRTUAL){ + int_array_append(intlist_find_err(fields.virtual_fields.indices, fields.virtual_fields.length, monomial.values[i]), &virtual_field_list); } else{ int_array_append(monomial.values[i], &base_monomial); } } - power=symbol_list.length; + power=virtual_field_list.length; // trivial case if(power==0){ mean(monomial, &mean_num, fields, propagator); polynomial_concat_noinit(mean_num, output); - free_Int_Array(symbol_list); + free_Int_Array(virtual_field_list); free_Int_Array(base_monomial); return(0); } else{ // initialize current term to a position that has no repetitions current_term=calloc(power,sizeof(int)); - exists_next=init_prod(current_term, symbol_list, fields, power, base_monomial)+1; + exists_next=init_prod(current_term, virtual_field_list, fields, power, base_monomial)+1; } - // loop over terms; the loop stops when all the pointers are at the end of the first symbol + // loop over terms; the loop stops when all the pointers are at the end of the first virtual field while(exists_next==1){ // construct monomial int_array_cpy(base_monomial, &tmp_monomial); tmp_num=number_one(); for(i=0;i<power;i++){ - int_array_concat(fields.symbols.expr[symbol_list.values[i]].monomials[current_term[i]], &tmp_monomial); - number_prod_chain(fields.symbols.expr[symbol_list.values[i]].nums[current_term[i]], &tmp_num); + int_array_concat(fields.virtual_fields.expr[virtual_field_list.values[i]].monomials[current_term[i]], &tmp_monomial); + number_prod_chain(fields.virtual_fields.expr[virtual_field_list.values[i]].nums[current_term[i]], &tmp_num); } // check whether the monomial vanishes if(check_monomial_match(tmp_monomial, fields)==1){ @@ -432,7 +502,7 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po free_Int_Array(tmp_monomial); // next term - exists_next=next_prod(current_term, symbol_list, fields, power, base_monomial)+1; + exists_next=next_prod(current_term, virtual_field_list, fields, power, base_monomial)+1; // simplfiy every 25 steps (improves both memory usage and performance) @@ -451,13 +521,13 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po // free memory free(current_term); - free_Int_Array(symbol_list); + free_Int_Array(virtual_field_list); free_Int_Array(base_monomial); return(0); } // first term in product with no repetitions -int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial){ +int init_prod(int* current_term, Int_Array virtual_field_list, Fields_Table fields, int power, Int_Array base_monomial){ // index we want to increment int move=0; // tmp monomial @@ -474,7 +544,7 @@ int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int // loop until move is out of range while(move>=0 && move<power){ // move - current_term[move]=next_term_norepeat(current_term[move], fields.symbols.expr[symbol_list.values[move]], &monomial, fields); + current_term[move]=next_term_norepeat(current_term[move], fields.virtual_fields.expr[virtual_field_list.values[move]], &monomial, fields); // if the next term does not exist, then move previous index if(current_term[move]==-1){ move--; @@ -495,7 +565,7 @@ int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int } // next term in product with no repetitions -int next_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial){ +int next_prod(int* current_term, Int_Array virtual_field_list, Fields_Table fields, int power, Int_Array base_monomial){ // index we want to increment int move=power-1; // tmp monomial @@ -506,13 +576,13 @@ int next_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int init_Int_Array(&monomial, base_monomial.length+5*power); int_array_cpy_noinit(base_monomial, &monomial); for(i=0;i<=move;i++){ - int_array_concat(fields.symbols.expr[symbol_list.values[i]].monomials[current_term[i]],&monomial); + int_array_concat(fields.virtual_fields.expr[virtual_field_list.values[i]].monomials[current_term[i]],&monomial); } // loop until move is out of range while(move>=0 && move<power){ // move - current_term[move]=next_term_norepeat(current_term[move], fields.symbols.expr[symbol_list.values[move]], &monomial, fields); + current_term[move]=next_term_norepeat(current_term[move], fields.virtual_fields.expr[virtual_field_list.values[move]], &monomial, fields); // if the next term does not exist, then move previous index if(current_term[move]==-1){ move--; @@ -623,13 +693,13 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol int group=-2; int next_group=-2; Polynomial tmp_poly; - int sign; + int sign=1; init_Polynomial(output, MONOMIAL_SIZE); - // check whether there are symbols - // IMPORTANT: the symbols must be at the end of the monomial - if(monomial.length==0 || field_type(monomial.values[monomial.length-1], fields)!=FIELD_SYMBOL){ + // check whether there are virtual fields + // IMPORTANT: the virtual fields must be at the end of the monomial + if(monomial.length==0 || field_type(monomial.values[monomial.length-1], fields)!=FIELD_VIRTUAL){ // mean mean(monomial, &num_mean, fields, propagator); // add to output @@ -650,7 +720,7 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol } // if group changes, take mean if((i>0 && next_group!=group) || i==monomial.length){ - mean_symbols(tmp_monomial, &tmp_poly, fields, propagator, groups, computed); + mean_virtual_fields(tmp_monomial, &tmp_poly, fields, propagator, groups, computed); // if zero if(polynomial_is_zero(tmp_poly)==1){ // set output to 0 @@ -702,10 +772,11 @@ struct mean_args{ Fields_Table fields; Polynomial_Matrix propagator; Groups groups; + int print_progress; }; // multithreaded -int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads){ +int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, int print_progress){ int i; Polynomial thread_polys[threads]; pthread_t thread_ids[threads]; @@ -720,11 +791,15 @@ int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Pol args[i].fields=fields; args[i].propagator=propagator; args[i].groups=groups; + args[i].print_progress=print_progress; } // split polynomial + // randomly choose the thread + // see random number generator + srand(time(NULL)); for(i=0;i<len;i++){ - polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], thread_polys+(i % threads)); + polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], thread_polys+(rand() % threads)); } // start threads @@ -750,12 +825,12 @@ int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Pol // mean for one of the threads void* polynomial_mean_thread(void* mean_args){ struct mean_args *args=mean_args; - polynomial_mean((*args).polynomial,(*args).fields,(*args).propagator,(*args).groups); + polynomial_mean((*args).polynomial,(*args).fields,(*args).propagator,(*args).groups, (*args).print_progress); return(NULL); } // single threaded version -int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups){ +int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int print_progress){ int i,j; Polynomial output; Polynomial tmp_poly; @@ -769,7 +844,9 @@ int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matr // mean of each monomial for(i=0;i<(*polynomial).length;i++){ - fprintf(stderr,"computing %d of %d means\n",i,(*polynomial).length-1); + if(print_progress==1){ + fprintf(stderr,"computing %d of %d means\n",i,(*polynomial).length-1); + } mean_groups((*polynomial).monomials[i], &tmp_poly, fields, propagator, groups, &computed); // write factors |