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 /// implement exp for mpf_class using MPFR ... I do not claim to know what issues the rounding presents here
exp(mpf_class x)25 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 }
39 /// implement pow for mpf_class using MPFR ... I do not claim to know what issues the rounding presents here
pow(mpf_class x,int a)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 }
53 /// this is needed to avoid ambiguity in pow(2.0, 2) ... the above pow competes with standard double pow(double, double)
pow(double x,int a)54 inline double pow(double x, int a) {
55 return std::pow(x, static_cast<double>(a));
56 }
57 /// implement erf for mpf_class using MPFR ... I do not claim to know what issues the rounding presents here
erf(mpf_class x)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 }
72 /// implement acos for mpf_class using MPFR ... I do not claim to know what issues the rounding presents here
acos(mpf_class x)73 inline mpf_class acos(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 acosx_r; mpfr_init2(acosx_r, prec);
79 mpfr_acos(acosx_r, x_r, MPFR_RNDN);
80
81 mpf_t acosx;
82 mpf_init2(acosx, prec);
83 mpfr_get_f(acosx, acosx_r, MPFR_RNDN);
84 mpf_class result(acosx, prec);
85 return result;
86 }
87 /// implement log for mpf_class using MPFR ... I do not claim to know what issues the rounding presents here
log(mpf_class x)88 inline mpf_class log(mpf_class x) {
89 const auto prec = x.get_prec();
90 mpfr_t x_r; mpfr_init2(x_r, prec);
91 mpfr_set_f(x_r, x.get_mpf_t(), MPFR_RNDN);
92
93 mpfr_t logx_r; mpfr_init2(logx_r, prec);
94 mpfr_log(logx_r, x_r, MPFR_RNDN);
95
96 mpf_t logx;
97 mpf_init2(logx, prec);
98 mpfr_get_f(logx, logx_r, MPFR_RNDN);
99 mpf_class result(logx, prec);
100 return result;
101 }
102 #endif
103
104 #ifdef LIBINT_HAS_MPFR
105 using LIBINT2_REF_REALTYPE = mpf_class;
106 #else
107 using LIBINT2_REF_REALTYPE = double;
108 #endif
109
110 namespace libint2 {
111 using value_type = LIBINT2_REALTYPE;
112 using scalar_type = libint2::vector_traits<value_type>::scalar_type;
113
114 template <typename Real>
115 inline Real get_epsilon(const Real& value);
116
117 #ifdef LIBINT_HAS_MPFR
118 template <>
get_epsilon(const mpf_class & value)119 inline mpf_class get_epsilon(const mpf_class& value) {
120 const auto nbits = value.get_prec();
121 return pow(mpf_class(2, nbits), -nbits);
122 };
123 #endif
124
125 template <typename Real>
get_epsilon(const Real & value)126 inline Real get_epsilon(const Real& value) {
127 return std::numeric_limits<Real>::epsilon();
128 }
129
130 template <typename Real>
131 inline int get_max_digits10(const Real& value);
132
133 #ifdef LIBINT_HAS_MPFR
134 template <>
get_max_digits10(const mpf_class & value)135 inline int get_max_digits10(const mpf_class& value) {
136 const auto nbits = value.get_prec();
137 return std::ceil(nbits * std::log10(2) + 1);
138 };
139 #endif
140
141 template <typename Real>
get_max_digits10(const Real & value)142 inline int get_max_digits10(const Real& value) {
143 return std::numeric_limits<Real>::max_digits10;
144 }
145
146 template <typename To, typename From>
147 typename std::enable_if<!std::is_same<typename std::decay<To>::type, typename std::decay<From>::type >::value,To>::type
sstream_convert(From && from)148 sstream_convert(From && from) {
149 std::stringstream ss;
150 ss << std::scientific << std::setprecision(get_max_digits10(from)) << from;
151 To to(ss.str().c_str());
152 return to;
153 }
154
155 template <typename To>
sstream_convert(const To & from)156 To sstream_convert(const To& from) {
157 return from;
158 }
159
160 }; // namespace libint2
161
162 #endif // _libint2_include_numeric_h_
163