propagate precision in mpfr functions

evaleev · Aug 10, 2018 · f36de87 · f36de87
1 parent 27bdaa5
commit f36de87
Showing 1 changed file with 18 additions and 14 deletions.
diff --git a/include/libint2/numeric.h b/include/libint2/numeric.h
@@ -21,57 +21,61 @@
 #if LIBINT_HAS_MPFR
  /// implement exp for mpf_class using MPFR ... I do not claim to know what issues the rounding presents here
  inline mpf_class exp(mpf_class x) {
-   mpfr_t x_r; mpfr_init(x_r);
+   const auto prec = x.get_prec();
+   mpfr_t x_r; mpfr_init2(x_r, prec);
    mpfr_set_f(x_r, x.get_mpf_t(), MPFR_RNDN);
 
-   mpfr_t expx_r; mpfr_init(expx_r);
+   mpfr_t expx_r; mpfr_init2(expx_r, prec);
    mpfr_exp(expx_r, x_r, MPFR_RNDN);
 
    mpf_t expx;
-   mpf_init(expx);
+   mpf_init2(expx, prec);
    mpfr_get_f(expx, expx_r, MPFR_RNDN);
-   mpf_class result(expx);
+   mpf_class result(expx, prec);
    return result;
  }
  /// implement pow for mpf_class using MPFR ... I do not claim to know what issues the rounding presents here
  inline mpf_class pow(mpf_class x, int a) {
+   const auto prec = x.get_prec();
    mpf_t x_to_a;
-   mpf_init(x_to_a);
+   mpf_init2(x_to_a, prec);
    if (a >= 0)
      mpf_pow_ui(x_to_a, x.get_mpf_t(), (unsigned int) a);
    else
      mpf_pow_ui(x_to_a, x.get_mpf_t(), (unsigned int) (-a));
-   mpf_class result(x_to_a);
+   mpf_class result(x_to_a, prec);
    if (a < 0)
      result = 1.0 / result;
    return result;
  }
  /// implement erf for mpf_class using MPFR ... I do not claim to know what issues the rounding presents here
  inline mpf_class erf(mpf_class x) {
-   mpfr_t x_r; mpfr_init(x_r);
+   const auto prec = x.get_prec();
+   mpfr_t x_r; mpfr_init2(x_r, prec);
    mpfr_set_f(x_r, x.get_mpf_t(), MPFR_RNDN);
 
-   mpfr_t erfx_r; mpfr_init(erfx_r);
+   mpfr_t erfx_r; mpfr_init2(erfx_r, prec);
    mpfr_erf(erfx_r, x_r, MPFR_RNDN);
 
    mpf_t erfx;
-   mpf_init(erfx);
+   mpf_init2(erfx, prec);
    mpfr_get_f(erfx, erfx_r, MPFR_RNDN);
-   mpf_class result(erfx);
+   mpf_class result(erfx, prec);
    return result;
  }
  /// implement log for mpf_class using MPFR ... I do not claim to know what issues the rounding presents here
  inline mpf_class log(mpf_class x) {
-   mpfr_t x_r; mpfr_init(x_r);
+   const auto prec = x.get_prec();
+   mpfr_t x_r; mpfr_init2(x_r, prec);
    mpfr_set_f(x_r, x.get_mpf_t(), MPFR_RNDN);
 
-   mpfr_t logx_r; mpfr_init(logx_r);
+   mpfr_t logx_r; mpfr_init2(logx_r, prec);
    mpfr_log(logx_r, x_r, MPFR_RNDN);
 
    mpf_t logx;
-   mpf_init(logx);
+   mpf_init2(logx, prec);
    mpfr_get_f(logx, logx_r, MPFR_RNDN);
-   mpf_class result(logx);
+   mpf_class result(logx, prec);
    return result;
  }
 #endif