diff options
| author | Ian Jauslin <ian@jauslin.org> | 2022-06-14 09:26:07 +0200 |
|---|---|---|
| committer | Ian Jauslin <ian@jauslin.org> | 2022-06-14 09:26:07 +0200 |
| commit | 167980ea437881ec56186332370afcc169f2e4dd (patch) | |
| tree | 723e0c8b02bdbf77702582c0e6a007c7af804b62 /src | |
| parent | 469bdc80712dbf9c12562059dc4594620b59a076 (diff) | |
Update to v1.5.v1.5
The update to version 1.5 is rather substantial, and introduces some minor
backward-incompatibilities:
* The header "#!symbols" has been replaced by "#!virtual_fields"
* Multiplying polynomials using the '*' symbol is no longer supported (or,
rather, the symbolic capabilities of meankondo were enhanced, and the
syntax has been changed).
* 'meantools exp' has been removed (its functionality is now handled by
other means)
* 'meantoolds derive' has been replaced by 'meantools differentiate'
* The symbolic capabilities were enhanced: polynomials can now be
multiplied, added, exponentiated, and their logarithms can be taken
directly in the configuration file.
* The flow equation can now be processed after being computed using the
various "#!postprocess_*" entries.
* Deprecated kondo_preprocess.
* Compute the mean using an LU decomposition if possible.
* More detailed checks for syntax errors in configuration file.
* Check that different '#!group' entries are indeed uncorrelated.
* New flags in meankondo: '-p' and '-A'.
* New tool: meantools expand.
* Improve conversion to LaTeX using meantools-convert
* Assign terms randomly to different threads.
* Multiple bug fixes
Diffstat (limited to 'src')
60 files changed, 3070 insertions, 682 deletions
diff --git a/src/array.c b/src/array.c index 11d14f7..7fb6bd8 100644 --- a/src/array.c +++ b/src/array.c @@ -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. @@ -72,6 +72,14 @@ int int_array_append(int val, Int_Array* output){ return(0); } +// add a value only if it is not already present +int int_array_append_unique(int val, Int_Array* output){ + if(int_array_find(val,*output)<0){ + int_array_append(val,output); + } + return(0); +} + // concatenate int int_array_concat(Int_Array input, Int_Array* output){ int i; @@ -87,6 +95,15 @@ int int_array_concat(Int_Array input, Int_Array* output){ return(0); } +// concat but only add values that are not already present in the array +int int_array_concat_unique(Int_Array input, Int_Array* output){ + int i; + for(i=0;i<input.length;i++){ + int_array_append_unique(input.values[i],output); + } + return(0); +} + // find (does not assume the array is sorted) int int_array_find(int val, Int_Array array){ int i; @@ -199,7 +216,12 @@ int int_array_print(Int_Array array){ for(i=0;i<array.length-1;i++){ printf("%d,",array.values[i]); } - printf("%d)",array.values[array.length-1]); + if(array.length>0){ + printf("%d)",array.values[array.length-1]); + } + else{ + printf(")"); + } return(0); } @@ -334,6 +356,19 @@ int char_array_concat(Char_Array input, Char_Array* output){ return(0); } +// substring +int char_array_substring(Char_Array str, int begin, int end, Char_Array* substr){ + int i; + if(begin>end || begin<0 || end>=str.length){ + fprintf(stderr,"error: cannot extract a substring [%d,%d] from a string of length %d\n", begin, end, str.length); + exit(-1); + } + init_Char_Array(substr,end-begin); + for(i=begin;i<=end;i++){ + char_array_append(str.str[i],substr); + } + return(0); +} // convert to char* @@ -343,7 +378,7 @@ int char_array_to_str(Char_Array input, char** output){ for(i=0;i<input.length;i++){ (*output)[i]=input.str[i]; } - if((*output)[input.length-1]!='\0'){ + if(input.length==0 || (*output)[input.length-1]!='\0'){ (*output)[input.length]='\0'; } return(0); @@ -371,6 +406,34 @@ int str_to_char_array(char* str, Char_Array* output){ return(0); } +// compare char_array's +int char_array_cmp(Char_Array char_array1, Char_Array char_array2){ + int j; + if(char_array1.length!=char_array2.length){ + return(0); + } + for(j=0;j<char_array1.length && j<char_array2.length;j++){ + if(char_array1.str[j]!=char_array2.str[j]){ + return(0); + } + } + return(1); +} + +// compare a char_array and a char* +int char_array_cmp_str(Char_Array char_array, char* str){ + int j; + for(j=0;j<char_array.length && str[j]!='\0';j++){ + if(char_array.str[j]!=str[j]){ + return(0); + } + } + if(j==char_array.length && str[j]=='\0'){ + return(1); + } + return(0); +} + // format strings int char_array_snprintf(Char_Array* output, char* fmt, ...){ diff --git a/src/array.h b/src/array.h index fb74e67..480b589 100644 --- a/src/array.h +++ b/src/array.h @@ -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. @@ -34,8 +34,12 @@ int int_array_resize(Int_Array* array, int newsize); // add a value int int_array_append(int val, Int_Array* output); +// add a value only if it is not already present +int int_array_append_unique(int val, Int_Array* output); // concatenate int int_array_concat(Int_Array input, Int_Array* output); +// concat but only add values that are not already present in the array +int int_array_concat_unique(Int_Array input, Int_Array* output); // find (does not assume the array is sorted) int int_array_find(int val, Int_Array array); @@ -75,6 +79,9 @@ int char_array_append_str(char* str, Char_Array* output); // concatenate int char_array_concat(Char_Array input, Char_Array* output); +// substring +int char_array_substring(Char_Array str, int begin, int end, Char_Array* substr); + // convert to char* int char_array_to_str(Char_Array input, char** output); // noinit (changes the size of input if needed) @@ -82,6 +89,11 @@ char* char_array_to_str_noinit(Char_Array* input); // convert from char* int str_to_char_array(char* str, Char_Array* output); +// compare char_array's +int char_array_cmp(Char_Array char_array1, Char_Array char_array2); +// compare a char_array and a char* +int char_array_cmp_str(Char_Array char_array, char* str); + // format strings int char_array_snprintf(Char_Array* output, char* fmt, ...); diff --git a/src/cli_parser.c b/src/cli_parser.c index bdca700..6012e8a 100644 --- a/src/cli_parser.c +++ b/src/cli_parser.c @@ -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. @@ -121,3 +121,23 @@ int find_str_arg(char* title, Str_Array str_args){ } } +// find a string argument with the specified title +// namespace support +int find_str_arg_ns(char* title, Char_Array namespace, Str_Array str_args){ + Char_Array buffer; + int ret; + + // append namespace to title + char_array_cpy(namespace,&buffer); + char_array_append(':',&buffer); + char_array_append_str(title,&buffer); + + // check whether the namespace entry exists + ret=find_str_arg(char_array_to_str_noinit(&buffer), str_args); + free_Char_Array(buffer); + // if not, use global entry + if(ret==-1){ + return(find_str_arg(title, str_args)); + } + return(ret); +} diff --git a/src/cli_parser.h b/src/cli_parser.h index b050e37..b66608c 100644 --- a/src/cli_parser.h +++ b/src/cli_parser.h @@ -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,8 @@ int read_config_file(Str_Array* str_args, const char* file, int read_from_stdin) int get_str_arg_title(Char_Array str_arg, Char_Array* out); // find a string argument with the specified title int find_str_arg(char* title, Str_Array str_args); +// namespace support +int find_str_arg_ns(char* title, Char_Array namespace, Str_Array str_args); #endif diff --git a/src/coefficient.c b/src/coefficient.c index d4cb6eb..1efe8fd 100644 --- a/src/coefficient.c +++ b/src/coefficient.c @@ -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. @@ -202,6 +202,235 @@ int coefficient_simplify(Coefficient* coefficient){ return(0); } +// put all terms under a common denominator and simplify the resulting fraction +int coefficient_simplify_rational(Coefficient constant, Coefficient* coefficient){ + int ret; + Coefficient remainder; + Coefficient quotient; + Coefficient quotient_prev; + Coefficient out; + int power; + int max_power; + + // common denominator + coefficient_common_denominator(constant, coefficient); + + // init + init_Coefficient(&out, COEF_SIZE); + + // simplify, one power at a time + // largest power (larger powers are at the end) + max_power=(*coefficient).denoms[(*coefficient).length-1].power; + + quotient_prev=*coefficient; + for(power=max_power;power>=1;power--){ + ret=coefficient_simplify_fraction(constant, quotient_prev, &remainder, "ient); + + // if fail to simplify, stop + if(ret<0){ + if(power<max_power){ + coefficient_concat_noinit(quotient_prev, &out); + } + else{ + coefficient_concat(quotient_prev, &out); + } + break; + } + + // add to output + coefficient_concat_noinit(remainder, &out); + } + // if the factorization always succeeded + if(max_power>=1 && power==0){ + coefficient_concat_noinit(quotient, &out); + } + + coefficient_simplify(&out); + + // set coefficient to out + free_Coefficient(*coefficient); + *coefficient=out; + + return 0; +} + +// put all terms under a common denominator +// only supports coefficients with only one constant +int coefficient_common_denominator(Coefficient constant, Coefficient* coefficient){ + int max_power; + int i,j; + Coefficient tmp; + Coefficient out; + Coefficient* C_n; + + init_Coefficient(&out, COEF_SIZE); + + // largest power (larger powers are at the end) + max_power=(*coefficient).denoms[(*coefficient).length-1].power; + + // store powers of the constant + C_n=calloc(sizeof(Coefficient), max_power-1); + for(i=0;i<max_power-1;i++){ + // start from previous product + if(i==0){ + coefficient_cpy(constant, C_n+i); + } + else{ + coefficient_cpy(C_n[i-1], C_n+i); + } + // multiply by constant + coefficient_prod_chain(constant, C_n+i); + } + + // multiply each term + for (i=0;i<(*coefficient).length;i++){ + init_Coefficient(&tmp, COEF_SIZE); + // start with numerator + coefficient_append_noinit((*coefficient).factors[i], (*coefficient).nums[i], (*coefficient).denoms[i], &tmp); + // multiply + if((*coefficient).denoms[i].power<max_power){ + if((*coefficient).denoms[i].power==max_power-1){ + coefficient_prod_chain(constant, &tmp); + } + else{ + coefficient_prod_chain(C_n[max_power-(*coefficient).denoms[i].power-2], &tmp); + } + } + + // set denom + for(j=0;j<tmp.length;j++){ + tmp.denoms[j].power=max_power; + } + + // add to out + coefficient_concat_noinit(tmp, &out); + } + + // free C_n + for(i=0;i<max_power-1;i++){ + free_Coefficient(C_n[i]); + } + free(C_n); + + // free coefficient vectors + free((*coefficient).factors); + free((*coefficient).nums); + free((*coefficient).denoms); + + // set output + *coefficient=out; + + coefficient_simplify(coefficient); + + return(0); +} + +// simplify coefficient / constant +// returns both the remainder and the quotient +// assumes both coefficient and constant are ordered with the highest order terms last +int coefficient_simplify_fraction(Coefficient constant, Coefficient coefficient, Coefficient* remainder, Coefficient* out){ + Coefficient tmp; + int step_counter=0; + int max_order; + int i,j,k; + Int_Array rfactors; + + if(constant.length==0){ + // nothing to do + return 0; + } + + coefficient_cpy(coefficient, remainder); + init_Coefficient(out, COEF_SIZE); + + // continue until (*remainder) is of lower order than constant + while((*remainder).length>0 && (*remainder).factors[(*remainder).length-1].length>=constant.factors[constant.length-1].length){ + step_counter++; + + // interrupt if too long + if(step_counter>=coefficient.length*100){ + free_Coefficient(*remainder); + free_Coefficient(*out); + return -1; + } + + // try to find a term in the constant that divides the last term of the (*remainder) + rfactors=(*remainder).factors[(*remainder).length-1]; + + // highest order in constant + max_order=constant.factors[constant.length-1].length; + + // start from one of the highest order term and look for a common factor + for(i=constant.length-1; i>=0; i--){ + // fail: no highest order terms have been matched + if(constant.factors[i].length<max_order){ + free_Coefficient(*remainder); + free_Coefficient(*out); + return -2; + } + + // check whether the term can be a factor of the last term of the (*remainder) + if(int_array_is_subarray_ordered(constant.factors[i], rfactors)==1){ + // extract the factors that are not in constant + init_Coefficient(&tmp, constant.length); + + // init with one term + tmp.length=1; + init_Int_Array(tmp.factors,MONOMIAL_SIZE); + + for(j=0,k=0;j<rfactors.length;j++){ + // check that index is not in constant + if(k<constant.factors[i].length){ + if(rfactors.values[j]!=constant.factors[i].values[k]){ + int_array_append(rfactors.values[j],tmp.factors); + } + else{ + // move to next term in constant + k++; + } + } + } + + // numerical prefactor: term in the (*remainder) / term in the constant + number_quot((*remainder).nums[(*remainder).length-1], constant.nums[i], tmp.nums); + // denominator (dummy) + tmp.denoms[0]=(*remainder).denoms[(*remainder).length-1]; + + // add to out + coefficient_concat(tmp, out); + + // multiply by -1 + Q minus_1; + minus_1.numerator=-1; + minus_1.denominator=1; + number_Qprod_chain(minus_1, tmp.nums); + + // multiply by constant + coefficient_prod_chain(constant, &tmp); + + // add to remainder + coefficient_concat(tmp, remainder); + + // free memory + free_Coefficient(tmp); + + // simplify + coefficient_simplify(remainder); + + break; + } + } + } + + // success! + // decrease power of constant + for(i=0;i<(*out).length;i++){ + (*out).denoms[i].power=(*out).denoms[i].power-1; + } + + return(0); +} + // sort the terms in an equation (quicksort algorithm) int sort_coefficient(Coefficient* coefficient, int begin, int end){ int i; @@ -251,7 +480,7 @@ int exchange_coefficient_terms(int i, int j, Coefficient* coefficient){ return(0); } -// derive a coefficient with respect to an index +// differentiate a coefficient with respect to an index int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){ int i,j; // temp list of indices @@ -325,7 +554,7 @@ int coefficient_deriv(Coefficient input, int index, Coefficient* output){ } /* -// derive a coefficient with respect to an index (as a polynomial) (does not derive the 1/(1+C)^p ) +// differentiate a coefficient with respect to an index (as a polynomial) (does not differentiate the 1/(1+C)^p ) int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){ int i; // temp list of indices @@ -365,7 +594,7 @@ int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){ return(0); } -// derive a monomial with respect to an index +// differentiate a monomial with respect to an index int monomial_deriv(Int_Array factor, int index, Int_Array* out_factor, int* match_count){ int j; @@ -414,14 +643,14 @@ int coefficient_prod(Coefficient coef1, Coefficient coef2, Coefficient* output){ int_array_concat(coef2.factors[j],&factor); // don't throw an error if the power is 0 - if(coef2.denoms[i].power==0){ - coef2.denoms[i].index=coef1.denoms[i].index; + if(coef2.denoms[j].power==0){ + coef2.denoms[j].index=coef1.denoms[i].index; } else if(coef1.denoms[i].power==0){ - coef1.denoms[i].index=coef2.denoms[i].index; + coef1.denoms[i].index=coef2.denoms[j].index; } if(coef1.denoms[i].index!=coef2.denoms[j].index){ - fprintf(stderr,"error: cannot multiply flow equations with different constants\n"); + fprintf(stderr,"error: cannot multiply flow equations with different constants: got %d and %d\n", coef1.denoms[i].index, coef2.denoms[j].index); exit(-1); } denom=coef1.denoms[i]; diff --git a/src/coefficient.h b/src/coefficient.h index 43dd313..f8a813a 100644 --- a/src/coefficient.h +++ b/src/coefficient.h @@ -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. @@ -48,12 +48,20 @@ int coefficient_concat_noinit(Coefficient input, Coefficient* output); // simplify a Coefficient int coefficient_simplify(Coefficient* coefficient); + +// put all terms under a common denominator and simplify the resulting fraction +int coefficient_simplify_rational(Coefficient constant, Coefficient* coefficient); +// put all terms under a common denominator +int coefficient_common_denominator(Coefficient constant, Coefficient* coefficient); +// simplify coefficient / constant +int coefficient_simplify_fraction(Coefficient constant, Coefficient coefficient, Coefficient* remainder, Coefficient* out); + // sort the terms in an equation (quicksort algorithm) int sort_coefficient(Coefficient* coefficient, int begin, int end); // exchange two terms (for the sorting algorithm) int exchange_coefficient_terms(int i, int j, Coefficient* coefficient); -// derive a coefficient with respect to an index (as a polynomial) (does not derive the 1/(1+C)^p ) +// differentiate a coefficient with respect to an index (as a polynomial) (does not differentiate the 1/(1+C)^p ) int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output); int coefficient_deriv(Coefficient input, int index, Coefficient* output); diff --git a/src/definitions.cpp b/src/definitions.cpp index 30f021f..db09e57 100644 --- a/src/definitions.cpp +++ b/src/definitions.cpp @@ -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. @@ -33,10 +33,16 @@ limitations under the License. #define EQUATION_SIZE 20 // number of fields #define FIELDS_SIZE 50 +// number of variables +#define VARIABLES_SIZE 10 // number of elements in numbers #define NUMBER_SIZE 5 // number of elements in a group #define GROUP_SIZE 5 +// number of children per node in a symbol_tree +#define SYMBOL_TREE_SIZE 2 +// size of character string in symbol tree +#define SYMBOL_TREE_LABEL_SIZE 10 // display options @@ -55,6 +61,6 @@ limitations under the License. #define FIELD_PARAMETER 1 #define FIELD_EXTERNAL 2 #define FIELD_INTERNAL 3 -#define FIELD_SYMBOL 4 +#define FIELD_VIRTUAL 4 #endif diff --git a/src/determinant.c b/src/determinant.c new file mode 100644 index 0000000..906c75f --- /dev/null +++ b/src/determinant.c @@ -0,0 +1,93 @@ +#include "determinant.h" + +#include "number.h" +#include "rational.h" +#include "definitions.cpp" + +// determinant of a matrix +// replaces the matrix by its LU decomposition +int determinant_inplace(Number_Matrix M, Number* out){ + int i; + int sign_correction; + + LU_dcmp_inplace(M, &sign_correction); + + if(sign_correction==0){ + *out=number_zero(); + return(0); + } + + *out=number_one(); + if(sign_correction==-1){ + number_Qprod_chain(quot(-1,1), out); + } + + for(i=0;i<M.length;i++){ + number_prod_chain(M.matrix[i][i], out); + } + + return(0); +} + +// LU decomposition +// uses pivoting to avoid dividing by 0 +// the sign_correction should be multiplied to the determinant to obtain the right value +// if dividing by 0 is unavoidable, then the determinant is 0, and sign_correction is set to 0 +int LU_dcmp_inplace(Number_Matrix M, int* sign_correction){ + int i,j,k,pivot; + Number tmp; + + *sign_correction=1; + + for(j=0;j<M.length;j++){ + for(i=0;i<=j;i++){ + for(k=0;k<i;k++){ + // -M[i][k]*M[k][j] + number_prod(M.matrix[i][k], M.matrix[k][j], &tmp); + number_Qprod_chain(quot(-1,1), &tmp); + number_add_chain(tmp, M.matrix[i]+j); + free_Number(tmp); + } + } + for(i=j+1;i<M.length;i++){ + for(k=0;k<j;k++){ + // -M[i][k]*M[k][j] + number_prod(M.matrix[i][k], M.matrix[k][j], &tmp); + number_Qprod_chain(quot(-1,1), &tmp); + number_add_chain(tmp, M.matrix[i]+j); + free_Number(tmp); + } + } + + // pivot if M[j][j]==0 + // find first M[j][j] that is not 0 + for(pivot=j;pivot<M.length && number_is_zero(M.matrix[pivot][j])==1;pivot++){} + + // no non-zero M[j][j] left: return + if(pivot>=M.length){ + *sign_correction=0; + return(0); + } + // pivot if needed + if(pivot!=j){ + for(k=0;k<M.length;k++){ + tmp=M.matrix[j][k]; + M.matrix[j][k]=M.matrix[pivot][k]; + M.matrix[pivot][k]=tmp; + } + *sign_correction*=-1; + + } + + for(i=j+1;i<M.length;i++){ + // do not use the inplace algorithm if M[j][j] has more than one terms, since it would be modified by the inplace function + if(M.matrix[j][j].length<=1){ + number_quot_inplace(M.matrix[i]+j, M.matrix[j]+j); + } + else{ + number_quot_chain(M.matrix[i]+j, M.matrix[j][j]); + } + } + } + return(0); +} diff --git a/src/determinant.h b/src/determinant.h new file mode 100644 index 0000000..613be87 --- /dev/null +++ b/src/determinant.h @@ -0,0 +1,17 @@ +/* + Compute the determinant of a number matrix +*/ + +#ifndef DETERMINANT_H +#define DETERMINANT_H + +#include "types.h" + +// determinant of a matrix +int determinant_inplace(Number_Matrix M, Number* out); + +// LU decomposition +int LU_dcmp_inplace(Number_Matrix M, int* sign_correction); + + +#endif diff --git a/src/expansions.c b/src/expansions.c deleted file mode 100644 index c889829..0000000 --- a/src/expansions.c +++ /dev/null @@ -1,28 +0,0 @@ -/* -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 "expansions.h" - -#include <stdio.h> -#include <stdlib.h> -#include "definitions.cpp" -#include "tools.h" -#include "array.h" -#include "polynomial.h" -#include "number.h" -#include "rational.h" - - diff --git a/src/expansions.h b/src/expansions.h deleted file mode 100644 index 85d6c2b..0000000 --- a/src/expansions.h +++ /dev/null @@ -1,34 +0,0 @@ -/* -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. -*/ - -/* -Compute exp(V) and log(1+W) -*/ - -#ifndef EXPANSIONS_H -#define EXPANSIONS_H - -#include "polynomial.h" -#include "fields.h" - -// exp(V) -int expand_exponential(Polynomial input_polynomial,Polynomial* output, Fields_Table fields); - -// log(1+W) -int expand_logarithm(Polynomial input_polynomial, Polynomial* output, Fields_Table fields); - - -#endif diff --git a/src/fields.c b/src/fields.c index bfd1f39..7be84df 100644 --- a/src/fields.c +++ b/src/fields.c @@ -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. @@ -23,6 +23,13 @@ limitations under the License. #include "polynomial.h" #include "array.h" #include "rational.h" +#include "tree.h" + +//--------------------------------------------------------------------- +// +// Fields_Table +// +//--------------------------------------------------------------------- // init and free for Fields_Table int init_Fields_Table(Fields_Table* fields){ @@ -30,7 +37,7 @@ int init_Fields_Table(Fields_Table* fields){ 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_Virtual_fields(&((*fields).virtual_fields), FIELDS_SIZE); init_Int_Array(&((*fields).fermions),FIELDS_SIZE); init_Int_Array(&((*fields).noncommuting),FIELDS_SIZE); return(0); @@ -40,7 +47,7 @@ int free_Fields_Table(Fields_Table fields){ free_Int_Array(fields.external); free_Int_Array(fields.internal); free_Identities(fields.ids); - free_Symbols(fields.symbols); + free_Virtual_fields(fields.virtual_fields); free_Int_Array(fields.fermions); free_Int_Array(fields.noncommuting); return(0); @@ -57,11 +64,11 @@ int field_type(int index, Fields_Table fields){ 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); + else if(intlist_find(fields.virtual_fields.indices, fields.virtual_fields.length, index)>=0){ + return(FIELD_VIRTUAL); } - fprintf(stderr,"error: index %d is neither a parameter nor an external or an internal field, nor a symbol\n",index); + fprintf(stderr,"error: index %d is neither a parameter nor an external or an internal field, nor a virtual field\n",index); exit(-1); } @@ -86,7 +93,11 @@ int is_noncommuting(int index, Fields_Table fields){ } -// ------------------ Identities -------------------- +//--------------------------------------------------------------------- +// +// Identities +// +//--------------------------------------------------------------------- // allocate memory int init_Identities(Identities* identities,int size){ @@ -314,6 +325,11 @@ int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fi int match_nc; int first=-1; + // cannot fit + if(test_array.length<input.length){ + return(-1); + } + // bound noncommuting elements while(is_noncommuting(input.values[post_nc], fields)==1){ post_nc++; @@ -324,7 +340,7 @@ int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fi match_nc=1; } for(j=1;j<post_nc;j++){ - if(test_array.values[i+j]!=input.values[j]){ + if(i+j>=test_array.length || test_array.values[i+j]!=input.values[j]){ match_nc=0; } } @@ -359,57 +375,61 @@ int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fi } -// ------------------ Symbols -------------------- +//--------------------------------------------------------------------- +// +// Virtual_fields +// +//--------------------------------------------------------------------- // 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; +int init_Virtual_fields(Virtual_fields* virtual_fields,int size){ + (*virtual_fields).indices=calloc(size,sizeof(int)); + (*virtual_fields).expr=calloc(size,sizeof(Polynomial)); + (*virtual_fields).length=0; + (*virtual_fields).memory=size; return(0); } // free memory -int free_Symbols(Symbols symbols){ +int free_Virtual_fields(Virtual_fields virtual_fields){ int i; - for(i=0;i<symbols.length;i++){ - free_Polynomial(symbols.expr[i]); + for(i=0;i<virtual_fields.length;i++){ + free_Polynomial(virtual_fields.expr[i]); } - free(symbols.indices); - free(symbols.expr); + free(virtual_fields.indices); + free(virtual_fields.expr); return(0); } // resize -int resize_symbols(Symbols* symbols,int new_size){ - Symbols new_symbols; +int resize_virtual_fields(Virtual_fields* virtual_fields,int new_size){ + Virtual_fields new_virtual_fields; 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]; + init_Virtual_fields(&new_virtual_fields,new_size); + for(i=0;i<(*virtual_fields).length;i++){ + new_virtual_fields.indices[i]=(*virtual_fields).indices[i]; + new_virtual_fields.expr[i]=(*virtual_fields).expr[i]; } - new_symbols.length=(*symbols).length; + new_virtual_fields.length=(*virtual_fields).length; - free((*symbols).indices); - free((*symbols).expr); + free((*virtual_fields).indices); + free((*virtual_fields).expr); - *symbols=new_symbols; + *virtual_fields=new_virtual_fields; return(0); } // copy -int symbols_cpy(Symbols input, Symbols* output){ - init_Symbols(output,input.length); - symbols_cpy_noinit(input,output); +int virtual_fields_cpy(Virtual_fields input, Virtual_fields* output){ + init_Virtual_fields(output,input.length); + virtual_fields_cpy_noinit(input,output); return(0); } -int symbols_cpy_noinit(Symbols input, Symbols* output){ +int virtual_fields_cpy_noinit(Virtual_fields input, Virtual_fields* 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); + fprintf(stderr,"error: trying to copy a virtual fields collection of length %d to another with memory %d\n",input.length,(*output).memory); exit(-1); } for(i=0;i<input.length;i++){ @@ -421,12 +441,12 @@ int symbols_cpy_noinit(Symbols input, Symbols* output){ return(0); } -// append an element to a symbols -int symbols_append(int index, Polynomial expr, Symbols* output){ +// append an element to a virtual_fields +int virtual_fields_append(int index, Polynomial expr, Virtual_fields* output){ int offset=(*output).length; if((*output).length>=(*output).memory){ - resize_symbols(output,2*(*output).memory+1); + resize_virtual_fields(output,2*(*output).memory+1); } // copy and allocate @@ -436,12 +456,12 @@ int symbols_append(int index, Polynomial expr, Symbols* output){ (*output).length++; return(0); } -// append an element to a symbols without allocating memory -int symbols_append_noinit(int index, Polynomial expr, Symbols* output){ +// append an element to a virtual_fields without allocating memory +int virtual_fields_append_noinit(int index, Polynomial expr, Virtual_fields* output){ int offset=(*output).length; if((*output).length>=(*output).memory){ - resize_symbols(output,2*(*output).memory+1); + resize_virtual_fields(output,2*(*output).memory+1); } // copy without allocating @@ -452,18 +472,22 @@ int symbols_append_noinit(int index, Polynomial expr, Symbols* output){ return(0); } -// concatenate two symbolss -int symbols_concat(Symbols input, Symbols* output){ +// concatenate two virtual_fields +int virtual_fields_concat(Virtual_fields input, Virtual_fields* output){ int i; for(i=0;i<input.length;i++){ - symbols_append(input.indices[i],input.expr[i],output); + virtual_fields_append(input.indices[i],input.expr[i],output); } return(0); } -// ------------------ Groups -------------------- +//--------------------------------------------------------------------- +// +// Groups +// +//--------------------------------------------------------------------- // allocate memory int init_Groups(Groups* groups,int size){ @@ -570,3 +594,164 @@ int find_group(int index, Groups groups){ } return(-1); } + + +//--------------------------------------------------------------------- +// +// Variables +// +//--------------------------------------------------------------------- + +// allocate memory +int init_Variables(Variables* variables,int size){ + (*variables).var_names=calloc(size,sizeof(Char_Array)); + (*variables).symbol_trees=calloc(size,sizeof(Tree)); + (*variables).length=0; + (*variables).memory=size; + return(0); +} + +// free memory +int free_Variables(Variables variables){ + int i; + for(i=0;i<variables.length;i++){ + free_Char_Array(variables.var_names[i]); + free_Tree(variables.symbol_trees[i]); + } + free(variables.var_names); + free(variables.symbol_trees); + return(0); +} + +// resize +int resize_variables(Variables* variables,int new_size){ + Variables new_variables; + int i; + + init_Variables(&new_variables,new_size); + for(i=0;i<(*variables).length;i++){ + new_variables.var_names[i]=(*variables).var_names[i]; + new_variables.symbol_trees[i]=(*variables).symbol_trees[i]; + } + new_variables.length=(*variables).length; + + free((*variables).var_names); + free((*variables).symbol_trees); + + *variables=new_variables; + return(0); +} + +// copy +int variables_cpy(Variables input, Variables* output){ + init_Variables(output,input.length); + variables_cpy_noinit(input,output); + return(0); +} +int variables_cpy_noinit(Variables input, Variables* output){ + int i; + if((*output).memory<input.length){ + fprintf(stderr,"error: trying to copy a variables collection of length %d to another with memory %d\n",input.length,(*output).memory); + exit(-1); + } + for(i=0;i<input.length;i++){ + char_array_cpy(input.var_names[i], (*output).var_names+i); + tree_cpy(input.symbol_trees[i],(*output).symbol_trees+i); + } + (*output).length=input.length; + + return(0); +} + +// append an element to a variables collection +int variables_append(Char_Array var_name, Tree symbol_tree, Variables* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_variables(output,2*(*output).memory+1); + } + + // copy and allocate + char_array_cpy(var_name,(*output).var_names+offset); + tree_cpy(symbol_tree,(*output).symbol_trees+offset); + // increment length + (*output).length++; + return(0); +} +// append an element to a variables collection without allocating memory +int variables_append_noinit(Char_Array var_name, Tree symbol_tree, Variables* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_variables(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).var_names[offset]=var_name; + (*output).symbol_trees[offset]=symbol_tree; + // increment length + (*output).length++; + return(0); +} + +// concatenate two variables collections +int variables_concat(Variables input, Variables* output){ + int i; + for(i=0;i<input.length;i++){ + variables_append(input.var_names[i], input.symbol_trees[i], output); + } + return(0); +} + + +// find a variable matching a var_name +int variables_find_var(Char_Array name, Variables variables, Tree* output){ + Char_Array varname; + int i; + + // drop inital '$' + char_array_substring(name, 1, name.length-1, &varname); + for(i=0;i<variables.length;i++){ + if(char_array_cmp(varname, variables.var_names[i])==1){ + tree_cpy(variables.symbol_trees[i], output); + break; + } + } + + // error if no variable was found + if(i==variables.length){ + fprintf(stderr, "error: variable '$%s' not found\n",char_array_to_str_noinit(&varname)); + exit(-1); + } + free_Char_Array(varname); + + return(0); +} + + +// add a polynomials as a new named variable +int add_polynomial_to_variables(char* name, Polynomial polynomial, Variables* variables){ + // save polynomial to string (to convert it to a variable, it must first be a string) + Char_Array poly_str; + Char_Array out_name; + Tree out_tree; + + init_Char_Array(&poly_str, STR_SIZE); + polynomial_sprint(polynomial, &poly_str); + + // convert name to Char_Array + init_Char_Array(&out_name,STR_SIZE); + char_array_append_str(name, &out_name); + + // trivial tree containing the polynomial + init_Tree(&out_tree,0,poly_str.length); + tree_set_label(poly_str, &out_tree); + free_Char_Array(poly_str); + + // add variable + variables_append(out_name, out_tree, variables); + free_Tree(out_tree); + free_Char_Array(out_name); + + return 0; +} diff --git a/src/fields.h b/src/fields.h index b6f05c2..7b7a07a 100644 --- a/src/fields.h +++ b/src/fields.h @@ -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. @@ -21,6 +21,9 @@ limitations under the License. #include "types.h" + +// Fields_Table + // init int init_Fields_Table(Fields_Table* fields); int free_Fields_Table(Fields_Table fields); @@ -33,6 +36,8 @@ int is_fermion(int index, Fields_Table fields); int is_noncommuting(int index, Fields_Table fields); +// Identities + // init int init_Identities(Identities* identities,int size); int free_Identities(Identities identities); @@ -57,25 +62,29 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields); int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fields_Table fields); +// Virtual_fields + // init -int init_Symbols(Symbols* symbols,int size); -int free_Symbols(Symbols symbols); +int init_Virtual_fields(Virtual_fields* virtual_fields,int size); +int free_Virtual_fields(Virtual_fields virtual_fields); // resize -int resize_symbols(Symbols* symbols,int new_size); +int resize_virtual_fields(Virtual_fields* virtual_fields,int new_size); // copy -int symbols_cpy(Symbols input, Symbols* output); -int symbols_cpy_noinit(Symbols input, Symbols* output); +int virtual_fields_cpy(Virtual_fields input, Virtual_fields* output); +int virtual_fields_cpy_noinit(Virtual_fields input, Virtual_fields* output); -// append an element to a symbols -int symbols_append(int index, Polynomial expr, Symbols* output); -int symbols_append_noinit(int index, Polynomial expr, Symbols* output); +// append an element to a virtual_fields +int virtual_fields_append(int index, Polynomial expr, Virtual_fields* output); +int virtual_fields_append_noinit(int index, Polynomial expr, Virtual_fields* output); -// concatenate two symbolss -int symbols_concat(Symbols input, Symbols* output); +// concatenate two virtual_fields +int virtual_fields_concat(Virtual_fields input, Virtual_fields* output); +// Groups + // init int init_Groups(Groups* groups,int size); int free_Groups(Groups groups); @@ -98,5 +107,32 @@ int groups_concat(Groups input, Groups* output); int find_group(int index, Groups groups); +// Variables + +// allocate memory +int init_Variables(Variables* variables,int size); +// free memory +int free_Variables(Variables variables); + +// resize +int resize_variables(Variables* variables,int new_size); + +// copy +int variables_cpy(Variables input, Variables* output); +int variables_cpy_noinit(Variables input, Variables* output); + +// append an element to a variables collection +int variables_append(Char_Array var_name, Tree symbol_tree, Variables* output); +int variables_append_noinit(Char_Array var_name, Tree symbol_tree, Variables* output); + +// concatenate two variables collections +int variables_concat(Variables input, Variables* output); + +// find a variable matching a var_name +int variables_find_var(Char_Array name, Variables variables, Tree* output); + +// add a polynomials as a new named variable +int add_polynomial_to_variables(char* name, Polynomial polynomial, Variables* variables); + #define FIELDS_H_DONE #endif @@ -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. @@ -25,14 +25,20 @@ limitations under the License. #include "array.h" #include "coefficient.h" #include "rcc.h" +#include "grouped_polynomial.h" // compute flow numerically, no exponentials -int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, int niter, int display_mode){ +int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Grouped_Polynomial postprocess_flow_equation, Labels labels, int niter, int display_mode){ // running coupling contants RCC rccs=init; + // for printing + RCC rcc_print; int i,j; + // init printing rcc + init_RCC(&rcc_print, rccs.length); + if(display_mode==DISPLAY_NUMERICAL){ // print labels printf("%5s ","n"); @@ -42,9 +48,25 @@ int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, in printf("\n\n"); // print initial values + RCC_cpy_noinit(rccs,&rcc_print); + if(postprocess_flow_equation.length>0){ + // ignore constants + for(j=0;j<rcc_print.length;j++){ + if(rcc_print.indices[j]<0){ + rcc_print.values[j]=1.; + } + } + evaleq(rcc_print, rcc_print, postprocess_flow_equation); + } printf("%5d ",0); - for(j=0;j<rccs.length;j++){ - printf("% 14.7Le ",rccs.values[j]); + for(j=0;j<rcc_print.length;j++){ + // use constants from rcc + if(rcc_print.indices[j]<0){ + printf("% 14.7Le ",rccs.values[j]); + } + else{ + printf("% 14.7Le ",rcc_print.values[j]); + } } printf("\n"); } @@ -52,12 +74,32 @@ int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, in for(i=0;i<niter;i++){ // compute a single step step_flow(&rccs, flow_equation); - // convert ls to alphas + + // print + if(postprocess_flow_equation.length>0){ + RCC_cpy_noinit(rccs,&rcc_print); + // ignore constants + for(j=0;j<rcc_print.length;j++){ + if(rcc_print.indices[j]<0){ + rcc_print.values[j]=1.; + } + } + evaleq(rcc_print, rcc_print, postprocess_flow_equation); + } + else{ + RCC_cpy_noinit(rccs,&rcc_print); + } if(display_mode==DISPLAY_NUMERICAL){ // print the result printf("%5d ",i+1); - for(j=0;j<rccs.length;j++){ - printf("% 14.7Le ",rccs.values[j]); + for(j=0;j<rcc_print.length;j++){ + // use constants from rcc + if(rcc_print.indices[j]<0){ + printf("% 14.7Le ",rccs.values[j]); + } + else{ + printf("% 14.7Le ",rcc_print.values[j]); + } } printf("\n"); } @@ -74,9 +116,24 @@ int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, in } if(display_mode==DISPLAY_FINAL){ - RCC_print(rccs); + if(postprocess_flow_equation.length>0){ + RCC_cpy_noinit(rccs,&rcc_print); + // ignore constants + for(j=0;j<rcc_print.length;j++){ + if(rcc_print.indices[j]<0){ + rcc_print.values[j]=1.; + } + } + evaleq(rcc_print, rcc_print, postprocess_flow_equation); + } + else{ + RCC_cpy_noinit(rccs,&rcc_print); + } + RCC_print(rcc_print); } + free_RCC(rcc_print); + return(0); } @@ -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. @@ -24,7 +24,7 @@ Compute flow numerically #include "types.h" // compute flow -int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, int niter, int display_mode); +int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Grouped_Polynomial postprocess_flow_equation, Labels labels, int niter, int display_mode); // single step int step_flow(RCC* rccs, Grouped_Polynomial flow_equation); diff --git a/src/flow_mpfr.c b/src/flow_mpfr.c index e03c130..1701b15 100644 --- a/src/flow_mpfr.c +++ b/src/flow_mpfr.c @@ -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. @@ -32,14 +32,21 @@ limitations under the License. #include "coefficient.h" #include "flow.h" #include "rcc_mpfr.h" +#include "grouped_polynomial.h" // compute flow numerically -int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels labels, int niter, int display_mode){ +int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Grouped_Polynomial postprocess_flow_equation, Labels labels, int niter, int display_mode){ // running coupling contants RCC_mpfr rccs=init; int i,j; + // for printing + RCC_mpfr rcc_print; + + + // init printing rcc + init_RCC_mpfr(&rcc_print, rccs.length); if(display_mode==DISPLAY_NUMERICAL){ // print labels @@ -50,9 +57,25 @@ int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels printf("\n\n"); // print initial values + RCC_mpfr_cpy_noinit(rccs,&rcc_print); + if(postprocess_flow_equation.length>0){ + // ignore constants + for(j=0;j<rcc_print.length;j++){ + if(rcc_print.indices[j]<0){ + mpfr_set_ui(rcc_print.values[j], 1, MPFR_RNDN); + } + } + evaleq_mpfr(rcc_print, rccs, postprocess_flow_equation); + } printf("%5d ",0); - for(j=0;j<rccs.length;j++){ - mpfr_printf("% 14.7Re ",rccs.values[j]); + for(j=0;j<rcc_print.length;j++){ + // use constants from rcc + if(rcc_print.indices[j]<0){ + mpfr_printf("% 14.7Re ",rccs.values[j]); + } + else{ + mpfr_printf("% 14.7Re ",rcc_print.values[j]); + } } printf("\n"); } @@ -60,12 +83,29 @@ int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels for(i=0;i<niter;i++){ // compute a single step step_flow_mpfr(&rccs, flow_equation); - // convert ls to alphas + + // print + RCC_mpfr_cpy_noinit(rccs,&rcc_print); + if(postprocess_flow_equation.length>0){ + // ignore constants + for(j=0;j<rcc_print.length;j++){ + if(rcc_print.indices[j]<0){ + mpfr_set_ui(rcc_print.values[j], 1, MPFR_RNDN); + } + } + evaleq_mpfr(rcc_print, rccs, postprocess_flow_equation); + } if(display_mode==DISPLAY_NUMERICAL){ // print the result printf("%5d ",i+1); - for(j=0;j<rccs.length;j++){ - mpfr_printf("% 14.7Re ",rccs.values[j]); + for(j=0;j<rcc_print.length;j++){ + // use constants from rcc + if(rcc_print.indices[j]<0){ + mpfr_printf("% 14.7Re ",rccs.values[j]); + } + else{ + mpfr_printf("% 14.7Re ",rcc_print.values[j]); + } } printf("\n"); } @@ -82,9 +122,16 @@ int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels } if(display_mode==DISPLAY_FINAL){ - RCC_mpfr_print(rccs); + if(postprocess_flow_equation.length>0){ + evaleq_mpfr(rcc_print, rccs, postprocess_flow_equation); + } + else{ + rcc_print=rccs; + } + RCC_mpfr_print(rcc_print); } + free_RCC_mpfr(rcc_print); return(0); } diff --git a/src/flow_mpfr.h b/src/flow_mpfr.h index e06cbd8..203770d 100644 --- a/src/flow_mpfr.h +++ b/src/flow_mpfr.h @@ -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. @@ -25,7 +25,7 @@ Compute flow numerically #include "types.h" // compute flow -int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels labels, int niter, int display_mode); +int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Grouped_Polynomial postprocess_flow_equation, Labels labels, int niter, int display_mode); // single step int step_flow_mpfr(RCC_mpfr* rccs, Grouped_Polynomial flow_equation); diff --git a/src/grouped_polynomial.c b/src/grouped_polynomial.c index 0d4d75d..9e245a8 100644 --- a/src/grouped_polynomial.c +++ b/src/grouped_polynomial.c @@ -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. @@ -205,9 +205,9 @@ int group_polynomial(Polynomial polynomial, Grouped_Polynomial* grouped_polynomi if(index==-2){ fprintf(stderr,"error: monomial ("); - for(j=0;j<polynomial.monomials[i].length;j++){ - fprintf(stderr,"%d", polynomial.monomials[i].values[j]); - if(j<polynomial.monomials[i].length-1){ + for(j=0;j<remainder.monomials[i].length;j++){ + fprintf(stderr,"%d", remainder.monomials[i].values[j]); + if(j<remainder.monomials[i].length-1){ fprintf(stderr,","); } } @@ -436,7 +436,7 @@ int simplify_grouped_polynomial(Grouped_Polynomial* polynomial){ } -// derive a flow equation with respect to an unknown variable +// differentiate a flow equation with respect to an unknown variable // equivalent to DB.dl where dl are symbols for the derivatives of the indices in the flow equation with respect to the unknown variable // indices specifies the list of indices that depend on the variable int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Grouped_Polynomial* dflow){ @@ -487,7 +487,7 @@ int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Gr /* -// derive a flow equation with respect to an index +// differentiate a flow equation with respect to an index int flow_equation_deriv(Grouped_Polynomial flow_equation, int index, Grouped_Polynomial* output){ int i,k; // temp list of indices @@ -603,7 +603,7 @@ int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre init_Char_Array(&buffer, STR_SIZE); coefficient_sprint(grouped_polynomial.coefs[i],&buffer,9,rhs_pre); if(buffer.length>0){ - printf("%s",buffer.str); + printf("%s",char_array_to_str_noinit(&buffer)); } free_Char_Array(buffer); @@ -732,29 +732,33 @@ int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output) } -// eValuate an equation on a vector -int evaleq(RCC* rccs, Grouped_Polynomial poly){ +// evaluate an equation on a vector +int evaleq(RCC out, RCC in, Grouped_Polynomial poly){ int i; - long double* res=calloc((*rccs).length,sizeof(long double)); + long double* res=calloc(out.length,sizeof(long double)); - if((*rccs).length!=poly.length){ - fprintf(stderr, "error: trying to evaluate an flow equation with %d components on an rcc with %d\n",poly.length,(*rccs).length); + if(in.length!=poly.length){ + fprintf(stderr, "error: trying to evaluate a flow equation with %d components on an rcc with %d\n",poly.length,in.length); + exit(-1); + } + if(out.length!=poly.length){ + fprintf(stderr, "error: trying to write the output of a flow equation with %d components on an rcc with %d\n",poly.length,out.length); exit(-1); } - // initialize vectors to 0 - for(i=0;i<(*rccs).length;i++){ + // initialize vectors to 0 in an auxiliary vector (to allow for out=in without interference) + for(i=0;i<in.length;i++){ res[i]=0.; } // for each equation for(i=0;i<poly.length;i++){ - evalcoef(*rccs, poly.coefs[i], res+i); + evalcoef(in, poly.coefs[i], res+i); } // copy res to rccs - for(i=0;i<(*rccs).length;i++){ - (*rccs).values[i]=res[i]; + for(i=0;i<out.length;i++){ + out.values[i]=res[i]; } // free memory @@ -763,25 +767,29 @@ int evaleq(RCC* rccs, Grouped_Polynomial poly){ } // evaluate an equation on a vector (using mpfr floats) -int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly){ +int evaleq_mpfr(RCC_mpfr out, RCC_mpfr in, Grouped_Polynomial poly){ int i; mpfr_t* res; - if((*rccs).length!=poly.length){ - fprintf(stderr, "error: trying to evaluate an flow equation with %d components on an rcc with %d\n",poly.length,(*rccs).length); + if(in.length!=poly.length){ + fprintf(stderr, "error: trying to evaluate a flow equation with %d components on an rcc with %d\n",poly.length,in.length); + exit(-1); + } + if(out.length!=poly.length){ + fprintf(stderr, "error: trying to write the output of a flow equation with %d components on an rcc with %d\n",poly.length,out.length); exit(-1); } - res=calloc((*rccs).length,sizeof(mpfr_t)); + res=calloc(out.length,sizeof(mpfr_t)); // for each equation for(i=0;i<poly.length;i++){ - evalcoef_mpfr(*rccs, poly.coefs[i], res[i]); + evalcoef_mpfr(in, poly.coefs[i], res[i]); } // copy res to rccs - for(i=0;i<(*rccs).length;i++){ - mpfr_set((*rccs).values[i], res[i], MPFR_RNDN); + for(i=0;i<out.length;i++){ + mpfr_set(out.values[i], res[i], MPFR_RNDN); mpfr_clear(res[i]); } @@ -791,3 +799,85 @@ int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly){ } +// compose two flow equations (replace the rcc's of flow1 by the right hand side of flow2) +int compose_flow_equations(Grouped_Polynomial flow1, Grouped_Polynomial flow2, Grouped_Polynomial* out){ + if(flow1.length!=flow2.length){ + fprintf(stderr, "error: trying to compose two flow equations of different size\n"); + exit(-1); + } + int i,j,k; + Coefficient constant; + + // init + init_Grouped_Polynomial(out, flow1.length); + (*out).length=flow1.length; + + // init constant (so we can tell when the constant was not found) + constant.length=0; + + // loop over rcc's + for(i=0;i<flow1.length;i++){ + // set indices + (*out).indices[i]=flow1.indices[i]; + + // passthrough constant terms + if((*out).indices[i]<0){ + int index=intlist_find_err(flow2.indices,flow2.length,(*out).indices[i]); + coefficient_cpy(flow2.coefs[index], (*out).coefs+i); + constant=flow2.coefs[index]; + continue; + } + + // init + init_Coefficient((*out).coefs+i, COEF_SIZE); + + // loop over terms + for(j=0;j<flow1.coefs[i].length;j++){ + Coefficient tmp_coef; + + // init + init_Coefficient(&tmp_coef, COEF_SIZE); + // init factor + Int_Array tmp_factor; + init_Int_Array(&tmp_factor, MONOMIAL_SIZE); + // init denom + coef_denom denom; + // index should be that appearing in flow2 + if(flow2.coefs[i].length<1){ + fprintf(stderr,"error: composing two flow equations: the %d-th term in the flow equation is empty\n",flow1.indices[i]); + exit(-1); + } + denom.index=flow2.coefs[i].denoms[0].index; + denom.power=0; + // init num + Number tmp_num; + number_cpy(flow1.coefs[i].nums[j], &tmp_num); + + // init coefficient with numerical prefactor + coefficient_append_noinit(tmp_factor, tmp_num, denom, &tmp_coef); + + // loop over factors + for(k=0;k<flow1.coefs[i].factors[j].length;k++){ + // multiply factors together + coefficient_prod_chain(flow2.coefs[intlist_find_err(flow2.indices,flow2.length,flow1.coefs[i].factors[j].values[k])], &tmp_coef); + } + + // add to out + coefficient_concat_noinit(tmp_coef, (*out).coefs+i); + } + + } + + // simplify fractions + if(constant.length!=0){ + for(i=0;i<(*out).length;i++){ + if((*out).indices[i]>=0){ + // reduce them to a common denominator (not much is gained from trying to simplify them) + coefficient_common_denominator(constant, (*out).coefs+i); + //coefficient_simplify_rational(constant, (*out).coefs+i); + } + } + } + + return(0); +} diff --git a/src/grouped_polynomial.h b/src/grouped_polynomial.h index e369721..7d95a2e 100644 --- a/src/grouped_polynomial.h +++ b/src/grouped_polynomial.h @@ -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. @@ -59,7 +59,7 @@ int find_id(Int_Array monomial, Id_Table idtable, int start); // simplify grouped polynomial int simplify_grouped_polynomial(Grouped_Polynomial* polynomial); -// derive a flow equation with respect to an unknown variable +// differentiate a flow equation with respect to an unknown variable int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Grouped_Polynomial* dflow); // print a grouped polynomial @@ -69,8 +69,10 @@ int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output); // evaluate an equation on an RCC -int evaleq(RCC* rccs, Grouped_Polynomial poly); +int evaleq(RCC out, RCC in, Grouped_Polynomial poly); // evaluate an equation on a vector (using mpfr floats) -int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly); +int evaleq_mpfr(RCC_mpfr out, RCC_mpfr in, Grouped_Polynomial poly); +// compose two flow equations (replace the rcc's of flow1 by the right hand side of flow2) +int compose_flow_equations(Grouped_Polynomial flow1, Grouped_Polynomial flow2, Grouped_Polynomial* out); #endif diff --git a/src/idtable.c b/src/idtable.c index fe409c2..fc72f56 100644 --- a/src/idtable.c +++ b/src/idtable.c @@ -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. diff --git a/src/idtable.h b/src/idtable.h index 38ab249..aa99dca 100644 --- a/src/idtable.h +++ b/src/idtable.h @@ -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. diff --git a/src/istring.c b/src/istring.c index 663c06a..d78f67b 100644 --- a/src/istring.c +++ b/src/istring.c @@ -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. diff --git a/src/istring.h b/src/istring.h index 5b47b9b..5027981 100644 --- a/src/istring.h +++ b/src/istring.h @@ -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. diff --git a/src/kondo.c b/src/kondo.c index 79b7c56..93318d5 100644 --- a/src/kondo.c +++ b/src/kondo.c @@ -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. @@ -35,7 +35,7 @@ limitations under the License. #define KONDO_SPIN 2 // offsets for indices of A, B, h and t -// order matters for symbols table +// order matters for virtual_fields table #define KONDO_A_OFFSET 1 #define KONDO_B_OFFSET 2 #define KONDO_H_OFFSET 3 @@ -72,6 +72,7 @@ limitations under the License. int kondo_generate_conf(Str_Array* str_args, int box_count){ Str_Array new_args; Fields_Table fields; + Variables variables; Char_Array tmp_str; int arg_index; int i; @@ -83,16 +84,19 @@ int kondo_generate_conf(Str_Array* str_args, int box_count){ kondo_fields_table(box_count, &tmp_str, &fields); str_array_append_noinit(tmp_str, &new_args); - // symbols - kondo_symbols(&tmp_str, box_count, &fields); - arg_index=find_str_arg("symbols", *str_args); + // dummy variables + init_Variables(&variables,1); + + // virtual fields + kondo_virtual_fields(&tmp_str, box_count, &fields); + arg_index=find_str_arg("virtual_fields", *str_args); if(arg_index>=0){ if(tmp_str.length>0){ char_array_snprintf(&tmp_str,",\n"); } char_array_concat((*str_args).strs[arg_index], &tmp_str); } - parse_input_symbols(tmp_str, &fields); + parse_input_virtual_fields(tmp_str, &fields, variables); str_array_append_noinit(tmp_str, &new_args); // identities @@ -104,7 +108,7 @@ int kondo_generate_conf(Str_Array* str_args, int box_count){ } char_array_concat((*str_args).strs[arg_index], &tmp_str); } - parse_input_identities(tmp_str, &fields); + parse_input_identities(tmp_str, &fields, variables); str_array_append_noinit(tmp_str, &new_args); // groups @@ -136,7 +140,7 @@ int kondo_generate_conf(Str_Array* str_args, int box_count){ // copy remaining entries for(i=0;i<(*str_args).length;i++){ get_str_arg_title((*str_args).strs[i], &title); - if(str_cmp(title.str, "symbols")==0 &&\ + if(str_cmp(title.str, "virtual_fields")==0 &&\ str_cmp(title.str, "identities")==0 &&\ str_cmp(title.str, "propagator")==0 &&\ str_cmp(title.str, "input_polynomial")==0 &&\ @@ -149,6 +153,7 @@ int kondo_generate_conf(Str_Array* str_args, int box_count){ } free_Fields_Table(fields); + free_Variables(variables); free_Str_Array(*str_args); *str_args=new_args; @@ -250,15 +255,15 @@ int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fiel } -// generate Kondo symbols -int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){ +// generate Kondo virtual_fields +int kondo_virtual_fields(Char_Array* str_virtual_fields, int box_count, Fields_Table* fields){ int i,j,k,l; Char_Array tmp_str; Polynomial poly; char letters[3]={'A','B','h'}; - init_Char_Array(str_symbols, STR_SIZE); - char_array_snprintf(str_symbols, "#!symbols\n"); + init_Char_Array(str_virtual_fields, STR_SIZE); + char_array_snprintf(str_virtual_fields, "#!virtual_fields\n"); // loop over box index for(i=1;i<=box_count;i++){ @@ -267,7 +272,7 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){ // loop over space dimension for(k=0;k<KONDO_DIM;k++){ // write index - char_array_snprintf(str_symbols, "%d=", 100*(10*(KONDO_A_OFFSET+j)+k)+i); + char_array_snprintf(str_virtual_fields, "%d=", 100*(10*(KONDO_A_OFFSET+j)+k)+i); // write the name of the scalar product init_Char_Array(&tmp_str, 6); char_array_snprintf(&tmp_str, "%c%d%d", letters[j], k, i); @@ -275,10 +280,10 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){ kondo_resolve_ABht(tmp_str.str, &poly, *fields); free_Char_Array(tmp_str); // write to output - polynomial_sprint(poly, str_symbols); + polynomial_sprint(poly, str_virtual_fields); free_Polynomial(poly); // add , - char_array_snprintf(str_symbols,",\n"); + char_array_snprintf(str_virtual_fields,",\n"); } } } @@ -290,46 +295,46 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){ for(j=0;j<3;j++){ for(k=0;k<3;k++){ // write index - char_array_snprintf(str_symbols, "%d=", 1000*(10*(KONDO_A_OFFSET+j)+KONDO_A_OFFSET+k)+i); + char_array_snprintf(str_virtual_fields, "%d=", 1000*(10*(KONDO_A_OFFSET+j)+KONDO_A_OFFSET+k)+i); for(l=0;l<KONDO_DIM;l++){ - char_array_snprintf(str_symbols, "(1)"); + char_array_snprintf(str_virtual_fields, "(1)"); if(j<2){ - char_array_snprintf(str_symbols,"[f%d]", 100*(10*(KONDO_A_OFFSET+j)+l)+i); + char_array_snprintf(str_virtual_fields,"[f%d]", 100*(10*(KONDO_A_OFFSET+j)+l)+i); } else{ - char_array_snprintf(str_symbols,"[f%d]", 10*(10*(KONDO_A_OFFSET+j)+l)); + char_array_snprintf(str_virtual_fields,"[f%d]", 10*(10*(KONDO_A_OFFSET+j)+l)); } if(k<2){ - char_array_snprintf(str_symbols,"[f%d]", 100*(10*(KONDO_A_OFFSET+k)+l)+i); + char_array_snprintf(str_virtual_fields,"[f%d]", 100*(10*(KONDO_A_OFFSET+k)+l)+i); } else{ - char_array_snprintf(str_symbols,"[f%d]", 10*(10*(KONDO_A_OFFSET+k)+l)); + char_array_snprintf(str_virtual_fields,"[f%d]", 10*(10*(KONDO_A_OFFSET+k)+l)); } if(l<KONDO_DIM-1){ - char_array_append('+',str_symbols); + char_array_append('+',str_virtual_fields); } } // add , - char_array_snprintf(str_symbols,",\n"); + char_array_snprintf(str_virtual_fields,",\n"); } } } // vector products for(i=1;i<=box_count;i++){ - char_array_snprintf(str_symbols, "%d=", 100*(100*(KONDO_A_OFFSET)+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+i); + char_array_snprintf(str_virtual_fields, "%d=", 100*(100*(KONDO_A_OFFSET)+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+i); for(l=0;l<KONDO_DIM;l++){ // remember (-1 %3 = -1) - char_array_snprintf(str_symbols, "(1)[f%d][f%d][f%d]+(-1)[f%d][f%d][f%d]", 100*(10*KONDO_A_OFFSET+((l+1)%KONDO_DIM))+i, 100*(10*KONDO_B_OFFSET+((l+2)%KONDO_DIM))+i, 10*(10*KONDO_H_OFFSET+l), 100*(10*KONDO_A_OFFSET+((l+2)%KONDO_DIM))+i, 100*(10*KONDO_B_OFFSET+((l+1)%KONDO_DIM))+i, 10*(10*KONDO_H_OFFSET+l)); + char_array_snprintf(str_virtual_fields, "(1)[f%d][f%d][f%d]+(-1)[f%d][f%d][f%d]", 100*(10*KONDO_A_OFFSET+((l+1)%KONDO_DIM))+i, 100*(10*KONDO_B_OFFSET+((l+2)%KONDO_DIM))+i, 10*(10*KONDO_H_OFFSET+l), 100*(10*KONDO_A_OFFSET+((l+2)%KONDO_DIM))+i, 100*(10*KONDO_B_OFFSET+((l+1)%KONDO_DIM))+i, 10*(10*KONDO_H_OFFSET+l)); if(l<KONDO_DIM-1){ - char_array_append('+',str_symbols); + char_array_append('+',str_virtual_fields); } } // add , if(i<box_count){ - char_array_snprintf(str_symbols,",\n"); + char_array_snprintf(str_virtual_fields,",\n"); } } @@ -778,13 +783,13 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_ // if polynomial exists, add to each monomial if(tmp_poly.length>0){ for(i=0;i<tmp_poly.length;i++){ - int_array_append(get_symbol_index(buffer), tmp_poly.monomials+i); + int_array_append(get_virtual_field_index(buffer), tmp_poly.monomials+i); } } // if not, create a new term in the polynomial else{ init_Int_Array(&tmp_monomial, MONOMIAL_SIZE); - int_array_append(get_symbol_index(buffer), &tmp_monomial); + int_array_append(get_virtual_field_index(buffer), &tmp_monomial); init_Int_Array(&dummy_factor, 1); polynomial_append_noinit(tmp_monomial, dummy_factor, number_one(), &tmp_poly); } @@ -1387,8 +1392,8 @@ int get_offsets_index(char* str, int* offset1, int* offset2, int* index){ return(0); } -// get the index of the symbol corresponding to a given string -int get_symbol_index(char* str){ +// get the index of the virtual_field corresponding to a given string +int get_virtual_field_index(char* str){ char* ptr; int offset=-1; int index=0; @@ -1439,7 +1444,7 @@ int get_symbol_index(char* str){ if(offset==-1){ return(-1); } - // no symbol for h or t + // no virtual field for h or t if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){ return(10*(10*offset+dim)); } diff --git a/src/kondo.h b/src/kondo.h index c534145..18488ae 100644 --- a/src/kondo.h +++ b/src/kondo.h @@ -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. @@ -27,10 +27,10 @@ int kondo_generate_conf(Str_Array* str_args, int box_count); // generate the Kondo fields table int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fields); -// generate Kondo symbols -int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields); -// generate Kondo symbols (older method: one symbol for each scalar product) -int kondo_symbols_scalarprod(Char_Array* str_symbols, int box_count, Fields_Table* fields); +// generate Kondo virtual_fields +int kondo_virtual_fields(Char_Array* str_virtual_fields, int box_count, Fields_Table* fields); +// generate Kondo virtual_fields (older method: one virtual_field for each scalar product) +int kondo_virtual_fields_scalarprod(Char_Array* str_virtual_fields, int box_count, Fields_Table* fields); // generate Kondo groups (groups of independent variables) int kondo_groups(Char_Array* str_groups, int box_count); @@ -71,6 +71,6 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output); int get_offset_index(char* str, int* offset, int* index); // get the offsets and index of a scalar product int get_offsets_index(char* str, int* offset1, int* offset2, int* index); -// get the index of the symbol corresponding to a given string -int get_symbol_index(char* str); +// get the index of the virtual_field corresponding to a given string +int get_virtual_field_index(char* str); #endif diff --git a/src/kondo_preprocess.c b/src/kondo_preprocess.c index 3371014..c472e9d 100644 --- a/src/kondo_preprocess.c +++ b/src/kondo_preprocess.c @@ -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. @@ -104,6 +104,11 @@ for(i=1;i<argc;i++){ // number of dimensions else if (flag==CP_FLAG_DIMENSION){ sscanf(argv[i],"%d",&((*opts).dimension)); + // check value of the dimension + if((*opts).dimension<=0 || (*opts).dimension>=4){ + fprintf(stderr,"error: kondo_preprocess only supports dimensions 1, 2 and 3 (got %d)\n",(*opts).dimension); + exit(-1); + } flag=0; } // read file name from command-line @@ -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 @@ -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. @@ -28,6 +28,12 @@ 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); + +// compute the mean of a monomial by computing a determinant +int mean_determinant(Int_Array internal_plus, Int_Array internal_minus, Number* out, Polynomial_Matrix propagator,Fields_Table fields); + +// 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); // 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); // next pairing with a non-vanishing propagator @@ -43,12 +49,12 @@ int mean_internal_slow(Int_Array internal_plus, Int_Array internal_minus, Number // 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); -// compute the mean of a monomial containing symbolic expressions -int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed); +// compute the mean of a monomial containing virtual fields +int mean_virtual_fields(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed); // 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); // 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); // 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); @@ -62,9 +68,9 @@ int sort_fermions(int* array, int begin, int end, int* sign); int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed); // compute the mean of a polynomial -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); // 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); // single thread mean void* polynomial_mean_thread(void* mean_args); #endif diff --git a/src/meankondo.c b/src/meankondo.c index 1c976d4..8b5c197 100644 --- a/src/meankondo.c +++ b/src/meankondo.c @@ -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. @@ -17,7 +17,7 @@ limitations under the License. /* meankondo -A simple tool to compute the renormalization group flow for Fermionic hierarchical models +A tool to compute the renormalization group flow for Fermionic hierarchical models */ @@ -49,15 +49,19 @@ A simple tool to compute the renormalization group flow for Fermionic hierarchic #include "mean.h" // various string operations #include "istring.h" +// symbolic trees +# include "tree.h" // read cli arguments int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meankondo_Options* opts); // print usage message int print_usage_meankondo(); +// check consistency of options +int check_meankondo_opts(Meankondo_Options opts); // compute flow int compute_flow(Str_Array str_args, Meankondo_Options opts); -// compute the flow equation -int compute_flow_equation(Polynomial init_poly, Id_Table idtable, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, Grouped_Polynomial* flow_equation); +// compute average +int compute_average(Polynomial init_poly, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, int print_progress, Polynomial* exp_poly); int main (int argc, const char* argv[]){ @@ -69,6 +73,9 @@ int main (int argc, const char* argv[]){ // read command-line arguments read_args_meankondo(argc,argv,&str_args,&opts); + // check command-line arguments + check_meankondo_opts(opts); + // warning message if representing rational numbers as floats #ifdef RATIONAL_AS_FLOAT fprintf(stderr,"info: representing rational numbers using floats\n"); @@ -101,6 +108,10 @@ int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meanko (*opts).threads=1; // do not chain (*opts).chain=0; + // do not print progress + (*opts).print_progress=0; + // print the flow equation + (*opts).group_poly=1; // loop over arguments for(i=1;i<argc;i++){ @@ -116,6 +127,13 @@ int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meanko case 'C': (*opts).chain=1; break; + // print progress + case 'p': + (*opts).print_progress=1; + break; + case 'A': + (*opts).group_poly=0; + break; // print version case 'v': printf("meankondo " VERSION "\n"); @@ -147,7 +165,16 @@ int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meanko // print usage message int print_usage_meankondo(){ - printf("\nusage:\n meankondo [-t threads] [-C] <filename>\n\n"); + printf("\nusage:\n meankondo [-t threads] [-C] [-p] [-A] <filename>\n\n"); + return(0); +} + +// check consistency of options +int check_meankondo_opts(Meankondo_Options opts){ + if(opts.chain==1 && opts.group_poly==0){ + fprintf(stderr,"aborting: the '-C' and '-A' options are incompatible\n"); + exit(-1); + } return(0); } @@ -163,6 +190,8 @@ int compute_flow(Str_Array str_args, Meankondo_Options opts){ Fields_Table fields; // their propagator Polynomial_Matrix propagator; + // preprocessor variables + Variables variables; // initial polynomial Polynomial init_poly; // list of rccs @@ -171,6 +200,8 @@ int compute_flow(Str_Array str_args, Meankondo_Options opts){ Groups groups; // flow equation Grouped_Polynomial flow_equation; + // polynomial produced by the averaging operation + Polynomial exp_poly; // parse fields @@ -183,29 +214,38 @@ int compute_flow(Str_Array str_args, Meankondo_Options opts){ parse_input_fields(str_args.strs[arg_index],&fields); } - // parse id table - arg_index=find_str_arg("id_table", str_args); - if(arg_index<0){ - fprintf(stderr,"error: no id table entry in the configuration file\n"); - exit(-1); + // parse variables + // must precede id_table, virtual_fields, identities and input_polynomial + arg_index=find_str_arg("preprocessor_variables", str_args); + if(arg_index>=0){ + parse_input_variables(str_args.strs[arg_index],&variables); } else{ - parse_input_id_table(str_args.strs[arg_index],&idtable, fields); + init_Variables(&variables,1); } - // parse symbols - arg_index=find_str_arg("symbols", str_args); - if(arg_index>=0){ - parse_input_symbols(str_args.strs[arg_index],&fields); + // parse id table + if(opts.group_poly==1){ + arg_index=find_str_arg("id_table", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no id table entry in the configuration file\n"); + exit(-1); + } + else{ + parse_input_id_table(str_args.strs[arg_index],&idtable, fields, variables); + } } - else{ - init_Symbols(&(fields.symbols),1); + + // parse virtual_fields + arg_index=find_str_arg("virtual_fields", str_args); + if(arg_index>=0){ + parse_input_virtual_fields(str_args.strs[arg_index], &fields, variables); } // parse input polynomial arg_index=find_str_arg("input_polynomial", str_args); if(arg_index>=0){ - parse_input_polynomial(str_args.strs[arg_index],&init_poly, fields); + parse_input_polynomial(str_args.strs[arg_index],&init_poly, fields, variables); } else{ fprintf(stderr,"error: no input polynomial entry in the configuration file\n"); @@ -225,41 +265,90 @@ int compute_flow(Str_Array str_args, Meankondo_Options opts){ // parse identities arg_index=find_str_arg("identities", str_args); if(arg_index>=0){ - parse_input_identities(str_args.strs[arg_index],&fields); - } - else{ - init_Identities(&(fields.ids),1); + parse_input_identities(str_args.strs[arg_index],&fields, variables); } - // parse groups + // parse groups (must come after virtual_fields and propagator) arg_index=find_str_arg("groups", str_args); if(arg_index>=0){ - parse_input_groups(str_args.strs[arg_index],&groups); + parse_input_groups(str_args.strs[arg_index],&groups, propagator, fields); } else{ init_Groups(&groups, 1); } - // flow equation - compute_flow_equation(init_poly, idtable, fields, propagator, groups, opts.threads, &flow_equation); + // compute the average + compute_average(init_poly, fields, propagator, groups, opts.threads, opts.print_progress, &exp_poly); + free_Polynomial(init_poly); free_Polynomial_Matrix(propagator); - free_Fields_Table(fields); free_Groups(groups); + // parse postprocessing entry + arg_index=find_str_arg("postprocess_operation", str_args); + if(arg_index>=0){ + add_polynomial_to_variables("OUT", exp_poly, &variables); + // parse postprocess entry + Polynomial postprocess_operation; + parse_input_polynomial(str_args.strs[arg_index], &postprocess_operation, fields, variables); + // replace exp_poly + free_Polynomial(exp_poly); + exp_poly=postprocess_operation; + } + + if(opts.group_poly==1){ + // flow equation + group_polynomial(exp_poly, &flow_equation, idtable, fields); + } + + // postprocess flow equation + arg_index=find_str_arg("postprocess_flow_equation", str_args); + if(arg_index>=0){ + Polynomial flow_polynomial; + // polynomial made of the rcc's multiplied by the corresponding fields (parsed from idtable) + idtable_to_polynomial(idtable, &flow_polynomial); + + // add to variables + add_polynomial_to_variables("FLOW", flow_polynomial, &variables); + free_Polynomial(flow_polynomial); + + // parse postprocess entry + Polynomial postprocess_polynomial; + parse_input_polynomial(str_args.strs[arg_index], &postprocess_polynomial, fields, variables); + + // convert to flow equation + Grouped_Polynomial postprocess_flow_equation; + group_polynomial(postprocess_polynomial, &postprocess_flow_equation, idtable, fields); + free_Polynomial(postprocess_polynomial); + + // apply postprocessing to flow equation + Grouped_Polynomial new_flow; + compose_flow_equations(postprocess_flow_equation, flow_equation, &new_flow); + free_Grouped_Polynomial(postprocess_flow_equation); + + // replace flow_equation + free_Grouped_Polynomial(flow_equation); + flow_equation=new_flow; + } + + // if chain then print config file if(opts.chain==1){ for(i=0;i<str_args.length;i++){ // check whether to print the str_arg get_str_arg_title(str_args.strs[i], &arg_header); if (\ - str_cmp(arg_header.str, "symbols")==0 &&\ + str_cmp(arg_header.str, "variables")==0 &&\ + str_cmp(arg_header.str, "virtual_fields")==0 &&\ str_cmp(arg_header.str, "groups")==0 &&\ str_cmp(arg_header.str, "fields")==0 &&\ str_cmp(arg_header.str, "identities")==0 &&\ str_cmp(arg_header.str, "propagator")==0 &&\ str_cmp(arg_header.str, "input_polynomial")==0 &&\ - str_cmp(arg_header.str, "id_table")==0 ){ + str_cmp(arg_header.str, "id_table")==0 &&\ + str_cmp(arg_header.str, "postprocess_operation")==0 &&\ + str_cmp(arg_header.str, "numerical_postprocess_operation")==0 + ){ printf("%s\n&\n",str_args.strs[i].str); } free_Char_Array(arg_header); @@ -268,34 +357,66 @@ int compute_flow(Str_Array str_args, Meankondo_Options opts){ printf("#!flow_equation\n"); } - // print flow equation - grouped_polynomial_print(flow_equation,'%','%'); + // print result + if(opts.group_poly==1){ + grouped_polynomial_print(flow_equation,'%','%'); + + // free memory + free_Grouped_Polynomial(flow_equation); + } + else{ + polynomial_print(exp_poly); + } // free memory - free_Id_Table(idtable); - free_Grouped_Polynomial(flow_equation); + free_Polynomial(exp_poly); + + // parse numerical_postprocessing entry + arg_index=find_str_arg("numerical_postprocess_operation", str_args); + if(arg_index>=0){ + Polynomial rcc_polynomial; + // polynomial made of the rcc's multiplied by the corresponding fields (parsed from idtable) + idtable_to_polynomial(idtable, &rcc_polynomial); + + // add to variables + add_polynomial_to_variables("RCC", rcc_polynomial, &variables); + free_Polynomial(rcc_polynomial); + + // parse postprocess entry + Polynomial numerical_postprocess_operation; + parse_input_polynomial(str_args.strs[arg_index], &numerical_postprocess_operation, fields, variables); + + // convert to flow equation + Grouped_Polynomial numerical_postprocess_flow_equation; + group_polynomial(numerical_postprocess_operation, &numerical_postprocess_flow_equation, idtable, fields); + free_Polynomial(numerical_postprocess_operation); + + // print postprocessing flow equation + printf("\n&\n#!postprocess_operation\n"); + grouped_polynomial_print(numerical_postprocess_flow_equation,'%','%'); + free_Grouped_Polynomial(numerical_postprocess_flow_equation); + } + + if(opts.group_poly==1){ + free_Id_Table(idtable); + } + free_Fields_Table(fields); + free_Variables(variables); + return(0); } -// compute the flow equation -int compute_flow_equation(Polynomial init_poly, Id_Table idtable, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, Grouped_Polynomial* flow_equation){ - // expectation - Polynomial exp_poly; - - polynomial_cpy(init_poly,&exp_poly); +// compute average +int compute_average(Polynomial init_poly, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, int print_progress, Polynomial* exp_poly){ + polynomial_cpy(init_poly,exp_poly); // average if(threads>1){ - polynomial_mean_multithread(&exp_poly, fields, propagator, groups, threads); + polynomial_mean_multithread(exp_poly, fields, propagator, groups, threads, print_progress); } else{ - polynomial_mean(&exp_poly, fields, propagator, groups); + polynomial_mean(exp_poly, fields, propagator, groups, print_progress); } - - // grouped representation of expanded_poly - group_polynomial(exp_poly,flow_equation,idtable, fields); - free_Polynomial(exp_poly); - return(0); } diff --git a/src/meantools.c b/src/meantools.c index e3e7247..64cf3f1 100644 --- a/src/meantools.c +++ b/src/meantools.c @@ -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. @@ -39,13 +39,13 @@ Utility to perform various operations on flow equations // string functions #include "istring.h" // tools -#include "meantools_exp.h" #include "meantools_deriv.h" #include "meantools_eval.h" +#include "meantools_expand.h" -#define EXP_COMMAND 1 -#define DERIV_COMMAND 2 -#define EVAL_COMMAND 3 +#define DERIV_COMMAND 1 +#define EVAL_COMMAND 2 +#define EXPAND_COMMAND 3 // read cli arguments int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meantools_Options* opts); @@ -63,12 +63,12 @@ int main (int argc, const char* argv[]){ read_args_meantools(argc,argv,&str_args, &opts); switch(opts.command){ - case EXP_COMMAND: tool_exp(str_args); - break; case DERIV_COMMAND: tool_deriv(str_args,opts); break; case EVAL_COMMAND: tool_eval(str_args,opts); break; + case EXPAND_COMMAND: tool_expand(str_args,opts); + break; } //free memory @@ -86,11 +86,7 @@ int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meanto exit(-1); } - if(str_cmp((char*)argv[1],"exp")==1){ - (*opts).command=EXP_COMMAND; - tool_exp_read_args(argc, argv, str_args); - } - else if(str_cmp((char*)argv[1],"derive")==1){ + if(str_cmp((char*)argv[1],"differentiate")==1){ (*opts).command=DERIV_COMMAND; tool_deriv_read_args(argc, argv, str_args, opts); } @@ -98,6 +94,10 @@ int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meanto (*opts).command=EVAL_COMMAND; tool_eval_read_args(argc, argv, str_args, opts); } + else if(str_cmp((char*)argv[1],"expand")==1){ + (*opts).command=EXPAND_COMMAND; + tool_expand_read_args(argc, argv, str_args, opts); + } else{ print_usage_meantools(); exit(-1); @@ -108,9 +108,6 @@ int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meanto // print usage message int print_usage_meantools(){ - printf("\nusage:\n meantools exp [config_file]\n meantools derive [-d derivatives] [-V variables] [-C] [config_file]\n meantools eval [-R values] [-P precision] [-E max_exponent] [config_file]\n\n"); + printf("\nusage:\n meantools differentiate [-d derivatives] [-V variables] [-C] [config_file]\n meantools eval [-R values] [-P precision] [-E max_exponent] [config_file]\n meantools expand [-N namespace] [config_file]\n\n"); return(0); } - - - diff --git a/src/meantools_deriv.c b/src/meantools_deriv.c index b096380..3530808 100644 --- a/src/meantools_deriv.c +++ b/src/meantools_deriv.c @@ -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. @@ -42,9 +42,9 @@ int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Mean char* ptr; // defaults - // derive once + // differentiate once (*opts).deriv_derivs=1; - // derive with respect to all variables + // differentiate with respect to all variables (*opts).deriv_vars.length=-1; // do not chain (*opts).chain=0; @@ -77,7 +77,7 @@ int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Mean } // variables else if(flag==CP_FLAG_VARS){ - // if the argument is "all" then derive wrt all variables + // if the argument is "all" then differentiate wrt all variables if(str_cmp((char*)argv[i],"all")){ (*opts).deriv_vars.length=-2; } @@ -101,7 +101,7 @@ int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Mean } -// derive a flow equation +// differentiate a flow equation int tool_deriv(Str_Array str_args, Meantools_Options opts){ // index of the entry in the input file int arg_index; @@ -158,7 +158,7 @@ int tool_deriv(Str_Array str_args, Meantools_Options opts){ get_str_arg_title(str_args.strs[i], &arg_header); if (\ str_cmp(arg_header.str, "flow_equation")==0 &&\ - str_cmp(arg_header.str, "symbols")==0 &&\ + str_cmp(arg_header.str, "virtual_fields")==0 &&\ str_cmp(arg_header.str, "groups")==0 &&\ str_cmp(arg_header.str, "fields")==0 &&\ str_cmp(arg_header.str, "identities")==0 &&\ @@ -208,7 +208,7 @@ int flow_equation_derivative(int n, Int_Array variables, Grouped_Polynomial flow // free tmpflow free_Grouped_Polynomial(tmpflow); - // add the derived indices as variables for the next derivative + // add the differentiated indices as variables for the next derivative for(i=0;i<variables.length;i++){ if(variables.values[i]>=0){ int_array_append((j+1)*DOFFSET+variables.values[i], &indices); diff --git a/src/meantools_deriv.h b/src/meantools_deriv.h index 4ea36a9..7425fd1 100644 --- a/src/meantools_deriv.h +++ b/src/meantools_deriv.h @@ -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. @@ -21,7 +21,7 @@ limitations under the License. // read arguments int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts); -// derive a flow equation +// differentiate a flow equation int tool_deriv(Str_Array str_args, Meantools_Options opts); // n first derivatives of a flow equation wrt to variables int flow_equation_derivative(int n, Int_Array variables, Grouped_Polynomial flow_equation, Grouped_Polynomial* flow_equation_derivs); diff --git a/src/meantools_eval.c b/src/meantools_eval.c index 4f166cd..027d18d 100644 --- a/src/meantools_eval.c +++ b/src/meantools_eval.c @@ -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. @@ -163,12 +163,12 @@ int tool_eval(Str_Array str_args, Meantools_Options opts){ // evaluate if(mpfr_flag==0){ - evaleq(&rccs, flow_equation); + evaleq(rccs, rccs, flow_equation); RCC_print(rccs); free_RCC(rccs); } else{ - evaleq_mpfr(&rccs_mpfr, flow_equation); + evaleq_mpfr(rccs_mpfr, rccs_mpfr, flow_equation); RCC_mpfr_print(rccs_mpfr); free_RCC_mpfr(rccs_mpfr); } diff --git a/src/meantools_eval.h b/src/meantools_eval.h index 518049d..54f6e51 100644 --- a/src/meantools_eval.h +++ b/src/meantools_eval.h @@ -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. diff --git a/src/meantools_exp.c b/src/meantools_exp.c deleted file mode 100644 index 1aca928..0000000 --- a/src/meantools_exp.c +++ /dev/null @@ -1,130 +0,0 @@ -/* -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 "meantools_exp.h" - -#include <stdio.h> -#include <stdlib.h> -#include "parse_file.h" -#include "cli_parser.h" -#include "polynomial.h" -#include "fields.h" -#include "grouped_polynomial.h" -#include "idtable.h" - -// read command line arguments -int tool_exp_read_args(int argc, const char* argv[], Str_Array* str_args){ - // file to read the polynomial from in flow mode - const char* file=""; - // whether a file was specified on the command-line - int exists_file=0; - - if(argc>=3){ - file=argv[2]; - exists_file=1; - } - read_config_file(str_args, file, 1-exists_file); - - return(0); -} - - -// compute the exponential of the input polynomial -int tool_exp(Str_Array str_args){ - // index of the entry in the input file - int arg_index; - // list of fields - Fields_Table fields; - // input polynomial - Polynomial poly; - // exp as a polynomial - Polynomial exp_poly; - // list of rccs - Id_Table idtable; - // exp - Grouped_Polynomial exp; - int i,j; - - // parse fields - arg_index=find_str_arg("fields", str_args); - if(arg_index<0){ - fprintf(stderr,"error: no fields entry in the configuration file\n"); - exit(-1); - } - else{ - parse_input_fields(str_args.strs[arg_index],&fields); - } - - // parse id table - arg_index=find_str_arg("id_table", str_args); - if(arg_index<0){ - fprintf(stderr,"error: no id table entry in the configuration file\n"); - exit(-1); - } - else{ - parse_input_id_table(str_args.strs[arg_index],&idtable, fields); - } - - // parse input polynomial - arg_index=find_str_arg("input_polynomial", str_args); - if(arg_index>=0){ - parse_input_polynomial(str_args.strs[arg_index],&poly, fields); - } - else{ - fprintf(stderr,"error: no input polynomial entry in the configuration file\n"); - exit(-1); - } - - // parse symbols - arg_index=find_str_arg("symbols", str_args); - if(arg_index>=0){ - parse_input_symbols(str_args.strs[arg_index],&fields); - } - else{ - init_Symbols(&(fields.symbols),1); - } - - // parse identities - arg_index=find_str_arg("identities", str_args); - if(arg_index>=0){ - parse_input_identities(str_args.strs[arg_index],&fields); - } - else{ - init_Identities(&(fields.ids),1); - } - - // exp(V) - polynomial_exponential(poly,&exp_poly, fields); - // grouped representation - group_polynomial(exp_poly, &exp, idtable, fields); - free_Polynomial(exp_poly); - free_Polynomial(poly); - - // no denominators - for(i=0;i<exp.length;i++){ - for(j=0;j<exp.coefs[i].length;j++){ - exp.coefs[i].denoms[j].power=0; - } - } - - grouped_polynomial_print(exp,'%','%'); - - // free memory - free_Fields_Table(fields); - free_Id_Table(idtable); - free_Grouped_Polynomial(exp); - return(0); -} diff --git a/src/meantools_expand.c b/src/meantools_expand.c new file mode 100644 index 0000000..2586ff5 --- /dev/null +++ b/src/meantools_expand.c @@ -0,0 +1,166 @@ +/* +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. +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 "meantools_expand.h" + +#include <stdio.h> +#include <stdlib.h> +#include "cli_parser.h" +#include "parse_file.h" +#include "polynomial.h" +#include "array.h" +#include "fields.h" + + +#define CP_FLAG_NAMESPACE 1 +// read command line arguments +int tool_expand_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts){ + // file to read the polynomial from in flow mode + const char* file=""; + // whether a file was specified on the command-line + int exists_file=0; + // flag + int flag=0; + int i; + char* ptr; + + // defaults + // no namespace + (*opts).namespace.length=-1; + + + // loop over arguments + for(i=2;i<argc;i++){ + // flag + if(argv[i][0]=='-'){ + for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){ + switch(*ptr){ + // number of derivatives + case 'N': + flag=CP_FLAG_NAMESPACE; + break; + } + } + } + // namespace + else if(flag==CP_FLAG_NAMESPACE){ + str_to_char_array((char*)argv[i], &(opts->namespace)); + flag=0; + } + // read file name from command-line + else{ + file=argv[i]; + exists_file=1; + } + } + + read_config_file(str_args, file, 1-exists_file); + + return(0); +} + + +// expand a polynomial +int tool_expand(Str_Array str_args, Meantools_Options opts){ + // index of the entry in the input file + int arg_index; + // input polynomial + Polynomial polynomial; + // fields table + Fields_Table fields; + // preprocessor variables + Variables variables; + + // parse fields + if(opts.namespace.length>=0){ + arg_index=find_str_arg_ns("fields", opts.namespace, str_args); + } + else{ + arg_index=find_str_arg("fields", str_args); + } + if(arg_index<0){ + fprintf(stderr,"error: no fields entry in the configuration file\n"); + exit(-1); + } + else{ + parse_input_fields(str_args.strs[arg_index],&fields); + } + + // parse variables + // must precede id_table, virtual_fields, identities and input_polynomial + if(opts.namespace.length>=0){ + arg_index=find_str_arg_ns("preprocessor_variables", opts.namespace, str_args); + } + else{ + arg_index=find_str_arg("preprocessor_variables", str_args); + } + if(arg_index>=0){ + parse_input_variables(str_args.strs[arg_index],&variables); + } + else{ + init_Variables(&variables,1); + } + + // parse virtual_fields + if(opts.namespace.length>=0){ + arg_index=find_str_arg_ns("virtual_fields", opts.namespace, str_args); + } + else{ + arg_index=find_str_arg("virtual_fields", str_args); + } + if(arg_index>=0){ + parse_input_virtual_fields(str_args.strs[arg_index], &fields, variables); + } + + // parse identities + if(opts.namespace.length>=0){ + arg_index=find_str_arg_ns("identities", opts.namespace, str_args); + } + else{ + arg_index=find_str_arg("identities", str_args); + } + if(arg_index>=0){ + parse_input_identities(str_args.strs[arg_index],&fields, variables); + } + + // parse input polynomial + if(opts.namespace.length>=0){ + arg_index=find_str_arg_ns("input_polynomial", opts.namespace, str_args); + } + else{ + arg_index=find_str_arg("input_polynomial", str_args); + } + if(arg_index>=0){ + parse_input_polynomial(str_args.strs[arg_index], &polynomial, fields, variables); + } + else{ + fprintf(stderr,"error: no input polynomial entry in the configuration file\n"); + exit(-1); + } + + // print polynomial + polynomial_print(polynomial); + + // free memory + free_Variables(variables); + free_Fields_Table(fields); + free_Polynomial(polynomial); + + if(opts.namespace.length>=0){ + free_Char_Array(opts.namespace); + } + return(0); +} diff --git a/src/meantools_exp.h b/src/meantools_expand.h index f3f6666..9d988b1 100644 --- a/src/meantools_exp.h +++ b/src/meantools_expand.h @@ -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. @@ -14,14 +14,16 @@ See the License for the specific language governing permissions and limitations under the License. */ -#ifndef MEANTOOLS_EXP_H -#define MEANTOOLS_EXP_H +#ifndef MEANTOOLS_EXPAND_H +#define MEANTOOLS_EXPAND_H #include "types.h" // read arguments -int tool_exp_read_args(int argc, const char* argv[], Str_Array* str_args); -// compute the exponential of the input polynomial -int tool_exp(Str_Array str_args); +int tool_expand_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts); +// expand a flow equation +int tool_expand(Str_Array str_args, Meantools_Options opts); #endif + + diff --git a/src/number.c b/src/number.c index e63427f..5bc87de 100644 --- a/src/number.c +++ b/src/number.c @@ -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. @@ -27,6 +27,7 @@ limitations under the License. #include "tools.h" #include "rational.h" #include "array.h" +#include "parse_file.h" // init int init_Number(Number* number, int memory){ @@ -271,50 +272,107 @@ Number number_Qprod_ret(Q q, Number x){ return(ret); } -// inverse -int number_inverse_inplace(Number* inout){ - int i; - for(i=0;i<(*inout).length;i++){ - if((*inout).base[i]>0){ - (*inout).scalars[i]=Q_inverse((*inout).scalars[i]); - (*inout).scalars[i].denominator*=(*inout).base[i]; +// quotient of two numbers +// recursively simplify numerator/denominator until denominator only has one term +// the output is set to 'numerator' +// both numerator and denominator may be changed by this function +// this algorithm is not optimal in cases where denominator has several terms, in particular if their bases are large +// it is optimal if denominator only has one term +int number_quot_inplace(Number* numerator, Number* denominator){ + Number tmp; + int i,factor; + + switch((*denominator).length){ + // error + case 0: + fprintf(stderr,"error: attempting to invert 0\n"); + exit(-1); + break; + + // trivial case + case 1: + // trivial base + if((*denominator).base[0]==1){ + number_Qprod_chain(Q_inverse((*denominator).scalars[0]), numerator); } - else if((*inout).base[i]<0){ - (*inout).scalars[i]=Q_inverse((*inout).scalars[i]); - (*inout).scalars[i].denominator*=-(*inout).base[i]; - (*inout).scalars[i].numerator*=-1; + // non-trivial base + else{ + // set tmp=1/denominator + tmp=number_one(); + tmp.base[0]=(*denominator).base[0]; + tmp.scalars[0].denominator=(*denominator).scalars[0].numerator*abs((*denominator).base[0]); + if((*denominator).base[0]>0){ + tmp.scalars[0].numerator=(*denominator).scalars[0].denominator; + } + else if((*denominator).base[0]<0){ + tmp.scalars[0].numerator=-(*denominator).scalars[0].denominator; + } + else{ + fprintf(stderr,"error: attempting to invert 0\n"); + exit(-1); + } + number_prod_chain(tmp, numerator); + } + break; + + default: + // find non-trivial basis + for(i=0;(*denominator).base[i]==1;i++){} + + // smallest prime factor of the base + if((*denominator).base[i]<0){ + factor=-1; + } + else if((*denominator).base[i]==2){ + factor=2; } else{ - fprintf(stderr,"error: attempting to invert 0\n"); - exit(-1); + // COMMENT: this is not optimal, but, provided the basis is not too large, this should not be a problem + for(factor=3;is_factor(factor, (*denominator).base[i])==0;factor=factor+2){} + } + + // tmp is set to ((terms that k does not divide)-(terms that k divides)) + // tmp will be multiplied to numerator and denominator + init_Number(&tmp,(*denominator).length); + + // find all terms whose base is a multiple of k + for(i=0;i<(*denominator).length;i++){ + if(is_factor(factor, (*denominator).base[i])){ + // add to tmp with a - sign + number_add_elem(quot(-(*denominator).scalars[i].numerator,(*denominator).scalars[i].denominator), (*denominator).base[i], &tmp); + } + else{ + // add to tmp + number_add_elem((*denominator).scalars[i], (*denominator).base[i], &tmp); + } } + + number_prod_chain(tmp, numerator); + number_prod_chain(tmp, denominator); + free_Number(tmp); + + // recurse + number_quot_inplace(numerator, denominator); } + return(0); } -// write to output -int number_inverse(Number input, Number* output){ - number_cpy(input,output); - number_inverse_inplace(output); - return(0); -} -// return result -Number number_inverse_ret(Number x){ - Number ret; - number_inverse(x,&ret); - return(ret); -} -// quotient +// not inplace int number_quot(Number x1, Number x2, Number* output){ - Number inv; - number_inverse(x2, &inv); - number_prod(x1, inv, output); - free_Number(inv); + Number numerator, denominator; + number_cpy(x1, &numerator); + number_cpy(x2, &denominator); + number_quot_inplace(&numerator, &denominator); + *output=numerator; + free_Number(denominator); return(0); } -int number_quot_chain(Number x1, Number* inout){ - number_inverse_inplace(inout); - number_prod_chain(x1, inout); +int number_quot_chain(Number* inout, Number x2){ + Number tmp; + number_quot(*inout,x2,&tmp); + free_Number(*inout); + *inout=tmp; return(0); } Number number_quot_ret(Number x1, Number x2){ @@ -419,7 +477,7 @@ int number_print(Number number){ Char_Array buffer; init_Char_Array(&buffer,5*number.length); number_sprint(number, &buffer); - printf("%s",buffer.str); + printf("%s",char_array_to_str_noinit(&buffer)); return(0); } @@ -434,18 +492,54 @@ int str_to_Number(char* str, Number* number){ char* buffer_ptr=buffer; Q num; int base; - // whether there are parentheses in the string - int exist_parenthesis=0; + int ret; + int j; + char* aux_str; + int aux_free=0; init_Number(number, NUMBER_SIZE); + // check whether the string is blank (return 0 in that case) + for(ptr=str;*ptr!='\0';ptr++){ + if(*ptr!=' ' && *ptr!='\n'){ + break; + } + } + // blank string + if(*ptr=='\0'){ + number_add_elem(quot(0,1), 1, number); + free(buffer); + return(0); + } + // init num and base - // init to 0 so that if str is empty, then the number is set to 0 - num=quot(0,1); + num=quot(1,1); base=1; + // check whether the str only contains a rational number, and add parentheses + // keep rtack of the length of str + for(j=0,ptr=str;*ptr!='\0';j++,ptr++){ + if((*ptr<'0' || *ptr>'9') && *ptr!='-' && *ptr!='/'){ + break; + } + } + // only rational + if(*ptr=='\0'){ + aux_str=calloc(j+3,sizeof(char)); + aux_str[0]='('; + for(j=0,ptr=str;*ptr!='\0';ptr++,j++){ + aux_str[j+1]=*ptr; + } + aux_str[j+1]=')'; + aux_str[j+2]='\0'; + aux_free=1; + } + else{ + aux_str=str; + } + mode=PP_NULL_MODE; - for(ptr=str;*ptr!='\0';ptr++){ + for(ptr=aux_str;*ptr!='\0';ptr++){ switch(*ptr){ // read number case '(': @@ -453,7 +547,10 @@ int str_to_Number(char* str, Number* number){ // init base base=1; mode=PP_NUM_MODE; - exist_parenthesis=1; + } + else{ + fprintf(stderr,"syntax error: misplaced '(' in number '%s'\n",str); + exit(-1); } break; case ')': @@ -463,6 +560,10 @@ int str_to_Number(char* str, Number* number){ *buffer_ptr='\0'; mode=PP_NULL_MODE; } + else{ + fprintf(stderr,"syntax error: mismatched ')' in number '%s'\n",str); + exit(-1); + } break; // read sqrt @@ -474,16 +575,26 @@ int str_to_Number(char* str, Number* number){ if(mode==PP_NULL_MODE){ mode=PP_SQRT_MODE; } - // if there is a square root, then do not read a fraction (see end of loop) - exist_parenthesis=1; + else{ + fprintf(stderr,"syntax error: misplaced '{' in number '%s'\n",str); + exit(-1); + } break; case '}': if(mode==PP_SQRT_MODE){ - sscanf(buffer,"%d",&base); + ret=read_int(buffer,&base); + if(ret<0){ + fprintf(stderr,"syntax error: number base should be an integer, got '%s' in '%s'",buffer,str); + exit(-1); + } buffer_ptr=buffer; *buffer_ptr='\0'; mode=PP_NULL_MODE; } + else{ + fprintf(stderr,"syntax error: mismatched '}' in number '%s'\n",str); + exit(-1); + } break; // write num @@ -494,24 +605,40 @@ int str_to_Number(char* str, Number* number){ num=quot(0,1); base=1; } + else{ + fprintf(stderr,"syntax error: misplaced '+' in number '%s'\n",str); + exit(-1); + } break; + // ignore 's', ' ' and '\n' + case 's':break; + case ' ':break; + case '\n':break; default: if(mode!=PP_NULL_MODE){ buffer_ptr=str_addchar(buffer_ptr,*ptr); } + else{ + fprintf(stderr,"syntax error: unrecognized character '%c' in number '%s'\n",*ptr,str); + exit(-1); + } break; } } - // last step if(mode==PP_NULL_MODE){ - if(exist_parenthesis==0){ - str_to_Q(str, &num); - } number_add_elem(num, base, number); } + else{ + fprintf(stderr,"syntax error: mismatched '(' in number '%s'\n",str); + exit(-1); + } + + if(aux_free==1){ + free(aux_str); + } free(buffer); return(0); } diff --git a/src/number.h b/src/number.h index a33d425..5eacdaa 100644 --- a/src/number.h +++ b/src/number.h @@ -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. @@ -78,16 +78,11 @@ int number_Qprod(Q q, Number x, Number* inout); // return result Number number_Qprod_ret(Q q, Number x); -// inverse -int number_inverse_inplace(Number* inout); -// write to output -int number_inverse(Number input, Number* output); -// return result -Number number_inverse_ret(Number x); - -// quotient +// quotient of two numbers +int number_quot_inplace(Number* numerator, Number* denominator); +// not inplace int number_quot(Number x1, Number x2, Number* output); -int number_quot_chain(Number x1, Number* inout); +int number_quot_chain(Number* inout, Number x2); Number number_quot_ret(Number x1, Number x2); // remove 0's diff --git a/src/numkondo.c b/src/numkondo.c index 0568861..91e6be2 100644 --- a/src/numkondo.c +++ b/src/numkondo.c @@ -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. @@ -86,12 +86,6 @@ int read_args_numkondo(int argc,const char* argv[], Str_Array* str_args, Numkond // whether a file was specified on the command-line int exists_file=0; - // if there are no arguments - if(argc==1){ - print_usage_numkondo(); - exit(-1); - } - // defaults // display entire flow (*opts).display_mode=DISPLAY_NUMERICAL; @@ -193,6 +187,8 @@ int numflow(Str_Array str_args, Numkondo_Options opts){ Grouped_Polynomial flow_equation; // whether or not to use mpfr floats int mpfr_flag=0; + // postprocess flow equation + Grouped_Polynomial postprocess_flow_equation; // set mpfr defaults if(opts.mpfr_prec!=0){ @@ -205,7 +201,7 @@ int numflow(Str_Array str_args, Numkondo_Options opts){ } - // parse id table + // parse labels arg_index=find_str_arg("labels", str_args); if(arg_index<0){ fprintf(stderr,"error: no labels entry in the configuration file\n"); @@ -225,6 +221,16 @@ int numflow(Str_Array str_args, Numkondo_Options opts){ char_array_to_Grouped_Polynomial(str_args.strs[arg_index], &flow_equation); } + // parse postprocess operation + arg_index=find_str_arg("postprocess_operation", str_args); + if(arg_index>=0){ + char_array_to_Grouped_Polynomial(str_args.strs[arg_index], &postprocess_flow_equation); + } + else{ + init_Grouped_Polynomial(&postprocess_flow_equation,1); + } + + // initial conditions // check they were not specified on the command line if(opts.eval_rccstring.length==-1){ @@ -252,17 +258,18 @@ int numflow(Str_Array str_args, Numkondo_Options opts){ } if(mpfr_flag==0){ - numerical_flow(flow_equation, init_cd, labels, opts.niter, opts.display_mode); + numerical_flow(flow_equation, init_cd, postprocess_flow_equation, labels, opts.niter, opts.display_mode); free_RCC(init_cd); } else{ - numerical_flow_mpfr(flow_equation, init_cd_mpfr, labels, opts.niter, opts.display_mode); + numerical_flow_mpfr(flow_equation, init_cd_mpfr, postprocess_flow_equation, labels, opts.niter, opts.display_mode); free_RCC_mpfr(init_cd_mpfr); } // free memory free_Labels(labels); + free_Grouped_Polynomial(postprocess_flow_equation); free_Grouped_Polynomial(flow_equation); return(0); } diff --git a/src/parse_file.c b/src/parse_file.c index ae3a9af..e998803 100644 --- a/src/parse_file.c +++ b/src/parse_file.c @@ -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,8 @@ limitations under the License. #include "istring.h" #include "tools.h" #include "idtable.h" +#include "tree.h" +#include "symbolic_parser.h" // parsing modes @@ -59,6 +61,70 @@ limitations under the License. #define PP_GROUP_MODE 16 +// read a positive integer from a string +int read_positive_int(char* str, int* out){ + char* ptr; + int res; + for(ptr=str;*ptr!='\0';ptr++){ + if(*ptr-'0'>=10 || *ptr-'0'<0){ + return(-1); + } + } + res=sscanf(str,"%d",out); + if(res!=1){ + return(-2); + } + return(0); +} +// read an integer from a string +int read_int(char* str, int* out){ + char* ptr; + int res; + for(ptr=str;*ptr!='\0';ptr++){ + if((*ptr-'0'>=10 || *ptr-'0'<0) && (*ptr!='-')){ + return(-1); + } + } + res=sscanf(str,"%d",out); + if(res!=1){ + return(-2); + } + return(0); +} +// read an long int from a string +int read_long_int(char* str, long int* out){ + char* ptr; + int res; + for(ptr=str;*ptr!='\0';ptr++){ + if((*ptr-'0'>=10 || *ptr-'0'<0) && (*ptr!='-')){ + return(-1); + } + } + res=sscanf(str,"%ld",out); + if(res!=1){ + return(-2); + } + return(0); +} + +// read a long double +int read_long_double(char* str, long double* out){ + char* ptr; + int res; + for(ptr=str;*ptr!='\0';ptr++){ + if((*ptr-'0'>=10 || *ptr-'0'<0) && (*ptr!='-') && (*ptr!='e') && (*ptr!='E') && (*ptr!='.')){ + return(-1); + } + } + res=sscanf(str,"%Lf",out); + if(res!=1){ + return(-2); + } + return(0); +} + + + // parse fields list int parse_input_fields(Char_Array str_fields, Fields_Table* fields){ // buffer @@ -67,6 +133,8 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){ int i,j; int mode; int comment=0; + int res; + int expect_colon=0; // allocate memory init_Fields_Table(fields); @@ -86,39 +154,100 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){ if(mode==PP_NULL_MODE){ mode=PP_PARAMETER_MODE; } + else{ + fprintf(stderr,"syntax error: fields entry (%d): 'h' should be at the beginning of the line\n",j); + exit(-1); + } + // flag: colon is expected immediately after letter + expect_colon=1; break; // external fields case 'x': if(mode==PP_NULL_MODE){ mode=PP_EXTERNAL_MODE; } + else{ + fprintf(stderr,"syntax error: fields entry (%d): 'x' should be at the beginning of the line\n",j); + exit(-1); + } + // flag: colon is expected immediately after letter + expect_colon=1; break; // internal fields case 'i': if(mode==PP_NULL_MODE){ mode=PP_INTERNAL_MODE; } + else{ + fprintf(stderr,"syntax error: fields entry (%d): 'i' should be at the beginning of the line\n",j); + exit(-1); + } + // flag: colon is expected immediately after letter + expect_colon=1; break; case 'f': if(mode==PP_NULL_MODE){ mode=PP_FERMIONS_MODE; } + else{ + fprintf(stderr,"syntax error: fields entry (%d): 'f' should be at the beginning of the line\n",j); + exit(-1); + } + // flag: colon is expected immediately after letter + expect_colon=1; break; case 'a': if(mode==PP_NULL_MODE){ mode=PP_NONCOMMUTING_MODE; } + else{ + fprintf(stderr,"syntax error: fields entry (%d): 'a' should be at the beginning of the line\n",j); + exit(-1); + } + // flag: colon is expected immediately after letter + expect_colon=1; break; // reset buffer case ':': + // check mode + if(mode==PP_NULL_MODE || expect_colon==0){ + fprintf(stderr,"syntax error: fields entry (%d): ':' must be preceded by 'h', 'x', 'i', 'f', or 'a' on the same line\n",j); + exit(-1); + } + expect_colon=0; buffer_ptr=buffer; *buffer_ptr='\0'; break; // write to fields case ',': - sscanf(buffer,"%d",&i); + res=read_positive_int(buffer,&i); + + // check i + if(res<0){ + fprintf(stderr,"syntax error: fields entry (%d): fields must be non-negative integers, got '%s'\n",j,buffer); + exit(-1); + } + if(mode==PP_PARAMETER_MODE || mode==PP_EXTERNAL_MODE || mode==PP_INTERNAL_MODE){ + if(int_array_find(i,fields->parameter)>=0 || int_array_find(i,fields->external)>=0 || int_array_find(i,fields->internal)>=0){ + fprintf(stderr,"error: fields entry (%d): field %d is already present in the fields table\n",j,i); + exit(-1); + } + } + else if(mode==PP_FERMIONS_MODE){ + if(int_array_find(i,fields->fermions)>=0){ + fprintf(stderr,"error: fields entry (%d): field %d is already present in the Fermions table\n",j,i); + exit(-1); + } + } + else if(mode==PP_NONCOMMUTING_MODE){ + if(int_array_find(i,fields->noncommuting)>=0){ + fprintf(stderr,"error: fields entry (%d): field %d is already present in the noncommuting table\n",j,i); + exit(-1); + } + } + if(mode==PP_PARAMETER_MODE){ int_array_append(i,&((*fields).parameter)); } @@ -134,31 +263,65 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){ else if(mode==PP_NONCOMMUTING_MODE){ int_array_append(i,&((*fields).noncommuting)); } + else{ + fprintf(stderr,"syntax error: fields entry (%d): misplaced ','\n",j); + exit(-1); + } buffer_ptr=buffer; *buffer_ptr='\0'; break; // back to null mode case '\n': - sscanf(buffer,"%d",&i); - if(mode==PP_PARAMETER_MODE){ - int_array_append(i,&((*fields).parameter)); - } - else if(mode==PP_EXTERNAL_MODE){ - int_array_append(i,&((*fields).external)); - } - else if(mode==PP_INTERNAL_MODE){ - int_array_append(i,&((*fields).internal)); - } - else if(mode==PP_FERMIONS_MODE){ - int_array_append(i,&((*fields).fermions)); - } - else if(mode==PP_NONCOMMUTING_MODE){ - int_array_append(i,&((*fields).noncommuting)); + // ignore if in NULL_MODE + if(mode!=PP_NULL_MODE){ + res=read_positive_int(buffer,&i); + // check i + if(res<0){ + fprintf(stderr,"syntax error: fields entry (%d): fields must be non-negative integers, got '%s'\n",j,buffer); + exit(-1); + } + if(mode==PP_PARAMETER_MODE || mode==PP_EXTERNAL_MODE || mode==PP_INTERNAL_MODE){ + if(int_array_find(i,fields->parameter)>=0 || int_array_find(i,fields->external)>=0 || int_array_find(i,fields->internal)>=0){ + fprintf(stderr,"error: fields entry (%d): field %d is already present in the fields table\n",j,i); + exit(-1); + } + } + else if(mode==PP_FERMIONS_MODE){ + if(int_array_find(i,fields->fermions)>=0){ + fprintf(stderr,"error: fields entry (%d): field %d is already present in the Fermions table\n",j,i); + exit(-1); + } + } + else if(mode==PP_NONCOMMUTING_MODE){ + if(int_array_find(i,fields->noncommuting)>=0){ + fprintf(stderr,"error: fields entry (%d): field %d is already present in the noncommuting table\n",j,i); + exit(-1); + } + } + + if(mode==PP_PARAMETER_MODE){ + int_array_append(i,&((*fields).parameter)); + } + else if(mode==PP_EXTERNAL_MODE){ + int_array_append(i,&((*fields).external)); + } + else if(mode==PP_INTERNAL_MODE){ + int_array_append(i,&((*fields).internal)); + } + else if(mode==PP_FERMIONS_MODE){ + int_array_append(i,&((*fields).fermions)); + } + else if(mode==PP_NONCOMMUTING_MODE){ + int_array_append(i,&((*fields).noncommuting)); + } + mode=PP_NULL_MODE; } - mode=PP_NULL_MODE; break; + + // ignore spaces + case ' ':break; // comment case '#': comment=1; @@ -168,57 +331,85 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){ if(mode!=PP_NULL_MODE){ buffer_ptr=str_addchar(buffer_ptr,str_fields.str[j]); } + else{ + fprintf(stderr,"syntax error: fields entry (%d): unrecognized character '%c'\n",j,str_fields.str[j]); + exit(-1); + } + break; + break; } } } free(buffer); + + // check whether there are non-Fermionic internal fields + for(i=0;i<(*fields).internal.length;i++){ + if(is_fermion((*fields).internal.values[i], *fields)==0){ + fprintf(stderr,"warning: some internal fields are not Fermions: means may be computed using sums over permutations rather than the more efficient LU decomposition\n"); + break; + } + } + return(0); } -// parse symbols list +// parse virtual_fields list // write result to fields -int parse_input_symbols(Char_Array str_symbols, Fields_Table* fields){ +int parse_input_virtual_fields(Char_Array str_virtual_fields, Fields_Table* fields, Variables variables){ // buffer - char* buffer=calloc(str_symbols.length+1,sizeof(char)); + char* buffer=calloc(str_virtual_fields.length+1,sizeof(char)); char* buffer_ptr=buffer; Polynomial polynomial; int index; int i,j; int mode; int comment=0; + int res; // loop over input mode=PP_INDEX_MODE; - for(j=0;j<str_symbols.length;j++){ + for(j=0;j<str_virtual_fields.length;j++){ if(comment==1){ - if(str_symbols.str[j]=='\n'){ + if(str_virtual_fields.str[j]=='\n'){ comment=0; } } // stay in polynomial mode until ',' else if(mode==PP_POLYNOMIAL_MODE){ - if(str_symbols.str[j]==','){ + if(str_virtual_fields.str[j]==','){ // parse polynomial - str_to_Polynomial(buffer, &polynomial); + parse_symbolic_expression_str(buffer, *fields, variables, &polynomial); // write index and polynomial - symbols_append_noinit(index, polynomial, &((*fields).symbols)); + virtual_fields_append_noinit(index, polynomial, &((*fields).virtual_fields)); mode=PP_INDEX_MODE; buffer_ptr=buffer; *buffer_ptr='\0'; } + // check there are no '=' signs in polynomial (forgotten ',') + else if(str_virtual_fields.str[j]=='='){ + fprintf(stderr,"syntax error: virtual fields entry (%d): '=' in polynomial '%s'\n",j, buffer); + exit(-1); + } + else if(str_virtual_fields.str[j]=='#'){ + comment=1; + } else{ - buffer_ptr=str_addchar(buffer_ptr,str_symbols.str[j]); + buffer_ptr=str_addchar(buffer_ptr,str_virtual_fields.str[j]); } } else{ - switch(str_symbols.str[j]){ + switch(str_virtual_fields.str[j]){ // polynomial mode case '=': if(mode==PP_INDEX_MODE){ // read index - sscanf(buffer,"%d",&index); + res=read_positive_int(buffer, &index); + if(res<0){ + fprintf(stderr,"syntax error: virtual fields entry (%d): virtual field index must be a non-negative integer, got '%s'\n",j,buffer); + exit(-1); + } // reset buffer buffer_ptr=buffer; *buffer_ptr='\0'; @@ -226,43 +417,53 @@ int parse_input_symbols(Char_Array str_symbols, Fields_Table* fields){ } break; + // misplaced ',' + case ',': + fprintf(stderr,"syntax error: virtual fields entry (%d): misplaced ',' or forgotten virtual field index\n",j); + exit(-1); + break; + + // ignore spaces and newlines + case ' ':break; + case '\n':break; + // comment case '#': comment=1; break; default: - buffer_ptr=str_addchar(buffer_ptr,str_symbols.str[j]); + buffer_ptr=str_addchar(buffer_ptr,str_virtual_fields.str[j]); break; } } } // last step - if(polynomial.length>0){ - str_to_Polynomial(buffer, &polynomial); - symbols_append_noinit(index, polynomial, &((*fields).symbols)); - } + parse_symbolic_expression_str(buffer, *fields, variables, &polynomial); + virtual_fields_append_noinit(index, polynomial, &((*fields).virtual_fields)); // simplify - for(i=0;i<(*fields).symbols.length;i++){ - polynomial_simplify((*fields).symbols.expr+i, *fields); + for(i=0;i<(*fields).virtual_fields.length;i++){ + polynomial_simplify((*fields).virtual_fields.expr+i, *fields); } free(buffer); return(0); } + // parse groups of independent fields -int parse_input_groups(Char_Array str_groups, Groups* groups){ +int parse_input_groups(Char_Array str_groups, Groups* groups, Polynomial_Matrix propagator, Fields_Table fields){ // buffer char* buffer=calloc(str_groups.length+1,sizeof(char)); char* buffer_ptr=buffer; int index; - int j; + int i,j; Int_Array group; int mode; int comment=0; + int res; // alloc init_Groups(groups, GROUP_SIZE); @@ -287,26 +488,60 @@ int parse_input_groups(Char_Array str_groups, Groups* groups){ init_Int_Array(&group, GROUP_SIZE); mode=PP_GROUP_MODE; } + else{ + fprintf(stderr,"syntax error: groups entry (%d): nested parentheses are not allowed in this entry\n",j); + exit(-1); + } break; case')': if(mode==PP_GROUP_MODE){ - sscanf(buffer,"%d",&index); + res=read_positive_int(buffer,&index); + if(res<0){ + fprintf(stderr,"syntax error: groups entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer); + exit(-1); + } + for(i=0;i<groups->length;i++){ + if(int_array_find(index,groups->indices[i])>=0){ + fprintf(stderr,"syntax error: groups entry (%d): index %d is already present in group %d\n",j,index,i); + exit(-1); + } + } + int_array_append(index, &group); groups_append_noinit(group, groups); mode=PP_NULL_MODE; } + else{ + fprintf(stderr,"syntax error: groups entry (%d): extra ')'\n",j); + exit(-1); + } break; // read index case',': if(mode==PP_GROUP_MODE){ - sscanf(buffer,"%d",&index); + res=read_positive_int(buffer,&index); + if(res<0){ + fprintf(stderr,"syntax errorL groups entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer); + exit(-1); + } + for(i=0;i<groups->length;i++){ + if(int_array_find(index,groups->indices[i])>=0){ + fprintf(stderr,"syntax error: groups entry (%d): index %d is already present in group %d\n",j,index,i); + exit(-1); + } + } + int_array_append(index, &group); buffer_ptr=buffer; *buffer_ptr='\0'; } break; + // skip spaces and newlines + case ' ':break; + case '\n':break; + // comment case '#': comment=1; @@ -316,11 +551,204 @@ int parse_input_groups(Char_Array str_groups, Groups* groups){ if(mode!=PP_NULL_MODE){ buffer_ptr=str_addchar(buffer_ptr,str_groups.str[j]); } + else{ + fprintf(stderr,"syntax error: groups entry (%d): unrecognized character '%c'\n",j,str_groups.str[j]); + exit(-1); + } break; } } } + // check fields in groups + check_groups(*groups, propagator, fields); + + free(buffer); + return(0); +} + +// check that the members of groups are independent (assuming the virtual fields and propagator were already parsed) +int check_groups(Groups groups, Polynomial_Matrix propagator, Fields_Table fields){ + Int_Array* group_fields=calloc(groups.length,sizeof(Int_Array)); + int i,j,k,l,a,b; + + // get the lists of fields involved in each group + for(i=0;i<groups.length;i++){ + // expand virtual_fields + fields_in_virtual_field_list(groups.indices[i], fields, group_fields+i); + } + + // check whether pairs of fields are present in the propagator + for(i=0;i<groups.length;i++){ + for(j=0;j<group_fields[i].length;j++){ + if(field_type(group_fields[i].values[j], fields)==FIELD_INTERNAL){ + for(k=0;k<groups.length;k++){ + if(k!=i){ + for(l=0;l<group_fields[k].length;l++){ + if(field_type(group_fields[k].values[l], fields)==FIELD_INTERNAL){ + if(group_fields[i].values[j]==group_fields[k].values[l]){ + fprintf(stderr,"error: groups %d and %d both contain internal field %d\n",i,k,group_fields[i].values[j]); + exit(-1); + } + a=intlist_find(propagator.indices, propagator.length, group_fields[i].values[j]); + b=intlist_find(propagator.indices, propagator.length, group_fields[k].values[l]); + if(a>=0 && b>=0 && polynomial_is_zero(propagator.matrix[a][b])==0){ + fprintf(stderr,"error: groups %d and %d are not independent: they contain %d and %d which have a non-vanishing propagator\n",i,k,group_fields[i].values[j], group_fields[k].values[l]); + exit(-1); + } + } + } + } + } + } + } + } + + for(i=0;i<groups.length;i++){ + free_Int_Array(group_fields[i]); + } + free(group_fields); + + return(0); +} + +// list of fields involved in a list of virtual_fields +int fields_in_virtual_field_list(Int_Array indices, Fields_Table fields, Int_Array* output){ + int i,j,k; + Int_Array tmp_array; + + init_Int_Array(output,FIELDS_SIZE); + for(k=0;k<indices.length;k++){ + if(field_type(indices.values[k],fields)!=FIELD_VIRTUAL){ + int_array_append(abs(indices.values[k]), output); + } + else{ + i=intlist_find_err(fields.virtual_fields.indices, fields.virtual_fields.length, indices.values[k]); + for(j=0;j<fields.virtual_fields.expr[i].length;j++){ + fields_in_virtual_field_list(fields.virtual_fields.expr[i].monomials[j], fields, &tmp_array); + int_array_concat_unique(tmp_array, output); + free_Int_Array(tmp_array); + } + } + } + + return(0); +} + + + +// parse variables list +int parse_input_variables(Char_Array str_variables, Variables* variables){ + // buffer + char* buffer=calloc(str_variables.length+1,sizeof(char)); + char* buffer_ptr=buffer; + Tree symbol_tree; + Char_Array varname; + int j; + int mode; + int comment=0; + char* ptr; + + init_Variables(variables,VARIABLES_SIZE); + + // loop over input + mode=PP_INDEX_MODE; + for(j=0;j<str_variables.length;j++){ + if(comment==1){ + if(str_variables.str[j]=='\n'){ + comment=0; + } + } + // stay in polynomial mode until ',' + else if(mode==PP_POLYNOMIAL_MODE){ + if(str_variables.str[j]==','){ + // parse symbol tree + str_to_symbol_tree(buffer, &symbol_tree); + variables_append_noinit(varname , symbol_tree, variables); + mode=PP_INDEX_MODE; + buffer_ptr=buffer; + *buffer_ptr='\0'; + } + // check there are no '=' signs in polynomial (forgotten ',') + else if(str_variables.str[j]=='='){ + fprintf(stderr,"syntax error: preprocessor_variables entry (%d): '=' in polynomial '%s'\n",j, buffer); + exit(-1); + } + else if(str_variables.str[j]=='#'){ + comment=1; + } + else{ + buffer_ptr=str_addchar(buffer_ptr,str_variables.str[j]); + } + } + else{ + switch(str_variables.str[j]){ + // polynomial mode + case '=': + if(mode==PP_INDEX_MODE){ + if(*buffer=='\0'){ + fprintf(stderr,"syntax error: preprocessor_variables entry (%d): empty variable name\n",j); + exit(-1); + } + for(ptr=buffer;*ptr!='\0';ptr++){ + if(*ptr=='<' || *ptr=='>'){ + fprintf(stderr,"syntax error: preprocessor_variables entry (%d): forbidden '%c' character in variable name '%s'\n",j,*ptr,buffer); + exit(-1); + } + } + str_to_char_array(buffer,&varname); + // check that the variable name is not 'OUT' or 'FLOW' or 'RCC' + if(char_array_cmp_str(varname, "OUT")==1){ + fprintf(stderr,"error: preprocessor_variables entry (%d): variable name cannot be 'OUT' (this name is reserved)\n",j); + exit(-1); + } + if(char_array_cmp_str(varname, "FLOW")==1){ + fprintf(stderr,"error: preprocessor_variables entry (%d): variable name cannot be 'FLOW' (this name is reserved)\n",j); + exit(-1); + } + if(char_array_cmp_str(varname, "RCC")==1){ + fprintf(stderr,"error: preprocessor_variables entry (%d): variable name cannot be 'RCC' (this name is reserved)\n",j); + exit(-1); + } + // reset buffer + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_POLYNOMIAL_MODE; + } + break; + + // misplaced ',' + case ',': + fprintf(stderr,"syntax error: preprocessor_variables entry (%d): misplaced ',' or forgotten variable name\n",j); + exit(-1); + break; + + // ignore characters + case ' ':break; + case '\n': break; + + // comment + case '#': + comment=1; + break; + + default: + buffer_ptr=str_addchar(buffer_ptr,str_variables.str[j]); + break; + } + } + } + + // last step + if(mode==PP_POLYNOMIAL_MODE){ + str_to_symbol_tree(buffer, &symbol_tree); + variables_append_noinit(varname , symbol_tree, variables); + } + else if(*buffer!='\0'){ + fprintf(stderr,"syntax error: preprocessor_variables entry (%d): mismatched variable/polynomial pair, got '%s'\n",j,buffer); + exit(-1); + } + free(buffer); return(0); } @@ -328,7 +756,7 @@ int parse_input_groups(Char_Array str_groups, Groups* groups){ // parse identities between fields // write result to fields -int parse_input_identities(Char_Array str_identities, Fields_Table* fields){ +int parse_input_identities(Char_Array str_identities, Fields_Table* fields, Variables variables){ // buffer char* buffer=calloc(str_identities.length+1,sizeof(char)); char* buffer_ptr=buffer; @@ -339,6 +767,7 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){ int tmp; int mode; int comment=0; + int ret; init_Int_Array(&monomial, MONOMIAL_SIZE); @@ -354,13 +783,21 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){ else if(mode==PP_POLYNOMIAL_MODE){ if(str_identities.str[j]==','){ // parse polynomial - str_to_Polynomial(buffer, &polynomial); + parse_symbolic_expression_str(buffer, *fields, variables, &polynomial); // write monomial and polynomial identities_append_noinit(monomial, polynomial, &((*fields).ids)); // realloc init_Int_Array(&monomial, MONOMIAL_SIZE); mode=PP_NULL_MODE; } + // check there are no '=' signs in polynomial (forgotten ',') + else if(str_identities.str[j]=='='){ + fprintf(stderr,"syntax error: identities entry (%d): '=' in polynomial '%s'\n",j, buffer); + exit(-1); + } + else if(str_identities.str[j]=='#'){ + comment=1; + } else{ buffer_ptr=str_addchar(buffer_ptr,str_identities.str[j]); } @@ -377,20 +814,36 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){ buffer_ptr=buffer; *buffer_ptr='\0'; } + else{ + fprintf(stderr,"syntax error: identities entry (%d): misplaced '['\n",j); + exit(-1); + } break; // for notational homogeneity case 'f': if(mode==PP_BRACKET_MODE){ mode=PP_FIELD_MODE; } + else{ + fprintf(stderr,"syntax error: identities entry (%d): misplaced 'f'\n",j); + exit(-1); + } break; // write monomial term case ']': if(mode==PP_FIELD_MODE){ - sscanf(buffer,"%d",&i); + ret=read_int(buffer,&i); + if(ret<0){ + fprintf(stderr,"syntax error: identities entry (%d): field indices must be integers, got '%s'\n",j,buffer); + exit(-1); + } int_array_append(i,&monomial); mode=PP_INDEX_MODE; } + else{ + fprintf(stderr,"syntax error: identities entry (%d): mismatched ']'\n",j); + exit(-1); + } break; // polynomial mode @@ -400,8 +853,22 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){ *buffer_ptr='\0'; mode=PP_POLYNOMIAL_MODE; } + else{ + fprintf(stderr,"syntax error: identities entry (%d): misplaced '='\n",j); + exit(-1); + } break; + // misplaced ',' + case ',': + fprintf(stderr,"syntax error: identities entry (%d): misplaced ',' or forgotten left hand side\n",j); + exit(-1); + break; + + // ignore spaces and newlines + case ' ':break; + case '\n':break; + // comment case '#': comment=1; @@ -411,6 +878,10 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){ if(mode!=PP_NULL_MODE){ buffer_ptr=str_addchar(buffer_ptr,str_identities.str[j]); } + else{ + fprintf(stderr,"syntax error: identities entry (%d): unrecognized character '%c'\n",j,str_identities.str[j]); + exit(-1); + } break; } } @@ -418,9 +889,13 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){ // last step if(mode==PP_POLYNOMIAL_MODE){ - str_to_Polynomial(buffer, &polynomial); + parse_symbolic_expression_str(buffer, *fields, variables, &polynomial); identities_append_noinit(monomial, polynomial, &((*fields).ids)); } + else if(*buffer!='\0'){ + fprintf(stderr,"syntax error: identities entry (%d): mismatched left/right hand side, got '%s'\n",j,buffer); + exit(-1); + } else{ free_Int_Array(monomial); } @@ -451,6 +926,7 @@ int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propaga int index2=-1; int mode; int comment=0; + int ret; // allocate memory init_Polynomial_Matrix(propagator, fields.internal.length); @@ -473,20 +949,36 @@ int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propaga // indices case ';': if(mode==PP_INDEX_MODE){ - sscanf(buffer,"%d",&i); + ret=read_positive_int(buffer,&i); + if(ret<0){ + fprintf(stderr,"syntax error: propagator entry (%d): indices must be non-negative integers, got '%s'\n",j,buffer); + exit(-1); + } index1=intlist_find_err((*propagator).indices, (*propagator).length, i); buffer_ptr=buffer; *buffer_ptr='\0'; } + else{ + fprintf(stderr,"syntax error: propagator entry (%d): misplaced ';'\n",j); + exit(-1); + } break; case ':': if(mode==PP_INDEX_MODE){ - sscanf(buffer,"%d",&i); + ret=read_positive_int(buffer,&i); + if(ret<0){ + fprintf(stderr,"syntax error: propagator entry (%d): indices must be non-negative integers, got '%s'\n",j,buffer); + exit(-1); + } index2=intlist_find_err((*propagator).indices, (*propagator).length, i); buffer_ptr=buffer; *buffer_ptr='\0'; mode=PP_POLYNOMIAL_MODE; } + else{ + fprintf(stderr,"syntax error: propagator entry (%d): misplaced ':'\n",j); + exit(-1); + } break; // num @@ -494,12 +986,22 @@ int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propaga if(mode==PP_POLYNOMIAL_MODE && index1>=0 && index2>=0){ free_Polynomial((*propagator).matrix[index1][index2]); str_to_Polynomial(buffer,(*propagator).matrix[index1]+index2); + index1=-1; + index2=-1; buffer_ptr=buffer; *buffer_ptr='\0'; mode=PP_INDEX_MODE; } + else{ + fprintf(stderr,"syntax error: propagator entry (%d): mismatched index/polynomial pair\n",j); + exit(-1); + } break; + // ignore ' 'and newline + case ' ':break; + case '\n':break; + // comment case '#': comment=1; @@ -517,6 +1019,15 @@ int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propaga free_Polynomial((*propagator).matrix[index1][index2]); str_to_Polynomial(buffer,(*propagator).matrix[index1]+index2); } + else if (*buffer!='\0'){ + fprintf(stderr,"syntax error: propagator entry (%d): mismatched index/polynomial pair, got '%s'\n",j,buffer); + exit(-1); + } + + // check whether the polynomial is numeric + if(polynomial_matrix_is_numeric(*propagator)==0){ + fprintf(stderr,"warning: the propagator contains polynomial entries: means may be computed using sums over permutations rather than the more efficient LU decomposition\n"); + } free(buffer); return(0); @@ -524,49 +1035,9 @@ int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propaga // parse input polynomial -int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields){ - int j; - // buffer - char* buffer=calloc(str_polynomial.length+1,sizeof(char)); - char* buffer_ptr=buffer; - Polynomial tmp_poly; - - // allocate memory - init_Polynomial(output,POLY_SIZE); - - for(j=0;j<str_polynomial.length;j++){ - switch(str_polynomial.str[j]){ - case '*': - str_to_Polynomial(buffer, &tmp_poly); - if((*output).length==0){ - polynomial_concat(tmp_poly, output); - } - else{ - polynomial_prod_chain(tmp_poly, output, fields); - } - free_Polynomial(tmp_poly); - - buffer_ptr=buffer; - *buffer_ptr='\0'; - break; - - default: - buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]); - break; - } - } - - //last step - str_to_Polynomial(buffer, &tmp_poly); - if((*output).length==0){ - polynomial_concat(tmp_poly, output); - } - else{ - polynomial_prod_chain(tmp_poly, output, fields); - } - free_Polynomial(tmp_poly); - - free(buffer); +int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields, Variables variables){ + parse_symbolic_expression(str_polynomial, fields, variables, output); + polynomial_simplify(output, fields); return(0); } @@ -574,7 +1045,7 @@ int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields // parse id table // fields argument for sorting -int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields){ +int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields, Variables variables){ // buffer char* buffer=calloc(str_idtable.length+1,sizeof(char)); char* buffer_ptr=buffer; @@ -583,6 +1054,7 @@ int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table int j; int mode; int comment=0; + int ret; // allocate memory init_Id_Table(idtable,EQUATION_SIZE); @@ -601,24 +1073,40 @@ int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table case ',': // write polynomial if(mode==PP_POLYNOMIAL_MODE){ - str_to_Polynomial(buffer,&polynomial); + parse_symbolic_expression_str(buffer, fields, variables, &polynomial); // add to idtable idtable_append_noinit(index,polynomial,idtable); mode=PP_INDEX_MODE; buffer_ptr=buffer; *buffer_ptr='\0'; } + else{ + fprintf(stderr,"syntax error: id_table entry (%d): misplaced ','\n",j); + exit(-1); + } break; case ':': if(mode==PP_INDEX_MODE){ - sscanf(buffer,"%d",&index); + ret=read_positive_int(buffer,&index); + if(ret<0){ + fprintf(stderr,"syntax error: id_table entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer); + exit(-1); + } mode=PP_POLYNOMIAL_MODE; buffer_ptr=buffer; *buffer_ptr='\0'; } + else{ + fprintf(stderr,"syntax error: id_table entry (%d): misplaced ':'\n",j); + exit(-1); + } break; + // ignore ' ' and '\n' + case ' ':break; + case '\n':break; + // comment case '#': comment=1; @@ -628,6 +1116,10 @@ int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table if(mode!=PP_NULL_MODE){ buffer_ptr=str_addchar(buffer_ptr,str_idtable.str[j]); } + else{ + fprintf(stderr,"syntax error: id_table entry (%d): unrecognized character '%c'\n",j, str_idtable.str[j]); + exit(-1); + } break; } } @@ -635,9 +1127,13 @@ int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table //last step if(mode==PP_POLYNOMIAL_MODE){ - str_to_Polynomial(buffer,&polynomial); + parse_symbolic_expression_str(buffer, fields, variables, &polynomial); idtable_append_noinit(index,polynomial,idtable); } + else if(*buffer!='\0'){ + fprintf(stderr,"syntax error: id_table entry (%d): mismatched index/polynomial pair, got '%s'\n",j,buffer); + exit(-1); + } // sort for(j=0;j<(*idtable).length;j++){ @@ -658,6 +1154,7 @@ int parse_labels(Char_Array str_labels, Labels* labels){ int j; int mode; int comment=0; + int ret; // allocate memory init_Labels(labels,EQUATION_SIZE); @@ -686,18 +1183,30 @@ int parse_labels(Char_Array str_labels, Labels* labels){ buffer_ptr=buffer; *buffer_ptr='\0'; } + else{ + fprintf(stderr,"syntax error: labels entry (%d): mismatched '\"'\n",j); + exit(-1); + } break; case ':': // write if(mode==PP_INDEX_MODE){ sscanf(buffer,"%d",&index); + ret=read_positive_int(buffer,&index); + if(ret<0){ + fprintf(stderr,"syntax error: labels entry (%d): index must be a non-negative integer, got '%s'\n",j,buffer); + exit(-1); + } + } + else{ + fprintf(stderr,"syntax error: labels entry (%d): misplaced ':'\n",j); + exit(-1); } break; // characters to ignore case ' ':break; - case '&':break; case '\n':break; case ',':break; @@ -726,6 +1235,7 @@ int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_f int i,j; int comment_mode=0; int dcount=0; + int ret; *buffer_ptr='\0'; // loop over the input @@ -741,7 +1251,11 @@ int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_f case ',': // write init if(mpfr_flag==0){ - sscanf(buffer,"%Lf",(*init).values+intlist_find_err((*init).indices,(*init).length,index)); + ret=read_long_double(buffer,(*init).values+intlist_find_err((*init).indices,(*init).length,index)); + if(ret<0){ + fprintf(stderr,"syntax error: initial_condition entry (%d): failed to extract long double from '%s'\n",j,buffer); + exit(-1); + } } else{ mpfr_strtofr((*init_mpfr).values[intlist_find_err((*init_mpfr).indices,(*init_mpfr).length,index)], buffer, &buffer_ptr, 10, MPFR_RNDN); @@ -755,7 +1269,11 @@ int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_f // separator case ':': // write index - sscanf(buffer,"%d",&i); + ret=read_int(buffer,&i); + if(ret<0){ + fprintf(stderr,"syntax error: initial_condition entry (%d): index must be an integer, got '%s'\n",j,buffer); + exit(-1); + } if(i<0){ index=i-dcount*DOFFSET; } @@ -772,7 +1290,6 @@ int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_f // characters to ignore case ' ':break; - case '&':break; case '\n':break; // comments @@ -790,7 +1307,11 @@ int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_f // write init if(mpfr_flag==0){ - sscanf(buffer,"%Lf",(*init).values+intlist_find_err((*init).indices,(*init).length,index)); + ret=read_long_double(buffer,(*init).values+intlist_find_err((*init).indices,(*init).length,index)); + if(ret<0){ + fprintf(stderr,"syntax error: initial_condition entry (%d): failed to extract long double from '%s'\n",j,buffer); + exit(-1); + } } else{ mpfr_strtofr((*init_mpfr).values[intlist_find_err((*init_mpfr).indices,(*init_mpfr).length,index)], buffer, &buffer_ptr, 10, MPFR_RNDN); diff --git a/src/parse_file.h b/src/parse_file.h index 0430763..e22393e 100644 --- a/src/parse_file.h +++ b/src/parse_file.h @@ -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. @@ -23,26 +23,43 @@ Parse the input file #include "types.h" +// read a positive integer from a string +int read_positive_int(char* str, int* out); +// read an integer from a string +int read_int(char* str, int* out); +// read an long int from a string +int read_long_int(char* str, long int* out); +// read a long double +int read_long_double(char* str, long double* out); + + // parse fields list int parse_input_fields(Char_Array str_fields, Fields_Table* fields); -// parse symbols list -int parse_input_symbols(Char_Array str_symbols, Fields_Table* fields); +// parse virtual_fields list +int parse_input_virtual_fields(Char_Array str_virtual_fields, Fields_Table* fields, Variables variables); // parse groups of independent fields -int parse_input_groups(Char_Array str_groups, Groups* groups); +int parse_input_groups(Char_Array str_groups, Groups* groups, Polynomial_Matrix propagator, Fields_Table fields); +// check that the members of groups are independent (assuming the virtual_fields and propagator were already parsed) +int check_groups(Groups groups, Polynomial_Matrix propagator, Fields_Table fields); +// list of fields involved in a list of virtual_fields +int fields_in_virtual_field_list(Int_Array indices, Fields_Table fields, Int_Array* output); + +// parse variables list +int parse_input_variables(Char_Array str_variables, Variables* variables); // parse identities between fields -int parse_input_identities(Char_Array str_identities, Fields_Table* fields); +int parse_input_identities(Char_Array str_identities, Fields_Table* fields, Variables variables); // parse propagator int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propagator, Fields_Table fields); // parse input polynomial -int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields); +int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields, Variables variables); // parse id table -int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields); +int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields, Variables variables); // parse a list of labels int parse_labels(Char_Array str_labels, Labels* labels); diff --git a/src/polynomial.c b/src/polynomial.c index d082b60..fae0c08 100644 --- a/src/polynomial.c +++ b/src/polynomial.c @@ -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. @@ -27,6 +27,7 @@ limitations under the License. #include "array.h" #include "number.h" #include "fields.h" +#include "parse_file.h" // allocate memory @@ -368,7 +369,7 @@ int polynomial_prod_chain_nosimplify(Polynomial input, Polynomial* inout, Fields Int_Array out_monomial; Int_Array out_factor; Number out_num; - // save length of inout (which changes during the loop + // 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; @@ -630,8 +631,8 @@ int remove_unmatched_plusminus(Polynomial* polynomial, Fields_Table fields){ match_internals--; } } - // don't remove a term containing symbols - else if(type==FIELD_SYMBOL){ + // don't remove a term containing virtual_field + else if(type==FIELD_VIRTUAL){ match_internals=0; break; } @@ -1093,7 +1094,7 @@ int polynomial_print(Polynomial polynomial){ Char_Array buffer; init_Char_Array(&buffer, STR_SIZE); polynomial_sprint(polynomial, &buffer); - printf("%s",buffer.str); + printf("%s",char_array_to_str_noinit(&buffer)); free_Char_Array(buffer); return(0); } @@ -1104,7 +1105,7 @@ int polynomial_print(Polynomial polynomial){ #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){ +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; @@ -1115,6 +1116,7 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){ int comment=0; int i,j; int parenthesis_count=0; + int ret; // allocate memory init_Polynomial(output,POLY_SIZE); @@ -1174,6 +1176,10 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){ init_Int_Array(&factor, MONOMIAL_SIZE); num=number_one(); } + else{ + fprintf(stderr,"syntax error: misplaced '+' in polynomial\n"); + exit(-1); + } break; // enter monomial or factor mode @@ -1181,6 +1187,10 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){ if(mode==PP_NULL_MODE){ mode=PP_BRACKET_MODE; } + else{ + fprintf(stderr,"syntax error: misplaced '[' in polynomial\n"); + exit(-1); + } break; // factor mode case 'l': @@ -1189,6 +1199,10 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){ buffer_ptr=buffer; *buffer_ptr='\0'; } + else{ + fprintf(stderr,"syntax error: misplaced 'l' in polynomial\n"); + exit(-1); + } break; // monomial mode case 'f': @@ -1197,16 +1211,30 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){ buffer_ptr=buffer; *buffer_ptr='\0'; } + else{ + fprintf(stderr,"syntax error: misplaced 'j' in polynomial\n"); + exit(-1); + } break; // read monomial or factor case ']': sscanf(buffer,"%d",&i); + ret=read_int(buffer,&i); + if(ret<0){ + fprintf(stderr,"syntax error: in polynomial, expected integer field or factor index, got '%s'\n",buffer); + exit(-1); + } + if(mode==PP_FACTOR_MODE){ int_array_append(i,&factor); } else if(mode==PP_MONOMIAL_MODE){ int_array_append(i,&monomial); } + else{ + fprintf(stderr,"syntax error: mismatched ']' in polynomial\n"); + exit(-1); + } // switch back to null mode mode=PP_NULL_MODE; break; @@ -1224,6 +1252,10 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){ parenthesis_count++; buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]); } + else{ + fprintf(stderr,"syntax error: misplaced '(' in polynomial\n"); + exit(-1); + } break; case ')': if(mode==PP_NUMBER_MODE){ @@ -1240,11 +1272,14 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){ buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]); } } + else{ + fprintf(stderr,"syntax error: mismatched ')' in polynomial\n"); + exit(-1); + } break; // characters to ignore case ' ':break; - case '&':break; case '\n':break; // comments @@ -1257,6 +1292,10 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){ // write to buffer buffer_ptr=str_addchar(buffer_ptr,str_polynomial.str[j]); } + else{ + fprintf(stderr,"syntax error: in polynomial, unrecognized character '%c'\n",str_polynomial.str[j]); + exit(-1); + } break; } } @@ -1270,7 +1309,7 @@ int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output){ } // with str input -int str_to_Polynomial(char* str_polynomial,Polynomial* output){ +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); @@ -1278,6 +1317,15 @@ int str_to_Polynomial(char* str_polynomial,Polynomial* output){ return(0); } +// check whether the polynomial is a constant +int polynomial_is_number(Polynomial poly){ + if(poly.length==0 || (poly.length==1 && poly.monomials[0].length==0 && poly.factors[0].length==0)){ + return(1); + } + else{ + return(0); + } +} // -------------------- Polynomial_Matrix --------------------- @@ -1307,3 +1355,16 @@ int free_Polynomial_Matrix(Polynomial_Matrix matrix){ free(matrix.indices); return(0); } + +// check whether the entries are numbers +int polynomial_matrix_is_numeric(Polynomial_Matrix matrix){ + int i,j; + for(i=0;i<matrix.length;i++){ + for(j=0;j<matrix.length;j++){ + if(polynomial_is_number(matrix.matrix[i][j])==0){ + return(0); + } + } + } + return(1); +} diff --git a/src/polynomial.h b/src/polynomial.h index 6aff76d..f3488f5 100644 --- a/src/polynomial.h +++ b/src/polynomial.h @@ -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. @@ -124,12 +124,18 @@ int replace_factors(Grouped_Polynomial equations, Polynomial* polynomial); int polynomial_sprint(Polynomial polynomial, Char_Array* output); int polynomial_print(Polynomial polynomial); // read a polynomial -int Char_Array_to_Polynomial(Char_Array str_polynomial,Polynomial* output); -int str_to_Polynomial(char* str_polynomial,Polynomial* output); +int Char_Array_to_Polynomial(Char_Array str_polynomial, Polynomial* output); +int str_to_Polynomial(char* str_polynomial, Polynomial* output); + +// check whether the polynomial is a constant +int polynomial_is_number(Polynomial poly); //------------------------ Polynomial_Matrix -------------------------- // init int init_Polynomial_Matrix(Polynomial_Matrix* matrix, int length); int free_Polynomial_Matrix(Polynomial_Matrix matrix); +// check whether the entries are numbers +int polynomial_matrix_is_numeric(Polynomial_Matrix matrix); + #endif diff --git a/src/rational.h b/src/rational.h index 5beab8d..aea0625 100644 --- a/src/rational.h +++ b/src/rational.h @@ -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. diff --git a/src/rational_float.c b/src/rational_float.c index d6a6a4c..2538d41 100644 --- a/src/rational_float.c +++ b/src/rational_float.c @@ -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. diff --git a/src/rational_float.h b/src/rational_float.h index 02c4a22..ef5541b 100644 --- a/src/rational_float.h +++ b/src/rational_float.h @@ -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. diff --git a/src/rational_int.c b/src/rational_int.c index 7aab561..f3e5475 100644 --- a/src/rational_int.c +++ b/src/rational_int.c @@ -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. @@ -25,6 +25,7 @@ limitations under the License. #include <mpfr.h> #include "istring.h" #include "array.h" +#include "parse_file.h" Q quot(long int p, long int q){ Q ret; @@ -171,13 +172,18 @@ int str_to_Q(char* str, Q* num){ int mode; char* buffer=calloc(str_len(str)+1,sizeof(char)); char* buffer_ptr=buffer; + int ret; *num=quot(0,1); mode=PP_NUMERATOR_MODE; for(ptr=str;*ptr!='\0';ptr++){ if(*ptr=='/'){ - sscanf(buffer,"%ld",&((*num).numerator)); + ret=read_long_int(buffer,&((*num).numerator)); + if(ret<0){ + fprintf(stderr,"syntax error: numerator should be an integer, got '%s' in '%s'\n",buffer,str); + exit(-1); + } buffer_ptr=buffer; *buffer_ptr='\0'; mode=PP_DENOMINATOR_MODE; @@ -189,10 +195,18 @@ int str_to_Q(char* str, Q* num){ // last step if(mode==PP_NUMERATOR_MODE){ - sscanf(buffer,"%ld",&((*num).numerator)); + ret=read_long_int(buffer,&((*num).numerator)); + if(ret<0){ + fprintf(stderr,"syntax error: numerator should be an integer, got '%s' in '%s'\n",buffer,str); + exit(-1); + } } else if(mode==PP_DENOMINATOR_MODE){ - sscanf(buffer,"%ld",&((*num).denominator)); + ret=read_long_int(buffer,&((*num).denominator)); + if(ret<0){ + fprintf(stderr,"syntax error: numerator should be an integer, got '%s' in '%s'\n",buffer,str); + exit(-1); + } } free(buffer); diff --git a/src/rational_int.h b/src/rational_int.h index a4dd000..1aba15d 100644 --- a/src/rational_int.h +++ b/src/rational_int.h @@ -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. @@ -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. @@ -49,6 +49,13 @@ int RCC_cpy(RCC input,RCC* output){ } return(0); } +int RCC_cpy_noinit(RCC input,RCC* output){ + int i; + for(i=0;i<input.length;i++){ + RCC_set_elem(input.values[i], input.indices[i], output, i); + } + return(0); +} // concatenate rccs int RCC_concat(RCC rccs1, RCC rccs2, RCC* output){ @@ -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. @@ -28,6 +28,7 @@ int free_RCC(RCC rccs); int RCC_set_elem(long double value, int index, RCC* rcc, int pos); // copy int RCC_cpy(RCC input,RCC* output); +int RCC_cpy_noinit(RCC input,RCC* output); // concatenate 2 rccs int RCC_concat(RCC rccs1, RCC rccs2, RCC* output); // append an rcc to another diff --git a/src/rcc_mpfr.c b/src/rcc_mpfr.c index 8a362e3..7703b1e 100644 --- a/src/rcc_mpfr.c +++ b/src/rcc_mpfr.c @@ -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. @@ -65,6 +65,13 @@ int RCC_mpfr_cpy(RCC_mpfr input,RCC_mpfr* output){ } return(0); } +int RCC_mpfr_cpy_noinit(RCC_mpfr input,RCC_mpfr* output){ + int i; + for(i=0;i<input.length;i++){ + RCC_mpfr_set_elem(input.values[i], input.indices[i], output, i); + } + return(0); +} // concatenate rcc_mpfr int RCC_mpfr_concat(RCC_mpfr rcc_mpfr1, RCC_mpfr rcc_mpfr2, RCC_mpfr* output){ diff --git a/src/rcc_mpfr.h b/src/rcc_mpfr.h index 18ebed3..940fbb1 100644 --- a/src/rcc_mpfr.h +++ b/src/rcc_mpfr.h @@ -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. @@ -32,6 +32,7 @@ int free_RCC_mpfr(RCC_mpfr rcc_mpfr); int RCC_mpfr_set_elem(mpfr_t value, int index, RCC_mpfr* rcc_mpfr, int pos); // copy int RCC_mpfr_cpy(RCC_mpfr input,RCC_mpfr* output); +int RCC_mpfr_cpy_noinit(RCC_mpfr input,RCC_mpfr* output); // concatenate 2 rcc_mpfr_mpfr int RCC_mpfr_concat(RCC_mpfr rcc_mpfr_mpfr1, RCC_mpfr rcc_mpfr_mpfr2, RCC_mpfr* output); // append an rcc_mpfr to another diff --git a/src/symbolic_parser.c b/src/symbolic_parser.c new file mode 100644 index 0000000..08cae88 --- /dev/null +++ b/src/symbolic_parser.c @@ -0,0 +1,330 @@ +/* +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. +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 "symbolic_parser.h" +#include <stdlib.h> +#include <stdio.h> +#include "tree.h" +#include "definitions.cpp" +#include "array.h" +#include "istring.h" +#include "fields.h" +#include "polynomial.h" + +#define SP_NULL_MODE 0 +#define SP_FUNCTION_MODE 1 + +// parse a symbolic expression from a char_array +int parse_symbolic_expression(Char_Array str, Fields_Table fields, Variables variables, Polynomial* polynomial){ + Tree symbol_tree; + char_array_to_symbol_tree(str, &symbol_tree); + resolve_symbol_tree(symbol_tree, fields, variables, polynomial); + free_Tree(symbol_tree); + return(0); +} +// from char* +int parse_symbolic_expression_str(char* str, Fields_Table fields, Variables variables, Polynomial* polynomial){ + Char_Array char_array; + str_to_char_array(str,&char_array); + parse_symbolic_expression(char_array, fields, variables, polynomial); + free_Char_Array(char_array); + return(0); +} + +// compute the symbol tree from a string +int char_array_to_symbol_tree(Char_Array str, Tree* symbol_tree){ + // buffer + char* buffer=calloc(str.length+1,sizeof(char)); + char* buffer_ptr=buffer; + Tree child; + int match; + Char_Array nodestr; + Char_Array label; + Char_Array str_clean; + int mode; + int comment=0; + int j; + int gotanode=0; + + // allocate memory + init_Tree(symbol_tree,SYMBOL_TREE_SIZE, SYMBOL_TREE_LABEL_SIZE); + + // remove comments, ' ' and '\n' + init_Char_Array(&str_clean,str.length); + for(j=0;j<str.length;j++){ + if(comment==1){ + if(str.str[j]=='\n'){ + comment=0; + } + } + else{ + switch(str.str[j]){ + case ' ':break; + case '\n':break; + // comments + case '#': + comment=1; + break; + default: + char_array_append(str.str[j],&str_clean); + break; + } + } + } + + // if the string contains no '<', then trivial tree + for(j=0;j<str_clean.length;j++){ + if(str_clean.str[j]=='<'){ + break; + } + } + // no '<': trivial tree + if(j==str_clean.length){ + tree_set_label(str_clean, symbol_tree); + free(buffer); + free_Char_Array(str_clean); + return(0); + } + + *buffer_ptr='\0'; + // loop over the input string + // start in null mode + mode=SP_NULL_MODE; + for(j=0;j<str_clean.length;j++){ + switch(str_clean.str[j]){ + // new node + case '<': + // find matching bracket + match=matching_bracket(str_clean,j); + // check whether it exists + if(match<0){ + fprintf(stderr,"syntax error: unmatched brackets in %s\n",char_array_to_str_noinit(&str)); + exit(-1); + } + // extract substring until bracket + char_array_substring(str_clean,j+1,match-1,&nodestr); + + // check whether node is trivial + if(j==0 && match==str_clean.length-1){ + free_Tree(*symbol_tree); + char_array_to_symbol_tree(nodestr, symbol_tree); + free_Char_Array(nodestr); + j=match; + break; + } + + // parse subexpression + char_array_to_symbol_tree(nodestr, &child); + free_Char_Array(nodestr); + // add child to tree + tree_append_child_noinit(child, symbol_tree); + // boolean indicating a node has been found + gotanode=1; + // set next position after the node + j=match; + + // if function mode, then check that the match is at the end of the node + if(mode==SP_FUNCTION_MODE){ + if(match<str_clean.length-1){ + fprintf(stderr,"syntax error: functions must occupy an entire node (e.g. <%%exp<...>>), got %s\n",char_array_to_str_noinit(&str)); + exit(-1); + } + else{ + // set label + str_to_char_array(buffer,&label); + tree_set_label(label,symbol_tree); + free_Char_Array(label); + } + } + break; + + // function + case '%': + if(j>0){ + fprintf(stderr,"syntax error: functions must occupy an entire node (e.g. <%%exp<...>>), got %s\n",char_array_to_str_noinit(&str)); + exit(-1); + } + mode=SP_FUNCTION_MODE; + break; + + // product + case '*': + if(gotanode==0){ + fprintf(stderr,"syntax error: '*' is not preceded by a node in %s\n",char_array_to_str_noinit(&str)); + exit(-1); + } + if(j>=str_clean.length-1){ + fprintf(stderr,"syntax error: '*' cannot be at the end of an expression, got %s\n",char_array_to_str_noinit(&str)); + exit(-1); + } + // set label + init_Char_Array(&label,1); + char_array_append('*',&label); + tree_set_label(label,symbol_tree); + free_Char_Array(label); + // next child + char_array_substring(str_clean,j+1,str_clean.length-1,&nodestr); + // parse subexpression + char_array_to_symbol_tree(nodestr, &child); + free_Char_Array(nodestr); + // append next child + tree_append_child_noinit(child, symbol_tree); + + // make it stop + j=str_clean.length-1; + break; + + // sum + case '+': + if(gotanode==0){ + fprintf(stderr,"syntax error: '+' is not preceded by a node in %s\n",char_array_to_str_noinit(&str)); + exit(-1); + } + if(j>=str_clean.length-1){ + fprintf(stderr,"syntax error: '+' cannot be at the end of an expression, got %s\n",char_array_to_str_noinit(&str)); + exit(-1); + } + // set label + init_Char_Array(&label,1); + char_array_append('+',&label); + tree_set_label(label,symbol_tree); + free_Char_Array(label); + // next child + char_array_substring(str_clean,j+1,str_clean.length-1,&nodestr); + // parse subexpression + char_array_to_symbol_tree(nodestr, &child); + free_Char_Array(nodestr); + // append next child + tree_append_child_noinit(child, symbol_tree); + + // make it stop + j=str_clean.length-1; + break; + + default: + if(mode!=SP_NULL_MODE){ + // write to buffer + buffer_ptr=str_addchar(buffer_ptr,str_clean.str[j]); + } + break; + } + } + + free_Char_Array(str_clean); + free(buffer); + return(0); +} +// from char* +int str_to_symbol_tree(char* str, Tree* symbol_tree){ + Char_Array char_array; + str_to_char_array(str,&char_array); + char_array_to_symbol_tree(char_array,symbol_tree); + free_Char_Array(char_array); + return(0); +} + + +// find matching '<' and '>' +int matching_bracket(Char_Array str, int start){ + int bracket_count=0; + int i; + for(i=start;i<str.length;i++){ + if(str.str[i]=='<'){ + bracket_count++; + } + else if(str.str[i]=='>'){ + bracket_count--; + if(bracket_count==0){ + return(i); + } + } + } + // if the function has not returned, then no matching bracket + return(-1); +} + + +// resolve a symbol tree to its corresponding polynomial +int resolve_symbol_tree(Tree symbol_tree, Fields_Table fields, Variables variables, Polynomial* output){ + Polynomial poly; + Tree variable_tree; + + // trivial tree + if(symbol_tree.length==0){ + // variable + if(symbol_tree.root_label.length>0 && symbol_tree.root_label.str[0]=='$'){ + variables_find_var(symbol_tree.root_label, variables, &variable_tree); + resolve_symbol_tree(variable_tree, fields, variables, output); + free_Tree(variable_tree); + } + //polynomial + else{ + Char_Array_to_Polynomial(symbol_tree.root_label, output); + } + } + + // exp + else if (char_array_cmp_str(symbol_tree.root_label,"exp")==1){ + if(symbol_tree.length!=1){ + fprintf(stderr,"syntax error: exp must have 1 argument\n"); + exit(-1); + } + resolve_symbol_tree(symbol_tree.children[0], fields, variables, &poly); + polynomial_exponential(poly, output, fields); + free_Polynomial(poly); + } + + // log + else if (char_array_cmp_str(symbol_tree.root_label,"log_1")==1){ + if(symbol_tree.length!=1){ + fprintf(stderr,"syntax error: log_1 must have 1 argument\n"); + exit(-1); + } + resolve_symbol_tree(symbol_tree.children[0], fields, variables, &poly); + polynomial_logarithm(poly, output, fields); + free_Polynomial(poly); + } + + // product + else if (char_array_cmp_str(symbol_tree.root_label,"*")==1){ + if(symbol_tree.length!=2){ + fprintf(stderr,"syntax error: '*' must have 2 arguments\n"); + exit(-1); + } + resolve_symbol_tree(symbol_tree.children[0], fields, variables, output); + resolve_symbol_tree(symbol_tree.children[1], fields, variables, &poly); + polynomial_prod_chain(poly, output, fields); + free_Polynomial(poly); + } + + // sum + else if (char_array_cmp_str(symbol_tree.root_label,"+")==1){ + if(symbol_tree.length!=2){ + fprintf(stderr,"syntax error: '+' must have 2 arguments\n"); + exit(-1); + } + resolve_symbol_tree(symbol_tree.children[0], fields, variables, output); + resolve_symbol_tree(symbol_tree.children[1], fields, variables, &poly); + polynomial_add_chain_noinit(poly, output, fields); + } + + else{ + fprintf(stderr,"syntax error: unrecognized operation '%s'\n",char_array_to_str_noinit(&(symbol_tree.root_label))); + exit(-1); + } + + return(0); +} diff --git a/src/symbolic_parser.h b/src/symbolic_parser.h new file mode 100644 index 0000000..e91c902 --- /dev/null +++ b/src/symbolic_parser.h @@ -0,0 +1,42 @@ +/* +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. +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. +*/ + +/* + parse symbolic expressions +*/ + +#ifndef SYMBOLIC_PARSER_H +#define SYMBOLIC_PARSER_H + +#include "types.h" + +// parse a symbolic expression from a char_array +int parse_symbolic_expression(Char_Array str, Fields_Table fields, Variables variables, Polynomial* polynomial); +// from char* +int parse_symbolic_expression_str(char* str, Fields_Table fields, Variables variables, Polynomial* polynomial); + +// compute the symbol tree from a string +int char_array_to_symbol_tree(Char_Array str, Tree* symbol_tree); +// from char* +int str_to_symbol_tree(char* str, Tree* symbol_tree); + +// find matching '<' and '>' +int matching_bracket(Char_Array str, int start); + +// resolve a symbol tree to its corresponding polynomial +int resolve_symbol_tree(Tree symbol_tree, Fields_Table fields, Variables variables, Polynomial* output); + +#endif diff --git a/src/tools.c b/src/tools.c index 2b7e85d..5912819 100644 --- a/src/tools.c +++ b/src/tools.c @@ -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. @@ -140,3 +140,13 @@ int min(int x1, int x2){ return(x2); } } + +// check whether a divides b +int is_factor(int a, int b){ + if(b-a*(b/a)==0){ + return(1); + } + else{ + return(0); + } +} diff --git a/src/tools.h b/src/tools.h index e5f319f..8f29cf8 100644 --- a/src/tools.h +++ b/src/tools.h @@ -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. @@ -46,4 +46,7 @@ int intlist_find_err(int* list, int size, int x); int max(int x1, int x2); int min(int x1, int x2); +// check whether a divides b +int is_factor(int a, int b); + #endif diff --git a/src/tree.c b/src/tree.c new file mode 100644 index 0000000..0edf558 --- /dev/null +++ b/src/tree.c @@ -0,0 +1,117 @@ +/* +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. +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 "tree.h" +#include <stdlib.h> +#include <stdio.h> +#include "array.h" + +// init +int init_Tree(Tree* tree, int memory_children, int memory_label){ + init_Char_Array(&(tree->root_label),memory_label); + (*tree).children=calloc(memory_children,sizeof(Tree)); + (*tree).memory=memory_children; + (*tree).length=0; + return(0); +} +int free_Tree(Tree tree){ + int i; + free_Char_Array(tree.root_label); + for(i=0;i<tree.length;i++){ + free_Tree(tree.children[i]); + } + free(tree.children); + return(0); +} + +// copy +int tree_cpy(Tree input, Tree* output){ + init_Tree(output,input.length, input.root_label.length); + tree_cpy_noinit(input,output); + return(0); +} +int tree_cpy_noinit(Tree input, Tree* output){ + int i; + if((*output).memory<input.length){ + fprintf(stderr,"error: trying to copy a tree of length %d to another with memory %d\n",input.length, (*output).memory); + exit(-1); + } + char_array_cpy_noinit(input.root_label,&(output->root_label)); + for(i=0;i<input.length;i++){ + tree_cpy(input.children[i],(*output).children+i); + } + + (*output).length=input.length; + return(0); +} + +// resize memory +int tree_resize(Tree* tree, int newsize){ + Tree new_tree; + init_Tree(&new_tree,newsize,tree->root_label.memory); + tree_cpy_noinit(*tree,&new_tree); + free_Tree(*tree); + *tree=new_tree; + return(0); +} + +// set label +int tree_set_label(Char_Array label, Tree* tree){ + if(label.length > tree->root_label.memory){ + char_array_resize(&(tree->root_label),label.length); + } + char_array_cpy_noinit(label,&(tree->root_label)); + return(0); +} + +// add a child to a tree +int tree_append_child(Tree child, Tree* output){ + if((*output).length>=(*output).memory){ + tree_resize(output,2*(*output).memory+1); + } + tree_cpy(child,(*output).children+(*output).length); + (*output).length++; + return(0); +} +// add a child to a tree without allocating memory for the new child +int tree_append_child_noinit(Tree child, Tree* output){ + if((*output).length>=(*output).memory){ + tree_resize(output,2*(*output).memory+1); + } + (*output).children[(*output).length]=child; + (*output).length++; + return(0); +} + +// concatenate the children of two trees +int tree_concat_children(Tree input, Tree* output){ + int i; + for(i=0;i<input.length;i++){ + tree_append_child(input.children[i],output); + } + return(0); +} +// noinit +int tree_concat_children_noinit(Tree input, Tree* output){ + int i; + for(i=0;i<input.length;i++){ + tree_append_child_noinit(input.children[i],output); + } + // free input array + free(input.children); + return(0); +} + diff --git a/src/tree.h b/src/tree.h new file mode 100644 index 0000000..60b222c --- /dev/null +++ b/src/tree.h @@ -0,0 +1,48 @@ +/* +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. +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. +*/ + +/* Trees */ + +#ifndef TREE_H +#define TREE_H + +#include "types.h" + +int init_Tree(Tree* tree, int memory_children, int memory_label); +int free_Tree(Tree tree); + +// copy +int tree_cpy(Tree input, Tree* output); +int tree_cpy_noinit(Tree input, Tree* output); + +// resize memory +int tree_resize(Tree* tree, int newsize); + +// set label +int tree_set_label(Char_Array label, Tree* tree); + +// add a child to a tree +int tree_append_child(Tree tree, Tree* output); +// add a child to a tree without allocating memory for the new child +int tree_append_child_noinit(Tree tree, Tree* output); + +// concatenate the children of two trees +int tree_concat_children(Tree input, Tree* output); +// noinit +int tree_concat_children_noinit(Tree input, Tree* output); + +#endif + diff --git a/src/types.h b/src/types.h index d30840e..00a34ac 100644 --- a/src/types.h +++ b/src/types.h @@ -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. @@ -71,6 +71,15 @@ typedef struct Str_Array{ int memory; } Str_Array; +// tree +typedef struct Tree Tree; +struct Tree{ + Char_Array root_label; + Tree* children; + int length; + int memory; +}; + // polynomial typedef struct Polynomial{ Int_Array* monomials; @@ -133,13 +142,21 @@ typedef struct Identities{ int memory; } Identities; -// symbolic expressions -typedef struct Symbols{ +// virtual_fields +typedef struct Virtual_fields{ int* indices; Polynomial* expr; int length; int memory; -} Symbols; +} Virtual_fields; + +// variables used in symbolic expressions +typedef struct Variables{ + Char_Array* var_names; + Tree* symbol_trees; + int length; + int memory; +} Variables; // groups of independent fields typedef struct Groups{ @@ -159,10 +176,10 @@ typedef struct Fields_Table{ // identities between fields Identities ids; // symbolic expressions (commuting) - Symbols symbols; - // list of anti-commuting variables (fields or symbols) + Virtual_fields virtual_fields; + // list of anti-commuting variables Int_Array fermions; - // list of non-commuting variables (fields or symbols) + // list of non-commuting variables Int_Array noncommuting; } Fields_Table; @@ -187,6 +204,8 @@ typedef struct Id_Table{ typedef struct Meankondo_Options{ int threads; int chain; + int print_progress; + int group_poly; } Meankondo_Options; typedef struct Numkondo_Options{ @@ -203,6 +222,7 @@ typedef struct Meantools_Options{ Int_Array deriv_vars; Char_Array eval_rccstring; int chain; + Char_Array namespace; mpfr_prec_t mpfr_prec; mpfr_exp_t mpfr_emax; } Meantools_Options; |