From 3f0510629e422e979b57d3f93791937912a4183a Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Tue, 14 Jun 2022 09:26:07 +0200 Subject: Update to 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. * Created vim files to implement syntax highlighting for configuration files. * Multiple bug fixes --- src/array.c | 69 ++++- src/array.h | 14 +- src/cli_parser.c | 22 +- src/cli_parser.h | 4 +- src/coefficient.c | 245 +++++++++++++++- src/coefficient.h | 12 +- src/definitions.cpp | 10 +- src/determinant.c | 93 ++++++ src/determinant.h | 17 ++ src/expansions.c | 28 -- src/expansions.h | 34 --- src/fields.c | 273 +++++++++++++++--- src/fields.h | 58 +++- src/flow.c | 73 ++++- src/flow.h | 4 +- src/flow_mpfr.c | 63 ++++- src/flow_mpfr.h | 4 +- src/grouped_polynomial.c | 138 +++++++-- src/grouped_polynomial.h | 10 +- src/idtable.c | 2 +- src/idtable.h | 2 +- src/istring.c | 2 +- src/istring.h | 2 +- src/kondo.c | 69 ++--- src/kondo.h | 14 +- src/kondo_preprocess.c | 7 +- src/mean.c | 141 +++++++--- src/mean.h | 20 +- src/meankondo.c | 213 +++++++++++--- src/meantools.c | 29 +- src/meantools_deriv.c | 14 +- src/meantools_deriv.h | 4 +- src/meantools_eval.c | 6 +- src/meantools_eval.h | 2 +- src/meantools_exp.c | 130 --------- src/meantools_exp.h | 27 -- src/meantools_expand.c | 166 +++++++++++ src/meantools_expand.h | 29 ++ src/number.c | 223 +++++++++++---- src/number.h | 15 +- src/numkondo.c | 27 +- src/parse_file.c | 717 ++++++++++++++++++++++++++++++++++++++++------- src/parse_file.h | 31 +- src/polynomial.c | 77 ++++- src/polynomial.h | 12 +- src/rational.h | 2 +- src/rational_float.c | 2 +- src/rational_float.h | 2 +- src/rational_int.c | 22 +- src/rational_int.h | 2 +- src/rcc.c | 9 +- src/rcc.h | 3 +- src/rcc_mpfr.c | 9 +- src/rcc_mpfr.h | 3 +- src/symbolic_parser.c | 330 ++++++++++++++++++++++ src/symbolic_parser.h | 42 +++ src/tools.c | 12 +- src/tools.h | 5 +- src/tree.c | 117 ++++++++ src/tree.h | 48 ++++ src/types.h | 34 ++- 61 files changed, 3091 insertions(+), 703 deletions(-) create mode 100644 src/determinant.c create mode 100644 src/determinant.h delete mode 100644 src/expansions.c delete mode 100644 src/expansions.h delete mode 100644 src/meantools_exp.c delete mode 100644 src/meantools_exp.h create mode 100644 src/meantools_expand.c create mode 100644 src/meantools_expand.h create mode 100644 src/symbolic_parser.c create mode 100644 src/symbolic_parser.h create mode 100644 src/tree.c create mode 100644 src/tree.h (limited to 'src') 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;i0){ + 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=1;power--){ + ret=coefficient_simplify_fraction(constant, quotient_prev, &remainder, "ient); + + // if fail to simplify, stop + if(ret<0){ + if(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;i0 && (*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=M.length){ + *sign_correction=0; + return(0); + } + // pivot if needed + if(pivot!=j){ + for(k=0;k -#include -#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=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=(*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=(*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;i0){ + // ignore constants + for(j=0;j0){ + RCC_cpy_noinit(rccs,&rcc_print); + // ignore constants + for(j=0;j0){ + RCC_cpy_noinit(rccs,&rcc_print); + // ignore constants + for(j=0;j0){ + // ignore constants + for(j=0;j0){ + // ignore constants + for(j=0;j0){ + 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;j0){ - 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=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;k0){ for(i=0;i=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 diff --git a/src/mean.c b/src/mean.c index 39ece56..0c5707e 100644 --- a/src/mean.c +++ b/src/mean.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,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=0 && move=0 && move0 && 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\n\n"); + printf("\nusage:\n meankondo [-t threads] [-C] [-p] [-A] \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=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=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 -#include -#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 +#include +#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;inamespace)); + 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_expand.h b/src/meantools_expand.h new file mode 100644 index 0000000..9d988b1 --- /dev/null +++ b/src/meantools_expand.h @@ -0,0 +1,29 @@ +/* +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. +*/ + +#ifndef MEANTOOLS_EXPAND_H +#define MEANTOOLS_EXPAND_H + +#include "types.h" + +// read arguments +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;j0){ - 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;ilength;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;ilength;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=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'){ + 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 #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. diff --git a/src/rcc.c b/src/rcc.c index 484e20a..89d2c96 100644 --- a/src/rcc.c +++ b/src/rcc.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. @@ -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 +#include +#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>), 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'){ + 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 +#include +#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;iroot_label)); + for(i=0;iroot_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