From 469bdc80712dbf9c12562059dc4594620b59a076 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Wed, 7 Oct 2015 12:51:41 +0000 Subject: Support MPFR floats in numkondo Remove '-D' option (error tolerance) in numkondo --- src/coefficient.c | 47 ++++++++++++++++- src/coefficient.h | 2 + src/definitions.cpp | 2 +- src/flow.c | 24 ++------- src/flow.h | 8 ++- src/flow_mpfr.c | 128 +++++++++++++++++++++++++++++++++++++++++++++++ src/flow_mpfr.h | 32 ++++++++++++ src/grouped_polynomial.c | 30 ++++++++++- src/grouped_polynomial.h | 2 + src/meantools.c | 2 +- src/meantools_eval.c | 64 +++++++++++++++++++++--- src/number.c | 28 +++++++++++ src/number.h | 2 + src/numkondo.c | 80 ++++++++++++++++++++++------- src/parse_file.c | 37 ++++++++++++-- src/parse_file.h | 6 ++- src/rational_float.c | 25 +++++---- src/rational_float.h | 2 + src/rational_int.c | 14 ++++++ src/rational_int.h | 2 + src/rcc_mpfr.c | 124 +++++++++++++++++++++++++++++++++++++++++++++ src/rcc_mpfr.h | 43 ++++++++++++++++ src/types.h | 32 +++++------- 23 files changed, 645 insertions(+), 91 deletions(-) create mode 100644 src/flow_mpfr.c create mode 100644 src/flow_mpfr.h create mode 100644 src/rcc_mpfr.c create mode 100644 src/rcc_mpfr.h (limited to 'src') diff --git a/src/coefficient.c b/src/coefficient.c index c26b668..d4cb6eb 100644 --- a/src/coefficient.c +++ b/src/coefficient.c @@ -18,6 +18,10 @@ limitations under the License. #include #include +#include +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include #include "definitions.cpp" #include "rational.h" #include "istring.h" @@ -721,7 +725,7 @@ int evalcoef(RCC rccs, Coefficient coef, long double* out){ int i,j; long double num_factor; - *out=0; + *out=0.; // for each monomial for(i=0;i0){ + mpfr_pow_si(y, rccs.values[intlist_find_err(rccs.indices,rccs.length,coef.denoms[i].index)], coef.denoms[i].power, MPFR_RNDN); + mpfr_div(x, num_factor, y, MPFR_RNDN); + mpfr_set(num_factor, x, MPFR_RNDN); + } + + number_mpfr_val(Z, coef.nums[i]); + mpfr_mul(x, num_factor, Z, MPFR_RNDN); + mpfr_add(y, x, out, MPFR_RNDN); + mpfr_set(out, y, MPFR_RNDN); + mpfr_clear(Z); + } + + // free numbers + mpfr_clears(num_factor, x, y, (mpfr_ptr)NULL); + + return(0); +} diff --git a/src/coefficient.h b/src/coefficient.h index a7a151c..43dd313 100644 --- a/src/coefficient.h +++ b/src/coefficient.h @@ -74,5 +74,7 @@ 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); +// evaluate a coefficient on a vector (using mpfr floats) +int evalcoef_mpfr(RCC_mpfr rccs, Coefficient coef, mpfr_t out); #endif diff --git a/src/definitions.cpp b/src/definitions.cpp index 25715e5..30f021f 100644 --- a/src/definitions.cpp +++ b/src/definitions.cpp @@ -17,7 +17,7 @@ limitations under the License. #ifndef DEFINITIONS_GCC #define DEFINITIONS_GCC -#define VERSION "1.3.1" +#define VERSION "1.4" // number of entries in a configuration file #define ARG_COUNT 10 diff --git a/src/flow.c b/src/flow.c index d271e44..b294bb0 100644 --- a/src/flow.c +++ b/src/flow.c @@ -24,13 +24,11 @@ limitations under the License. #include "number.h" #include "array.h" #include "coefficient.h" - +#include "rcc.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){ +int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, int niter, int display_mode){ // running coupling contants RCC rccs=init; int i,j; @@ -53,7 +51,7 @@ int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, in for(i=0;i=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 +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +// define MPFR_USE_FILE to enable the use of mpfr_printf +#define MPFR_USE_FILE +#include +#include "tools.h" +#include "math.h" +#include "definitions.cpp" +#include "number.h" +#include "array.h" +#include "coefficient.h" +#include "flow.h" +#include "rcc_mpfr.h" + + + +// compute flow numerically +int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels labels, int niter, int display_mode){ + // running coupling contants + RCC_mpfr rccs=init; + int i,j; + + if(display_mode==DISPLAY_NUMERICAL){ + // print labels + printf("%5s ","n"); + for(j=0;j=0){ + evalcoef_mpfr(*rccs, flow_equation.coefs[i], res[i]); + } + } + + // set new rccs + for(i=0;i\n meantools derive [-d derivatives] -V \n meantools eval -R \n\n"); + printf("\nusage:\n meantools exp [config_file]\n meantools derive [-d derivatives] [-V variables] [-C] [config_file]\n meantools eval [-R values] [-P precision] [-E max_exponent] [config_file]\n\n"); return(0); } diff --git a/src/meantools_eval.c b/src/meantools_eval.c index 50fff5b..4f166cd 100644 --- a/src/meantools_eval.c +++ b/src/meantools_eval.c @@ -18,16 +18,22 @@ limitations under the License. #include #include +#include #include "parse_file.h" #include "cli_parser.h" #include "grouped_polynomial.h" #include "array.h" #include "rcc.h" +#include "rcc_mpfr.h" #define CP_FLAG_RCCS 1 +#define CP_FLAG_MPFR_PREC 2 +#define CP_FLAG_MPFR_EXP 3 // read command line arguments int tool_eval_read_args(int argc, const char* argv[], Str_Array* str_args, Meantools_Options* opts){ + // temporary long int + long int tmp_lint; // file to read the polynomial from in flow mode const char* file=""; // whether a file was specified on the command-line @@ -40,6 +46,9 @@ int tool_eval_read_args(int argc, const char* argv[], Str_Array* str_args, Meant // defaults // mark rccstring so that it can be recognized whether it has been set or not (*opts).eval_rccstring.length=-1; + // no mpfr + (*opts).mpfr_prec=0; + (*opts).mpfr_emax=0; // loop over arguments for(i=2;i #include #include +#include +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include #include "istring.h" #include "definitions.cpp" #include "tools.h" @@ -358,6 +362,30 @@ long double number_double_val(Number x){ } return(ret); } +// approximate numerical expression (as mpfr float) +int number_mpfr_val(mpfr_t out, Number x){ + int i; + // auxiliary variables (do not initialize A) + mpfr_t A,b,c; + mpfr_inits(b,c, (mpfr_ptr)NULL); + + mpfr_init(out); + mpfr_set_zero(out,1); + + for(i=0;i\n\n"); + printf("\nusage:\n numkondo [-F] [-N niter] [-I initial_condition] [-P precision] [-E exponent_range] \n\n"); return(0); } @@ -171,8 +188,22 @@ int numflow(Str_Array str_args, Numkondo_Options opts){ Labels labels; // initial condition RCC init_cd; + RCC_mpfr init_cd_mpfr; // flow equation Grouped_Polynomial flow_equation; + // whether or not to use mpfr floats + int mpfr_flag=0; + + // set mpfr defaults + if(opts.mpfr_prec!=0){ + mpfr_set_default_prec(opts.mpfr_prec); + mpfr_flag=1; + } + if(opts.mpfr_emax!=0){ + mpfr_set_emax(opts.mpfr_emax); + mpfr_flag=1; + } + // parse id table arg_index=find_str_arg("labels", str_args); @@ -207,20 +238,31 @@ int numflow(Str_Array str_args, Numkondo_Options opts){ } } // initialize the rccs - prepare_init(flow_equation.indices,flow_equation.length,&init_cd); + if(mpfr_flag==0){ + prepare_init(flow_equation.indices,flow_equation.length,&init_cd); + } + else{ + prepare_init_mpfr(flow_equation.indices,flow_equation.length,&init_cd_mpfr); + } + // read rccs from string if(opts.eval_rccstring.length!=-1){ - parse_init_cd(opts.eval_rccstring, &init_cd); + parse_init_cd(opts.eval_rccstring, &init_cd, &init_cd_mpfr, mpfr_flag); free_Char_Array(opts.eval_rccstring); } - numerical_flow(flow_equation, init_cd, labels, opts.niter, opts.tol, opts.display_mode); + if(mpfr_flag==0){ + numerical_flow(flow_equation, init_cd, labels, opts.niter, opts.display_mode); + free_RCC(init_cd); + } + else{ + numerical_flow_mpfr(flow_equation, init_cd_mpfr, labels, opts.niter, opts.display_mode); + free_RCC_mpfr(init_cd_mpfr); + } - 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 index 19b91c6..ae3a9af 100644 --- a/src/parse_file.c +++ b/src/parse_file.c @@ -18,12 +18,14 @@ limitations under the License. #include #include +#include #include "array.h" #include "fields.h" #include "rational.h" #include "number.h" #include "polynomial.h" #include "rcc.h" +#include "rcc_mpfr.h" #include "definitions.cpp" #include "istring.h" #include "tools.h" @@ -716,8 +718,8 @@ int parse_labels(Char_Array str_labels, Labels* labels){ } -// read initial condition for numerical computation -int parse_init_cd(Char_Array init_cd, RCC* init){ +// read initial condition for numerical computation (using either RCC or RCC_mpfr, as specified by mpfr_flag) +int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_flag){ char* buffer=calloc(init_cd.length+1,sizeof(char)); char* buffer_ptr=buffer; int index=0; @@ -738,7 +740,12 @@ int parse_init_cd(Char_Array init_cd, RCC* init){ // new term case ',': // write init - sscanf(buffer,"%Lf",(*init).values+intlist_find_err((*init).indices,(*init).length,index)); + if(mpfr_flag==0){ + sscanf(buffer,"%Lf",(*init).values+intlist_find_err((*init).indices,(*init).length,index)); + } + else{ + mpfr_strtofr((*init_mpfr).values[intlist_find_err((*init_mpfr).indices,(*init_mpfr).length,index)], buffer, &buffer_ptr, 10, MPFR_RNDN); + } // reset buffer buffer_ptr=buffer; *buffer_ptr='\0'; @@ -782,7 +789,12 @@ int parse_init_cd(Char_Array init_cd, RCC* init){ } // write init - sscanf(buffer,"%Lf",(*init).values+unlist_find((*init).indices,(*init).length,index)); + if(mpfr_flag==0){ + sscanf(buffer,"%Lf",(*init).values+intlist_find_err((*init).indices,(*init).length,index)); + } + else{ + mpfr_strtofr((*init_mpfr).values[intlist_find_err((*init_mpfr).indices,(*init_mpfr).length,index)], buffer, &buffer_ptr, 10, MPFR_RNDN); + } free(buffer); return(0); @@ -806,3 +818,20 @@ int prepare_init(int* indices, int length, RCC* init){ } return(0); } +// set indices and length of init for RCC_mpfr +int prepare_init_mpfr(int* indices, int length, RCC_mpfr* init){ + int i; + init_RCC_mpfr(init, length); + for(i=0;i-DOFFSET){ + mpfr_set_ui((*init).values[i],1,MPFR_RNDN); + } + else{ + mpfr_set_zero((*init).values[i],1); + } + + } + return(0); +} diff --git a/src/parse_file.h b/src/parse_file.h index b51103a..0430763 100644 --- a/src/parse_file.h +++ b/src/parse_file.h @@ -47,10 +47,12 @@ int parse_input_id_table(Char_Array str_idtable, Id_Table* idtable, Fields_Table // 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); +// read initial condition for numerical computation (using either RCC or RCC_mpfr, as specified by mpfr_flag) +int parse_init_cd(Char_Array init_cd, RCC* init, RCC_mpfr* init_mpfr, int mpfr_flag); // set indices and length of init int prepare_init(int* indices, int length, RCC* init); +// set indices and length of init for RCC_mpfr +int prepare_init_mpfr(int* indices, int length, RCC_mpfr* init); #endif diff --git a/src/rational_float.c b/src/rational_float.c index eebc4f4..d6a6a4c 100644 --- a/src/rational_float.c +++ b/src/rational_float.c @@ -19,21 +19,14 @@ limitations under the License. #include "rational_float.h" #include #include +#include +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#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; @@ -141,6 +134,16 @@ long double lcm(long double x,long double y){ double Q_double_value(Q q){ return(1.0*q.numerator/q.denominator); } +// approximate value as mpfr float +int Q_mpfr_value(mpfr_t out, Q q){ + mpfr_t x; + mpfr_init(out); + mpfr_init(x); + mpfr_set_ld(x, q.denominator, MPFR_RNDN); + mpfr_ld_div(out, q.numerator, x, MPFR_RNDN); + mpfr_clear(x); + return(0); +} // print to string diff --git a/src/rational_float.h b/src/rational_float.h index 931b5ec..02c4a22 100644 --- a/src/rational_float.h +++ b/src/rational_float.h @@ -54,6 +54,8 @@ long double lcm(long double x,long double y); // approximate value as double double Q_double_value(Q q); +// approximate value as mpfr float +int Q_mpfr_value(mpfr_t out, Q q); // print to string int Q_sprint(Q num, Char_Array* out); diff --git a/src/rational_int.c b/src/rational_int.c index ec601a8..7aab561 100644 --- a/src/rational_int.c +++ b/src/rational_int.c @@ -19,6 +19,10 @@ limitations under the License. #include "rational_int.h" #include #include +#include +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include #include "istring.h" #include "array.h" @@ -135,6 +139,16 @@ long int lcm(long int x,long int y){ double Q_double_value(Q q){ return(1.0*q.numerator/q.denominator); } +// approximate value as mpfr float +int Q_mpfr_value(mpfr_t out, Q q){ + mpfr_t x; + mpfr_init(out); + mpfr_init(x); + mpfr_set_si(x, q.denominator, MPFR_RNDN); + mpfr_si_div(out, q.numerator, x, MPFR_RNDN); + mpfr_clear(x); + return(0); +} // print to string diff --git a/src/rational_int.h b/src/rational_int.h index a9f164d..a4dd000 100644 --- a/src/rational_int.h +++ b/src/rational_int.h @@ -49,6 +49,8 @@ long int lcm(long int x,long int y); // approximate value as double double Q_double_value(Q q); +// approximate value as mpfr float +int Q_mpfr_value(mpfr_t out, Q q); // print to string int Q_sprint(Q num, Char_Array* out); diff --git a/src/rcc_mpfr.c b/src/rcc_mpfr.c new file mode 100644 index 0000000..8a362e3 --- /dev/null +++ b/src/rcc_mpfr.c @@ -0,0 +1,124 @@ +/* +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_mpfr.h" +#include +#include +#include +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +// define MPFR_USE_FILE to enable the use of mpfr_printf +#define MPFR_USE_FILE +#include +#include +#include "array.h" + +// init +int init_RCC_mpfr(RCC_mpfr* rcc_mpfr, int size){ + int i; + (*rcc_mpfr).values=calloc(size,sizeof(mpfr_t)); + (*rcc_mpfr).indices=calloc(size,sizeof(int)); + (*rcc_mpfr).length=size; + for(i=0;i + // rational number typedef struct Q{ @@ -114,6 +116,12 @@ typedef struct RCC{ int* indices; int length; } RCC; +// rcc using mpfr floats +typedef struct RCC_mpfr{ + mpfr_t* values; + int* indices; + int length; +} RCC_mpfr; // identities between fields typedef struct Identities{ @@ -174,25 +182,6 @@ typedef struct Id_Table{ 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{ @@ -203,8 +192,9 @@ typedef struct Meankondo_Options{ typedef struct Numkondo_Options{ int display_mode; int niter; - long double tol; Char_Array eval_rccstring; + mpfr_prec_t mpfr_prec; + mpfr_exp_t mpfr_emax; } Numkondo_Options; typedef struct Meantools_Options{ @@ -213,6 +203,8 @@ typedef struct Meantools_Options{ Int_Array deriv_vars; Char_Array eval_rccstring; int chain; + mpfr_prec_t mpfr_prec; + mpfr_exp_t mpfr_emax; } Meantools_Options; typedef struct Kondopp_Options{ -- cgit v1.2.3-54-g00ecf