LIBINT  2.6.0
numeric.h
1 //
2 // Created by Eduard Valeyev on 8/7/18.
3 //
4 
5 #ifndef _libint2_include_numeric_h_
6 #define _libint2_include_numeric_h_
7 
8 #include <libint2/config.h>
9 #include <iomanip>
10 #include <limits>
11 #include <sstream>
12 #include <type_traits>
13 
14 #if LIBINT_HAS_MPFR
15 # include <cstddef>
16 # include <gmpxx.h>
17 # include <mpfr.h>
18 #endif
19 
20 #include <libint2/util/generated/libint2_params.h>
21 #include <libint2/util/type_traits.h>
22 
23 #if LIBINT_HAS_MPFR
24  inline mpf_class exp(mpf_class x) {
26  const auto prec = x.get_prec();
27  mpfr_t x_r; mpfr_init2(x_r, prec);
28  mpfr_set_f(x_r, x.get_mpf_t(), MPFR_RNDN);
29 
30  mpfr_t expx_r; mpfr_init2(expx_r, prec);
31  mpfr_exp(expx_r, x_r, MPFR_RNDN);
32 
33  mpf_t expx;
34  mpf_init2(expx, prec);
35  mpfr_get_f(expx, expx_r, MPFR_RNDN);
36  mpf_class result(expx, prec);
37  return result;
38  }
40  inline mpf_class pow(mpf_class x, int a) {
41  const auto prec = x.get_prec();
42  mpf_t x_to_a;
43  mpf_init2(x_to_a, prec);
44  if (a >= 0)
45  mpf_pow_ui(x_to_a, x.get_mpf_t(), (unsigned int) a);
46  else
47  mpf_pow_ui(x_to_a, x.get_mpf_t(), (unsigned int) (-a));
48  mpf_class result(x_to_a, prec);
49  if (a < 0)
50  result = 1.0 / result;
51  return result;
52  }
54  inline double pow(double x, int a) {
55  return std::pow(x, static_cast<double>(a));
56  }
58  inline mpf_class erf(mpf_class x) {
59  const auto prec = x.get_prec();
60  mpfr_t x_r; mpfr_init2(x_r, prec);
61  mpfr_set_f(x_r, x.get_mpf_t(), MPFR_RNDN);
62 
63  mpfr_t erfx_r; mpfr_init2(erfx_r, prec);
64  mpfr_erf(erfx_r, x_r, MPFR_RNDN);
65 
66  mpf_t erfx;
67  mpf_init2(erfx, prec);
68  mpfr_get_f(erfx, erfx_r, MPFR_RNDN);
69  mpf_class result(erfx, prec);
70  return result;
71  }
73  inline mpf_class log(mpf_class x) {
74  const auto prec = x.get_prec();
75  mpfr_t x_r; mpfr_init2(x_r, prec);
76  mpfr_set_f(x_r, x.get_mpf_t(), MPFR_RNDN);
77 
78  mpfr_t logx_r; mpfr_init2(logx_r, prec);
79  mpfr_log(logx_r, x_r, MPFR_RNDN);
80 
81  mpf_t logx;
82  mpf_init2(logx, prec);
83  mpfr_get_f(logx, logx_r, MPFR_RNDN);
84  mpf_class result(logx, prec);
85  return result;
86  }
87 #endif
88 
89 #ifdef LIBINT_HAS_MPFR
90 using LIBINT2_REF_REALTYPE = mpf_class;
91 #else
92 using LIBINT2_REF_REALTYPE = double;
93 #endif
94 
95 namespace libint2 {
96  using value_type = LIBINT2_REALTYPE;
97  using scalar_type = libint2::vector_traits<value_type>::scalar_type;
98 
99  template <typename Real>
100  inline Real get_epsilon(const Real& value);
101 
102 #ifdef LIBINT_HAS_MPFR
103  template <>
104  inline mpf_class get_epsilon(const mpf_class& value) {
105  const auto nbits = value.get_prec();
106  return pow(mpf_class(2, nbits), -nbits);
107  };
108 #endif
109 
110  template <typename Real>
111  inline Real get_epsilon(const Real& value) {
112  return std::numeric_limits<Real>::epsilon();
113  }
114 
115 template <typename Real>
116 inline int get_max_digits10(const Real& value);
117 
118 #ifdef LIBINT_HAS_MPFR
119 template <>
120  inline int get_max_digits10(const mpf_class& value) {
121  const auto nbits = value.get_prec();
122  return std::ceil(nbits * std::log10(2) + 1);
123  };
124 #endif
125 
126 template <typename Real>
127 inline int get_max_digits10(const Real& value) {
128  return std::numeric_limits<Real>::max_digits10;
129 }
130 
131 template <typename To, typename From>
132 typename std::enable_if<!std::is_same<typename std::decay<To>::type, typename std::decay<From>::type >::value,To>::type
133 sstream_convert(From && from) {
134  std::stringstream ss;
135  ss << std::scientific << std::setprecision(get_max_digits10(from)) << from;
136  To to(ss.str().c_str());
137  return to;
138 }
139 
140 template <typename To>
141 To sstream_convert(const To& from) {
142  return from;
143 }
144 
145 }; // namespace libint2
146 
147 #endif // _libint2_include_numeric_h_
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24