Support MPFR floats in numkondov1.4

Remove '-D' option (error tolerance) in numkondo

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{