diff options
Diffstat (limited to 'src/polynomial.c')
| -rw-r--r-- | src/polynomial.c | 1263 |
1 files changed, 1263 insertions, 0 deletions
diff --git a/src/polynomial.c b/src/polynomial.c new file mode 100644 index 0000000..639728a --- /dev/null +++ b/src/polynomial.c @@ -0,0 +1,1263 @@ +/* +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 "polynomial.h" + +#include <stdio.h> +#include <stdlib.h> +#include "definitions.cpp" +#include "rational.h" +#include "tools.h" +#include "mean.h" +#include "coefficient.h" +#include "istring.h" +#include "array.h" +#include "number.h" +#include "fields.h" + + +// allocate memory +int init_Polynomial(Polynomial* polynomial,int size){ + (*polynomial).monomials=calloc(size,sizeof(Int_Array)); + (*polynomial).factors=calloc(size,sizeof(Int_Array)); + (*polynomial).nums=calloc(size,sizeof(Number)); + (*polynomial).length=0; + (*polynomial).memory=size; + return(0); +} + +// free memory +int free_Polynomial(Polynomial polynomial){ + int i; + for(i=0;i<polynomial.length;i++){ + free_Int_Array(polynomial.monomials[i]); + free_Int_Array(polynomial.factors[i]); + free_Number(polynomial.nums[i]); + } + free(polynomial.monomials); + free(polynomial.factors); + free(polynomial.nums); + + return(0); +} + +// resize the memory allocated to a polynomial +int resize_polynomial(Polynomial* polynomial,int new_size){ + Polynomial new_poly; + int i; + + init_Polynomial(&new_poly,new_size); + for(i=0;i<(*polynomial).length;i++){ + new_poly.monomials[i]=(*polynomial).monomials[i]; + new_poly.factors[i]=(*polynomial).factors[i]; + new_poly.nums[i]=(*polynomial).nums[i]; + } + new_poly.length=(*polynomial).length; + + free((*polynomial).monomials); + free((*polynomial).factors); + free((*polynomial).nums); + + *polynomial=new_poly; + return(0); +} + +// copy a polynomial +int polynomial_cpy(Polynomial input, Polynomial* output){ + init_Polynomial(output,input.length); + polynomial_cpy_noinit(input,output); + return(0); +} +int polynomial_cpy_noinit(Polynomial input, Polynomial* output){ + int i; + if((*output).memory<input.length){ + fprintf(stderr,"error: trying to copy a polynomial 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.monomials[i],(*output).monomials+i); + int_array_cpy(input.factors[i],(*output).factors+i); + number_cpy(input.nums[i],(*output).nums+i); + } + (*output).length=input.length; + + return(0); +} + +// append an element to a polynomial +int polynomial_append(Int_Array monomial, Int_Array factor, Number num, Polynomial* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_polynomial(output,2*(*output).memory+1); + } + + // copy and allocate + int_array_cpy(monomial,(*output).monomials+offset); + int_array_cpy(factor,(*output).factors+offset); + number_cpy(num,(*output).nums+offset); + //increment length + (*output).length++; + + return(0); +} +// append an element to a polynomial without allocating memory +int polynomial_append_noinit(Int_Array monomial, Int_Array factor, Number num, Polynomial* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_polynomial(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).monomials[offset]=monomial; + (*output).factors[offset]=factor; + (*output).nums[offset]=num; + // increment length + (*output).length++; + return(0); +} +// noinit, and if there already exists an element with the same monomial and factor, then just add numbers +int polynomial_append_noinit_inplace(Int_Array monomial, Int_Array factor, Number num, Polynomial* output){ + int i; + int foundit=0; + for(i=0;i<(*output).length;i++){ + if(int_array_cmp(monomial, (*output).monomials[i])==0 && int_array_cmp(factor, (*output).factors[i])==0){ + number_add_chain(num,(*output).nums+i); + foundit=1; + free_Number(num); + free_Int_Array(monomial); + free_Int_Array(factor); + break; + } + } + if(foundit==0){ + polynomial_append_noinit(monomial, factor, num, output); + } + return(0); +} + +// concatenate two polynomials +int polynomial_concat(Polynomial input, Polynomial* output){ + int i; + for(i=0;i<input.length;i++){ + polynomial_append(input.monomials[i],input.factors[i],input.nums[i],output); + } + return(0); +} +// noinit +int polynomial_concat_noinit(Polynomial input, Polynomial* output){ + int i; + for(i=0;i<input.length;i++){ + polynomial_append_noinit(input.monomials[i],input.factors[i],input.nums[i],output); + } + + // free input arrays + free(input.monomials); + free(input.factors); + free(input.nums); + return(0); +} +// noinit, inplace +int polynomial_concat_noinit_inplace(Polynomial input, Polynomial* output){ + int i; + for(i=0;i<input.length;i++){ + polynomial_append_noinit_inplace(input.monomials[i],input.factors[i],input.nums[i],output); + } + + // free input arrays + free(input.monomials); + free(input.factors); + free(input.nums); + return(0); +} + +// add polynomials +int polynomial_add_chain(Polynomial input, Polynomial* inout, Fields_Table fields){ + polynomial_concat(input,inout); + polynomial_simplify(inout, fields); + return(0); +} +// add polynomials (noinit) +int polynomial_add_chain_noinit(Polynomial input, Polynomial* inout, Fields_Table fields){ + polynomial_concat_noinit(input,inout); + polynomial_simplify(inout, fields); + return(0); +} + +// multiply a polynomial by a scalar +int polynomial_multiply_scalar(Polynomial polynomial, Number num){ + int i; + for(i=0;i<polynomial.length;i++){ + number_prod_chain(num,polynomial.nums+i); + } + return(0); +} +// multiply a polynomial by a rational number +int polynomial_multiply_Qscalar(Polynomial polynomial, Q q){ + int i; + for(i=0;i<polynomial.length;i++){ + number_Qprod_chain(q,polynomial.nums+i); + } + return(0); +} + +// change the sign of the monomials in a polynomial +int polynomial_conjugate(Polynomial polynomial){ + int i,j; + for(i=0;i<polynomial.length;i++){ + for(j=0;j<polynomial.monomials[i].length;j++){ + polynomial.monomials[i].values[j]*=-1; + } + } + return(0); +} + + +// returns an initialized polynomial, equal to 1 +Polynomial polynomial_one(){ + Polynomial ret; + Int_Array dummy_monomial; + Int_Array dummy_factor; + + init_Polynomial(&ret,1); + init_Int_Array(&dummy_monomial,1); + init_Int_Array(&dummy_factor,1); + polynomial_append_noinit(dummy_monomial,dummy_factor,number_one(),&ret); + + return(ret); +} + +// returns an initialized polynomial, equal to 0 +Polynomial polynomial_zero(){ + Polynomial ret; + Int_Array dummy_monomial; + Int_Array dummy_factor; + + init_Polynomial(&ret,1); + init_Int_Array(&dummy_monomial,1); + init_Int_Array(&dummy_factor,1); + polynomial_append_noinit(dummy_monomial,dummy_factor,number_zero(),&ret); + + return(ret); +} + +// check whether a polynomial is 0 +int polynomial_is_zero(Polynomial poly){ + int i; + for(i=0;i<poly.length;i++){ + if(number_is_zero(poly.nums[i])==0){ + return(0); + } + } + return(1); +} + +// compute V^p +int polynomial_power(Polynomial input_polynomial, Polynomial* output, int power, Fields_Table fields){ + int i,j; + int* current_term; + // out buffers: since we first check whether the monomial vanishes before adding it to the output, + // it is more efficient to concatenate the monomials in a separate variable, and then add it to the ouput + // instead of incrementally adding terms to the output + Int_Array out_monomial; + int monomial_size; + Int_Array out_factor; + int factor_size; + Number out_num; + + init_Polynomial(output, POLY_SIZE); + + // trivial case + if(power==0){ + init_Int_Array(&out_monomial, 1); + init_Int_Array(&out_factor, 1); + polynomial_append_noinit(out_monomial, out_factor, number_one(), output); + return(0); + } + + // initialize current term + current_term=calloc(power,sizeof(int)); + for(i=0;i<power;i++){ + current_term[i]=0; + } + + // loop over terms; the loop stops when all the pointers in the list + // mentioned above are at the end of the input polynomial + while(current_term[0]<input_polynomial.length){ + // compute the amount of memory to allocate + for(i=0,monomial_size=0,factor_size=0;i<power;i++){ + monomial_size+=input_polynomial.monomials[current_term[i]].length; + factor_size+=input_polynomial.factors[current_term[i]].length; + } + // allocate + init_Int_Array(&out_monomial, monomial_size); + init_Int_Array(&out_factor, factor_size); + out_num=number_one(); + + // concatenate monomial and factor + for(i=0;i<power;i++){ + int_array_concat(input_polynomial.monomials[current_term[i]],&out_monomial); + int_array_concat(input_polynomial.factors[current_term[i]],&out_factor); + number_prod_chain(input_polynomial.nums[current_term[i]], &out_num); + } + + // check whether the monomial vanishes + if(check_monomial(out_monomial, fields)==1){ + // multinomial combinatorial factor n!/(m1!m2!m3!...) + number_Qprod_chain(quot(factorial(power),multinomial(power,current_term)),&out_num); + // write monomials and factors + polynomial_append_noinit(out_monomial, out_factor, out_num, output); + } + else{ + free_Int_Array(out_monomial); + free_Int_Array(out_factor); + free_Number(out_num); + } + + // move to the next term of V^p by advancing the last pointer if possible, + // the next to last if not, and so forth, until all the pointers have + // reached the end of the input polynomial + (current_term[power-1])++; + // loop until the last pointer is in an adequate place or the first + // pointer has reached the end + for(i=power-2;current_term[0]<input_polynomial.length && current_term[power-1]>=input_polynomial.length;i--){ + // try to advance pointer + (current_term[i])++; + // if the step fails, keep on looping, if not, set all the following + // pointers to the latest position + if(current_term[i]<input_polynomial.length){ + for(j=i+1;j<power;j++){ + current_term[j]=current_term[i]; + } + } + } + } + + polynomial_simplify(output, fields); + + // free memory + free(current_term); + return(0); +} + + +// compute V*W +int polynomial_prod(Polynomial input1, Polynomial input2, Polynomial* output, Fields_Table fields){ + polynomial_cpy(input2,output); + polynomial_prod_chain(input1,output, fields); + return(0); +} +// chain +int polynomial_prod_chain_nosimplify(Polynomial input, Polynomial* inout, Fields_Table fields){ + // position in V and W + int pos1, pos2; + Int_Array out_monomial; + Int_Array out_factor; + Number out_num; + // save length of inout (which changes during the loop + int inout_length=(*inout).length; + // first position in input which can multiply a term of inout without vanishing + int firstpos; + + // loop over terms + for(pos1=0;pos1<inout_length;pos1++){ + // multiply first term of input to inout[pos1] + // search for the first possible product + firstpos=-1; + for(pos2=0; pos2<input.length; pos2++){ + if(check_monomial_willnot_vanish((*inout).monomials[pos1],input.monomials[pos2],fields)==1){ + firstpos=pos2; + pos2++; + break; + } + } + + // if there was no possible product + if(firstpos==-1){ + number_Qprod_chain(quot(0,1),(*inout).nums+pos1); + } + else{ + // add other terms at the end of inout + for(;pos2<input.length;pos2++){ + // check whether the term will vanish + if(check_monomial_willnot_vanish((*inout).monomials[pos1],input.monomials[pos2],fields)==1){ + // allocate + init_Int_Array(&out_monomial, (*inout).monomials[pos1].length+input.monomials[pos2].length); + init_Int_Array(&out_factor, (*inout).factors[pos1].length+input.factors[pos2].length); + + // concatenate monomial and factor + int_array_concat((*inout).monomials[pos1],&out_monomial); + int_array_concat(input.monomials[pos2],&out_monomial); + int_array_concat((*inout).factors[pos1],&out_factor); + int_array_concat(input.factors[pos2],&out_factor); + number_prod((*inout).nums[pos1],input.nums[pos2],&out_num); + + // write monomials and factors + polynomial_append_noinit(out_monomial, out_factor, out_num, inout); + } + } + // first term + int_array_concat(input.monomials[firstpos],(*inout).monomials+pos1); + int_array_concat(input.factors[firstpos],(*inout).factors+pos1); + number_prod_chain(input.nums[firstpos],(*inout).nums+pos1); + } + } + + return(0); +} +// simplify +int polynomial_prod_chain(Polynomial input, Polynomial* inout, Fields_Table fields){ + polynomial_prod_chain_nosimplify(input, inout, fields); + polynomial_simplify(inout, fields); + return(0); +} + +// exp(V) +int polynomial_exponential(Polynomial input_polynomial, Polynomial* output, Fields_Table fields){ + // a buffer for the result of a given power + Polynomial tmp_poly; + // a buffer for the previous power + Polynomial previous_power; + // power + int power=1; + Int_Array out_monomial; + Int_Array out_factor; + + // allocate memory + init_Polynomial(output,POLY_SIZE); + + // 1 + init_Int_Array(&out_monomial, 1); + init_Int_Array(&out_factor, 1); + polynomial_append_noinit(out_monomial, out_factor, number_one(), output); + + while(1){ + if(power>1){ + if(power>33){ + fprintf(stderr,"error: trying to take a power of a polynomial that is too high (>33)\n"); + exit(-1); + } + // next power + polynomial_prod(input_polynomial,previous_power,&tmp_poly, fields); + + // free + free_Polynomial(previous_power); + } + else{ + polynomial_cpy(input_polynomial,&tmp_poly); + } + + // if the power is high enough that V^p=0, then stop + if(tmp_poly.length==0){ + free_Polynomial(tmp_poly); + break; + } + + // copy for next power + polynomial_cpy(tmp_poly,&previous_power); + + // 1/p! + polynomial_multiply_Qscalar(tmp_poly,quot(1,factorial(power))); + // append tmp to the output + polynomial_concat_noinit(tmp_poly,output); + + // increase power + power++; + } + + return(0); +} + + +// log(1+W) +int polynomial_logarithm(Polynomial input_polynomial,Polynomial* output, Fields_Table fields){ + // a buffer for the result of a given power + Polynomial tmp_poly; + // a buffer for the previous power + Polynomial previous_power; + // power + int power=1; + + // allocate memory + init_Polynomial(output,POLY_SIZE); + + while(1){ + if(power>1){ + // next power + polynomial_prod(input_polynomial,previous_power,&tmp_poly, fields); + + // free + free_Polynomial(previous_power); + } + else{ + polynomial_cpy(input_polynomial,&tmp_poly); + } + + // if the power is high enough that V^p=0, then stop + if(tmp_poly.length==0){ + free_Polynomial(tmp_poly); + break; + } + + // copy for next power + polynomial_cpy(tmp_poly,&previous_power); + + // (-1)^{p-1}/p + polynomial_multiply_Qscalar(tmp_poly,quot(ipower(-1,power-1),power)); + // append tmp to the output + polynomial_concat_noinit(tmp_poly,output); + + // increase power + power++; + } + + return(0); +} + +// check whether a monomial vanishes +int check_monomial(Int_Array monomial, Fields_Table fields){ + int i,j; + for(i=0;i<monomial.length;i++){ + // check for repetitions + if(is_fermion(monomial.values[i], fields)==1){ + for(j=i+1;j<monomial.length;j++){ + if(monomial.values[j]==monomial.values[i]){ + return(0); + } + } + } + } + return(1); +} +// check whether the product of two monomials will vanish +int check_monomial_willnot_vanish(Int_Array monomial1, Int_Array monomial2, Fields_Table fields){ + int i,j; + for(i=0;i<monomial1.length;i++){ + // check for repetitions + if(is_fermion(monomial1.values[i], fields)==1){ + for(j=0;j<monomial2.length;j++){ + if(monomial2.values[j]==monomial1.values[i]){ + return(0); + } + } + } + } + return(1); +} + +// check whether one can add a term to a monomial without creating repetitions +int check_monomial_addterm(Int_Array monomial, Int_Array term, Fields_Table fields){ + int i,j; + for(i=0;i<term.length;i++){ + // check for repetitions + if(is_fermion(term.values[i], fields)==1){ + for(j=0;j<monomial.length;j++){ + if(monomial.values[j]==term.values[i]){ + return(0); + } + } + } + } + return(1); +} + +// check whether a monomial vanishes or has unmatched +/- fields +int check_monomial_match(Int_Array monomial, Fields_Table fields){ + int i,j; + int match=0; + for(i=0;i<monomial.length;i++){ + // count match + if(field_type(monomial.values[i], fields)==FIELD_INTERNAL){ + if(monomial.values[i]>0){ + match++; + } + else if(monomial.values[i]<0){ + match--; + } + } + + // check for repetitions + if(is_fermion(monomial.values[i], fields)==1){ + for(j=i+1;j<monomial.length;j++){ + if(monomial.values[j]==monomial.values[i]){ + return(0); + } + } + } + } + if(match==0){ + return(1); + } + else{ + // different return codes depending on why the monomial was rejected + return(-1); + } +} + +// remove terms with more plus internal fields than minus ones +int remove_unmatched_plusminus(Polynomial* polynomial, Fields_Table fields){ + int i,j; + int match_internals; + int type; + Polynomial output; + + init_Polynomial(&output, (*polynomial).length); + + for(i=0;i<(*polynomial).length;i++){ + match_internals=0; + for(j=0;j<(*polynomial).monomials[i].length;j++){ + // check for unmatched internal fields + type=field_type((*polynomial).monomials[i].values[j],fields); + if(type==FIELD_INTERNAL){ + if((*polynomial).monomials[i].values[j]>0){ + match_internals++; + } + else if((*polynomial).monomials[i].values[j]<0){ + match_internals--; + } + } + // don't remove a term containing symbols + else if(type==FIELD_SYMBOL){ + match_internals=0; + break; + } + } + if(match_internals==0){ + polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], &output); + } + } + + free_Polynomial(*polynomial); + *polynomial=output; + return(0); +} + + +// denominator of a multinomal factor: m1!m2!... +// requires terms to be sorted +int multinomial(int power,int* terms){ + int multiple=1; + int ret=1; + int i; + // the number of numbers to be multiplied in the multinomial is + // equal to power-1 (the first is 1) + for(i=1;i<power;i++){ + // if there is a degeneracy, then increment the multiple by + // which the multinomial is multiplied + if(terms[i-1]==terms[i]){ + multiple++; + } + // if there is no degeneracy, reset it to 1 + else{ + multiple=1; + } + // multiply the result by the multiple + ret*=multiple; + } + return(ret); +} + + +// simplify a Polynomial +// the fields table is there to compute the sign coming from re-arranging monomials +int polynomial_simplify(Polynomial* polynomial, Fields_Table fields){ + int i; + int monomial_cmp; + int factor_cmp; + int sign; + Polynomial output; + init_Polynomial(&output,(*polynomial).length); + // the combination of numerical factors + Number new_num; + init_Number(&new_num,NUMBER_SIZE); + + // sort monomials and factors + for(i=0;i<(*polynomial).length;i++){ + sign=1; + monomial_sort((*polynomial).monomials[i],0,(*polynomial).monomials[i].length-1,fields,&sign); + number_Qprod_chain(quot(sign,1),(*polynomial).nums+i); + int_array_sort((*polynomial).factors[i],0,(*polynomial).factors[i].length-1); + } + + // resolve the identities specified in the fields table + resolve_ids(polynomial, fields); + + // in order to perform a simplification, the list of terms must be + // sorted (so that terms that are proportional are next to each other) + polynomial_sort(polynomial,0,(*polynomial).length-1); + + for(i=0;i<(*polynomial).length;i++){ + // if the term actually exists + if(number_is_zero((*polynomial).nums[i])!=1){ + // combine numerical factors + number_add_chain((*polynomial).nums[i], &new_num); + } + // if the numerator is 0, the previous terms that may have the same factors should still be added, hence the 'if' ends here + + // if either the monomial or the factor is different from the next then add term + if(i<(*polynomial).length-1){ + monomial_cmp=int_array_cmp((*polynomial).monomials[i],(*polynomial).monomials[i+1]); + factor_cmp=int_array_cmp((*polynomial).factors[i],(*polynomial).factors[i+1]); + } + if(i>=(*polynomial).length-1 || monomial_cmp!=0 || factor_cmp!=0 ){ + // check that the polynomial is not trivial + if(number_is_zero(new_num)!=1){ + polynomial_append((*polynomial).monomials[i],(*polynomial).factors[i],new_num,&output); + } + + // reset new numerical factor + free_Number(new_num); + init_Number(&new_num,NUMBER_SIZE); + } + } + + free_Number(new_num); + free_Polynomial(*polynomial); + *polynomial=output; + return(0); +} + +// sort a polynomial +// requires the monomials and factors to be sorted +int polynomial_sort(Polynomial* polynomial, int begin, int end){ + int i; + int index; + int monomial_cmp; + int factor_cmp; + + // the pivot: middle of the array + int pivot=(begin+end)/2; + // if the array is non trivial + if(begin<end){ + // send pivot to the end + exchange_polynomial_terms(pivot,end,polynomial); + // loop over the others + for(i=begin, index=begin;i<end;i++){ + // compare with pivot + monomial_cmp=int_array_cmp((*polynomial).monomials[i],(*polynomial).monomials[end]); + factor_cmp=int_array_cmp((*polynomial).factors[i],(*polynomial).factors[end]); + if(monomial_cmp<0 || (monomial_cmp==0 && factor_cmp<0)){ + // if smaller, exchange with reference index + exchange_polynomial_terms(i,index,polynomial); + // move reference index + index++; + } + } + // put pivot (which we had sent to the end) in the right place + exchange_polynomial_terms(index,end,polynomial); + // recurse + polynomial_sort(polynomial, begin, index-1); + polynomial_sort(polynomial, index+1, end); + } + return(0); +} + +// exchange two terms (for the sorting algorithm) +int exchange_polynomial_terms(int i, int j, Polynomial* polynomial){ + Int_Array ptmp; + Number tmp; + + ptmp=(*polynomial).monomials[j]; + (*polynomial).monomials[j]=(*polynomial).monomials[i]; + (*polynomial).monomials[i]=ptmp; + + ptmp=(*polynomial).factors[j]; + (*polynomial).factors[j]=(*polynomial).factors[i]; + (*polynomial).factors[i]=ptmp; + + tmp=(*polynomial).nums[j]; + (*polynomial).nums[j]=(*polynomial).nums[i]; + (*polynomial).nums[i]=tmp; + + return(0); +} + +// sort a monomial (with sign coming from exchanging two Fermions) +int monomial_sort(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign){ + int i; + 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 + exchange_monomial_terms(monomial, pivot, end, fields, sign); + + // loop over the others + for(i=begin, index=begin;i<end;i++){ + // compare with pivot + if(compare_monomial_terms(monomial, i, end, fields)<0){ + // if smaller, exchange with reference index + exchange_monomial_terms(monomial, index, i, fields, sign); + // move reference index + index++; + } + } + // put pivot (which we had sent to the end) in the right place + exchange_monomial_terms(monomial, index, end, fields, sign); + + // recurse + monomial_sort(monomial, begin, index-1, fields, sign); + monomial_sort(monomial, index+1, end, fields, sign); + } + return(0); +} + +// order fields: parameter, external, internal +int compare_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields){ + int type1=field_type(monomial.values[pos1], fields); + int type2=field_type(monomial.values[pos2],fields); + + if(type1 < type2){ + return(-1); + } + else if(type1 > type2){ + return(1); + } + + if(monomial.values[pos1] < monomial.values[pos2]){ + return(-1); + } + else if (monomial.values[pos1] > monomial.values[pos2]){ + return(1); + } + + return(0); +} + +// exchange two fields (with sign) +int exchange_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields, int* sign){ + int tmp=monomial.values[pos1]; + int f1,f2; + int i; + + monomial.values[pos1]=monomial.values[pos2]; + monomial.values[pos2]=tmp; + + // sign change + // only if the exchange is not trivial + if(pos1!=pos2){ + f1=is_fermion(monomial.values[pos1],fields); + f2=is_fermion(monomial.values[pos2],fields); + // if both Fermions then sign + if(f1==1 && f2==1){ + *sign*=-1; + } + // if only one of them is a Fermion, then count the number of Fermions between them + else if(f1==1 || f2==1){ + for(i=min(pos1,pos2)+1;i<max(pos1,pos2);i++){ + if(is_fermion(monomial.values[i],fields)==1){ + *sign*=-1; + } + } + } + } + return(0); +} + + +// sort a monomial by putting each group together +int monomial_sort_groups(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign){ + int i; + 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 + exchange_monomial_terms(monomial, pivot, end, fields, sign); + + // loop over the others + for(i=begin, index=begin;i<end;i++){ + // compare with pivot + if(compare_monomial_terms_groups(monomial, i, end, fields, groups)<0){ + // if smaller, exchange with reference index + exchange_monomial_terms(monomial, index, i, fields, sign); + // move reference index + index++; + } + } + // put pivot (which we had sent to the end) in the right place + exchange_monomial_terms(monomial, index, end, fields, sign); + + // recurse + monomial_sort(monomial, begin, index-1, fields, sign); + monomial_sort(monomial, index+1, end, fields, sign); + } + return(0); +} + +// order fields: group, then parameter, external, internal +int compare_monomial_terms_groups(Int_Array monomial, int pos1, int pos2, Fields_Table fields, Groups groups){ + int group1=find_group(monomial.values[pos1], groups); + int group2=find_group(monomial.values[pos2], groups); + + if(group1 < group2){ + return(-1); + } + else if(group1 > group2){ + return(1); + } + + int type1=field_type(monomial.values[pos1], fields); + int type2=field_type(monomial.values[pos2],fields); + + if(type1 < type2){ + return(-1); + } + else if(type1 > type2){ + return(1); + } + + if(monomial.values[pos1] < monomial.values[pos2]){ + return(-1); + } + else if (monomial.values[pos1] > monomial.values[pos2]){ + return(1); + } + + return(0); +} + + +// convert and idtable to a polynomial +int idtable_to_polynomial(Id_Table idtable, Polynomial* polynomial){ + int i,j; + int start=0; + + // allocate memory + init_Polynomial(polynomial,POLY_SIZE); + + for(i=0;i<idtable.length;i++){ + polynomial_concat(idtable.polynomials[i],polynomial); + // set factors + for(j=start;j<(*polynomial).length;j++){ + int_array_append(idtable.indices[i],(*polynomial).factors+j); + } + // shift start point + start+=idtable.polynomials[i].length; + } + + return(0); +} + + +// replace the factors in a polynomial as prescribed by an equation in the form of a Grouped_Polynomial +int replace_factors(Grouped_Polynomial equations, Polynomial* polynomial){ + int i,j; + int index; + Polynomial output; + Coefficient coef; + + init_Polynomial(&output, POLY_SIZE); + + // loop over monomials + for(i=0;i<(*polynomial).length;i++){ + if((*polynomial).factors[i].length>0){ + // initialize coef to store the product of factors + init_Coefficient(&coef, POLY_SIZE); + + // first term + index=intlist_find_err(equations.indices, equations.length, (*polynomial).factors[i].values[0]); + if(index>=0){ + coefficient_cpy_noinit(equations.coefs[index],&coef); + } + + // other terms + for(j=1;j<(*polynomial).factors[i].length;j++){ + index=intlist_find_err(equations.indices, equations.length, (*polynomial).factors[i].values[j]); + if(index>=0){ + coefficient_prod_chain(equations.coefs[index],&coef); + } + } + + // new polynomial terms + for(j=0;j<coef.length;j++){ + // add to output + polynomial_append((*polynomial).monomials[i], coef.factors[j], number_prod_ret((*polynomial).nums[i],coef.nums[j]), &output); + } + + // free memory + free_Coefficient(coef); + } + // if no factors + else{ + polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], &output); + } + } + + // replace output + free_Polynomial(*polynomial); + *polynomial=output; + + return(0); +} + + +// print a polynomial to a string +int polynomial_sprint(Polynomial polynomial, Char_Array* output){ + int i,j; + + // for each monomial + for(i=0;i<polynomial.length;i++){ + if(i==0){ + char_array_snprintf(output, " "); + } + else{ + char_array_snprintf(output, "+ ",i); + } + + // print num + char_array_append('(',output); + number_sprint(polynomial.nums[i],output); + char_array_append(')',output); + + // print factors + for(j=0;j<polynomial.factors[i].length;j++){ + char_array_snprintf(output,"[l%d]",polynomial.factors[i].values[j]); + } + + // print monomials + for(j=0;j<polynomial.monomials[i].length;j++){ + char_array_snprintf(output,"[f%d]",polynomial.monomials[i].values[j]); + } + + char_array_append('\n',output); + } + return(0); +} +// print +int polynomial_print(Polynomial polynomial){ + Char_Array buffer; + init_Char_Array(&buffer, STR_SIZE); + polynomial_sprint(polynomial, &buffer); + printf("%s",buffer.str); + free_Char_Array(buffer); + return(0); +} + +// read a polynomial from a Char_Array +#define PP_NULL_MODE 0 +#define PP_BRACKET_MODE 1 +#define PP_MONOMIAL_MODE 2 +#define PP_FACTOR_MODE 3 +#define PP_NUMBER_MODE 4 +int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){ + // buffer + char* buffer=calloc(str_polynomial.length+1,sizeof(char)); + char* buffer_ptr=buffer; + Int_Array monomial; + Int_Array factor; + Number num, tmp1_num; + int mode; + int comment=0; + int i,j; + int parenthesis_count=0; + + // allocate memory + init_Polynomial(output,POLY_SIZE); + + // init + init_Int_Array(&monomial, MONOMIAL_SIZE); + init_Int_Array(&factor, MONOMIAL_SIZE); + num=number_one(); + + // if the string contains no '[', then read a number + for(j=0;j<str_polynomial.length;j++){ + // ignore comments + if(comment==1){ + if(str_polynomial.str[j]=='\n'){ + comment=0; + } + } + else{ + if(str_polynomial.str[j]=='['){ + break; + } + if(str_polynomial.str[j]=='#'){ + comment=1; + } + } + } + // no '[': read a number + if(j==str_polynomial.length){ + free_Number(num); + char_array_to_Number(str_polynomial,&num); + polynomial_append_noinit(monomial, factor, num, output); + free(buffer); + return(0); + } + + // reset comment flag + comment=0; + + *buffer_ptr='\0'; + // loop over the input polynomial + // start in null mode + mode=PP_NULL_MODE; + for(j=0;j<str_polynomial.length;j++){ + if(comment==1){ + if(str_polynomial.str[j]=='\n'){ + comment=0; + } + } + else{ + switch(str_polynomial.str[j]){ + // new monomial + case '+': + if(mode==PP_NULL_MODE){ + polynomial_append_noinit(monomial, factor, num, output); + // reset monomial, factor, num + init_Int_Array(&monomial, MONOMIAL_SIZE); + init_Int_Array(&factor, MONOMIAL_SIZE); + num=number_one(); + } + break; + + // enter monomial or factor mode + case '[': + if(mode==PP_NULL_MODE){ + mode=PP_BRACKET_MODE; + } + break; + // factor mode + case 'l': + if(mode==PP_BRACKET_MODE){ + mode=PP_FACTOR_MODE; + buffer_ptr=buffer; + *buffer_ptr='\0'; + } + break; + // monomial mode + case 'f': + if(mode==PP_BRACKET_MODE){ + mode=PP_MONOMIAL_MODE; + buffer_ptr=buffer; + *buffer_ptr='\0'; + } + break; + // read monomial or factor + case ']': + sscanf(buffer,"%d",&i); + if(mode==PP_FACTOR_MODE){ + int_array_append(i,&factor); + } + else if(mode==PP_MONOMIAL_MODE){ + int_array_append(i,&monomial); + } + // switch back to null mode + mode=PP_NULL_MODE; + break; + + // numerical pre-factor + case '(': + if(mode==PP_NULL_MODE){ + mode=PP_NUMBER_MODE; + parenthesis_count=0; + buffer_ptr=buffer; + *buffer_ptr='\0'; + } + else if(mode==PP_NUMBER_MODE){ + // match parentheses + parenthesis_count++; + buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]); + } + break; + case ')': + if(mode==PP_NUMBER_MODE){ + if(parenthesis_count==0){ + // write num + str_to_Number(buffer,&tmp1_num); + number_prod_chain(tmp1_num,&num); + free_Number(tmp1_num); + // back to null mode + mode=PP_NULL_MODE; + } + else{ + parenthesis_count--; + buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]); + } + } + break; + + // characters to ignore + case ' ':break; + case '&':break; + case '\n':break; + + // comments + case '#': + comment=1; + break; + + default: + if(mode!=PP_NULL_MODE){ + // write to buffer + buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]); + } + break; + } + } + } + + // last term + polynomial_append_noinit(monomial, factor, num, output); + + free(buffer); + return(0); +} + +// with str input +int str_to_Polynomial(char* str_polynomial,Polynomial* output){ + Char_Array buffer; + str_to_char_array(str_polynomial, &buffer); + Char_Array_to_Polynomial(buffer, output); + free_Char_Array(buffer); + return(0); +} + + +// -------------------- Polynomial_Matrix --------------------- + +// init +int init_Polynomial_Matrix(Polynomial_Matrix* matrix, int length){ + int i,j; + (*matrix).matrix=calloc(length,sizeof(Polynomial*)); + (*matrix).indices=calloc(length,sizeof(int)); + for(i=0;i<length;i++){ + (*matrix).matrix[i]=calloc(length,sizeof(Polynomial)); + for(j=0;j<length;j++){ + (*matrix).matrix[i][j]=polynomial_zero(); + } + } + (*matrix).length=length; + return(0); +} +int free_Polynomial_Matrix(Polynomial_Matrix matrix){ + int i,j; + for(i=0;i<matrix.length;i++){ + for(j=0;j<matrix.length;j++){ + free_Polynomial(matrix.matrix[i][j]); + } + free(matrix.matrix[i]); + } + free(matrix.matrix); + free(matrix.indices); + return(0); +} |