Support MPFR floats in numkondov1.4

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

1 files changed, 46 insertions, 1 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);
+}