From aa0f3ae2988d372b190b9bde2e75a6d17e744e93 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Sun, 14 Jun 2015 00:52:45 +0000 Subject: Initial commit --- src/array.c | 628 ++++++++++++++++++++ src/array.h | 134 +++++ src/cli_parser.c | 123 ++++ src/cli_parser.h | 35 ++ src/coefficient.c | 739 +++++++++++++++++++++++ src/coefficient.h | 78 +++ src/definitions.cpp | 60 ++ src/expansions.c | 28 + src/expansions.h | 34 ++ src/fields.c | 489 ++++++++++++++++ src/fields.h | 98 ++++ src/flow.c | 167 ++++++ src/flow.h | 35 ++ src/grouped_polynomial.c | 765 ++++++++++++++++++++++++ src/grouped_polynomial.h | 74 +++ src/idtable.c | 122 ++++ src/idtable.h | 44 ++ src/istring.c | 99 ++++ src/istring.h | 40 ++ src/kondo.c | 1449 ++++++++++++++++++++++++++++++++++++++++++++++ src/kondo.h | 76 +++ src/kondo_preprocess.c | 126 ++++ src/mean.c | 793 +++++++++++++++++++++++++ src/mean.h | 70 +++ src/meankondo.c | 301 ++++++++++ src/meantools.c | 116 ++++ src/meantools_deriv.c | 195 +++++++ src/meantools_deriv.h | 30 + src/meantools_eval.c | 129 +++++ src/meantools_eval.h | 29 + src/meantools_exp.c | 130 +++++ src/meantools_exp.h | 27 + src/number.c | 551 ++++++++++++++++++ src/number.h | 120 ++++ src/numkondo.c | 226 ++++++++ src/parse_file.c | 796 +++++++++++++++++++++++++ src/parse_file.h | 56 ++ src/polynomial.c | 1263 ++++++++++++++++++++++++++++++++++++++++ src/polynomial.h | 131 +++++ src/rational.h | 23 + src/rational_float.c | 196 +++++++ src/rational_float.h | 64 ++ src/rational_int.c | 190 ++++++ src/rational_int.h | 59 ++ src/rcc.c | 95 +++ src/rcc.h | 39 ++ src/tools.c | 142 +++++ src/tools.h | 49 ++ src/types.h | 219 +++++++ 49 files changed, 11482 insertions(+) create mode 100644 src/array.c create mode 100644 src/array.h create mode 100644 src/cli_parser.c create mode 100644 src/cli_parser.h create mode 100644 src/coefficient.c create mode 100644 src/coefficient.h create mode 100644 src/definitions.cpp create mode 100644 src/expansions.c create mode 100644 src/expansions.h create mode 100644 src/fields.c create mode 100644 src/fields.h create mode 100644 src/flow.c create mode 100644 src/flow.h create mode 100644 src/grouped_polynomial.c create mode 100644 src/grouped_polynomial.h create mode 100644 src/idtable.c create mode 100644 src/idtable.h create mode 100644 src/istring.c create mode 100644 src/istring.h create mode 100644 src/kondo.c create mode 100644 src/kondo.h create mode 100644 src/kondo_preprocess.c create mode 100644 src/mean.c create mode 100644 src/mean.h create mode 100644 src/meankondo.c create mode 100644 src/meantools.c create mode 100644 src/meantools_deriv.c create mode 100644 src/meantools_deriv.h create mode 100644 src/meantools_eval.c create mode 100644 src/meantools_eval.h create mode 100644 src/meantools_exp.c create mode 100644 src/meantools_exp.h create mode 100644 src/number.c create mode 100644 src/number.h create mode 100644 src/numkondo.c create mode 100644 src/parse_file.c create mode 100644 src/parse_file.h create mode 100644 src/polynomial.c create mode 100644 src/polynomial.h create mode 100644 src/rational.h create mode 100644 src/rational_float.c create mode 100644 src/rational_float.h create mode 100644 src/rational_int.c create mode 100644 src/rational_int.h create mode 100644 src/rcc.c create mode 100644 src/rcc.h create mode 100644 src/tools.c create mode 100644 src/tools.h create mode 100644 src/types.h (limited to 'src') diff --git a/src/array.c b/src/array.c new file mode 100644 index 0000000..11d14f7 --- /dev/null +++ b/src/array.c @@ -0,0 +1,628 @@ +/* +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 "array.h" +#include +#include +#include +#include "istring.h" +#include "definitions.cpp" + +// init +int init_Int_Array(Int_Array* array, int memory){ + (*array).values=calloc(memory,sizeof(int)); + (*array).memory=memory; + (*array).length=0; + return(0); +} +int free_Int_Array(Int_Array array){ + free(array.values); + return(0); +} + +// copy +int int_array_cpy(Int_Array input, Int_Array* output){ + init_Int_Array(output,input.length); + int_array_cpy_noinit(input,output); + return(0); +} +int int_array_cpy_noinit(Int_Array input, Int_Array* output){ + int i; + if((*output).memory=(*output).memory){ + int_array_resize(output,2*(*output).memory+1); + } + (*output).values[(*output).length]=val; + (*output).length++; + return(0); +} + +// concatenate +int int_array_concat(Int_Array input, Int_Array* output){ + int i; + int offset=(*output).length; + if((*output).length+input.length>(*output).memory){ + // make it longer than needed by (*output).length (for speed) + int_array_resize(output,2*(*output).length+input.length); + } + for(i=0;iarray2.length){ + return(1); + } + + // compare terms + for(i=0;iarray2.values[i]){ + return(1); + } + } + + // if equal + return(0); +} + +// check whether an array is a sub-array of another +// allows for the elements to be separated by others, but the order must be respected +int int_array_is_subarray_ordered(Int_Array input, Int_Array test_array){ + int i; + int matches=0; + + for(i=0;i=(*output).memory){ + char_array_resize(output,2*(*output).memory+1); + } + (*output).str[(*output).length]=val; + (*output).length++; + return(0); +} +// append a string +int char_array_append_str(char* str, Char_Array* output){ + char* ptr; + for (ptr=str;*ptr!='\0';ptr++){ + char_array_append(*ptr, output); + } + return(0); +} + +// concatenate +int char_array_concat(Char_Array input, Char_Array* output){ + int i; + int offset=(*output).length; + if((*output).length+input.length>(*output).memory){ + // make it longer than needed by (*output).length (for speed) + char_array_resize(output,2*(*output).length+input.length); + } + for(i=0;isize){ + // resize + free(out_str); + // +1 for '\0' + size=extra_size+1; + out_str=calloc(size,sizeof(char)); + // read format again + va_start(vaptr, fmt); + vsnprintf(out_str,size,fmt,vaptr); + va_end(vaptr); + } + + // write to char array + for(ptr=out_str;*ptr!='\0';ptr++){ + char_array_append(*ptr, output); + } + + free(out_str); + + return(0); +} + + +// replace '%' with given character +int replace_star(char c, Char_Array str, Char_Array* out){ + int i; + init_Char_Array(out, str.length); + for(i=0;i=(*output).memory){ + str_array_resize(output,2*(*output).memory+1); + } + char_array_cpy(val, (*output).strs+(*output).length); + (*output).length++; + return(0); +} +int str_array_append_noinit(Char_Array val, Str_Array* output){ + if((*output).length>=(*output).memory){ + str_array_resize(output,2*(*output).memory+1); + } + (*output).strs[(*output).length]=val; + (*output).length++; + return(0); +} + +// concatenate +int str_array_concat(Str_Array input, Str_Array* output){ + int i; + int offset=(*output).length; + if((*output).length+input.length>(*output).memory){ + // make it longer than needed by (*output).length (for speed) + str_array_resize(output,2*(*output).length+input.length); + } + for(i=0;i(*output).memory){ + // make it longer than needed by (*output).length (for speed) + str_array_resize(output,2*(*output).length+input.length); + } + for(i=0;i=(*output).memory){ + resize_labels(output,2*(*output).memory); + } + + // copy and allocate + char_array_cpy(label,(*output).labels+offset); + (*output).indices[offset]=index; + // increment length + (*output).length++; + return(0); +} +// append an element to a labels without allocating memory +int labels_append_noinit(Char_Array label, int index, Labels* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_labels(output,2*(*output).memory); + } + + // copy without allocating + (*output).labels[offset]=label; + (*output).indices[offset]=index; + // increment length + (*output).length++; + return(0); +} + +// concatenate two labels tables +int labels_concat(Labels input, Labels* output){ + int i; + for(i=0;i +#include +#include "definitions.cpp" +#include "array.h" + +// read a config file and write the output to str_args +int read_config_file(Str_Array* str_args, const char* file, int read_from_stdin){ + FILE* infile; + char c; + Char_Array buffer; + + // allocate memory + init_Str_Array(str_args, ARG_COUNT); + init_Char_Array(&buffer, STR_SIZE); + + // read from file + if(read_from_stdin==0){ + infile=fopen(file,"r"); + if(infile==NULL){ + fprintf(stderr,"error: can't open file %s\n",file); + exit(-1); + } + } + else{ + infile=stdin; + } + + while(1){ + c=fgetc(infile); + // add to str_args + if(c=='&'){ + str_array_append_noinit(buffer,str_args); + init_Char_Array(&buffer, STR_SIZE); + } + else if(c==EOF){ + str_array_append_noinit(buffer,str_args); + break; + } + else{ + char_array_append(c,&buffer); + } + } + fclose(infile); + + return(0); +} + + +// get the title from a string argument +int get_str_arg_title(Char_Array str_arg, Char_Array* out){ + int j,k; + + init_Char_Array(out,STR_SIZE); + + // find "#!" at beginning of line + for(j=0;j=str_arg.length-2){ + fprintf(stderr, "error: an entry in the configuration file does not have a title (which should be preceeded by '#!' at the beginning of a line)\n"); + exit(-1); + } + + // get title until end of line + for(k=j+2;k +#include +#include "definitions.cpp" +#include "rational.h" +#include "istring.h" +#include "array.h" +#include "number.h" +#include "tools.h" + + +// allocate memory +int init_Coefficient(Coefficient* coef,int size){ + (*coef).factors=calloc(size,sizeof(Int_Array)); + (*coef).nums=calloc(size,sizeof(Number)); + (*coef).denoms=calloc(size,sizeof(coef_denom)); + (*coef).length=0; + (*coef).memory=size; + + return(0); +} + +// free memory +int free_Coefficient(Coefficient coef){ + int i; + for(i=0;i(*output).memory){ + resize_Coefficient(output,input.length); + } + + for(i=0;i=(*output).memory){ + resize_Coefficient(output,2*(*output).memory+1); + } + // copy and allocate + int_array_cpy(factor,(*output).factors+offset); + number_cpy(num,(*output).nums+offset); + (*output).denoms[offset]=denom; + // increment length + (*output).length++; + return(0); +} +// append an element to a coefficient without allocating memory +int coefficient_append_noinit(Int_Array factor, Number num, coef_denom denom, Coefficient* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_Coefficient(output,2*(*output).memory+1); + } + // copy without allocating + (*output).factors[offset]=factor; + (*output).nums[offset]=num; + (*output).denoms[offset]=denom; + // increment length + (*output).length++; + return(0); +} + +// concatenate coefficients and simplify result +int coefficient_concat(Coefficient input, Coefficient* output){ + int i; + for(i=0;i0){ + at_least_one=1; + coefficient_append_noinit(factor,number_Qprod_ret(quot(match_count,1),input.nums[i]), input.denoms[i],output); + } + else{ + free_Int_Array(factor); + } + } + + if(at_least_one>0){ + coefficient_simplify(output); + } + else{ + // add a trivial element to the coefficient + init_Int_Array(&factor,0); + denom.index=-1; + denom.power=0; + coefficient_append_noinit(factor,number_zero(),denom,output); + } + + return(0); +} + +// derive a monomial with respect to an index +int monomial_deriv(Int_Array factor, int index, Int_Array* out_factor, int* match_count){ + int j; + + init_Int_Array(out_factor,factor.length); + // init match count + *match_count=0; + + // loop over indices + for(j=0;jdenom2.index){ + return(-1); + } + + if(denom1.powerdenom2.power){ + return(1); + } + + return(0); +} + + +// evaluate a coefficient on a vector +int evalcoef(RCC rccs, Coefficient coef, long double* out){ + int i,j; + long double num_factor; + + *out=0; + + // for each monomial + for(i=0;i0){ + num_factor/=dpower(rccs.values[intlist_find_err(rccs.indices,rccs.length,coef.denoms[i].index)],coef.denoms[i].power); + } + *out+=num_factor*number_double_val(coef.nums[i]); + } + return(0); +} diff --git a/src/coefficient.h b/src/coefficient.h new file mode 100644 index 0000000..a7a151c --- /dev/null +++ b/src/coefficient.h @@ -0,0 +1,78 @@ +/* +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. +*/ + +/* +Coefficients represent the collection of terms which appear + on the right hand side of a single flow equation. + +They are used as an elementary building block of grouped polynomials. +*/ + +#ifndef COEFFICIENT_H +#define COEFFICIENT_H + +#include "types.h" + +// allocate memory +int init_Coefficient(Coefficient* coef,int size); +// free memory +int free_Coefficient(Coefficient coef); + +// copy a coefficient +int coefficient_cpy(Coefficient input, Coefficient* output); +int coefficient_cpy_noinit(Coefficient input, Coefficient* output); + +// resize the memory allocated to a coefficient +int resize_Coefficient(Coefficient* coefficient,int new_size); + +// append an element to a coefficient +int coefficient_append(Int_Array factor,Number num, coef_denom denom, Coefficient* output); +int coefficient_append_noinit(Int_Array factor, Number num, coef_denom denom, Coefficient* output); + +// concatenate coefficients and simplify result +int coefficient_concat(Coefficient input, Coefficient* output); +int coefficient_concat_noinit(Coefficient input, Coefficient* output); + +// simplify a Coefficient +int coefficient_simplify(Coefficient* coefficient); +// sort the terms in an equation (quicksort algorithm) +int sort_coefficient(Coefficient* coefficient, int begin, int end); +// exchange two terms (for the sorting algorithm) +int exchange_coefficient_terms(int i, int j, Coefficient* coefficient); + +// derive a coefficient with respect to an index (as a polynomial) (does not derive the 1/(1+C)^p ) +int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output); +int coefficient_deriv(Coefficient input, int index, Coefficient* output); + +// product of two coefficients +int coefficient_prod(Coefficient coef1, Coefficient coef2, Coefficient* output); +// product of coefficients, output replaces the second coefficient +int coefficient_prod_chain(Coefficient in, Coefficient* inout); + +// print coefficient +int coefficient_sprint(Coefficient coef, Char_Array* output, int offset, char index_pre); + +// read from a string +int char_array_to_Coefficient(Char_Array str_coef, Coefficient* output); +int str_to_Coefficient(char* str, Coefficient* output); + +// compare coefficient denominators +int coef_denom_cmp(coef_denom denom1, coef_denom denom2); + +// evaluate a coefficient on a vector +int evalcoef(RCC rccs, Coefficient coef, long double* out); + +#endif diff --git a/src/definitions.cpp b/src/definitions.cpp new file mode 100644 index 0000000..d32537d --- /dev/null +++ b/src/definitions.cpp @@ -0,0 +1,60 @@ +/* +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. +*/ + +#ifndef DEFINITIONS_GCC +#define DEFINITIONS_GCC + +#define VERSION "1.2" + +// number of entries in a configuration file +#define ARG_COUNT 10 +// size of string representing a monomial +#define MONOMIAL_SIZE 20 +// size of various strings +#define STR_SIZE 100 +// number of terms in coefficients +#define COEF_SIZE 100 +// number of terms in polynomials +#define POLY_SIZE 100 +// number of equations +#define EQUATION_SIZE 20 +// number of fields +#define FIELDS_SIZE 50 +// number of elements in numbers +#define NUMBER_SIZE 5 +// number of elements in a group +#define GROUP_SIZE 5 + + +// display options +#define DISPLAY_EQUATION 1 +#define DISPLAY_NUMERICAL 2 +#define DISPLAY_EXPLOG 3 +#define DISPLAY_FINAL 4 + +// available preprocessors +#define PREPROCESSOR_KONDO 1 + +// offset derivative indices +#define DOFFSET 1000000 + +// types of fields (the order matters) +#define FIELD_PARAMETER 1 +#define FIELD_EXTERNAL 2 +#define FIELD_INTERNAL 3 +#define FIELD_SYMBOL 4 + +#endif diff --git a/src/expansions.c b/src/expansions.c new file mode 100644 index 0000000..c889829 --- /dev/null +++ b/src/expansions.c @@ -0,0 +1,28 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "expansions.h" + +#include +#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 new file mode 100644 index 0000000..85d6c2b --- /dev/null +++ b/src/expansions.h @@ -0,0 +1,34 @@ +/* +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 new file mode 100644 index 0000000..1b221f2 --- /dev/null +++ b/src/fields.c @@ -0,0 +1,489 @@ +/* +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 "fields.h" +#include "definitions.cpp" +#include +#include +#include "number.h" +#include "tools.h" +#include "polynomial.h" +#include "array.h" +#include "rational.h" + +// init and free for Fields_Table +int init_Fields_Table(Fields_Table* fields){ + init_Int_Array(&((*fields).parameter),FIELDS_SIZE); + 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_Int_Array(&((*fields).fermions),FIELDS_SIZE); + return(0); +} +int free_Fields_Table(Fields_Table fields){ + free_Int_Array(fields.parameter); + free_Int_Array(fields.external); + free_Int_Array(fields.internal); + free_Identities(fields.ids); + free_Symbols(fields.symbols); + free_Int_Array(fields.fermions); + return(0); +} + +// determine field type +int field_type(int index, Fields_Table fields){ + if(int_array_find(abs(index), fields.parameter)>=0){ + return(FIELD_PARAMETER); + } + else if(int_array_find(abs(index), fields.external)>=0){ + return(FIELD_EXTERNAL); + } + 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); + } + + fprintf(stderr,"error: index %d is neither a parameter nor an external or an internal field, nor a symbol\n",index); + exit(-1); +} + +// check whether a field anticommutes +int is_fermion(int index, Fields_Table fields){ + if(int_array_find(abs(index), fields.fermions)>=0){ + return(1); + } + else{ + return(0); + } +} + + +// ------------------ Identities -------------------- + +// allocate memory +int init_Identities(Identities* identities,int size){ + (*identities).lhs=calloc(size,sizeof(Int_Array)); + (*identities).rhs=calloc(size,sizeof(Polynomial)); + (*identities).length=0; + (*identities).memory=size; + return(0); +} + +// free memory +int free_Identities(Identities identities){ + int i; + for(i=0;i=(*output).memory){ + resize_identities(output,2*(*output).memory+1); + } + + // copy and allocate + int_array_cpy(lhs,(*output).lhs+offset); + polynomial_cpy(rhs,(*output).rhs+offset); + // increment length + (*output).length++; + return(0); +} +// append an element to a identities without allocating memory +int identities_append_noinit(Int_Array lhs, Polynomial rhs, Identities* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_identities(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).lhs[offset]=lhs; + (*output).rhs[offset]=rhs; + // increment length + (*output).length++; + return(0); +} + +// concatenate two identitiess +int identities_concat(Identities input, Identities* output){ + int i; + for(i=0;i0){ + at_least_one=0; + + // prevent infinite loops + security++; + if(security>1000000){ + fprintf(stderr,"error: 1000000 iterations reached when trying to simplify a monomial\n"); + exit(-1); + } + + // loop over ids + for(j=0;j=fields.ids.lhs[j].length || (*polynomial).monomials[i].values[k]!=fields.ids.lhs[j].values[l]){ + int_array_append((*polynomial).monomials[i].values[k],&monomial); + // sign correction + if(fermion_count % 2 ==1 && is_fermion((*polynomial).monomials[i].values[k], fields)){ + sign*=-1; + } + } + else{ + // increment fermion_count + if(is_fermion(fields.ids.lhs[j].values[l],fields)){ + fermion_count++; + } + // increment "current" field in the id + l++; + } + } + + num=number_Qprod_ret(quot(sign,1),(*polynomial).nums[i]); + // add extra monomials (if there are more than 1) + for(k=1;k=(*output).memory){ + resize_symbols(output,2*(*output).memory+1); + } + + // copy and allocate + (*output).indices[offset]=index; + polynomial_cpy(expr,(*output).expr+offset); + // increment length + (*output).length++; + return(0); +} +// append an element to a symbols without allocating memory +int symbols_append_noinit(int index, Polynomial expr, Symbols* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_symbols(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).indices[offset]=index; + (*output).expr[offset]=expr; + // increment length + (*output).length++; + return(0); +} + +// concatenate two symbolss +int symbols_concat(Symbols input, Symbols* output){ + int i; + for(i=0;i=(*output).memory){ + resize_groups(output,2*(*output).memory+1); + } + + // copy and allocate + int_array_cpy(indices,(*output).indices+offset); + // increment length + (*output).length++; + return(0); +} +// append an element to a groups without allocating memory +int groups_append_noinit(Int_Array indices, Groups* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_groups(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).indices[offset]=indices; + // increment length + (*output).length++; + return(0); +} + +// concatenate two groupss +int groups_concat(Groups input, Groups* output){ + int i; + for(i=0;i +#include +#include "tools.h" +#include "math.h" +#include "definitions.cpp" +#include "number.h" +#include "array.h" +#include "coefficient.h" + + + +// compute flow numerically, no exponentials +// inputs: flow_equation +// init, niter, tol (the allowed error at each step), ls (whether to display the results in terms of ls), display_mode (what to print) +int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, int niter, long double tol, int display_mode){ + // running coupling contants + RCC rccs=init; + int i,j; + + if(display_mode==DISPLAY_NUMERICAL){ + // print labels + printf("%5s ","n"); + for(j=0;j=0){ + evalcoef(*rccs, flow_equation.coefs[i], new_rccs+i); + // if the new rcc is too small, then ignore it + if(fabs(new_rccs[i]) +#include +#include "definitions.cpp" +#include "rational.h" +#include "istring.h" +#include "coefficient.h" +#include "polynomial.h" +#include "array.h" +#include "number.h" +#include "tools.h" + + +// allocate memory +int init_Grouped_Polynomial(Grouped_Polynomial* gpolynomial, int size){ + (*gpolynomial).coefs=calloc(size,sizeof(Coefficient)); + (*gpolynomial).indices=calloc(size,sizeof(int)); + (*gpolynomial).length=0; + (*gpolynomial).memory=size; + + return(0); +} + +// free memory +int free_Grouped_Polynomial(Grouped_Polynomial gpolynomial){ + int i; + for(i=0;i=(*output).memory){ + resize_grouped_polynomial(output,2*(*output).memory+1); + } + + // copy and allocate + (*output).indices[offset]=index; + coefficient_cpy(coef, (*output).coefs+offset); + //increment length + (*output).length++; + + return(0); +} +// append an element to a grouped_polynomial without allocating memory +int grouped_polynomial_append_noinit(int index, Coefficient coef, Grouped_Polynomial* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_grouped_polynomial(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).indices[offset]=index; + (*output).coefs[offset]=coef; + // increment length + (*output).length++; + return(0); +} + +// concatenate two grouped_polynomials +int grouped_polynomial_concat(Grouped_Polynomial input, Grouped_Polynomial* output){ + int i; + for(i=0;i0){ + // stop if the number of iterations exceeds 100 times the length of the polynomial + if(security >= 100*polynomial.length){ + fprintf(stderr,"error: polynomial could not be grouped in less than %d groupings\n", 100*polynomial.length); + exit(-1); + } + security++; + + // index of the last element + i=remainder.length-1; + + // find entry + if(remainder.monomials[i].length==0){ + // constant + index=-1; + } + else{ + // loop over entries + for(j=0,index=-2;j=0){ + ratio=number_quot_ret(remainder.nums[i],idtable.polynomials[index].nums[pos]); + factor=remainder.factors[i]; + // add to coefficient + denom.index=-1; + denom.power=1; + coefficient_append(factor, ratio, denom, (*grouped_polynomial).coefs+index+1); + + // remove from remainder + free_Int_Array(remainder.monomials[i]); + // do not free factor yet + free_Number(remainder.nums[i]); + remainder.length--; + + // add terms from idtable with minus sign + for(j=0;j=0){ + // a vector in which to store the indices that were masked + init_Int_Array(&mask_tmp_flips,idtable.polynomials[index].length); + + // loop over all monomials in that entry of the idtable + for(j=0;j=0){ + (*dflow).indices[i]=flow_equation.indices[i]+DOFFSET; + } + else{ + (*dflow).indices[i]=flow_equation.indices[i]-DOFFSET; + } + + init_Coefficient((*dflow).coefs+i, COEF_SIZE); + // for each index + for(j=0;j=0){ + int_array_append(DOFFSET + indices.values[j], tmp_coef.factors+k); + } + // constants are offset with -doffset (so that the derivatives of constants also have a negative index) + else{ + int_array_append(-DOFFSET + indices.values[j], tmp_coef.factors+k); + } + } + } + + // add to output + coefficient_concat_noinit(tmp_coef, (*dflow).coefs+i); + } + } + + (*dflow).length=flow_equation.length; + + return(0); +} + + +/* +// derive a flow equation with respect to an index +int flow_equation_deriv(Grouped_Polynomial flow_equation, int index, Grouped_Polynomial* output){ + int i,k; + // temp list of indices + Int_Array factor; + // number of times index was found + int match_count; + coef_denom denom; + // store the computation of the derivative of the constant + int previous_constant_index=0; + Coefficient dC; + Coefficient tmp_coef; + + init_Grouped_Polynomial(output, flow_equation.length); + + // loop over equations + for(k=0;k0){ + coefficient_append_noinit(factor,number_Qprod_ret(quot(match_count,1),flow_equation.coefs[k].nums[i]), flow_equation.coefs[k].denoms[i], (*output).coefs+k); + } + else{ + free_Int_Array(factor); + } + + // derivative of the denominator + if(flow_equation.coefs[k].denoms[i].power>0){ + // check whether the derivative was already computed + if(flow_equation.coefs[k].denoms[i].index!=previous_constant_index){ + // if not first, then free + if(previous_constant_index!=0){ + free_Coefficient(dC); + previous_constant_index=0; + } + init_Coefficient(&dC,COEF_SIZE); + coefficient_deriv_noinit(flow_equation.coefs[intlist_find_err(flow_equation.indices, flow_equation.length, flow_equation.coefs[k].denoms[i].index)], index, &dC); + previous_constant_index=flow_equation.coefs[k].denoms[i].index; + } + + init_Coefficient(&tmp_coef, dC.length); + coefficient_append(flow_equation.coefs[k].factors[i], number_Qprod_ret(quot(-flow_equation.coefs[k].denoms[i].power,1), flow_equation.coefs[k].nums[i]), flow_equation.coefs[k].denoms[i], &tmp_coef); + (tmp_coef.denoms[0].power)++; + + coefficient_prod_chain(dC, &tmp_coef); + + coefficient_concat_noinit(tmp_coef, (*output).coefs+k); + } + } + + // memory safe + if((*output).coefs[k].length>0){ + coefficient_simplify((*output).coefs+k); + } + else{ + // add a trivial element to the coefficient + init_Int_Array(&factor,0); + denom.index=-1; + denom.power=0; + coefficient_append_noinit(factor,number_zero(),denom,(*output).coefs+k); + } + } + + free_Coefficient(dC); + return(0); +} +*/ + + +// print a grouped polynomial +// prepend the indices on the left side with lhs_pre, and those on the right by rhs_pre +int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre, char rhs_pre){ + int i,j; + Char_Array buffer; + int dcount; + + // for each equation + for(i=0;i +#include +#include "array.h" +#include "polynomial.h" + +// allocate memory +int init_Id_Table(Id_Table* idtable,int size){ + (*idtable).indices=calloc(size,sizeof(int)); + (*idtable).polynomials=calloc(size,sizeof(Polynomial)); + (*idtable).length=0; + (*idtable).memory=size; + return(0); +} + +// free memory +int free_Id_Table(Id_Table idtable){ + int i; + for(i=0;i=(*output).memory){ + resize_idtable(output,2*(*output).memory); + } + + // copy and allocate + polynomial_cpy(polynomial,(*output).polynomials+offset); + (*output).indices[offset]=index; + // increment length + (*output).length++; + return(0); +} +// append an element to a idtable without allocating memory +int idtable_append_noinit(int index, Polynomial polynomial, Id_Table* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_idtable(output,2*(*output).memory); + } + + // copy without allocating + (*output).indices[offset]=index; + (*output).polynomials[offset]=polynomial; + // increment length + (*output).length++; + return(0); +} + +// concatenate two idtables +int idtable_concat(Id_Table input, Id_Table* output){ + int i; + for(i=0;i +#include "istring.h" + +char* str_concat(char* ptr, const char* str){ + char* str_ptr=(char*)str; + + while(*str_ptr!='\0'){ + *ptr=*str_ptr; + ptr++; + str_ptr++; + } + *ptr='\0'; + return(ptr); +} + +int str_concat_memorysafe(char** str_out, int pos, const char* str, int* memory){ + char* out_ptr; + + if(str_len((char*)str)+pos>=*memory){ + *memory=*memory+str_len((char*)str)+pos+1; + resize_str(str_out,*memory); + } + out_ptr=*str_out+pos; + str_concat(out_ptr,str); + return(0); +} + +int resize_str(char** out, int memory){ + char* tmp_str=calloc(memory,sizeof(char)); + char* tmp_ptr=tmp_str; + char* out_ptr=*out; + + for(;*out_ptr!='\0';out_ptr++,tmp_ptr++){ + *tmp_ptr=*out_ptr; + } + *tmp_ptr='\0'; + + free(*out); + *out=tmp_str; + return(0); +} + +char* str_addchar(char* ptr, const char c){ + *ptr=c; + ptr++; + *ptr='\0'; + return(ptr); +} + +int str_len(char* str){ + char* ptr=str; + int ret=0; + while(*ptr!='\0'){ret++;ptr++;} + return(ret); +} + +int str_cmp(char* str1, char* str2){ + char* ptr1=str1; + char* ptr2=str2; + while(*ptr1==*ptr2 && *ptr1!='\0' && *ptr2!='\0'){ + ptr1++; + ptr2++; + } + if(*ptr1=='\0' && *ptr2=='\0'){ + return(1); + } + else{ + return(0); + } +} + + +int str_cpy_noalloc(char* input, char* output){ + char* ptr; + char* out_ptr; + + for(ptr=input,out_ptr=output;*ptr!='\0';ptr++,out_ptr++){ + *out_ptr=*ptr; + } + *out_ptr='\0'; + return(0); +} + diff --git a/src/istring.h b/src/istring.h new file mode 100644 index 0000000..5b47b9b --- /dev/null +++ b/src/istring.h @@ -0,0 +1,40 @@ +/* +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. +*/ + +/* +String manipulation +*/ + +#ifndef ISTRING_H +#define ISTRING_H + +// concatenate str at the position pointed to by ptr, and return a pointer +// to the end of the new string +char* str_concat(char* ptr, const char* str); +// concatenate strings and resize them if necessary +int str_concat_memorysafe(char** str_out, int pos, const char* str, int* memory); +// resize a string +int resize_str(char** out, int memory); +// idem with a single character +char* str_addchar(char* ptr, const char c); +// string length +int str_len(char* str); +// compare strings +int str_cmp(char* str1, char* str2); +// copy a string to another without allocating memory +int str_cpy_noalloc(char* input, char* output); + +#endif diff --git a/src/kondo.c b/src/kondo.c new file mode 100644 index 0000000..c86b4b3 --- /dev/null +++ b/src/kondo.c @@ -0,0 +1,1449 @@ +/* +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 "kondo.h" +#include +#include + +#include "idtable.h" +#include "array.h" +#include "number.h" +#include "istring.h" +#include "cli_parser.h" +#include "fields.h" +#include "parse_file.h" +#include "polynomial.h" +#include "definitions.cpp" +#include "rational.h" + +// dimension of A, B and h (must be <10) +#define KONDO_DIM 3 +// number of spin components +#define KONDO_SPIN 2 + +// offsets for indices of A, B and h +// order matters for symbols table +#define KONDO_A_OFFSET 1 +#define KONDO_B_OFFSET 2 +#define KONDO_H_OFFSET 3 + +// parsing modes (from parse_file.c) +#define PP_NULL_MODE 0 +// when reading a factor +#define PP_FACTOR_MODE 1 +// reading a monomial +#define PP_MONOMIAL_MODE 2 +// reading a numerator and denominator +#define PP_NUMBER_MODE 3 +// types of fields +#define PP_FIELD_MODE 6 +#define PP_PARAMETER_MODE 7 +#define PP_EXTERNAL_MODE 8 +#define PP_INTERNAL_MODE 9 +// indices +#define PP_INDEX_MODE 10 +// factors or monomials +#define PP_BRACKET_MODE 11 +// labels +#define PP_LABEL_MODE 12 +// polynomial +#define PP_POLYNOMIAL_MODE 13 +// field scalar product +#define PP_FIELD_SCALAR_MODE 14 +#define PP_FIELD_VECTOR_PROD_MODE 15 + + + +// generate configuration file +int kondo_generate_conf(Str_Array* str_args, int box_count){ + Str_Array new_args; + Fields_Table fields; + Char_Array tmp_str; + int arg_index; + int i; + Char_Array title; + + init_Str_Array(&new_args,8); + + // fields + 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); + 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); + str_array_append_noinit(tmp_str, &new_args); + + // identities + kondo_identities(&tmp_str); + arg_index=find_str_arg("identities", *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_identities(tmp_str, &fields); + str_array_append_noinit(tmp_str, &new_args); + + // groups + kondo_groups(&tmp_str, box_count); + str_array_append_noinit(tmp_str, &new_args); + + + // propagator + arg_index=find_str_arg("propagator", *str_args); + if(arg_index>=0){ + kondo_propagator((*str_args).strs[arg_index], &tmp_str); + str_array_append_noinit(tmp_str, &new_args); + } + + // input polynomial + arg_index=find_str_arg("input_polynomial", *str_args); + if(arg_index>=0){ + kondo_input_polynomial((*str_args).strs[arg_index], &tmp_str, fields, box_count); + str_array_append_noinit(tmp_str, &new_args); + } + + // id table + arg_index=find_str_arg("id_table", *str_args); + if(arg_index>=0){ + kondo_idtable((*str_args).strs[arg_index], &tmp_str, fields); + str_array_append_noinit(tmp_str, &new_args); + } + + // 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 &&\ + str_cmp(title.str, "identities")==0 &&\ + str_cmp(title.str, "propagator")==0 &&\ + str_cmp(title.str, "input_polynomial")==0 &&\ + str_cmp(title.str, "id_table")==0 ){ + + char_array_cpy((*str_args).strs[i], &tmp_str); + str_array_append_noinit(tmp_str, &new_args); + } + free_Char_Array(title); + } + + free_Fields_Table(fields); + free_Str_Array(*str_args); + *str_args=new_args; + + return(0); +} + + +// generate the Kondo fields table +int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fields){ + int i,j; + + init_Char_Array(str_fields,STR_SIZE); + char_array_snprintf(str_fields, "#!fields\n"); + + // external fields + char_array_append_str("x:",str_fields); + for(i=0;i=0 && offset[1]>=0 && index[0]>=0 && index[1]>=0){ + // write indices and num + for(i=0;i0){ + // write num + polynomial_multiply_scalar(tmp_poly, tmp_num); + // replace factor + for(i=0;i0){ + for(i=0;i0){ + for(i=0;i': + if(mode==PP_MONOMIAL_MODE){ + // resolve scalar product + kondo_resolve_scalar_prod(buffer, &scalar_prod_poly, fields); + // add to tmp_poly + if(tmp_poly.length==0){ + polynomial_concat(scalar_prod_poly,&tmp_poly); + } + else{ + polynomial_prod_chain(scalar_prod_poly,&tmp_poly,fields); + } + free_Polynomial(scalar_prod_poly); + + mode=PP_NULL_MODE; + } + break; + + // characters to ignore + case ' ':break; + case '&':break; + case '\n':break; + + // comments + case '#': + comment=1; + break; + + default: + if(mode!=PP_NULL_MODE){ + // write to buffer + buffer_ptr=str_addchar(buffer_ptr,*polynomial_ptr); + } + break; + } + } + } + + // last term + if(tmp_poly.length>0){ + polynomial_multiply_scalar(tmp_poly,tmp_num); + for(i=0;i=0){ + index=*ptr-'0'; + } + else{ + dim=*ptr-'0'; + } + } + } + + // turn B3 into B2 and B4 into B1 + if(offset==KONDO_B_OFFSET){ + switch(index){ + case 3: + index=2; + break; + case 4: + index=1; + break; + } + } + + // h's + if(offset==KONDO_H_OFFSET){ + // external field + init_Int_Array(&monomial,1); + init_Int_Array(&factor,1); + int_array_append(10*(dim+10*offset), &monomial); + polynomial_append_noinit(monomial, factor, number_one(), output); + } + // psi's + else{ + // construct spin indices + for(i=0;i0){ + init_Int_Array(&monomial,1); + init_Int_Array(&factor,1); + + int_array_append(10*(i+10*offset)+index, &monomial); + polynomial_append_noinit(monomial, factor, number_one(), psi+i); + } + } + + // multiply by Pauli matrices + Pauli_matrix(dim+1,&pauli_mat); + for(a=0;a=0){ + kondo_polynomial_vector(offset, index, &poly_vect1, fields); + operation=K_VECT_PROD; + } + break; + + // index + default: + // char to int + index=*ptr-'0'; + } + } + + // final scalar product + if(operation==K_SCALAR_PROD){ + if(offset>=0){ + kondo_polynomial_vector(offset, index, &poly_vect2, fields); + kondo_polynomial_scalar_product(poly_vect1, poly_vect2, output, fields); + } + } + + // free memory + for(i=0;i0){ + init_Int_Array(&monomial,1); + init_Int_Array(&factor,1); + + int_array_append(10*(i+10*offset)+index, &monomial); + polynomial_append_noinit(monomial, factor, number_one(), psi+i); + } + } + + // multiply by Pauli matrices + for(i=0;i=10 || *ptr-'0'<0) && (*ptr!='-')){ + break; + } + } + if(*ptr=='\0'){ + sscanf(str,"%d",&index); + return(index); + } + + for(ptr=str;*ptr!='\0';ptr++){ + switch(*ptr){ + case 'A': + offset=KONDO_A_OFFSET; + break; + case 'a': + offset=KONDO_A_OFFSET; + break; + case 'B': + offset=KONDO_B_OFFSET; + break; + case 'b': + offset=KONDO_B_OFFSET; + break; + case 'h': + offset=KONDO_H_OFFSET; + break; + default: + // set index if dim was already set + if(dim>=0){ + index=*ptr-'0'; + } + else{ + dim=*ptr-'0'; + } + } + } + + if(offset==-1){ + return(-1); + } + // no symbol for h + if(offset==KONDO_H_OFFSET){ + return(10*(10*offset+dim)); + } + else{ + return(100*(10*offset+dim)+index); + } +} diff --git a/src/kondo.h b/src/kondo.h new file mode 100644 index 0000000..23756ad --- /dev/null +++ b/src/kondo.h @@ -0,0 +1,76 @@ +/* +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. +*/ + +/* Generate configuration files specific to the Kondo model */ + +#ifndef KONDO_H +#define KONDO_H + +#include "types.h" + +// generate configuration file +int kondo_generate_conf(Str_Array* str_args, int box_count); + +// generate the Kondo fields table +int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fields); + +// generate Kondo symbols +int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields); +// generate Kondo symbols (older method: one symbol for each scalar product) +int kondo_symbols_scalarprod(Char_Array* str_symbols, int box_count, Fields_Table* fields); + +// generate Kondo groups (groups of independent variables) +int kondo_groups(Char_Array* str_groups, int box_count); + +// generate Kondo identities +int kondo_identities(Char_Array* str_identities); + +// convert the Kondo propagator +int kondo_propagator(Char_Array str_kondo_propagator, Char_Array* str_propagator); + +// read a product of polynomials +int parse_kondo_polynomial_factors(Char_Array str_polynomial, Polynomial* output, Fields_Table fields); + +// convert Kondo input polynomial +int kondo_input_polynomial(Char_Array str_kondo_polynomial, Char_Array* str_polynomial, Fields_Table fields, int box_count); + +// convert the Kondo idtable +int kondo_idtable(Char_Array str_kondo_idtable, Char_Array* str_idtable, Fields_Table fields); + +// read a kondo polynomial and convert it to a polynomial expressed in terms of the fields in the fields table +int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_Table fields); +int parse_kondo_polynomial(Char_Array kondo_polynomial_str, Polynomial* polynomial, Fields_Table fields); + +// read Aij, Bij, hi where i is a space dimension and j is a box index +int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields); +// read a Kondo scalar product +int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields); +// compute a scalar product of polynomial vectors +int kondo_polynomial_scalar_product(Polynomial poly_vect1[3], Polynomial poly_vect2[3], Polynomial* output, Fields_Table fields); +// compute a vector product of polynomial vectors +int kondo_polynomial_vector_product(Polynomial (*poly_vect1)[3], Polynomial poly_vect2[3], Fields_Table fields); +// compute the 3 components of a kondo vector +int kondo_polynomial_vector(int offset, int index, Polynomial (*polys)[3], Fields_Table fields); +// read a scalar product of symbols +int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output); + +// get the offset and index of a monomial term +int get_offset_index(char* str, int* offset, int* index); +// get the offsets and index of a scalar product +int get_offsets_index(char* str, int* offset1, int* offset2, int* index); +// get the index of the symbol corresponding to a given string +int get_symbol_index(char* str); +#endif diff --git a/src/kondo_preprocess.c b/src/kondo_preprocess.c new file mode 100644 index 0000000..3371014 --- /dev/null +++ b/src/kondo_preprocess.c @@ -0,0 +1,126 @@ +/* +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. +*/ + +/* +kondo_preprocess + +Generate configuration files for the Kondo problem + +*/ + + +#include +#include + +#include "definitions.cpp" +#include "kondo.h" +#include "cli_parser.h" +#include "array.h" + + +// read cli arguments +int read_args_kondo_pp(int argc,const char* argv[], Str_Array* str_args, Kondopp_Options* opts); +// print usage message +int print_usage_kondo_pp(); + +int main (int argc, const char* argv[]){ + int i; + // string arguments + Str_Array str_args; + // options + Kondopp_Options opts; + + // read command-line arguments + read_args_kondo_pp(argc,argv,&str_args,&opts); + + kondo_generate_conf(&str_args, 2*opts.dimension); + for(i=0;i\n\n"); + return(0); +} + diff --git a/src/mean.c b/src/mean.c new file mode 100644 index 0000000..aad7a5a --- /dev/null +++ b/src/mean.c @@ -0,0 +1,793 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +As of version 1.0, the mean of a monomial is computed directly +*/ + +#include "mean.h" + +#include +#include +#include +#include "definitions.cpp" +#include "tools.h" +#include "polynomial.h" +#include "rational.h" +#include "array.h" +#include "fields.h" +#include "number.h" + +// mean of a monomial +int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Matrix propagator){ + int sign=1; + // +/- internal fields + Int_Array internal_plus; + Int_Array internal_minus; + + // init out + *out=polynomial_one(); + + // sort first + monomial_sort(monomial, 0, monomial.length-1, fields, &sign); + polynomial_multiply_Qscalar(*out, quot(sign,1)); + // get internals + // (*out).monomials is the first element of out but it only has 1 element + // first, free (*out).monomials[0] + free_Int_Array((*out).monomials[0]); + get_internals(monomial, &internal_plus, &internal_minus, (*out).monomials, fields); + + if(internal_plus.length>0 && internal_minus.length>0){ + mean_internal(internal_plus, internal_minus, out, propagator, fields); + } + + free_Int_Array(internal_plus); + free_Int_Array(internal_minus); + return(0); +} + +// compute the mean of a monomial of internal fields (with split + and -) +int mean_internal(Int_Array internal_plus, Int_Array internal_minus, Polynomial* out, Polynomial_Matrix propagator, Fields_Table fields){ + if(internal_plus.length!=internal_minus.length){ + fprintf(stderr,"error: monomial contains unmatched +/- fields\n"); + exit(-1); + } + int n=internal_minus.length; + // pairing as an array of positions + int* pairing=calloc(n,sizeof(int)); + // specifies which indices are available for pairing + int* mask=calloc(n,sizeof(int)); + // signature of the permutation + int pairing_sign; + // sign from mixing - and + together + int mixing_sign; + Polynomial num; + Polynomial num_summed=polynomial_zero(); + // propagator matrix indices + int* indices_minus=calloc(n,sizeof(int)); + int* indices_plus=calloc(n,sizeof(int)); + int i; + // whether the next pairing exists + int exists_next=0; + + // indices + for(i=0;i=0 && move=0 && move=0 && start=0;i--){ + // move i-th + mask[pairing[i]]=0; + // search for next possible position + for(j=pairing[i]+1;j=n){ + break; + } + } + } + } + // if the next position was found, then don't try to move the previous pairings + break; + } + // if no next position is found, store the pairing outside the array (so the algorithm stops when the first pairing is outside the array) + else{ + pairing[i]=n; + } + } + } + + number_prod_chain(num_summed,outnum); + free_Number(num_summed); + free(pairing); + free(mask); + return(0); +} +*/ + + +// get lists of internal fields from a monomial (separate + and -) +// requires the monomial to be sorted (for the sign to be correct) +int get_internals(Int_Array monomial, Int_Array* internal_plus, Int_Array* internal_minus, Int_Array* others, Fields_Table fields){ + int i; + init_Int_Array(internal_plus, monomial.length); + init_Int_Array(internal_minus, monomial.length); + init_Int_Array(others, monomial.length); + for (i=0;i=0){ + // split +/- fields + if(monomial.values[i]>0){ + int_array_append(monomial.values[i],internal_plus); + } + else{ + int_array_append(monomial.values[i],internal_minus); + } + } + else{ + int_array_append(monomial.values[i], others); + } + } + return(0); +} + + +// compute the mean of a monomial containing symbolic expressions +// keep track of which means were already computed +int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed){ + Int_Array symbol_list; + int i; + int power; + int* current_term; + Polynomial mean_num; + Int_Array tmp_monomial; + Number tmp_num; + Int_Array base_monomial; + int sign; + // whether or not the next term exists + int exists_next=0; + // simplify polynomial periodically + int simplify_freq=1; + Polynomial mean_poly; + + init_Polynomial(output, POLY_SIZE); + + // check whether the mean was already computed + for(i=0;i<(*computed).length;i++){ + if(int_array_cmp((*computed).lhs[i], monomial)==0){ + // write polynomial + polynomial_concat((*computed).rhs[i], output); + return(0); + } + } + + init_Int_Array(&symbol_list, monomial.length); + init_Int_Array(&base_monomial, monomial.length); + + // generate symbols list + for(i=0;i=0 && move=0 && move=0 && start0){ + sign=1; + monomial_sort_groups(monomial, 0, monomial.length-1, fields, groups, &sign); + } + // construct groups and take mean + init_Int_Array(&tmp_monomial, MONOMIAL_SIZE); + for(i=0;i<=monomial.length;i++){ + // new group + if(i0 && next_group!=group) || i==monomial.length){ + mean_symbols(tmp_monomial, &tmp_poly, fields, propagator, groups, computed); + // if zero + if(polynomial_is_zero(tmp_poly)==1){ + // set output to 0 + free_Polynomial(*output); + init_Polynomial(output, 1); + free_Polynomial(tmp_poly); + break; + } + // add to output + if(polynomial_is_zero(*output)==1){ + polynomial_concat(tmp_poly, output); + } + else{ + polynomial_prod_chain_nosimplify(tmp_poly, output, fields); + } + free_Polynomial(tmp_poly); + + // reset tmp_monomial + free_Int_Array(tmp_monomial); + init_Int_Array(&tmp_monomial, MONOMIAL_SIZE); + } + + // add to monomial + if(i +#include + +// pre-compiler definitions +#include "definitions.cpp" + +// various arrays +#include "array.h" + +// list of fields +#include "fields.h" +// numbers +#include "number.h" +// polynomials +#include "polynomial.h" +// list of rccs +#include "idtable.h" +// grouped representation of polynomials +#include "grouped_polynomial.h" +// command line parser +#include "cli_parser.h" +// parse input file +#include "parse_file.h" +// means +#include "mean.h" +// various string operations +#include "istring.h" + +// read cli arguments +int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meankondo_Options* opts); +// print usage message +int print_usage_meankondo(); +// compute flow +int compute_flow(Str_Array str_args, Meankondo_Options opts); +// compute the flow equation +int compute_flow_equation(Polynomial init_poly, Id_Table idtable, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, Grouped_Polynomial* flow_equation); + + +int main (int argc, const char* argv[]){ + // string arguments + Str_Array str_args; + // options + Meankondo_Options opts; + + // read command-line arguments + read_args_meankondo(argc,argv,&str_args,&opts); + + // warning message if representing rational numbers as floats +#ifdef RATIONAL_AS_FLOAT + fprintf(stderr,"info: representing rational numbers using floats\n"); +#endif + + compute_flow(str_args, opts); + + //free memory + free_Str_Array(str_args); + return(0); +} + + +// parse command-line arguments +#define CP_FLAG_THREADS 1 +int read_args_meankondo(int argc,const char* argv[], Str_Array* str_args, Meankondo_Options* opts){ + int i; + // pointers + char* ptr; + // file to read the polynomial from in flow mode + const char* file=""; + // flag that indicates what argument is being read + int flag=0; + // whether a file was specified on the command-line + int exists_file=0; + + + // defaults + // single thread + (*opts).threads=1; + // do not chain + (*opts).chain=0; + + // loop over arguments + for(i=1;i\n\n"); + return(0); +} + + +// compute the renormalization group flow +int compute_flow(Str_Array str_args, Meankondo_Options opts){ + int i; + // index of the entry in the input file + int arg_index; + // header of the entry + Char_Array arg_header; + // list of fields + Fields_Table fields; + // their propagator + Polynomial_Matrix propagator; + // initial polynomial + Polynomial init_poly; + // list of rccs + Id_Table idtable; + // groups of independent fields + Groups groups; + // flow equation + Grouped_Polynomial flow_equation; + + + // 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 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 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); + } + else{ + fprintf(stderr,"error: no input polynomial entry in the configuration file\n"); + exit(-1); + } + + // propagator + arg_index=find_str_arg("propagator", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no propagator entry in the configuration file\n"); + exit(-1); + } + else{ + parse_input_propagator(str_args.strs[arg_index],&propagator, fields); + } + + // 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 groups + arg_index=find_str_arg("groups", str_args); + if(arg_index>=0){ + parse_input_groups(str_args.strs[arg_index],&groups); + } + else{ + init_Groups(&groups, 1); + } + + // flow equation + compute_flow_equation(init_poly, idtable, fields, propagator, groups, opts.threads, &flow_equation); + free_Polynomial(init_poly); + free_Polynomial_Matrix(propagator); + free_Fields_Table(fields); + free_Groups(groups); + + // if chain then print config file + if(opts.chain==1){ + for(i=0;i1){ + polynomial_mean_multithread(&exp_poly, fields, propagator, groups, threads); + } + else{ + polynomial_mean(&exp_poly, fields, propagator, groups); + } + + // 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 new file mode 100644 index 0000000..f45c376 --- /dev/null +++ b/src/meantools.c @@ -0,0 +1,116 @@ +/* +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. +*/ + +/* +meantools + +Utility to perform various operations on flow equations + +*/ + + +#include +#include + +// pre-compiler definitions +#include "definitions.cpp" + +// grouped representation of polynomials +#include "grouped_polynomial.h" +// command line parser +#include "cli_parser.h" +// parse input file +#include "parse_file.h" +// arrays +#include "array.h" +// string functions +#include "istring.h" +// tools +#include "meantools_exp.h" +#include "meantools_deriv.h" +#include "meantools_eval.h" + +#define EXP_COMMAND 1 +#define DERIV_COMMAND 2 +#define EVAL_COMMAND 3 + +// read cli arguments +int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meantools_Options* opts); +// print usage message +int print_usage_meantools(); + + +int main (int argc, const char* argv[]){ + // string arguments + Str_Array str_args; + // options + Meantools_Options opts; + + // read command-line arguments + 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; + } + + //free memory + free_Str_Array(str_args); + return(0); +} + + +// parse command-line arguments +int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meantools_Options* opts){ + + // if there are no arguments + if(argc==1){ + print_usage_meantools(); + 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){ + (*opts).command=DERIV_COMMAND; + tool_deriv_read_args(argc, argv, str_args, opts); + } + else if(str_cmp((char*)argv[1],"eval")==1){ + (*opts).command=EVAL_COMMAND; + tool_eval_read_args(argc, argv, str_args, opts); + } + else{ + print_usage_meantools(); + exit(-1); + } + + return(0); +} + +// print usage message +int print_usage_meantools(){ + printf("\nusage:\n meantools exp \n meantools derive [-d derivatives] -V \n meantools eval -R \n\n"); + return(0); +} + + + diff --git a/src/meantools_deriv.c b/src/meantools_deriv.c new file mode 100644 index 0000000..28d8641 --- /dev/null +++ b/src/meantools_deriv.c @@ -0,0 +1,195 @@ +/* +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_deriv.h" + +#include +#include +#include "parse_file.h" +#include "cli_parser.h" +#include "istring.h" +#include "definitions.cpp" +#include "array.h" +#include "grouped_polynomial.h" + + +#define CP_FLAG_DERIVS 1 +#define CP_FLAG_VARS 2 +// read command line arguments +int tool_deriv_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; + // buffer in which to read the variables + Char_Array buffer; + int i; + char* ptr; + + // defaults + // derive once + (*opts).deriv_derivs=1; + // derive with respect to all variables + (*opts).deriv_vars.length=-1; + + + // loop over arguments + for(i=2;i=0){ + int_array_read(str_args.strs[arg_index],&(opts.deriv_vars)); + } + } + } + + // if variables length is negative then set the variables to all of the available ones + if(opts.deriv_vars.length<0){ + init_Int_Array(&(opts.deriv_vars), flow_equation.length); + for(i=0;i=0){ + int_array_append((j+1)*DOFFSET+variables.values[i], &indices); + } + // constants have a negative index + else{ + int_array_append(-(j+1)*DOFFSET+variables.values[i], &indices); + } + } + + // add to flow equation + grouped_polynomial_concat(dflow, flow_equation_derivs); + } + + if(n>0){ + free_Grouped_Polynomial(dflow); + } + free_Int_Array(indices); + return(0); +} diff --git a/src/meantools_deriv.h b/src/meantools_deriv.h new file mode 100644 index 0000000..4ea36a9 --- /dev/null +++ b/src/meantools_deriv.h @@ -0,0 +1,30 @@ +/* +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. +*/ + +#ifndef MEANTOOLS_DERIV_H +#define MEANTOOLS_DERIV_H + +#include "types.h" + +// read arguments +int tool_deriv_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts); +// derive 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); + +#endif + diff --git a/src/meantools_eval.c b/src/meantools_eval.c new file mode 100644 index 0000000..50fff5b --- /dev/null +++ b/src/meantools_eval.c @@ -0,0 +1,129 @@ +/* +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_eval.h" + +#include +#include +#include "parse_file.h" +#include "cli_parser.h" +#include "grouped_polynomial.h" +#include "array.h" +#include "rcc.h" + + +#define CP_FLAG_RCCS 1 +// read command line arguments +int tool_eval_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 + // mark rccstring so that it can be recognized whether it has been set or not + (*opts).eval_rccstring.length=-1; + + // loop over arguments + for(i=2;i=0){ + char_array_cpy(str_args.strs[arg_index],&(opts.eval_rccstring)); + } + } + } + + // initialize the rccs + prepare_init(flow_equation.indices,flow_equation.length,&rccs); + // read rccs from string + if(opts.eval_rccstring.length!=-1){ + parse_init_cd(opts.eval_rccstring, &rccs); + free_Char_Array(opts.eval_rccstring); + } + + // evaluate + evaleq(&rccs, flow_equation); + + // print + RCC_print(rccs); + + // free memory + free_Grouped_Polynomial(flow_equation); + free_RCC(rccs); + return(0); +} + diff --git a/src/meantools_eval.h b/src/meantools_eval.h new file mode 100644 index 0000000..518049d --- /dev/null +++ b/src/meantools_eval.h @@ -0,0 +1,29 @@ +/* +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. +*/ + +#ifndef MEANTOOLS_EVAL_H +#define MEANTOOLS_EVAL_H + +#include "types.h" + +// read arguments +int tool_eval_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts); +// evaluate a flow equation +int tool_eval(Str_Array str_args, Meantools_Options opts); + +#endif + + diff --git a/src/meantools_exp.c b/src/meantools_exp.c new file mode 100644 index 0000000..1aca928 --- /dev/null +++ b/src/meantools_exp.c @@ -0,0 +1,130 @@ +/* +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 +#include "istring.h" +#include "definitions.cpp" +#include "tools.h" +#include "rational.h" +#include "array.h" + +// init +int init_Number(Number* number, int memory){ + (*number).scalars=calloc(memory,sizeof(Q)); + (*number).base=calloc(memory,sizeof(int)); + (*number).memory=memory; + (*number).length=0; + return(0); +} +int free_Number(Number number){ + free(number.scalars); + free(number.base); + return(0); +} + +// copy +int number_cpy(Number input, Number* output){ + init_Number(output,input.length); + number_cpy_noinit(input,output); + return(0); +} +int number_cpy_noinit(Number input, Number* output){ + int i; + if((*output).memory=(*output).memory){ + number_resize(output,2*(*output).memory); + } + (*output).scalars[(*output).length]=scalar; + (*output).base[(*output).length]=base; + (*output).length++; + // not optimal + number_sort(*output,0,(*output).length-1); + return(0); +} + +// concatenate +int number_concat(Number input, Number* output){ + int i; + int offset=(*output).length; + if((*output).length+input.length>(*output).memory){ + // make it longer than needed by (*output).length (for speed) + number_resize(output,2*(*output).length+input.length); + } + for(i=0;i=0){ + (*inout).scalars[index]=Q_add((*inout).scalars[index], scalar); + } + else{ + number_append(scalar, base, inout); + } + return(0); +} +// create a new number +int number_add(Number x1, Number x2, Number* out){ + number_cpy(x1,out); + number_add_chain(x2,out); + return(0); +} +// return the number +Number number_add_ret(Number x1, Number x2){ + Number out; + number_add(x1,x2,&out); + return(out); +} + +// multiply +int number_prod(Number x1, Number x2, Number* out){ + int i,j; + int div; + Q new_scalar; + int new_base; + init_Number(out, x1.length); + for(i=0;i0){ + (*inout).scalars[i]=Q_inverse((*inout).scalars[i]); + (*inout).scalars[i].denominator*=(*inout).base[i]; + } + 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; + } + else{ + fprintf(stderr,"error: attempting to invert 0\n"); + exit(-1); + } + } + 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 +int number_quot(Number x1, Number x2, Number* output){ + Number inv; + number_inverse(x2, &inv); + number_prod(x1, inv, output); + free_Number(inv); + return(0); +} +int number_quot_chain(Number x1, Number* inout){ + number_inverse_inplace(inout); + number_prod_chain(x1, inout); + return(0); +} +Number number_quot_ret(Number x1, Number x2){ + Number ret; + number_quot(x1, x2, &ret); + return(ret); +} + + +// remove 0's +int number_simplify(Number in, Number* out){ + int i; + init_Number(out, in.length); + for(i=0;i0){ + char_array_snprintf(out," + "); + } + if(simp.length>1 || (simp.length==1 && simp.base[0]!=1)){ + char_array_append('(',out); + } + Q_sprint(simp.scalars[i], out); + if(simp.length>1 || (simp.length==1 && simp.base[0]!=1)){ + char_array_append(')',out); + } + + if(simp.base[i]!=1){ + char_array_snprintf(out,"s{%d}",simp.base[i]); + } + } + + free_Number(simp); + return(0); +} + +// print to stdout +int number_print(Number number){ + Char_Array buffer; + init_Char_Array(&buffer,5*number.length); + number_sprint(number, &buffer); + printf("%s",buffer.str); + return(0); +} + +#define PP_NULL_MODE 0 +#define PP_NUM_MODE 1 +#define PP_SQRT_MODE 2 +// read from a string +int str_to_Number(char* str, Number* number){ + char* ptr; + int mode; + char* buffer=calloc(str_len(str)+1,sizeof(char)); + char* buffer_ptr=buffer; + Q num; + int base; + // whether there are parentheses in the string + int exist_parenthesis=0; + + init_Number(number, NUMBER_SIZE); + + // init num and base + // init to 0 so that if str is empty, then the number is set to 0 + num=quot(0,1); + base=1; + + mode=PP_NULL_MODE; + for(ptr=str;*ptr!='\0';ptr++){ + switch(*ptr){ + // read number + case '(': + if(mode==PP_NULL_MODE){ + // init base + base=1; + mode=PP_NUM_MODE; + exist_parenthesis=1; + } + break; + case ')': + if(mode==PP_NUM_MODE){ + str_to_Q(buffer,&num); + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_NULL_MODE; + } + break; + + // read sqrt + case '{': + // init num + if(num.numerator==0){ + num=quot(1,1); + } + 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; + break; + case '}': + if(mode==PP_SQRT_MODE){ + sscanf(buffer,"%d",&base); + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_NULL_MODE; + } + break; + + // write num + case '+': + if(mode==PP_NULL_MODE){ + number_add_elem(num, base, number); + // re-init num and base + num=quot(0,1); + base=1; + } + break; + + default: + if(mode!=PP_NULL_MODE){ + buffer_ptr=str_addchar(buffer_ptr,*ptr); + } + break; + } + } + + // last step + if(mode==PP_NULL_MODE){ + if(exist_parenthesis==0){ + str_to_Q(str, &num); + } + number_add_elem(num, base, number); + } + + free(buffer); + return(0); +} + +// with Char_Array input +int char_array_to_Number(Char_Array cstr_num,Number* output){ + char* buffer; + char_array_to_str(cstr_num,&buffer); + str_to_Number(buffer, output); + free(buffer); + return(0); +} + + +// -------------------- Number_Matrix --------------------- + +// init +int init_Number_Matrix(Number_Matrix* matrix, int length){ + int i,j; + (*matrix).matrix=calloc(length,sizeof(Number*)); + (*matrix).indices=calloc(length,sizeof(int)); + for(i=0;i +#include + +// pre-compiler definitions +#include "definitions.cpp" + +// rccs +#include "rcc.h" +// grouped representation of polynomials +#include "grouped_polynomial.h" +// command line parser +#include "cli_parser.h" +// parse input file +#include "parse_file.h" +// numerical flow +#include "flow.h" +// arrays +#include "array.h" + +// read cli arguments +int read_args_numkondo(int argc,const char* argv[], Str_Array* str_args, Numkondo_Options* opts); +// print usage message +int print_usage_numkondo(); +// compute flow +int numflow(Str_Array str_args, Numkondo_Options opts); + + +int main (int argc, const char* argv[]){ + // string arguments + Str_Array str_args; + // options + Numkondo_Options opts; + + // read command-line arguments + read_args_numkondo(argc,argv,&str_args,&opts); + + numflow(str_args, opts); + + //free memory + free_Str_Array(str_args); + return(0); +} + + +// parse command-line arguments +#define CP_FLAG_NITER 1 +#define CP_FLAG_TOL 2 +#define CP_FLAG_RCCS 3 +int read_args_numkondo(int argc,const char* argv[], Str_Array* str_args, Numkondo_Options* opts){ + int i; + // pointers + char* ptr; + // file to read the polynomial from in flow mode + const char* file=""; + // flag that indicates what argument is being read + int flag=0; + // 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; + // default niter + (*opts).niter=100; + // default to 0 tolerance + (*opts).tol=0; + // mark rccstring so that it can be recognized whether it has been set or not + (*opts).eval_rccstring.length=-1; + +// loop over arguments +for(i=1;i\n\n"); + return(0); +} + + +// numerical computation of the flow +int numflow(Str_Array str_args, Numkondo_Options opts){ + // index of the entry in the input file + int arg_index; + // list of rccs + Labels labels; + // initial condition + RCC init_cd; + // flow equation + Grouped_Polynomial flow_equation; + + // parse id table + arg_index=find_str_arg("labels", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no labels entry in the configuration file\n"); + exit(-1); + } + else{ + parse_labels(str_args.strs[arg_index], &labels); + } + + // parse flow equation + arg_index=find_str_arg("flow_equation", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no flow equation entry in the configuration file\n"); + exit(-1); + } + else{ + char_array_to_Grouped_Polynomial(str_args.strs[arg_index], &flow_equation); + } + + // initial conditions + // check they were not specified on the command line + if(opts.eval_rccstring.length==-1){ + arg_index=find_str_arg("initial_condition", str_args); + if(arg_index<0){ + fprintf(stderr,"error: no initial condition in the configuration file or on the command line\n"); + exit(-1); + } + else{ + char_array_cpy(str_args.strs[arg_index],&(opts.eval_rccstring)); + } + } + // initialize the rccs + prepare_init(flow_equation.indices,flow_equation.length,&init_cd); + // read rccs from string + if(opts.eval_rccstring.length!=-1){ + parse_init_cd(opts.eval_rccstring, &init_cd); + free_Char_Array(opts.eval_rccstring); + } + + numerical_flow(flow_equation, init_cd, labels, opts.niter, opts.tol, opts.display_mode); + + free_RCC(init_cd); + + // free memory + free_Labels(labels); + free_Grouped_Polynomial(flow_equation); + return(0); +} + diff --git a/src/parse_file.c b/src/parse_file.c new file mode 100644 index 0000000..6054372 --- /dev/null +++ b/src/parse_file.c @@ -0,0 +1,796 @@ +/* +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 "parse_file.h" + +#include +#include +#include "array.h" +#include "fields.h" +#include "rational.h" +#include "number.h" +#include "polynomial.h" +#include "rcc.h" +#include "definitions.cpp" +#include "istring.h" +#include "tools.h" +#include "idtable.h" + + +// parsing modes +#define PP_NULL_MODE 0 +// when reading a factor +#define PP_FACTOR_MODE 1 +// reading a monomial +#define PP_MONOMIAL_MODE 2 +// reading a numerator and denominator +#define PP_NUMBER_MODE 3 +// types of fields +#define PP_FIELD_MODE 6 +#define PP_PARAMETER_MODE 7 +#define PP_EXTERNAL_MODE 8 +#define PP_INTERNAL_MODE 9 +#define PP_FERMIONS_MODE 10 +// indices +#define PP_INDEX_MODE 11 +// factors or monomials +#define PP_BRACKET_MODE 12 +// labels +#define PP_LABEL_MODE 13 +// polynomial +#define PP_POLYNOMIAL_MODE 14 +// group +#define PP_GROUP_MODE 15 + + +// parse fields list +int parse_input_fields(Char_Array str_fields, Fields_Table* fields){ + // buffer + char* buffer=calloc(str_fields.length+1,sizeof(char)); + char* buffer_ptr=buffer; + int i,j; + int mode; + int comment=0; + + // allocate memory + init_Fields_Table(fields); + + // loop over input + mode=PP_NULL_MODE; + for(j=0;j0){ + str_to_Polynomial(buffer, &polynomial); + symbols_append_noinit(index, polynomial, &((*fields).symbols)); + } + + // simplify + for(i=0;i<(*fields).symbols.length;i++){ + polynomial_simplify((*fields).symbols.expr+i, *fields); + } + + free(buffer); + return(0); +} + +// parse groups of independent fields +int parse_input_groups(Char_Array str_groups, Groups* groups){ + // buffer + char* buffer=calloc(str_groups.length+1,sizeof(char)); + char* buffer_ptr=buffer; + int index; + int j; + Int_Array group; + int mode; + int comment=0; + + // alloc + init_Groups(groups, GROUP_SIZE); + + // loop over input + mode=PP_NULL_MODE; + for(j=0;j=0 && index2>=0){ + free_Polynomial((*propagator).matrix[index1][index2]); + str_to_Polynomial(buffer,(*propagator).matrix[index1]+index2); + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_INDEX_MODE; + } + break; + + // comment + case '#': + comment=1; + break; + + default: + buffer_ptr=str_addchar(buffer_ptr,str_propagator.str[j]); + break; + } + } + } + + // last step + if(mode==PP_POLYNOMIAL_MODE){ + free_Polynomial((*propagator).matrix[index1][index2]); + str_to_Polynomial(buffer,(*propagator).matrix[index1]+index2); + } + + free(buffer); + return(0); +} + + +// 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-DOFFSET){ + (*init).values[i]=1.; + } + else{ + (*init).values[i]=0.; + } + + } + return(0); +} diff --git a/src/parse_file.h b/src/parse_file.h new file mode 100644 index 0000000..b51103a --- /dev/null +++ b/src/parse_file.h @@ -0,0 +1,56 @@ +/* +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. +*/ + +/* +Parse the input file +*/ + +#ifndef PARSE_FILE_H +#define PARSE_FILE_H + +#include "types.h" + +// parse fields list +int parse_input_fields(Char_Array str_fields, Fields_Table* fields); + +// parse symbols list +int parse_input_symbols(Char_Array str_symbols, Fields_Table* fields); + +// parse groups of independent fields +int parse_input_groups(Char_Array str_groups, Groups* groups); + +// parse identities between fields +int parse_input_identities(Char_Array str_identities, Fields_Table* fields); + +// parse propagator +int parse_input_propagator(Char_Array str_propagator, Polynomial_Matrix* propagator, Fields_Table fields); + +// parse input polynomial +int parse_input_polynomial(Char_Array str_polynomial, Polynomial* output, Fields_Table fields); + +// parse id table +int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table fields); + +// parse a list of labels +int parse_labels(Char_Array str_labels, Labels* labels); + +// parse the initial condition +int parse_init_cd(Char_Array init_cd, RCC* init); + +// set indices and length of init +int prepare_init(int* indices, int length, RCC* init); + +#endif diff --git a/src/polynomial.c b/src/polynomial.c new file mode 100644 index 0000000..639728a --- /dev/null +++ b/src/polynomial.c @@ -0,0 +1,1263 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "polynomial.h" + +#include +#include +#include "definitions.cpp" +#include "rational.h" +#include "tools.h" +#include "mean.h" +#include "coefficient.h" +#include "istring.h" +#include "array.h" +#include "number.h" +#include "fields.h" + + +// allocate memory +int init_Polynomial(Polynomial* polynomial,int size){ + (*polynomial).monomials=calloc(size,sizeof(Int_Array)); + (*polynomial).factors=calloc(size,sizeof(Int_Array)); + (*polynomial).nums=calloc(size,sizeof(Number)); + (*polynomial).length=0; + (*polynomial).memory=size; + return(0); +} + +// free memory +int free_Polynomial(Polynomial polynomial){ + int i; + for(i=0;i=(*output).memory){ + resize_polynomial(output,2*(*output).memory+1); + } + + // copy and allocate + int_array_cpy(monomial,(*output).monomials+offset); + int_array_cpy(factor,(*output).factors+offset); + number_cpy(num,(*output).nums+offset); + //increment length + (*output).length++; + + return(0); +} +// append an element to a polynomial without allocating memory +int polynomial_append_noinit(Int_Array monomial, Int_Array factor, Number num, Polynomial* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_polynomial(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).monomials[offset]=monomial; + (*output).factors[offset]=factor; + (*output).nums[offset]=num; + // increment length + (*output).length++; + return(0); +} +// noinit, and if there already exists an element with the same monomial and factor, then just add numbers +int polynomial_append_noinit_inplace(Int_Array monomial, Int_Array factor, Number num, Polynomial* output){ + int i; + int foundit=0; + for(i=0;i<(*output).length;i++){ + if(int_array_cmp(monomial, (*output).monomials[i])==0 && int_array_cmp(factor, (*output).factors[i])==0){ + number_add_chain(num,(*output).nums+i); + foundit=1; + free_Number(num); + free_Int_Array(monomial); + free_Int_Array(factor); + break; + } + } + if(foundit==0){ + polynomial_append_noinit(monomial, factor, num, output); + } + return(0); +} + +// concatenate two polynomials +int polynomial_concat(Polynomial input, Polynomial* output){ + int i; + for(i=0;i=input_polynomial.length;i--){ + // try to advance pointer + (current_term[i])++; + // if the step fails, keep on looping, if not, set all the following + // pointers to the latest position + if(current_term[i]1){ + if(power>33){ + fprintf(stderr,"error: trying to take a power of a polynomial that is too high (>33)\n"); + exit(-1); + } + // next power + polynomial_prod(input_polynomial,previous_power,&tmp_poly, fields); + + // free + free_Polynomial(previous_power); + } + else{ + polynomial_cpy(input_polynomial,&tmp_poly); + } + + // if the power is high enough that V^p=0, then stop + if(tmp_poly.length==0){ + free_Polynomial(tmp_poly); + break; + } + + // copy for next power + polynomial_cpy(tmp_poly,&previous_power); + + // 1/p! + polynomial_multiply_Qscalar(tmp_poly,quot(1,factorial(power))); + // append tmp to the output + polynomial_concat_noinit(tmp_poly,output); + + // increase power + power++; + } + + return(0); +} + + +// log(1+W) +int polynomial_logarithm(Polynomial input_polynomial,Polynomial* output, Fields_Table fields){ + // a buffer for the result of a given power + Polynomial tmp_poly; + // a buffer for the previous power + Polynomial previous_power; + // power + int power=1; + + // allocate memory + init_Polynomial(output,POLY_SIZE); + + while(1){ + if(power>1){ + // next power + polynomial_prod(input_polynomial,previous_power,&tmp_poly, fields); + + // free + free_Polynomial(previous_power); + } + else{ + polynomial_cpy(input_polynomial,&tmp_poly); + } + + // if the power is high enough that V^p=0, then stop + if(tmp_poly.length==0){ + free_Polynomial(tmp_poly); + break; + } + + // copy for next power + polynomial_cpy(tmp_poly,&previous_power); + + // (-1)^{p-1}/p + polynomial_multiply_Qscalar(tmp_poly,quot(ipower(-1,power-1),power)); + // append tmp to the output + polynomial_concat_noinit(tmp_poly,output); + + // increase power + power++; + } + + return(0); +} + +// check whether a monomial vanishes +int check_monomial(Int_Array monomial, Fields_Table fields){ + int i,j; + for(i=0;i0){ + match++; + } + else if(monomial.values[i]<0){ + match--; + } + } + + // check for repetitions + if(is_fermion(monomial.values[i], fields)==1){ + for(j=i+1;j0){ + match_internals++; + } + else if((*polynomial).monomials[i].values[j]<0){ + match_internals--; + } + } + // don't remove a term containing symbols + else if(type==FIELD_SYMBOL){ + match_internals=0; + break; + } + } + if(match_internals==0){ + polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], &output); + } + } + + free_Polynomial(*polynomial); + *polynomial=output; + return(0); +} + + +// denominator of a multinomal factor: m1!m2!... +// requires terms to be sorted +int multinomial(int power,int* terms){ + int multiple=1; + int ret=1; + int i; + // the number of numbers to be multiplied in the multinomial is + // equal to power-1 (the first is 1) + for(i=1;i=(*polynomial).length-1 || monomial_cmp!=0 || factor_cmp!=0 ){ + // check that the polynomial is not trivial + if(number_is_zero(new_num)!=1){ + polynomial_append((*polynomial).monomials[i],(*polynomial).factors[i],new_num,&output); + } + + // reset new numerical factor + free_Number(new_num); + init_Number(&new_num,NUMBER_SIZE); + } + } + + free_Number(new_num); + free_Polynomial(*polynomial); + *polynomial=output; + return(0); +} + +// sort a polynomial +// requires the monomials and factors to be sorted +int polynomial_sort(Polynomial* polynomial, int begin, int end){ + int i; + int index; + int monomial_cmp; + int factor_cmp; + + // the pivot: middle of the array + int pivot=(begin+end)/2; + // if the array is non trivial + if(begin type2){ + return(1); + } + + if(monomial.values[pos1] < monomial.values[pos2]){ + return(-1); + } + else if (monomial.values[pos1] > monomial.values[pos2]){ + return(1); + } + + return(0); +} + +// exchange two fields (with sign) +int exchange_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields, int* sign){ + int tmp=monomial.values[pos1]; + int f1,f2; + int i; + + monomial.values[pos1]=monomial.values[pos2]; + monomial.values[pos2]=tmp; + + // sign change + // only if the exchange is not trivial + if(pos1!=pos2){ + f1=is_fermion(monomial.values[pos1],fields); + f2=is_fermion(monomial.values[pos2],fields); + // if both Fermions then sign + if(f1==1 && f2==1){ + *sign*=-1; + } + // if only one of them is a Fermion, then count the number of Fermions between them + else if(f1==1 || f2==1){ + for(i=min(pos1,pos2)+1;i group2){ + return(1); + } + + int type1=field_type(monomial.values[pos1], fields); + int type2=field_type(monomial.values[pos2],fields); + + if(type1 < type2){ + return(-1); + } + else if(type1 > type2){ + return(1); + } + + if(monomial.values[pos1] < monomial.values[pos2]){ + return(-1); + } + else if (monomial.values[pos1] > monomial.values[pos2]){ + return(1); + } + + return(0); +} + + +// convert and idtable to a polynomial +int idtable_to_polynomial(Id_Table idtable, Polynomial* polynomial){ + int i,j; + int start=0; + + // allocate memory + init_Polynomial(polynomial,POLY_SIZE); + + for(i=0;i0){ + // initialize coef to store the product of factors + init_Coefficient(&coef, POLY_SIZE); + + // first term + index=intlist_find_err(equations.indices, equations.length, (*polynomial).factors[i].values[0]); + if(index>=0){ + coefficient_cpy_noinit(equations.coefs[index],&coef); + } + + // other terms + for(j=1;j<(*polynomial).factors[i].length;j++){ + index=intlist_find_err(equations.indices, equations.length, (*polynomial).factors[i].values[j]); + if(index>=0){ + coefficient_prod_chain(equations.coefs[index],&coef); + } + } + + // new polynomial terms + for(j=0;j +#include +#include "istring.h" +#include "array.h" +#include "math.h" + +Q quot(long double p, long double q){ + Q ret; + if(q==0){ + fprintf(stderr,"error: %Lf/%Lf is ill defined\n",p,q); + exit(-1); + } + ret.numerator=p; + ret.denominator=q; + return(ret); +} + +// add +Q Q_add(Q x1,Q x2){ + Q ret; + ret.denominator=lcm(x1.denominator,x2.denominator); + ret.numerator=x1.numerator*(ret.denominator/x1.denominator)+x2.numerator*(ret.denominator/x2.denominator); + return(Q_simplify(ret)); +} +//multiply +Q Q_prod(Q x1,Q x2){ + return(Q_simplify(quot(x1.numerator*x2.numerator,x1.denominator*x2.denominator))); +} +// inverse +Q Q_inverse(Q x1){ + if(x1.numerator>0){ + return(quot(x1.denominator,x1.numerator)); + } + else if(x1.numerator<0){ + return(quot(-x1.denominator,-x1.numerator)); + } + else{ + fprintf(stderr,"error: attempting to invert %Lf/%Lf\n",x1.numerator, x1.denominator); + exit(-1); + } + +} +// quotient +Q Q_quot(Q x1, Q x2){ + if(x2.numerator>0){ + return(Q_simplify(quot(x1.numerator*x2.denominator,x1.denominator*x2.numerator))); + } + else if(x2.numerator<0){ + return(Q_simplify(quot(-x1.numerator*x2.denominator,-x1.denominator*x2.numerator))); + } + else{ + fprintf(stderr,"error: attempting to invert %Lf/%Lf\n",x2.numerator, x2.denominator); + exit(-1); + } +} + +// compare +int Q_cmp(Q x1, Q x2){ + Q quo=Q_quot(x1,x2); + if(quo.numerator > quo.denominator){ + return(1); + } + else if(quo.numerator < quo.denominator){ + return(-1); + } + else{ + return(0); + } +} + +// simplify +Q Q_simplify(Q x){ + return(quot(x.numerator/gcd(x.numerator,x.denominator),x.denominator/gcd(x.numerator,x.denominator))); +} +//simplify in place +int Q_simplify_inplace(Q* x){ + Q ret=Q_simplify(*x); + *x=ret; + return(0); +} + +// greatest common divisor +long double gcd(long double x, long double y){ + long double ax=fabsl(x); + long double ay=fabsl(y); + int security=0; + while(ax!=0 && ay!=0){ + if(ax>ay){ax=modld(ax,ay);} + else{ay=modld(ay,ax);} + + security++; + if(security>1000000){ + fprintf(stderr,"error: could not compute gcd( %Lf , %Lf ) after %d tries\n",x,y,security); + exit(-1); + } + } + // if both are negative, gcd should be negative (useful for simplification of fractions and product of square roots) + if(x<0 && y<0){ + ax*=-1; + ay*=-1; + } + if(fabsl(ay)>fabsl(ax)){return(ay);} + else{return(ax);} +} + +long double modld(long double x, long double m){ + long double q=floorl(x/m); + return(x-q*m); +} + +// least common multiple +long double lcm(long double x,long double y){ + if(x!=0 || y!=0){ + return((x*y)/gcd(x,y)); + } + else{ + return(0); + } +} + +// approximate value as double +double Q_double_value(Q q){ + return(1.0*q.numerator/q.denominator); +} + + +// print to string +int Q_sprint(Q num, Char_Array* out){ + if(num.denominator!=1){ + char_array_snprintf(out,"%Lf/%Lf", num.numerator,num.denominator); + } + else{ + char_array_snprintf(out,"%Lf",num.numerator); + } + + return(0); +} + +#define PP_NUMERATOR_MODE 1 +#define PP_DENOMINATOR_MODE 2 +// read from a string +int str_to_Q(char* str, Q* num){ + char* ptr; + int mode; + char* buffer=calloc(str_len(str)+1,sizeof(char)); + char* buffer_ptr=buffer; + + *num=quot(0,1); + + mode=PP_NUMERATOR_MODE; + for(ptr=str;*ptr!='\0';ptr++){ + if(*ptr=='/'){ + sscanf(buffer,"%Lf",&((*num).numerator)); + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_DENOMINATOR_MODE; + } + else{ + buffer_ptr=str_addchar(buffer_ptr,*ptr); + } + } + + // last step + if(mode==PP_NUMERATOR_MODE){ + sscanf(buffer,"%Lf",&((*num).numerator)); + } + else if(mode==PP_DENOMINATOR_MODE){ + sscanf(buffer,"%Lf",&((*num).denominator)); + } + + free(buffer); + return(0); +} + +#else +int null_declaration_so_that_the_compiler_does_not_complain; +#endif diff --git a/src/rational_float.h b/src/rational_float.h new file mode 100644 index 0000000..931b5ec --- /dev/null +++ b/src/rational_float.h @@ -0,0 +1,64 @@ +/* +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. +*/ + +/* + Rational numbers + uses long doubles to represent integers (to avoid overflow) +*/ + +#ifdef RATIONAL_AS_FLOAT + +#ifndef RATIONAL_H +#define RATIONAL_H + +#include "types.h" + +// Q from int/int +Q quot(long double p, long double q); + +// add +Q Q_add(Q x1,Q x2); +//multiply +Q Q_prod(Q x1,Q x2); +// inverse +Q Q_inverse(Q x1); +// quotient +Q Q_quot(Q x1, Q x2); + +// compare +int Q_cmp(Q x1, Q x2); + +// simplify +Q Q_simplify(Q x); +//simplify in place +int Q_simplify_inplace(Q* x); + +// greatest common divisor +long double gcd(long double x,long double y); +long double modld(long double x, long double m); +// least common multiple +long double lcm(long double x,long double y); + +// approximate value as double +double Q_double_value(Q q); + +// print to string +int Q_sprint(Q num, Char_Array* out); +// read from a string +int str_to_Q(char* str, Q* num); + +#endif +#endif diff --git a/src/rational_int.c b/src/rational_int.c new file mode 100644 index 0000000..ec601a8 --- /dev/null +++ b/src/rational_int.c @@ -0,0 +1,190 @@ +/* +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. +*/ + +#ifndef RATIONAL_AS_FLOAT + +#include "rational_int.h" +#include +#include +#include "istring.h" +#include "array.h" + +Q quot(long int p, long int q){ + Q ret; + if(q==0){ + fprintf(stderr,"error: %ld/%ld is ill defined\n",p,q); + exit(-1); + } + ret.numerator=p; + ret.denominator=q; + return(ret); +} + +// add +Q Q_add(Q x1,Q x2){ + Q ret; + ret.denominator=lcm(x1.denominator,x2.denominator); + ret.numerator=x1.numerator*(ret.denominator/x1.denominator)+x2.numerator*(ret.denominator/x2.denominator); + return(Q_simplify(ret)); +} +//multiply +Q Q_prod(Q x1,Q x2){ + return(Q_simplify(quot(x1.numerator*x2.numerator,x1.denominator*x2.denominator))); +} +// inverse +Q Q_inverse(Q x1){ + if(x1.numerator>0){ + return(quot(x1.denominator,x1.numerator)); + } + else if(x1.numerator<0){ + return(quot(-x1.denominator,-x1.numerator)); + } + else{ + fprintf(stderr,"error: attempting to invert %ld/%ld\n",x1.numerator, x1.denominator); + exit(-1); + } + +} +// quotient +Q Q_quot(Q x1, Q x2){ + if(x2.numerator>0){ + return(Q_simplify(quot(x1.numerator*x2.denominator,x1.denominator*x2.numerator))); + } + else if(x2.numerator<0){ + return(Q_simplify(quot(-x1.numerator*x2.denominator,-x1.denominator*x2.numerator))); + } + else{ + fprintf(stderr,"error: attempting to invert %ld/%ld\n",x2.numerator, x2.denominator); + exit(-1); + } +} + +// compare +int Q_cmp(Q x1, Q x2){ + Q quo=Q_quot(x1,x2); + if(quo.numerator > quo.denominator){ + return(1); + } + else if(quo.numerator < quo.denominator){ + return(-1); + } + else{ + return(0); + } +} + +// simplify +Q Q_simplify(Q x){ + return(quot(x.numerator/gcd(x.numerator,x.denominator),x.denominator/gcd(x.numerator,x.denominator))); +} +//simplify in place +int Q_simplify_inplace(Q* x){ + Q ret=Q_simplify(*x); + *x=ret; + return(0); +} + +// greatest common divisor +long int gcd(long int x, long int y){ + long int ax=labs(x); + long int ay=labs(y); + int security=0; + while(ax!=0 && ay!=0){ + if(ax>ay){ax=ax%ay;} + else{ay=ay%ax;} + + security++; + if(security>1000000){ + fprintf(stderr,"error: could not compute gcd( %ld , %ld ) after %d tries\n",x,y,security); + exit(-1); + } + } + // if both are negative, gcd should be negative (useful for simplification of fractions and product of square roots) + if(x<0 && y<0){ + ax*=-1; + ay*=-1; + } + if(labs(ay)>labs(ax)){return(ay);} + else{return(ax);} +} + +// least common multiple +long int lcm(long int x,long int y){ + if(x!=0 || y!=0){ + return((x*y)/gcd(x,y)); + } + else{ + return(0); + } +} + +// approximate value as double +double Q_double_value(Q q){ + return(1.0*q.numerator/q.denominator); +} + + +// print to string +int Q_sprint(Q num, Char_Array* out){ + if(num.denominator!=1){ + char_array_snprintf(out,"%ld/%ld", num.numerator,num.denominator); + } + else{ + char_array_snprintf(out,"%ld",num.numerator); + } + + return(0); +} + +#define PP_NUMERATOR_MODE 1 +#define PP_DENOMINATOR_MODE 2 +// read from a string +int str_to_Q(char* str, Q* num){ + char* ptr; + int mode; + char* buffer=calloc(str_len(str)+1,sizeof(char)); + char* buffer_ptr=buffer; + + *num=quot(0,1); + + mode=PP_NUMERATOR_MODE; + for(ptr=str;*ptr!='\0';ptr++){ + if(*ptr=='/'){ + sscanf(buffer,"%ld",&((*num).numerator)); + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_DENOMINATOR_MODE; + } + else{ + buffer_ptr=str_addchar(buffer_ptr,*ptr); + } + } + + // last step + if(mode==PP_NUMERATOR_MODE){ + sscanf(buffer,"%ld",&((*num).numerator)); + } + else if(mode==PP_DENOMINATOR_MODE){ + sscanf(buffer,"%ld",&((*num).denominator)); + } + + free(buffer); + return(0); +} + +#else +int null_declaration_so_that_the_compiler_does_not_complain; +#endif diff --git a/src/rational_int.h b/src/rational_int.h new file mode 100644 index 0000000..a9f164d --- /dev/null +++ b/src/rational_int.h @@ -0,0 +1,59 @@ +/* +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. +*/ + +/* Rational numbers */ + +#ifndef RATIONAL_AS_FLOAT +#ifndef RATIONAL_H +#define RATIONAL_H + +#include "types.h" + +// Q from int/int +Q quot(long int p, long int q); + +// add +Q Q_add(Q x1,Q x2); +//multiply +Q Q_prod(Q x1,Q x2); +// inverse +Q Q_inverse(Q x1); +// quotient +Q Q_quot(Q x1, Q x2); + +// compare +int Q_cmp(Q x1, Q x2); + +// simplify +Q Q_simplify(Q x); +//simplify in place +int Q_simplify_inplace(Q* x); + +// greatest common divisor +long int gcd(long int x,long int y); +// least common multiple +long int lcm(long int x,long int y); + +// approximate value as double +double Q_double_value(Q q); + +// print to string +int Q_sprint(Q num, Char_Array* out); +// read from a string +int str_to_Q(char* str, Q* num); + +#endif +#endif diff --git a/src/rcc.c b/src/rcc.c new file mode 100644 index 0000000..484e20a --- /dev/null +++ b/src/rcc.c @@ -0,0 +1,95 @@ +/* +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 "rcc.h" +#include +#include +#include "array.h" + +// init +int init_RCC(RCC* rccs, int size){ + (*rccs).values=calloc(size,sizeof(long double)); + (*rccs).indices=calloc(size,sizeof(int)); + (*rccs).length=size; + return(0); +} + +int free_RCC(RCC rccs){ + free(rccs.values); + free(rccs.indices); + return(0); +} + +// set a given element of an rcc +int RCC_set_elem(long double value, int index, RCC* rcc, int pos){ + (*rcc).values[pos]=value; + (*rcc).indices[pos]=index; + return(0); +} + +int RCC_cpy(RCC input,RCC* output){ + int i; + + init_RCC(output,input.length); + for(i=0;i +#include + +int factorial(int n){ + int i; + int res=1; + for(i=2;i<=n;i++){ + res*=i; + } + return(res); +} + +int ipower(int n, int p){ + int i; + int res=1; + for(i=1;i<=p;i++){ + res*=n; + } + return(res); +} + +// power of a double +long double dpower(long double x, int p){ + int i; + long double res=1.; + if(p>=0){ + for(i=1;i<=p;i++){ + res*=x; + } + } + else{ + for(i=1;i<=-p;i++){ + res/=x; + } + } + + return(res); +} + + +// find an element in a list whose first element is its size +int list_find(int* list, int x){ + int i; + for(i=1;i<=*list;i++){ + if(list[i]==x){return(i);} + } + return(0); +} + +// find an element in a list whose first element is not its number of elements (skips first element) +int unlist_find(int* list, int size, int x){ + int i; + for(i=1;i=x2){ + return(x1); + } + else{ + return(x2); + } +} +int min(int x1, int x2){ + if(x1<=x2){ + return(x1); + } + else{ + return(x2); + } +} diff --git a/src/tools.h b/src/tools.h new file mode 100644 index 0000000..e5f319f --- /dev/null +++ b/src/tools.h @@ -0,0 +1,49 @@ +/* +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. +*/ + +/* +Various tools +*/ + +#ifndef TOOLS_H +#define TOOLS_H + +#include "types.h" + +int factorial(int n); +int ipower(int n, int p); +long double dpower(long double x, int p); + +// find an element in a list whose first element is its size +int list_find(int* list, int x); +// compare two lists +int list_cmp(int* list1, int* list2); + +// find an element in a list whose first element is not its number of elements (skips first element) +int unlist_find(int* list, int size, int x); +// no error +int unlist_find_noerr(int* list, int size, int x); + +// find an element in a list (checks first element) +int intlist_find(int* list, int size, int x); +// with error +int intlist_find_err(int* list, int size, int x); + +// max and min +int max(int x1, int x2); +int min(int x1, int x2); + +#endif diff --git a/src/types.h b/src/types.h new file mode 100644 index 0000000..0452ebc --- /dev/null +++ b/src/types.h @@ -0,0 +1,219 @@ +/* +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. +*/ + +/* + typedefs used by meankondo +*/ + +#ifndef TYPES_H +#define TYPES_H + + +// rational number +typedef struct Q{ +#ifndef RATIONAL_AS_FLOAT + long int numerator; + long int denominator; +#else + long double numerator; + long double denominator; +#endif +} Q; + +// numbers +typedef struct Number{ + Q* scalars; + int* base; + int length; + int memory; +} Number; + +// matrix +typedef struct Number_Matrix{ + Number** matrix; + int* indices; + int length; +} Number_Matrix; + +// list of integers +typedef struct Int_Array{ + int* values; + int length; + int memory; +} Int_Array; + +// string +typedef struct Char_Array{ + char* str; + int length; + int memory; +} Char_Array; + +// string array +typedef struct Str_Array{ + Char_Array* strs; + int length; + int memory; +} Str_Array; + +// polynomial +typedef struct Polynomial{ + Int_Array* monomials; + Int_Array* factors; + Number* nums; + int length; + int memory; +} Polynomial; + +// matrix of polynomial +typedef struct Polynomial_Matrix{ + Polynomial** matrix; + int* indices; + int length; +} Polynomial_Matrix; + +// denominators in coefficients (index of the denominator and the power to which it is raised) +typedef struct coef_denom{ + int index; + int power; +} coef_denom; + +// coefficient +typedef struct Coefficient{ + Int_Array* factors; + Number* nums; + coef_denom* denoms; + int length; + int memory; +} Coefficient; + +// grouped polynomial +typedef struct Grouped_Polynomial{ + Coefficient* coefs; + int* indices; + int length; + int memory; +} Grouped_Polynomial; + +// rcc +typedef struct RCC{ + long double* values; + int* indices; + int length; +} RCC; + +// identities between fields +typedef struct Identities{ + // left hand sides + Int_Array* lhs; + // right hand sides + Polynomial* rhs; + int length; + int memory; +} Identities; + +// symbolic expressions +typedef struct Symbols{ + int* indices; + Polynomial* expr; + int length; + int memory; +} Symbols; + +// groups of independent fields +typedef struct Groups{ + Int_Array* indices; + int length; + int memory; +} Groups; + +// collection of fields +typedef struct Fields_Table{ + // constant parameters + Int_Array parameter; + // anticommuting external fields + Int_Array external; + // anticommuting internal fields + Int_Array internal; + // identities between fields + Identities ids; + // symbolic expressions (commuting) + Symbols symbols; + // list of anti-commuting variables (fields or symbols) + Int_Array fermions; +} Fields_Table; + +// index labels +typedef struct Labels{ + Char_Array* labels; + int* indices; + int length; + int memory; +} Labels; + +// id table +typedef struct Id_Table{ + int* indices; + Polynomial* polynomials; + int length; + int memory; +} Id_Table; + +/* +// polynomial scalar and vectors +typedef struct Polynomial_Scalar{ + Coefficient coef; + int* indices; + int length; +} Polynomial_Scalar; +typedef struct Polynomial_Vector{ + Coefficient* coefv; + int* indices; + int length; +} Polynomial_Vector; +typedef struct Polynomial_Matrix{ + Coefficient** coefm; + int* indices; + int length; +} Polynomial_Matrix; +*/ + + +// command line options +typedef struct Meankondo_Options{ + int threads; + int chain; +} Meankondo_Options; + +typedef struct Numkondo_Options{ + int display_mode; + int niter; + long double tol; + Char_Array eval_rccstring; +} Numkondo_Options; + +typedef struct Meantools_Options{ + int command; + int deriv_derivs; + Int_Array deriv_vars; + Char_Array eval_rccstring; +} Meantools_Options; + +typedef struct Kondopp_Options{ + int dimension; +} Kondopp_Options; + +#endif -- cgit v1.2.3-54-g00ecf