Ian Jauslin
summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/coefficient.c47
-rw-r--r--src/coefficient.h2
-rw-r--r--src/definitions.cpp2
-rw-r--r--src/flow.c24
-rw-r--r--src/flow.h8
-rw-r--r--src/flow_mpfr.c128
-rw-r--r--src/flow_mpfr.h32
-rw-r--r--src/grouped_polynomial.c30
-rw-r--r--src/grouped_polynomial.h2
-rw-r--r--src/meantools.c2
-rw-r--r--src/meantools_eval.c64
-rw-r--r--src/number.c28
-rw-r--r--src/number.h2
-rw-r--r--src/numkondo.c80
-rw-r--r--src/parse_file.c37
-rw-r--r--src/parse_file.h6
-rw-r--r--src/rational_float.c25
-rw-r--r--src/rational_float.h2
-rw-r--r--src/rational_int.c14
-rw-r--r--src/rational_int.h2
-rw-r--r--src/rcc_mpfr.c124
-rw-r--r--src/rcc_mpfr.h43
-rw-r--r--src/types.h32
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
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<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);
}
diff --git a/src/flow.h b/src/flow.h
index baef992..9491323 100644
--- a/src/flow.h
+++ b/src/flow.h
@@ -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{