diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/coefficient.c | 47 | ||||
| -rw-r--r-- | src/coefficient.h | 2 | ||||
| -rw-r--r-- | src/definitions.cpp | 2 | ||||
| -rw-r--r-- | src/flow.c | 24 | ||||
| -rw-r--r-- | src/flow.h | 8 | ||||
| -rw-r--r-- | src/flow_mpfr.c | 128 | ||||
| -rw-r--r-- | src/flow_mpfr.h | 32 | ||||
| -rw-r--r-- | src/grouped_polynomial.c | 30 | ||||
| -rw-r--r-- | src/grouped_polynomial.h | 2 | ||||
| -rw-r--r-- | src/meantools.c | 2 | ||||
| -rw-r--r-- | src/meantools_eval.c | 64 | ||||
| -rw-r--r-- | src/number.c | 28 | ||||
| -rw-r--r-- | src/number.h | 2 | ||||
| -rw-r--r-- | src/numkondo.c | 80 | ||||
| -rw-r--r-- | src/parse_file.c | 37 | ||||
| -rw-r--r-- | src/parse_file.h | 6 | ||||
| -rw-r--r-- | src/rational_float.c | 25 | ||||
| -rw-r--r-- | src/rational_float.h | 2 | ||||
| -rw-r--r-- | src/rational_int.c | 14 | ||||
| -rw-r--r-- | src/rational_int.h | 2 | ||||
| -rw-r--r-- | src/rcc_mpfr.c | 124 | ||||
| -rw-r--r-- | src/rcc_mpfr.h | 43 | ||||
| -rw-r--r-- | src/types.h | 32 |
23 files changed, 645 insertions, 91 deletions
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 <stdio.h> #include <stdlib.h> +#include <stdarg.h> +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include <mpfr.h> #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;i<coef.length;i++){ @@ -737,3 +741,44 @@ int evalcoef(RCC rccs, Coefficient coef, long double* out){ } return(0); } + +// evaluate a coefficient on a vector (using mpfr floats) +int evalcoef_mpfr(RCC_mpfr rccs, Coefficient coef, mpfr_t out){ + int i,j; + mpfr_t num_factor; + // tmp number (do not initialize Z) + mpfr_t x, y, Z; + + // init numbers + mpfr_inits(num_factor, x, y, (mpfr_ptr) NULL); + + mpfr_init(out); + mpfr_set_zero(out, 1); + + // for each monomial + for(i=0;i<coef.length;i++){ + // product of factors + mpfr_set_flt(num_factor, 1., MPFR_RNDN); + for(j=0;j<coef.factors[i].length;j++){ + mpfr_mul(x,num_factor,rccs.values[intlist_find_err(rccs.indices,rccs.length,coef.factors[i].values[j])], MPFR_RNDN); + mpfr_set(num_factor,x, MPFR_RNDN); + } + // denominator + if(coef.denoms[i].power>0){ + 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 @@ -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<niter;i++){ // compute a single step - step_flow(&rccs, flow_equation, tol); + step_flow(&rccs, flow_equation); // convert ls to alphas if(display_mode==DISPLAY_NUMERICAL){ // print the result @@ -83,14 +81,9 @@ int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, in } // single step in the flow no exponentials -// inputs: flow_equation, tol -// input/outputs: rccs -int step_flow(RCC* rccs, Grouped_Polynomial flow_equation, long double tol){ +int step_flow(RCC* rccs, Grouped_Polynomial flow_equation){ int i; long double* new_rccs=calloc((*rccs).length,sizeof(long double)); - Int_Array computed; - - init_Int_Array(&computed, (*rccs).length); // initialize vectors to 0 for(i=0;i<(*rccs).length;i++){ @@ -101,10 +94,6 @@ int step_flow(RCC* rccs, Grouped_Polynomial flow_equation, long double tol){ for(i=0;i<flow_equation.length;i++){ if(flow_equation.indices[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])<tol){ - new_rccs[i]=0.; - } (*rccs).values[i]=new_rccs[i]; } } @@ -113,10 +102,6 @@ int step_flow(RCC* rccs, Grouped_Polynomial flow_equation, long double tol){ for(i=0;i<flow_equation.length;i++){ if(flow_equation.indices[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])<tol){ - new_rccs[i]=0.; - } } } @@ -126,7 +111,6 @@ int step_flow(RCC* rccs, Grouped_Polynomial flow_equation, long double tol){ } // free memory - free_Int_Array(computed); free(new_rccs); return(0); } @@ -21,14 +21,12 @@ Compute flow numerically #ifndef NUMERICAL_FLOW_H #define NUMERICAL_FLOW_H - -#include "grouped_polynomial.h" -#include "rcc.h" +#include "types.h" // compute flow -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); // single step -int step_flow(RCC* rccs, Grouped_Polynomial flow_equation, long double tol); +int step_flow(RCC* rccs, Grouped_Polynomial flow_equation); // print the label of an rcc (takes constants and derivatives into account) int print_label(int index, Labels labels); diff --git a/src/flow_mpfr.c b/src/flow_mpfr.c new file mode 100644 index 0000000..e03c130 --- /dev/null +++ b/src/flow_mpfr.c @@ -0,0 +1,128 @@ +/* +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 "flow_mpfr.h" + +#include <stdio.h> +#include <stdlib.h> +#include <stdarg.h> +// 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 <mpfr.h> +#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<rccs.length;j++){ + print_label(rccs.indices[j], labels); + } + printf("\n\n"); + + // print initial values + printf("%5d ",0); + for(j=0;j<rccs.length;j++){ + mpfr_printf("% 14.7Re ",rccs.values[j]); + } + printf("\n"); + } + + for(i=0;i<niter;i++){ + // compute a single step + step_flow_mpfr(&rccs, flow_equation); + // convert ls to alphas + if(display_mode==DISPLAY_NUMERICAL){ + // print the result + printf("%5d ",i+1); + for(j=0;j<rccs.length;j++){ + mpfr_printf("% 14.7Re ",rccs.values[j]); + } + printf("\n"); + } + } + + if(display_mode==DISPLAY_NUMERICAL){ + // print labels + printf("\n"); + printf("%5s ","n"); + for(j=0;j<rccs.length;j++){ + print_label(rccs.indices[j], labels); + } + printf("\n\n"); + } + + if(display_mode==DISPLAY_FINAL){ + RCC_mpfr_print(rccs); + } + + return(0); +} + +// single step in the flow +int step_flow_mpfr(RCC_mpfr* rccs, Grouped_Polynomial flow_equation){ + int i; + mpfr_t* res; + + // security: this function assumes that the length of the rcc and the flow_equation are the same + if((*rccs).length!=flow_equation.length){ + fprintf(stderr,"error: mismatch in the size of the flow equation and the rccs"); + exit(-1); + } + + res=calloc((*rccs).length,sizeof(mpfr_t)); + + // compute the constants first + for(i=0;i<flow_equation.length;i++){ + if(flow_equation.indices[i]<0){ + evalcoef_mpfr(*rccs, flow_equation.coefs[i], res[i]); + mpfr_set((*rccs).values[i], res[i], MPFR_RNDN); + } + } + + // for each equation + for(i=0;i<flow_equation.length;i++){ + if(flow_equation.indices[i]>=0){ + evalcoef_mpfr(*rccs, flow_equation.coefs[i], res[i]); + } + } + + // set new rccs + for(i=0;i<flow_equation.length;i++){ + mpfr_set((*rccs).values[i], res[i], MPFR_RNDN); + mpfr_clear(res[i]); + } + + // free memory + free(res); + return(0); +} diff --git a/src/flow_mpfr.h b/src/flow_mpfr.h new file mode 100644 index 0000000..e06cbd8 --- /dev/null +++ b/src/flow_mpfr.h @@ -0,0 +1,32 @@ +/* +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 flow numerically +*/ + +#ifndef NUMERICAL_FLOW_MPFR_H +#define NUMERICAL_FLOW_MPFR_H + + +#include "types.h" + +// compute flow +int numerical_flow_mpfr(Grouped_Polynomial flow_equation, RCC_mpfr init, Labels labels, int niter, int display_mode); +// single step +int step_flow_mpfr(RCC_mpfr* rccs, Grouped_Polynomial flow_equation); + +#endif diff --git a/src/grouped_polynomial.c b/src/grouped_polynomial.c index b7bea42..0d4d75d 100644 --- a/src/grouped_polynomial.c +++ b/src/grouped_polynomial.c @@ -732,7 +732,7 @@ int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output) } -// evaluate an equation on a vector +// eValuate an equation on a vector int evaleq(RCC* rccs, Grouped_Polynomial poly){ int i; long double* res=calloc((*rccs).length,sizeof(long double)); @@ -762,4 +762,32 @@ int evaleq(RCC* rccs, Grouped_Polynomial poly){ return(0); } +// evaluate an equation on a vector (using mpfr floats) +int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly){ + int i; + mpfr_t* res; + + if((*rccs).length!=poly.length){ + fprintf(stderr, "error: trying to evaluate an flow equation with %d components on an rcc with %d\n",poly.length,(*rccs).length); + exit(-1); + } + + res=calloc((*rccs).length,sizeof(mpfr_t)); + + // for each equation + for(i=0;i<poly.length;i++){ + evalcoef_mpfr(*rccs, poly.coefs[i], res[i]); + } + + // copy res to rccs + for(i=0;i<(*rccs).length;i++){ + mpfr_set((*rccs).values[i], res[i], MPFR_RNDN); + mpfr_clear(res[i]); + } + + // free memory + free(res); + return(0); +} + diff --git a/src/grouped_polynomial.h b/src/grouped_polynomial.h index 3c38f1b..e369721 100644 --- a/src/grouped_polynomial.h +++ b/src/grouped_polynomial.h @@ -70,5 +70,7 @@ int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output) // evaluate an equation on an RCC int evaleq(RCC* rccs, Grouped_Polynomial poly); +// evaluate an equation on a vector (using mpfr floats) +int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly); #endif diff --git a/src/meantools.c b/src/meantools.c index f45c376..e3e7247 100644 --- a/src/meantools.c +++ b/src/meantools.c @@ -108,7 +108,7 @@ int read_args_meantools(int argc,const char* argv[], Str_Array* str_args, Meanto // print usage message int print_usage_meantools(){ - printf("\nusage:\n meantools exp <filename>\n meantools derive [-d derivatives] -V <variables> <filename>\n meantools eval -R <rccs> <filename>\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 <stdio.h> #include <stdlib.h> +#include <mpfr.h> #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<argc;i++){ @@ -51,6 +60,14 @@ int tool_eval_read_args(int argc, const char* argv[], Str_Array* str_args, Meant case 'R': flag=CP_FLAG_RCCS; break; + // mpfr precision + case 'P': + flag=CP_FLAG_MPFR_PREC; + break; + // mpfr emax + case 'E': + flag=CP_FLAG_MPFR_EXP; + break; } } } @@ -59,6 +76,18 @@ int tool_eval_read_args(int argc, const char* argv[], Str_Array* str_args, Meant str_to_char_array((char*)argv[i], &((*opts).eval_rccstring)); flag=0; } + // mpfr precision + else if(flag==CP_FLAG_MPFR_PREC){ + sscanf(argv[i],"%ld",&tmp_lint); + (*opts).mpfr_prec=(mpfr_prec_t)tmp_lint; + flag=0; + } + // mpfr emax + else if(flag==CP_FLAG_MPFR_EXP){ + sscanf(argv[i],"%ld",&tmp_lint); + (*opts).mpfr_emax=(mpfr_exp_t)tmp_lint; + flag=0; + } // read file name from command-line else{ file=argv[i]; @@ -78,9 +107,21 @@ int tool_eval(Str_Array str_args, Meantools_Options opts){ int arg_index; // rccs RCC rccs; + RCC_mpfr rccs_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 flow equation // if there is a unique argument, assume it is the flow equation @@ -108,22 +149,33 @@ int tool_eval(Str_Array str_args, Meantools_Options opts){ } // initialize the rccs - prepare_init(flow_equation.indices,flow_equation.length,&rccs); + if(mpfr_flag==0){ + prepare_init(flow_equation.indices,flow_equation.length,&rccs); + } + else{ + prepare_init_mpfr(flow_equation.indices,flow_equation.length,&rccs_mpfr); + } // read rccs from string if(opts.eval_rccstring.length!=-1){ - parse_init_cd(opts.eval_rccstring, &rccs); + parse_init_cd(opts.eval_rccstring, &rccs, &rccs_mpfr, mpfr_flag); free_Char_Array(opts.eval_rccstring); } // evaluate - evaleq(&rccs, flow_equation); + if(mpfr_flag==0){ + evaleq(&rccs, flow_equation); + RCC_print(rccs); + free_RCC(rccs); + } + else{ + evaleq_mpfr(&rccs_mpfr, flow_equation); + RCC_mpfr_print(rccs_mpfr); + free_RCC_mpfr(rccs_mpfr); + } - // print - RCC_print(rccs); // free memory free_Grouped_Polynomial(flow_equation); - free_RCC(rccs); return(0); } diff --git a/src/number.c b/src/number.c index 5d4cd18..e63427f 100644 --- a/src/number.c +++ b/src/number.c @@ -18,6 +18,10 @@ limitations under the License. #include <stdlib.h> #include <stdio.h> #include <math.h> +#include <stdarg.h> +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include <mpfr.h> #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<x.length;i++){ + if(x.scalars[i].numerator!=0){ + mpfr_sqrt_ui(b, x.base[i], MPFR_RNDN); + Q_mpfr_value(A, x.scalars[i]); + mpfr_mul(c, A, b, MPFR_RNDN); + mpfr_add(b, out, c, MPFR_RNDN); + mpfr_set(out, b, MPFR_RNDN); + } + } + + mpfr_clears(A,b,c, (mpfr_ptr)NULL); + + return(0); +} // print to string diff --git a/src/number.h b/src/number.h index 4095f9c..a33d425 100644 --- a/src/number.h +++ b/src/number.h @@ -98,6 +98,8 @@ int number_is_zero(Number x); // approximate numerical expression long double number_double_val(Number x); +// approximate numerical expression (as mpfr float) +int number_mpfr_val(mpfr_t out, Number x); // print to string int number_sprint(Number number, Char_Array* out); diff --git a/src/numkondo.c b/src/numkondo.c index dce7de6..0568861 100644 --- a/src/numkondo.c +++ b/src/numkondo.c @@ -30,6 +30,7 @@ Compute the flow of a flow equation numerically // rccs #include "rcc.h" +#include "rcc_mpfr.h" // grouped representation of polynomials #include "grouped_polynomial.h" // command line parser @@ -38,6 +39,7 @@ Compute the flow of a flow equation numerically #include "parse_file.h" // numerical flow #include "flow.h" +#include "flow_mpfr.h" // arrays #include "array.h" @@ -68,10 +70,13 @@ int main (int argc, const char* argv[]){ // parse command-line arguments #define CP_FLAG_NITER 1 -#define CP_FLAG_TOL 2 -#define CP_FLAG_RCCS 3 +#define CP_FLAG_RCCS 2 +#define CP_FLAG_MPFR_PREC 3 +#define CP_FLAG_MPFR_EXP 4 int read_args_numkondo(int argc,const char* argv[], Str_Array* str_args, Numkondo_Options* opts){ int i; + // temporary long int + long int tmp_lint; // pointers char* ptr; // file to read the polynomial from in flow mode @@ -92,10 +97,11 @@ int read_args_numkondo(int argc,const char* argv[], Str_Array* str_args, Numkond (*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; + // no mpfr + (*opts).mpfr_prec=0; + (*opts).mpfr_emax=0; // loop over arguments for(i=1;i<argc;i++){ @@ -111,14 +117,18 @@ for(i=1;i<argc;i++){ case 'N': flag=CP_FLAG_NITER; break; - // tolerance - case 'D': - flag=CP_FLAG_TOL; - break; // initial condition case 'I': flag=CP_FLAG_RCCS; break; + // mpfr precision + case 'P': + flag=CP_FLAG_MPFR_PREC; + break; + // mpfr emax + case 'E': + flag=CP_FLAG_MPFR_EXP; + break; // print version case 'v': printf("numkondo " VERSION "\n"); @@ -134,16 +144,23 @@ for(i=1;i<argc;i++){ // reset flag flag=0; } - // tolerance - else if (flag==CP_FLAG_TOL){ - sscanf(argv[i],"%Lf",&((*opts).tol)); - flag=0; - } // init condition else if(flag==CP_FLAG_RCCS){ str_to_char_array((char*)argv[i], &((*opts).eval_rccstring)); flag=0; } + // mpfr precision + else if(flag==CP_FLAG_MPFR_PREC){ + sscanf(argv[i],"%ld",&tmp_lint); + (*opts).mpfr_prec=(mpfr_prec_t)tmp_lint; + flag=0; + } + // mpfr emax + else if(flag==CP_FLAG_MPFR_EXP){ + sscanf(argv[i],"%ld",&tmp_lint); + (*opts).mpfr_emax=(mpfr_exp_t)tmp_lint; + flag=0; + } // read file name from command-line else{ file=argv[i]; @@ -158,7 +175,7 @@ for(i=1;i<argc;i++){ // print usage message int print_usage_numkondo(){ - printf("\nusage:\n numkondo [-F] [-N niter] [-D tolerance] [-I initial_condition] <filename>\n\n"); + printf("\nusage:\n numkondo [-F] [-N niter] [-I initial_condition] [-P precision] [-E exponent_range] <filename>\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 <stdio.h> #include <stdlib.h> +#include <mpfr.h> #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<length;i++){ + (*init).indices[i]=indices[i]; + // set constants to 1 + if(indices[i]<0 && indices[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 <stdio.h> #include <stdlib.h> +#include <stdarg.h> +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include <mpfr.h> #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 <stdio.h> #include <stdlib.h> +#include <stdarg.h> +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include <mpfr.h> #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 <stdio.h> +#include <stdlib.h> +#include <stdarg.h> +// 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 <mpfr.h> +#include <math.h> +#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<size;i++){ + mpfr_init((*rcc_mpfr).values[i]); + } + + return(0); +} + +int free_RCC_mpfr(RCC_mpfr rcc_mpfr){ + int i; + for(i=0;i<rcc_mpfr.length;i++){ + mpfr_clear(rcc_mpfr.values[i]); + } + free(rcc_mpfr.values); + free(rcc_mpfr.indices); + return(0); +} + +// set a given element of an rcc_mpfr +int RCC_mpfr_set_elem(mpfr_t value, int index, RCC_mpfr* rcc_mpfr, int pos){ + mpfr_set((*rcc_mpfr).values[pos], value, MPFR_RNDN); + (*rcc_mpfr).indices[pos]=index; + return(0); +} + +int RCC_mpfr_cpy(RCC_mpfr input,RCC_mpfr* output){ + int i; + + init_RCC_mpfr(output,input.length); + for(i=0;i<input.length;i++){ + RCC_mpfr_set_elem(input.values[i], input.indices[i], output, i); + } + return(0); +} + +// concatenate rcc_mpfr +int RCC_mpfr_concat(RCC_mpfr rcc_mpfr1, RCC_mpfr rcc_mpfr2, RCC_mpfr* output){ + int i; + + init_RCC_mpfr(output,rcc_mpfr1.length+rcc_mpfr2.length); + + for(i=0;i<rcc_mpfr1.length;i++){ + RCC_mpfr_set_elem(rcc_mpfr1.values[i], rcc_mpfr1.indices[i], output, i); + } + + for(i=0;i<rcc_mpfr2.length;i++){ + RCC_mpfr_set_elem(rcc_mpfr2.values[i], rcc_mpfr2.indices[i], output, i+rcc_mpfr1.length); + } + + return(0); +} + +// append an rcc_mpfr at the end of another +int RCC_mpfr_append(RCC_mpfr input, RCC_mpfr* output){ + int i; + for(i=0;i<input.length;i++){ + RCC_mpfr_set_elem(input.values[i], input.indices[i], output, i+(*output).length); + } + (*output).length+=input.length; + return(0); +} + +// print an rcc_mpfr vector with maximal precision +int RCC_mpfr_print(RCC_mpfr rcc_mpfr){ + int j; + // the printf format + Char_Array printf_format; + // number of digits in output + int size; + + // compute size + // WARNING: assumes mpfr_default_prec is an int + size=mpfr_get_default_prec()*log10(2)-1; + + init_Char_Array(&printf_format,12); + char_array_snprintf(&printf_format,"%%d:%%.%dRe",size); + + for(j=0;j<rcc_mpfr.length;j++){ + mpfr_printf(printf_format.str,rcc_mpfr.indices[j],rcc_mpfr.values[j]); + if(j<rcc_mpfr.length-1){ + printf(",\n"); + } + else{ + printf("\n"); + } + } + + free_Char_Array(printf_format); + + return(0); +} diff --git a/src/rcc_mpfr.h b/src/rcc_mpfr.h new file mode 100644 index 0000000..18ebed3 --- /dev/null +++ b/src/rcc_mpfr.h @@ -0,0 +1,43 @@ +/* +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. +*/ + +/* + RCC_mpfr struct + + This data type is similar to RCC but the values of the rcc's are specified as mpfr floats +*/ + +#ifndef RCC_MPFR_H +#define RCC_MPFR_H + +#include "types.h" + +// init +int init_RCC_mpfr(RCC_mpfr* rcc_mpfr, int size); +int free_RCC_mpfr(RCC_mpfr rcc_mpfr); +// set an element of an rcc_mpfr +int RCC_mpfr_set_elem(mpfr_t value, int index, RCC_mpfr* rcc_mpfr, int pos); +// copy +int RCC_mpfr_cpy(RCC_mpfr input,RCC_mpfr* output); +// concatenate 2 rcc_mpfr_mpfr +int RCC_mpfr_concat(RCC_mpfr rcc_mpfr_mpfr1, RCC_mpfr rcc_mpfr_mpfr2, RCC_mpfr* output); +// append an rcc_mpfr to another +int RCC_mpfr_append(RCC_mpfr input, RCC_mpfr* output); + +// print an rcc_mpfr vector with maximal precision +int RCC_mpfr_print(RCC_mpfr rcc_mpfr_mpfr); + +#endif diff --git a/src/types.h b/src/types.h index a5d47b9..d30840e 100644 --- a/src/types.h +++ b/src/types.h @@ -21,6 +21,8 @@ limitations under the License. #ifndef TYPES_H #define TYPES_H +#include <mpfr.h> + // 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{ |