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/mean.c
Initial commitv1.2
Diffstat (limited to 'src/mean.c')
-rw-r--r--src/mean.c793
1 files changed, 793 insertions, 0 deletions
diff --git a/src/mean.c b/src/mean.c
new file mode 100644
index 0000000..aad7a5a
--- /dev/null
+++ b/src/mean.c
@@ -0,0 +1,793 @@
+/*
+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.
+*/
+
+/*
+As of version 1.0, the mean of a monomial is computed directly
+*/
+
+#include "mean.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <pthread.h>
+#include "definitions.cpp"
+#include "tools.h"
+#include "polynomial.h"
+#include "rational.h"
+#include "array.h"
+#include "fields.h"
+#include "number.h"
+
+// mean of a monomial
+int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Matrix propagator){
+ int sign=1;
+ // +/- internal fields
+ Int_Array internal_plus;
+ Int_Array internal_minus;
+
+ // init out
+ *out=polynomial_one();
+
+ // sort first
+ monomial_sort(monomial, 0, monomial.length-1, fields, &sign);
+ polynomial_multiply_Qscalar(*out, quot(sign,1));
+ // get internals
+ // (*out).monomials is the first element of out but it only has 1 element
+ // first, free (*out).monomials[0]
+ free_Int_Array((*out).monomials[0]);
+ get_internals(monomial, &internal_plus, &internal_minus, (*out).monomials, fields);
+
+ if(internal_plus.length>0 && internal_minus.length>0){
+ mean_internal(internal_plus, internal_minus, out, propagator, fields);
+ }
+
+ free_Int_Array(internal_plus);
+ free_Int_Array(internal_minus);
+ return(0);
+}
+
+// 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){
+ if(internal_plus.length!=internal_minus.length){
+ fprintf(stderr,"error: monomial contains unmatched +/- fields\n");
+ exit(-1);
+ }
+ int n=internal_minus.length;
+ // pairing as an array of positions
+ int* pairing=calloc(n,sizeof(int));
+ // specifies which indices are available for pairing
+ int* mask=calloc(n,sizeof(int));
+ // signature of the permutation
+ int pairing_sign;
+ // sign from mixing - and + together
+ int mixing_sign;
+ Polynomial num;
+ Polynomial num_summed=polynomial_zero();
+ // propagator matrix indices
+ int* indices_minus=calloc(n,sizeof(int));
+ int* indices_plus=calloc(n,sizeof(int));
+ int i;
+ // whether the next pairing exists
+ int exists_next=0;
+
+ // indices
+ for(i=0;i<n;i++){
+ indices_plus[i]=intlist_find_err(propagator.indices, propagator.length, internal_plus.values[i]);
+ indices_minus[i]=intlist_find_err(propagator.indices, propagator.length, -internal_minus.values[i]);
+ }
+
+ // init pairing and mask
+ exists_next=init_pairing(pairing, mask, n, propagator, indices_plus, indices_minus)+1;
+
+ // initial sign
+ pairing_sign=permutation_signature(pairing,n);
+
+ // mixing sign (from ordering psi+psi-): (-1)^{n(n+1)/2}
+ if((n*(n+1))/2 %2 ==0){
+ mixing_sign=1;
+ }
+ else{
+ mixing_sign=-1;
+ }
+
+ // loop over pairings
+ // loop ends when the first pairing leaves the array
+ while(exists_next==1){
+ num=polynomial_one();
+ // propagator product for the current pairing (only simplify after all pairings)
+ for(i=0;i<n;i++){
+ polynomial_prod_chain_nosimplify(propagator.matrix[indices_plus[i]][indices_minus[pairing[i]]],&num, fields);
+ }
+ polynomial_multiply_Qscalar(num,quot(mixing_sign*pairing_sign,1));
+ polynomial_concat_noinit_inplace(num,&num_summed);
+
+ exists_next=next_pairing(pairing, mask, n, propagator, indices_plus, indices_minus)+1;
+ pairing_sign=permutation_signature(pairing,n);
+ }
+
+ // only simplify in mean_symbols
+ polynomial_prod_chain_nosimplify(num_summed,out,fields);
+ free_Polynomial(num_summed);
+ free(pairing);
+ free(mask);
+ free(indices_plus);
+ free(indices_minus);
+ return(0);
+}
+
+// first pairing with a non-vanishing propagator
+int init_pairing(int* pairing, int* mask, int n, Polynomial_Matrix propagator, int* indices_plus, int* indices_minus){
+ // index we want to increment
+ int move=0;
+ int i;
+ for(i=0;i<n;i++){
+ pairing[i]=-1;
+ mask[i]=0;
+ }
+
+ // loop until move is out of range
+ while(move>=0 && move<n){
+ // move
+ pairing[move]=next_wick(move, pairing[move], mask, n, propagator, indices_plus, indices_minus);
+ // if the next term does not exist, then move previous index
+ if(pairing[move]==-1){
+ move--;
+ }
+ // else move next index
+ else{
+ move++;
+ }
+ }
+
+ // if move==-1, then there is no first term, return -1
+ if(move==-1){
+ return(-1);
+ }
+ // if the first term exists
+ return(0);
+}
+
+// next pairing with a non-vanishing propagator
+int next_pairing(int* pairing, int* mask, int n, Polynomial_Matrix propagator, int* indices_plus, int* indices_minus){
+ // index we want to increment
+ int move=n-1;
+
+ // last index
+ mask[pairing[n-1]]=0;
+
+ // loop until move is out of range
+ while(move>=0 && move<n){
+ // move
+ pairing[move]=next_wick(move, pairing[move], mask, n, propagator, indices_plus, indices_minus);
+ // if the next term does not exist, then move previous index
+ if(pairing[move]==-1){
+ move--;
+ }
+ // else move next index
+ else{
+ move++;
+ }
+ }
+
+ // if move==-1, then there is no next term, return -1
+ if(move==-1){
+ return(-1);
+ }
+ // if the next term exists
+ return(0);
+}
+
+// next term in the Wick expansion
+int next_wick(int index, int start, int* mask, int n, Polynomial_Matrix propagator, int* indices_plus, int* indices_minus){
+ int i;
+ // unset mask
+ if(start>=0 && start<n){
+ mask[start]=0;
+ }
+ // find next position
+ for(i=start+1;i<n;i++){
+ // if the propagator doesn't vanish
+ if(mask[i]==0 && polynomial_is_zero(propagator.matrix[indices_plus[index]][indices_minus[i]])==0){
+ mask[i]=1;
+ return(i);
+ }
+ }
+ // no next term
+ return(-1);
+}
+
+
+/* Older function: propagator as number
+// compute the mean of a monomial of internal fields (with split + and -)
+// compute all contractions
+int mean_internal_slow(Int_Array internal_plus, Int_Array internal_minus, Number* outnum, Polynomial_Matrix propagator){
+ if(internal_plus.length!=internal_minus.length){
+ fprintf(stderr,"error: monomial contains unmatched +/- fields\n");
+ exit(-1);
+ }
+ int n=internal_minus.length;
+ // pairing as an array of positions
+ int* pairing=calloc(n,sizeof(int));
+ // specifies which indices are available for pairing
+ int* mask=calloc(n,sizeof(int));
+ // signature of the permutation
+ int pairing_sign;
+ // sign from mixing - and + together
+ int mixing_sign;
+ Number num;
+ Number num_summed=number_zero();
+ // propagator matrix indices
+ int index1, index2;
+ int i,j,k,l;
+
+ // init pairing and mask
+ for(i=0;i<n-1;i++){
+ pairing[i]=i;
+ mask[i]=1;
+ }
+ pairing[n-1]=n-1;
+ pairing_sign=1;
+
+ // mixing sign: (-1)^{n(n+1)/2}
+ if((n*(n+1))/2 %2 ==0){
+ mixing_sign=1;
+ }
+ else{
+ mixing_sign=-1;
+ }
+
+ // loop over pairings
+ // loop ends when the first pairing leaves the array
+ while(pairing[0]<n){
+ num=number_one();
+ // propagator product for the current pairing
+ for(i=0;i<n;i++){
+ // indices within the propagator matrix
+ index1=intlist_find_err(propagator.indices, propagator.length, internal_plus.values[i]);
+ index2=intlist_find_err(propagator.indices, propagator.length, -internal_minus.values[pairing[i]]);
+
+ number_prod_chain(propagator.matrix[index1][index2],&num);
+ }
+ number_Qprod_chain(quot(mixing_sign*pairing_sign,1),&num);
+ number_add_chain(num,&num_summed);
+ free_Number(num);
+
+ // next pairing
+ // last element of the pairing that we can move
+ for(i=n-1;i>=0;i--){
+ // move i-th
+ mask[pairing[i]]=0;
+ // search for next possible position
+ for(j=pairing[i]+1;j<n;j++){
+ if(mask[j]==0){
+ break;
+ }
+ }
+
+ // if the next position exists
+ if(j<n){
+ // sign correction: change sign by (-1)^{1+(n-i)(n-i-1)/2}
+ // actually (-1)^{1+(n-1-i)(n-1-i-1)/2
+ if(((n-i-1)*(n-i-2))/2 % 2==0){
+ pairing_sign*=-1;
+ }
+
+ pairing[i]=j;
+ mask[j]=1;
+ // put the remaining pairings at their left-most possible values
+ if(i<n-1){
+ k=i+1;
+ for(l=0;l<n;l++){
+ if(mask[l]==0){
+ mask[l]=1;
+ pairing[k]=l;
+ k++;
+ // if exhausted all indices
+ if(k>=n){
+ break;
+ }
+ }
+ }
+ }
+ // if the next position was found, then don't try to move the previous pairings
+ break;
+ }
+ // if no next position is found, store the pairing outside the array (so the algorithm stops when the first pairing is outside the array)
+ else{
+ pairing[i]=n;
+ }
+ }
+ }
+
+ number_prod_chain(num_summed,outnum);
+ free_Number(num_summed);
+ free(pairing);
+ free(mask);
+ return(0);
+}
+*/
+
+
+// get lists of internal fields from a monomial (separate + and -)
+// requires the monomial to be sorted (for the sign to be correct)
+int get_internals(Int_Array monomial, Int_Array* internal_plus, Int_Array* internal_minus, Int_Array* others, Fields_Table fields){
+ int i;
+ init_Int_Array(internal_plus, monomial.length);
+ init_Int_Array(internal_minus, monomial.length);
+ init_Int_Array(others, monomial.length);
+ for (i=0;i<monomial.length;i++){
+ if(int_array_find(abs(monomial.values[i]),fields.internal)>=0){
+ // split +/- fields
+ if(monomial.values[i]>0){
+ int_array_append(monomial.values[i],internal_plus);
+ }
+ else{
+ int_array_append(monomial.values[i],internal_minus);
+ }
+ }
+ else{
+ int_array_append(monomial.values[i], others);
+ }
+ }
+ return(0);
+}
+
+
+// compute the mean of a monomial containing symbolic expressions
+// 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 i;
+ int power;
+ int* current_term;
+ Polynomial mean_num;
+ Int_Array tmp_monomial;
+ Number tmp_num;
+ Int_Array base_monomial;
+ int sign;
+ // whether or not the next term exists
+ int exists_next=0;
+ // simplify polynomial periodically
+ int simplify_freq=1;
+ Polynomial mean_poly;
+
+ init_Polynomial(output, POLY_SIZE);
+
+ // check whether the mean was already computed
+ for(i=0;i<(*computed).length;i++){
+ if(int_array_cmp((*computed).lhs[i], monomial)==0){
+ // write polynomial
+ polynomial_concat((*computed).rhs[i], output);
+ return(0);
+ }
+ }
+
+ init_Int_Array(&symbol_list, monomial.length);
+ init_Int_Array(&base_monomial, monomial.length);
+
+ // generate symbols 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);
+ }
+ else{
+ int_array_append(monomial.values[i], &base_monomial);
+ }
+ }
+ power=symbol_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(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;
+ }
+
+ // loop over terms; the loop stops when all the pointers are at the end of the first symbol
+ 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);
+ }
+ // check whether the monomial vanishes
+ if(check_monomial_match(tmp_monomial, fields)==1){
+ // sort monomial
+ sign=1;
+ monomial_sort(tmp_monomial, 0, tmp_monomial.length-1, fields, &sign);
+ number_Qprod_chain(quot(sign,1), &tmp_num);
+
+ // mean
+ mean_groups(tmp_monomial, &mean_poly, fields, propagator, groups, computed);
+
+ // write to output
+ polynomial_multiply_scalar(mean_poly,tmp_num);
+ polynomial_concat_noinit_inplace(mean_poly, output);
+ }
+
+ free_Number(tmp_num);
+ free_Int_Array(tmp_monomial);
+
+ // next term
+ exists_next=next_prod(current_term, symbol_list, fields, power, base_monomial)+1;
+
+
+ // simplfiy every 25 steps (improves both memory usage and performance)
+ if(simplify_freq %25 ==0){
+ polynomial_simplify(output, fields);
+ simplify_freq=0;
+ }
+ simplify_freq++;
+ }
+
+ // simplify
+ polynomial_simplify(output, fields);
+
+ // write computed
+ identities_append(monomial, *output, computed);
+
+ // free memory
+ free(current_term);
+ free_Int_Array(symbol_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){
+ // index we want to increment
+ int move=0;
+ // tmp monomial
+ Int_Array monomial;
+ int i;
+
+ init_Int_Array(&monomial, base_monomial.length+5*power);
+ int_array_cpy_noinit(base_monomial, &monomial);
+ // init current term
+ for(i=0;i<power;i++){
+ current_term[i]=-1;
+ }
+
+ // 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);
+ // if the next term does not exist, then move previous index
+ if(current_term[move]==-1){
+ move--;
+ }
+ // else move next index
+ else{
+ move++;
+ }
+ }
+
+ free_Int_Array(monomial);
+ // if move==-1, then there is no first term, return -1
+ if(move==-1){
+ return(-1);
+ }
+ // if the next term exists
+ return(0);
+}
+
+// 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){
+ // index we want to increment
+ int move=power-1;
+ // tmp monomial
+ Int_Array monomial;
+ int i;
+
+ // init monomial
+ 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);
+ }
+
+ // 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);
+ // if the next term does not exist, then move previous index
+ if(current_term[move]==-1){
+ move--;
+ }
+ // else move next index
+ else{
+ move++;
+ }
+ }
+
+ free_Int_Array(monomial);
+ // if move==-1, then there is no next term, return -1
+ if(move==-1){
+ return(-1);
+ }
+ // if the next term exists
+ return(0);
+}
+
+// find the next term in a polynomial that can be multiplied to a given monomial and add it to the monomial
+int next_term_norepeat(int start, Polynomial polynomial, Int_Array* monomial, Fields_Table fields){
+ int i;
+ // remove last term from monomial
+ if(start>=0 && start<polynomial.length){
+ (*monomial).length-=polynomial.monomials[start].length;
+ }
+ // find next position
+ for(i=start+1;i<polynomial.length;i++){
+ // if no repetitions
+ if(check_monomial_addterm(*monomial,polynomial.monomials[i],fields)==1){
+ // append to monomial
+ int_array_concat(polynomial.monomials[i], monomial);
+ return(i);
+ }
+ }
+ // no next term
+ return(-1);
+}
+
+
+// signature of a permutation
+int permutation_signature(int* permutation, int n){
+ int* tmp_array=calloc(n,sizeof(int));
+ int i;
+ int ret=1;
+ for(i=0;i<n;i++){
+ tmp_array[i]=permutation[i];
+ }
+ sort_fermions(tmp_array, 0, n-1, &ret);
+ free(tmp_array);
+ return(ret);
+}
+
+// sort a list of anti-commuting variables
+int sort_fermions(int* array, int begin, int end, int* sign){
+ int i;
+ int tmp;
+ int index;
+ // the pivot: middle of the monomial
+ int pivot=(begin+end)/2;
+
+ // if the monomial is non trivial
+ if(begin<end){
+ // send pivot to the end
+ if(pivot!=end){
+ tmp=array[end];
+ array[end]=array[pivot];
+ array[pivot]=tmp;
+ *sign*=-1;
+ }
+
+ // loop over the others
+ for(i=begin, index=begin;i<end;i++){
+ // compare with pivot
+ if(array[i]<array[end]){
+ // if smaller, exchange with reference index
+ if(i!=index){
+ tmp=array[i];
+ array[i]=array[index];
+ array[index]=tmp;
+ *sign*=-1;
+ }
+ // move reference index
+ index++;
+ }
+ }
+ // put pivot (which we had sent to the end) in the right place
+ if(end!=index){
+ tmp=array[end];
+ array[end]=array[index];
+ array[index]=tmp;
+ *sign*=-1;
+ }
+
+ // recurse
+ sort_fermions(array, begin, index-1, sign);
+ sort_fermions(array, index+1, end, sign);
+ }
+ return(0);
+}
+
+
+// mean while factoring groups out
+int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed){
+ Polynomial num_mean;
+ Int_Array tmp_monomial;
+ int i;
+ int group=-2;
+ int next_group=-2;
+ Polynomial tmp_poly;
+ int sign;
+
+ init_Polynomial(output, MONOMIAL_SIZE);
+
+ // check whether there are symbols
+ // requires the symbols to be at the end of the monomial
+ if(monomial.length==0 || field_type(monomial.values[monomial.length-1], fields)!=FIELD_SYMBOL){
+ // mean
+ mean(monomial, &num_mean, fields, propagator);
+ // add to output
+ polynomial_concat_noinit(num_mean,output);
+ }
+ else{
+ // sort into groups
+ if(groups.length>0){
+ sign=1;
+ monomial_sort_groups(monomial, 0, monomial.length-1, fields, groups, &sign);
+ }
+ // construct groups and take mean
+ init_Int_Array(&tmp_monomial, MONOMIAL_SIZE);
+ for(i=0;i<=monomial.length;i++){
+ // new group
+ if(i<monomial.length){
+ next_group=find_group(monomial.values[i], groups);
+ }
+ // if group changes, take mean
+ if((i>0 && next_group!=group) || i==monomial.length){
+ mean_symbols(tmp_monomial, &tmp_poly, fields, propagator, groups, computed);
+ // if zero
+ if(polynomial_is_zero(tmp_poly)==1){
+ // set output to 0
+ free_Polynomial(*output);
+ init_Polynomial(output, 1);
+ free_Polynomial(tmp_poly);
+ break;
+ }
+ // add to output
+ if(polynomial_is_zero(*output)==1){
+ polynomial_concat(tmp_poly, output);
+ }
+ else{
+ polynomial_prod_chain_nosimplify(tmp_poly, output, fields);
+ }
+ free_Polynomial(tmp_poly);
+
+ // reset tmp_monomial
+ free_Int_Array(tmp_monomial);
+ init_Int_Array(&tmp_monomial, MONOMIAL_SIZE);
+ }
+
+ // add to monomial
+ if(i<monomial.length){
+ int_array_append(monomial.values[i], &tmp_monomial);
+ }
+ group=next_group;
+ }
+
+ // sign correction
+ if(sign==-1){
+ polynomial_multiply_Qscalar(*output,quot(sign,1));
+ }
+
+ free_Int_Array(tmp_monomial);
+
+ }
+
+ return(0);
+}
+
+
+
+// mean of a polynomial
+
+// argument struct for multithreaded mean
+struct mean_args{
+ Polynomial* polynomial;
+ Fields_Table fields;
+ Polynomial_Matrix propagator;
+ Groups groups;
+};
+
+// multithreaded
+int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads){
+ int i;
+ Polynomial thread_polys[threads];
+ pthread_t thread_ids[threads];
+ struct mean_args args[threads];
+ int len=(*polynomial).length;
+
+
+ // alloc
+ for(i=0;i<threads;i++){
+ init_Polynomial(thread_polys+i,(*polynomial).length/threads+1);
+ // arguments
+ args[i].fields=fields;
+ args[i].propagator=propagator;
+ args[i].groups=groups;
+ }
+
+ // split polynomial
+ for(i=0;i<len;i++){
+ polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], thread_polys+(i % threads));
+ }
+
+ // start threads
+ for(i=0;i<threads;i++){
+ args[i].polynomial=thread_polys+i;
+ pthread_create(thread_ids+i, NULL, polynomial_mean_thread, (void*)(args+i));
+ }
+
+ free_Polynomial(*polynomial);
+ init_Polynomial(polynomial, len);
+
+ // wait for completion and join
+ for(i=0;i<threads;i++){
+ pthread_join(thread_ids[i], NULL);
+ polynomial_concat_noinit(thread_polys[i], polynomial);
+ }
+
+ polynomial_simplify(polynomial, fields);
+
+ return(0);
+}
+
+// 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);
+ return(NULL);
+}
+
+// single threaded version
+int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups){
+ int i,j;
+ Polynomial output;
+ Polynomial tmp_poly;
+ // a list of already computed means
+ Identities computed;
+
+ init_Polynomial(&output, (*polynomial).length);
+ init_Identities(&computed, EQUATION_SIZE);
+
+ remove_unmatched_plusminus(polynomial, fields);
+
+ // mean of each monomial
+ for(i=0;i<(*polynomial).length;i++){
+ 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
+ for(j=0;j<tmp_poly.length;j++){
+ int_array_concat((*polynomial).factors[i], tmp_poly.factors+j);
+ number_prod_chain((*polynomial).nums[i], tmp_poly.nums+j);
+ }
+
+ // add to output
+ polynomial_concat_noinit(tmp_poly, &output);
+ // simplify (simplify here in order to keep memory usage low)
+ polynomial_simplify(&output, fields);
+ }
+ free_Identities(computed);
+
+ free_Polynomial(*polynomial);
+ *polynomial=output;
+
+ return(0);
+}
+