1 // Copyright 2016-2021 Francesco Biscani (bluescarni@gmail.com)
2 //
3 // This file is part of the mp++ library.
4 //
5 // This Source Code Form is subject to the terms of the Mozilla
6 // Public License v. 2.0. If a copy of the MPL was not distributed
7 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
9 #ifndef MPPP_REAL_HPP
10 #define MPPP_REAL_HPP
11
12 #include <mp++/config.hpp>
13
14 #if defined(MPPP_WITH_MPFR)
15
16 #include <algorithm>
17 #include <array>
18 #include <cassert>
19 #include <cmath>
20 #include <complex>
21 #include <cstddef>
22 #include <cstdint>
23 #include <iostream>
24 #include <limits>
25 #include <stdexcept>
26 #include <string>
27 #include <type_traits>
28 #include <utility>
29 #include <vector>
30
31 #if defined(MPPP_HAVE_STRING_VIEW)
32 #include <string_view>
33 #endif
34
35 #if defined(MPPP_WITH_BOOST_S11N)
36
37 #include <boost/archive/binary_iarchive.hpp>
38 #include <boost/archive/binary_oarchive.hpp>
39 #include <boost/serialization/access.hpp>
40 #include <boost/serialization/split_member.hpp>
41 #include <boost/serialization/string.hpp>
42 #include <boost/serialization/tracking.hpp>
43
44 #endif
45
46 #include <mp++/concepts.hpp>
47 #include <mp++/detail/gmp.hpp>
48 #include <mp++/detail/mpfr.hpp>
49 #include <mp++/detail/type_traits.hpp>
50 #include <mp++/detail/utils.hpp>
51 #include <mp++/detail/visibility.hpp>
52 #include <mp++/fwd.hpp>
53 #include <mp++/integer.hpp>
54 #include <mp++/rational.hpp>
55 #include <mp++/type_name.hpp>
56
57 #if defined(MPPP_WITH_QUADMATH)
58 #include <mp++/real128.hpp>
59 #endif
60
61 namespace mppp
62 {
63
64 namespace detail
65 {
66
67 // Clamp the precision between the min and max allowed values. This is used in the generic constructor/assignment
68 // operator.
clamp_mpfr_prec(::mpfr_prec_t p)69 constexpr ::mpfr_prec_t clamp_mpfr_prec(::mpfr_prec_t p)
70 {
71 return real_prec_check(p) ? p : (p < real_prec_min() ? real_prec_min() : real_prec_max());
72 }
73
74 // Helper function to print an mpfr to stream in a given base.
75 MPPP_DLL_PUBLIC void mpfr_to_stream(const ::mpfr_t, std::ostream &, int);
76
77 // Helpers to deduce the precision when constructing/assigning a real via another type.
78 // NOTE: it is important that these helpers return a valid (i.e., clamped) precision
79 // value, because they are used in contexts in which the output precision
80 // is not checked for validity.
81 template <typename T, enable_if_t<is_integral<T>::value, int> = 0>
real_deduce_precision(const T &)82 inline ::mpfr_prec_t real_deduce_precision(const T &)
83 {
84 static_assert(nl_digits<T>() < nl_max<::mpfr_prec_t>(), "Overflow error.");
85 // NOTE: for signed integers, include the sign bit as well.
86 return clamp_mpfr_prec(static_cast<::mpfr_prec_t>(nl_digits<T>()) + is_signed<T>::value);
87 }
88
89 // Utility function to determine the number of base-2 digits of the significand
90 // of a floating-point type which is not in base-2.
91 template <typename T>
dig2mpfr_prec()92 inline ::mpfr_prec_t dig2mpfr_prec()
93 {
94 static_assert(std::is_floating_point<T>::value, "Invalid type.");
95 // NOTE: just do a raw cast for the time being, it's not like we have ways of testing
96 // this in any case. In the future we could consider switching to a compile-time implementation
97 // of the integral log2, and do everything as compile-time integral computations.
98 return static_cast<::mpfr_prec_t>(std::ceil(nl_digits<T>() * std::log2(std::numeric_limits<T>::radix)));
99 }
100
101 template <typename T, enable_if_t<std::is_floating_point<T>::value, int> = 0>
real_deduce_precision(const T &)102 inline ::mpfr_prec_t real_deduce_precision(const T &)
103 {
104 static_assert(nl_digits<T>() <= nl_max<::mpfr_prec_t>(), "Overflow error.");
105 return clamp_mpfr_prec(std::numeric_limits<T>::radix == 2 ? static_cast<::mpfr_prec_t>(nl_digits<T>())
106 : dig2mpfr_prec<T>());
107 }
108
109 template <std::size_t SSize>
real_deduce_precision(const integer<SSize> & n)110 inline ::mpfr_prec_t real_deduce_precision(const integer<SSize> &n)
111 {
112 // Infer the precision from the bit size of n.
113 const auto ls = n.size();
114 // Check that ls * GMP_NUMB_BITS is representable by mpfr_prec_t.
115 // LCOV_EXCL_START
116 if (mppp_unlikely(ls > make_unsigned(nl_max<::mpfr_prec_t>()) / unsigned(GMP_NUMB_BITS))) {
117 throw std::overflow_error("The deduced precision for a real from an integer is too large");
118 }
119 // LCOV_EXCL_STOP
120 return clamp_mpfr_prec(static_cast<::mpfr_prec_t>(static_cast<::mpfr_prec_t>(ls) * GMP_NUMB_BITS));
121 }
122
123 template <std::size_t SSize>
real_deduce_precision(const rational<SSize> & q)124 inline ::mpfr_prec_t real_deduce_precision(const rational<SSize> &q)
125 {
126 // Infer the precision from the bit size of num/den.
127 const auto n_size = q.get_num().size();
128 const auto d_size = q.get_den().size();
129 // Overflow checks.
130 // LCOV_EXCL_START
131 if (mppp_unlikely(
132 // Overflow in total size.
133 (n_size > nl_max<decltype(q.get_num().size())>() - d_size)
134 // Check that tot_size * GMP_NUMB_BITS is representable by mpfr_prec_t.
135 || ((n_size + d_size) > make_unsigned(nl_max<::mpfr_prec_t>()) / unsigned(GMP_NUMB_BITS)))) {
136 throw std::overflow_error("The deduced precision for a real from a rational is too large");
137 }
138 // LCOV_EXCL_STOP
139 return clamp_mpfr_prec(static_cast<::mpfr_prec_t>(static_cast<::mpfr_prec_t>(n_size + d_size) * GMP_NUMB_BITS));
140 }
141
142 #if defined(MPPP_WITH_QUADMATH)
143
real_deduce_precision(const real128 &)144 inline ::mpfr_prec_t real_deduce_precision(const real128 &)
145 {
146 // The significand precision in bits is 113 for real128. Let's double-check it.
147 static_assert(real128_sig_digits() == 113u, "Invalid number of digits.");
148 return clamp_mpfr_prec(113);
149 }
150
151 #endif
152
153 // Fwd declare for friendship.
154 template <bool, typename F, typename Arg0, typename... Args>
155 real &mpfr_nary_op_impl(::mpfr_prec_t, const F &, real &, Arg0 &&, Args &&...);
156
157 template <bool, typename F, typename Arg0, typename... Args>
158 real mpfr_nary_op_return_impl(::mpfr_prec_t, const F &, Arg0 &&, Args &&...);
159
160 template <typename F>
161 real real_constant(const F &, ::mpfr_prec_t);
162
163 // Wrapper for calling mpfr_lgamma().
164 MPPP_DLL_PUBLIC void real_lgamma_wrapper(::mpfr_t, const ::mpfr_t, ::mpfr_rnd_t);
165
166 // Wrapper for calling mpfr_li2().
167 MPPP_DLL_PUBLIC void real_li2_wrapper(::mpfr_t, const ::mpfr_t, ::mpfr_rnd_t);
168
169 // Wrappers for calling integer and remainder-related functions
170 // with NaN checking.
171 MPPP_DLL_PUBLIC void real_ceil_wrapper(::mpfr_t, const ::mpfr_t);
172 MPPP_DLL_PUBLIC void real_floor_wrapper(::mpfr_t, const ::mpfr_t);
173 MPPP_DLL_PUBLIC void real_round_wrapper(::mpfr_t, const ::mpfr_t);
174 #if defined(MPPP_MPFR_HAVE_MPFR_ROUNDEVEN)
175 MPPP_DLL_PUBLIC void real_roundeven_wrapper(::mpfr_t, const ::mpfr_t);
176 #endif
177 MPPP_DLL_PUBLIC void real_trunc_wrapper(::mpfr_t, const ::mpfr_t);
178 MPPP_DLL_PUBLIC void real_frac_wrapper(::mpfr_t, const ::mpfr_t);
179
180 #if defined(MPPP_WITH_ARB)
181
182 // The Arb MPFR wrappers.
183 MPPP_DLL_PUBLIC void arb_sqrt1pm1(::mpfr_t, const ::mpfr_t);
184
185 MPPP_DLL_PUBLIC void arb_log_hypot(::mpfr_t, const ::mpfr_t, const ::mpfr_t);
186 MPPP_DLL_PUBLIC void arb_log_base_ui(::mpfr_t, const ::mpfr_t, unsigned long);
187
188 MPPP_DLL_PUBLIC void arb_sin_pi(::mpfr_t, const ::mpfr_t);
189 MPPP_DLL_PUBLIC void arb_cos_pi(::mpfr_t, const ::mpfr_t);
190 MPPP_DLL_PUBLIC void arb_tan_pi(::mpfr_t, const ::mpfr_t);
191 MPPP_DLL_PUBLIC void arb_cot_pi(::mpfr_t, const ::mpfr_t);
192 MPPP_DLL_PUBLIC void arb_sinc(::mpfr_t, const ::mpfr_t);
193 MPPP_DLL_PUBLIC void arb_sinc_pi(::mpfr_t, const ::mpfr_t);
194
195 MPPP_DLL_PUBLIC void arb_hypgeom_bessel_j(::mpfr_t, const ::mpfr_t, const ::mpfr_t);
196 MPPP_DLL_PUBLIC void arb_hypgeom_bessel_y(::mpfr_t, const ::mpfr_t, const ::mpfr_t);
197
198 MPPP_DLL_PUBLIC void arb_lambert_w0(::mpfr_t, const ::mpfr_t);
199 MPPP_DLL_PUBLIC void arb_lambert_wm1(::mpfr_t, const ::mpfr_t);
200
201 MPPP_DLL_PUBLIC void arb_polylog_si(::mpfr_t, long, const ::mpfr_t);
202 MPPP_DLL_PUBLIC void arb_polylog(::mpfr_t, const ::mpfr_t, const ::mpfr_t);
203
204 #endif
205
206 } // namespace detail
207
208 // Fwd declare swap.
209 void swap(real &, real &) noexcept;
210
211 template <typename T>
212 using is_real_interoperable = detail::disjunction<is_cpp_arithmetic<T>, detail::is_integer<T>, detail::is_rational<T>
213 #if defined(MPPP_WITH_QUADMATH)
214 ,
215 std::is_same<T, real128>
216 #endif
217 >;
218
219 #if defined(MPPP_HAVE_CONCEPTS)
220
221 template <typename T>
222 MPPP_CONCEPT_DECL real_interoperable = is_real_interoperable<T>::value;
223
224 #endif
225
226 template <typename T>
227 using is_cvr_real = std::is_same<detail::uncvref_t<T>, real>;
228
229 #if defined(MPPP_HAVE_CONCEPTS)
230
231 template <typename T>
232 MPPP_CONCEPT_DECL cvr_real = is_cvr_real<T>::value;
233
234 #endif
235
236 template <typename... Args>
237 using cvr_real_enabler = detail::enable_if_t<detail::conjunction<is_cvr_real<Args>...>::value, int>;
238
239 // Special initialisation tags for real.
240 enum class real_kind : std::underlying_type<::mpfr_kind_t>::type {
241 nan = MPFR_NAN_KIND,
242 inf = MPFR_INF_KIND,
243 zero = MPFR_ZERO_KIND
244 };
245
246 // For the future:
247 // - construction from/conversion to interoperables can probably be improved performance wise, especially
248 // if we exploit the mpfr_t internals.
249 // - probably we should have a build in the CI against the latest MPFR, built with sanitizers on.
250 // - probably we should have MPFR as well in the 32bit coverage build.
251 // - it seems like we might be doing multiple roundings when cting from real128.
252 // See if we can implement with only a single rounding, perhaps via an integer
253 // and mpfr_set_z_2exp()?
254 // - Do we need real_equal_to() to work also on invalid reals, the way
255 // real_lt/gt() do?
256 // - Not sure what the caching situation is currently. The MPFR 4 changelog mentions it:
257 // https://www.mpfr.org/mpfr-4.0.0/#changes
258 // But experiments with valgrind and the real_alloc benchmark seem to indicate
259 // that no caching is done. Need to revisit this.
260
261 // Multiprecision floating-point class.
262 class MPPP_DLL_PUBLIC real
263 {
264 #if defined(MPPP_WITH_BOOST_S11N)
265 friend class boost::serialization::access;
266
267 template <typename Archive>
save(Archive & ar,unsigned) const268 void save(Archive &ar, unsigned) const
269 {
270 ar << get_prec();
271 ar << to_string();
272 }
273
274 template <typename Archive>
load(Archive & ar,unsigned)275 void load(Archive &ar, unsigned)
276 {
277 // NOLINTNEXTLINE(cppcoreguidelines-init-variables)
278 ::mpfr_prec_t p;
279 ar >> p;
280 std::string tmp;
281 ar >> tmp;
282
283 *this = real{tmp, p};
284 }
285
286 // Overloads for binary archives.
287 void save(boost::archive::binary_oarchive &, unsigned) const;
288 void load(boost::archive::binary_iarchive &, unsigned);
289
290 BOOST_SERIALIZATION_SPLIT_MEMBER()
291 #endif
292
293 // Make friends, for accessing the non-checking prec setting funcs.
294 template <bool, typename F, typename Arg0, typename... Args>
295 // NOLINTNEXTLINE(readability-redundant-declaration)
296 friend real &detail::mpfr_nary_op_impl(::mpfr_prec_t, const F &, real &, Arg0 &&, Args &&...);
297 template <bool, typename F, typename Arg0, typename... Args>
298 // NOLINTNEXTLINE(readability-redundant-declaration)
299 friend real detail::mpfr_nary_op_return_impl(::mpfr_prec_t, const F &, Arg0 &&, Args &&...);
300 template <typename F>
301 // NOLINTNEXTLINE(readability-redundant-declaration)
302 friend real detail::real_constant(const F &, ::mpfr_prec_t);
303 // Utility function to check the precision upon init.
check_init_prec(::mpfr_prec_t p)304 static ::mpfr_prec_t check_init_prec(::mpfr_prec_t p)
305 {
306 if (mppp_unlikely(!detail::real_prec_check(p))) {
307 throw std::invalid_argument("Cannot init a real with a precision of " + detail::to_string(p)
308 + ": the maximum allowed precision is " + detail::to_string(real_prec_max())
309 + ", the minimum allowed precision is " + detail::to_string(real_prec_min()));
310 }
311 return p;
312 }
313
314 #if defined(MPPP_WITH_MPC)
315 // NOTE: the complex class needs access to some
316 // private bits of real.
317 friend class complex;
318
319 // Shallow copy constructor from mpfr_t, used
320 // only by the complex class.
321 struct shallow_copy_t {
322 };
real(shallow_copy_t,const::mpfr_t r)323 explicit real(shallow_copy_t, const ::mpfr_t r) : m_mpfr(r[0]) {}
324 #endif
325
326 public:
327 // Default constructor.
328 real();
329
330 private:
331 // A tag to call private ctors.
332 struct ptag {
333 };
334 // Private ctor that sets to NaN with a certain precision,
335 // without checking the input precision value.
336 explicit real(const ptag &, ::mpfr_prec_t, bool);
337
338 public:
339 // Copy constructor.
340 real(const real &);
341 // Move constructor.
real(real && other)342 real(real &&other) noexcept
343 : // Shallow copy other.
344 m_mpfr(other.m_mpfr)
345 {
346 // Mark the other as moved-from.
347 other.m_mpfr._mpfr_d = nullptr;
348 }
349
350 // Copy constructor with custom precision.
351 explicit real(const real &, ::mpfr_prec_t);
352 // Move constructor with custom precision.
353 explicit real(real &&, ::mpfr_prec_t);
354
355 // Constructor from a special value, sign and precision.
356 explicit real(real_kind, int, ::mpfr_prec_t);
357 // Constructor from a special value and precision.
358 explicit real(real_kind, ::mpfr_prec_t);
359
360 // Constructors from n*2**e.
361 template <std::size_t SSize>
362 explicit real(const integer<SSize> &, ::mpfr_exp_t, ::mpfr_prec_t);
363 explicit real(unsigned long, ::mpfr_exp_t, ::mpfr_prec_t);
364 explicit real(long, ::mpfr_exp_t, ::mpfr_prec_t);
365
366 private:
367 // Construction from FPs.
368 template <typename Func, typename T>
369 MPPP_DLL_LOCAL void dispatch_fp_construction(const Func &, const T &);
370 void dispatch_construction(const float &);
371 void dispatch_construction(const double &);
372 void dispatch_construction(const long double &);
373
374 // Construction from integral types.
375 // Special casing for bool, otherwise MSVC warns if we fold this into the
376 // constructor from unsigned.
377 void dispatch_construction(const bool &);
378 template <typename T,
379 detail::enable_if_t<detail::conjunction<detail::is_integral<T>, detail::is_unsigned<T>>::value, int> = 0>
dispatch_construction(const T & n)380 void dispatch_construction(const T &n)
381 {
382 if (n <= detail::nl_max<unsigned long>()) {
383 mpfr_set_ui(&m_mpfr, static_cast<unsigned long>(n), MPFR_RNDN);
384 } else {
385 // NOTE: here and elsewhere let's use a 2-limb integer, in the hope
386 // of avoiding dynamic memory allocation.
387 ::mpfr_set_z(&m_mpfr, integer<2>(n).get_mpz_view(), MPFR_RNDN);
388 }
389 }
390 template <typename T,
391 detail::enable_if_t<detail::conjunction<detail::is_integral<T>, detail::is_signed<T>>::value, int> = 0>
dispatch_construction(const T & n)392 void dispatch_construction(const T &n)
393 {
394 if (n <= detail::nl_max<long>() && n >= detail::nl_min<long>()) {
395 mpfr_set_si(&m_mpfr, static_cast<long>(n), MPFR_RNDN);
396 } else {
397 ::mpfr_set_z(&m_mpfr, integer<2>(n).get_mpz_view(), MPFR_RNDN);
398 }
399 }
400
401 // Construction from mppp::integer.
402 void dispatch_mpz_construction(const ::mpz_t);
403 template <std::size_t SSize>
dispatch_construction(const integer<SSize> & n)404 void dispatch_construction(const integer<SSize> &n)
405 {
406 dispatch_mpz_construction(n.get_mpz_view());
407 }
408
409 // Construction from mppp::rational.
410 void dispatch_mpq_construction(const ::mpq_t);
411 template <std::size_t SSize>
dispatch_construction(const rational<SSize> & q)412 void dispatch_construction(const rational<SSize> &q)
413 {
414 // NOTE: get_mpq_view() returns an mpq_struct, whose
415 // address we then need to use.
416 const auto v = detail::get_mpq_view(q);
417 dispatch_mpq_construction(&v);
418 }
419
420 #if defined(MPPP_WITH_QUADMATH)
421 void dispatch_construction(const real128 &);
422 // NOTE: split this off from the dispatch_construction() overload, so we can re-use it in the
423 // generic assignment.
424 void assign_real128(const real128 &);
425 #endif
426
427 public:
428 // Generic constructors.
429 #if defined(MPPP_HAVE_CONCEPTS)
430 template <real_interoperable T>
431 #else
432 template <typename T, detail::enable_if_t<is_real_interoperable<T>::value, int> = 0>
433 #endif
434 // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init, hicpp-member-init)
real(const T & x)435 real(const T &x) : real(ptag{}, detail::real_deduce_precision(x), true)
436 {
437 dispatch_construction(x);
438 }
439 #if defined(MPPP_HAVE_CONCEPTS)
440 template <real_interoperable T>
441 #else
442 template <typename T, detail::enable_if_t<is_real_interoperable<T>::value, int> = 0>
443 #endif
real(const T & x,::mpfr_prec_t p)444 explicit real(const T &x, ::mpfr_prec_t p) : real(ptag{}, check_init_prec(p), true)
445 {
446 dispatch_construction(x);
447 }
448
449 // Constructors from std::complex.
450 #if defined(MPPP_HAVE_CONCEPTS)
451 template <cpp_complex T>
452 #else
453 template <typename T, detail::enable_if_t<is_cpp_complex<T>::value, int> = 0>
454 #endif
455 // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init, hicpp-member-init)
real(const T & c)456 explicit real(const T &c)
457 : real(c.imag() == 0 ? c.real()
458 : throw std::domain_error(
459 "Cannot construct a real from a complex C++ value with a non-zero imaginary part of "
460 + detail::to_string(c.imag())))
461 {
462 }
463 #if defined(MPPP_HAVE_CONCEPTS)
464 template <cpp_complex T>
465 #else
466 template <typename T, detail::enable_if_t<is_cpp_complex<T>::value, int> = 0>
467 #endif
468 // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init, hicpp-member-init)
real(const T & c,::mpfr_prec_t p)469 explicit real(const T &c, ::mpfr_prec_t p)
470 : real(c.imag() == 0 ? c.real()
471 : throw std::domain_error(
472 "Cannot construct a real from a complex C++ value with a non-zero imaginary part of "
473 + detail::to_string(c.imag())),
474 p)
475 {
476 }
477
478 private:
479 MPPP_DLL_LOCAL void construct_from_c_string(const char *, int, ::mpfr_prec_t);
480 explicit real(const ptag &, const char *, int, ::mpfr_prec_t);
481 explicit real(const ptag &, const std::string &, int, ::mpfr_prec_t);
482 #if defined(MPPP_HAVE_STRING_VIEW)
483 explicit real(const ptag &, const std::string_view &, int, ::mpfr_prec_t);
484 #endif
485
486 public:
487 // Constructor from string, base and precision.
488 #if defined(MPPP_HAVE_CONCEPTS)
489 template <string_type T>
490 #else
491 template <typename T, detail::enable_if_t<is_string_type<T>::value, int> = 0>
492 #endif
493 // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init, hicpp-member-init)
real(const T & s,int base,::mpfr_prec_t p)494 explicit real(const T &s, int base, ::mpfr_prec_t p) : real(ptag{}, s, base, p)
495 {
496 }
497 // Constructor from string and precision.
498 #if defined(MPPP_HAVE_CONCEPTS)
499 template <string_type T>
500 #else
501 template <typename T, detail::enable_if_t<is_string_type<T>::value, int> = 0>
502 #endif
503 // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init, hicpp-member-init)
real(const T & s,::mpfr_prec_t p)504 explicit real(const T &s, ::mpfr_prec_t p) : real(s, 10, p)
505 {
506 }
507 // Constructor from range of characters, base and precision.
508 explicit real(const char *, const char *, int, ::mpfr_prec_t);
509 // Constructor from range of characters and precision.
510 explicit real(const char *, const char *, ::mpfr_prec_t);
511
512 // Copy constructor from mpfr_t.
513 explicit real(const ::mpfr_t);
514 #if !defined(_MSC_VER) || defined(__clang__)
515 // Move constructor from mpfr_t.
real(::mpfr_t && x)516 explicit real(::mpfr_t &&x) : m_mpfr(*x) {}
517 #endif
518
519 // Destructor.
520 ~real();
521
522 // Copy assignment operator.
523 real &operator=(const real &);
524
525 // Move assignment operator.
operator =(real && other)526 real &operator=(real &&other) noexcept
527 {
528 // NOTE: for generic code, std::swap() is not a particularly good way of implementing
529 // the move assignment:
530 //
531 // https://stackoverflow.com/questions/6687388/why-do-some-people-use-swap-for-move-assignments
532 //
533 // Here however it is fine, as we know there are no side effects we need to maintain.
534 //
535 // NOTE: we use a raw std::swap() here (instead of mpfr_swap()) because we don't know in principle
536 // if mpfr_swap() relies on the operands not to be in a moved-from state (although it's unlikely).
537 std::swap(m_mpfr, other.m_mpfr);
538 return *this;
539 }
540
541 private:
542 // Assignment from FPs.
543 template <bool SetPrec, typename Func, typename T>
dispatch_fp_assignment(const Func & func,const T & x)544 void dispatch_fp_assignment(const Func &func, const T &x)
545 {
546 if (SetPrec) {
547 set_prec_impl<false>(detail::real_deduce_precision(x));
548 }
549 func(&m_mpfr, x, MPFR_RNDN);
550 }
551 template <bool SetPrec>
dispatch_assignment(const float & x)552 void dispatch_assignment(const float &x)
553 {
554 dispatch_fp_assignment<SetPrec>(::mpfr_set_flt, x);
555 }
556 template <bool SetPrec>
dispatch_assignment(const double & x)557 void dispatch_assignment(const double &x)
558 {
559 dispatch_fp_assignment<SetPrec>(::mpfr_set_d, x);
560 }
561 template <bool SetPrec>
dispatch_assignment(const long double & x)562 void dispatch_assignment(const long double &x)
563 {
564 dispatch_fp_assignment<SetPrec>(::mpfr_set_ld, x);
565 }
566
567 // Assignment from integral types.
568 template <bool SetPrec, typename T>
dispatch_integral_ass_prec(const T & n)569 void dispatch_integral_ass_prec(const T &n)
570 {
571 if (SetPrec) {
572 set_prec_impl<false>(detail::real_deduce_precision(n));
573 }
574 }
575 // Special casing for bool.
576 template <bool SetPrec>
dispatch_assignment(const bool & b)577 void dispatch_assignment(const bool &b)
578 {
579 dispatch_integral_ass_prec<SetPrec>(b);
580 mpfr_set_ui(&m_mpfr, static_cast<unsigned long>(b), MPFR_RNDN);
581 }
582 template <bool SetPrec, typename T,
583 detail::enable_if_t<detail::conjunction<detail::is_integral<T>, detail::is_unsigned<T>>::value, int> = 0>
dispatch_assignment(const T & n)584 void dispatch_assignment(const T &n)
585 {
586 dispatch_integral_ass_prec<SetPrec>(n);
587 if (n <= detail::nl_max<unsigned long>()) {
588 mpfr_set_ui(&m_mpfr, static_cast<unsigned long>(n), MPFR_RNDN);
589 } else {
590 ::mpfr_set_z(&m_mpfr, integer<2>(n).get_mpz_view(), MPFR_RNDN);
591 }
592 }
593 template <bool SetPrec, typename T,
594 detail::enable_if_t<detail::conjunction<detail::is_integral<T>, detail::is_signed<T>>::value, int> = 0>
dispatch_assignment(const T & n)595 void dispatch_assignment(const T &n)
596 {
597 dispatch_integral_ass_prec<SetPrec>(n);
598 if (n <= detail::nl_max<long>() && n >= detail::nl_min<long>()) {
599 mpfr_set_si(&m_mpfr, static_cast<long>(n), MPFR_RNDN);
600 } else {
601 ::mpfr_set_z(&m_mpfr, integer<2>(n).get_mpz_view(), MPFR_RNDN);
602 }
603 }
604
605 // Assignment from integer.
606 template <bool SetPrec, std::size_t SSize>
dispatch_assignment(const integer<SSize> & n)607 void dispatch_assignment(const integer<SSize> &n)
608 {
609 if (SetPrec) {
610 set_prec_impl<false>(detail::real_deduce_precision(n));
611 }
612 ::mpfr_set_z(&m_mpfr, n.get_mpz_view(), MPFR_RNDN);
613 }
614
615 // Assignment from rational.
616 template <bool SetPrec, std::size_t SSize>
dispatch_assignment(const rational<SSize> & q)617 void dispatch_assignment(const rational<SSize> &q)
618 {
619 if (SetPrec) {
620 set_prec_impl<false>(detail::real_deduce_precision(q));
621 }
622 const auto v = detail::get_mpq_view(q);
623 ::mpfr_set_q(&m_mpfr, &v, MPFR_RNDN);
624 }
625
626 #if defined(MPPP_WITH_QUADMATH)
627 // Assignment from real128.
628 template <bool SetPrec>
dispatch_assignment(const real128 & x)629 void dispatch_assignment(const real128 &x)
630 {
631 if (SetPrec) {
632 set_prec_impl<false>(detail::real_deduce_precision(x));
633 }
634 assign_real128(x);
635 }
636 #endif
637
638 public:
639 // Generic assignment operator.
640 #if defined(MPPP_HAVE_CONCEPTS)
641 template <real_interoperable T>
642 #else
643 template <typename T, detail::enable_if_t<is_real_interoperable<T>::value, int> = 0>
644 #endif
operator =(const T & x)645 real &operator=(const T &x)
646 {
647 dispatch_assignment<true>(x);
648 return *this;
649 }
650 // Assignment from std::complex.
651 #if defined(MPPP_HAVE_CONCEPTS)
652 template <cpp_complex T>
653 #else
654 template <typename T, detail::enable_if_t<is_cpp_complex<T>::value, int> = 0>
655 #endif
operator =(const T & c)656 real &operator=(const T &c)
657 {
658 // NOLINTNEXTLINE(cppcoreguidelines-c-copy-assignment-signature, misc-unconventional-assign-operator)
659 return *this = static_cast<real>(c);
660 }
661 #if defined(MPPP_WITH_QUADMATH)
662 real &operator=(const complex128 &);
663 #endif
664 #if defined(MPPP_WITH_MPC)
665 real &operator=(const complex &);
666 #endif
667
668 // Copy assignment from mpfr_t.
669 real &operator=(const ::mpfr_t);
670
671 #if !defined(_MSC_VER) || defined(__clang__)
672 // Move assignment from mpfr_t.
673 real &operator=(::mpfr_t &&);
674 #endif
675
676 // Check validity.
is_valid() const677 MPPP_NODISCARD bool is_valid() const noexcept
678 {
679 return m_mpfr._mpfr_d != nullptr;
680 }
681
682 // Set to another real.
683 real &set(const real &);
684
685 // Generic setter.
686 #if defined(MPPP_HAVE_CONCEPTS)
687 template <real_interoperable T>
688 #else
689 template <typename T, detail::enable_if_t<is_real_interoperable<T>::value, int> = 0>
690 #endif
set(const T & x)691 real &set(const T &x)
692 {
693 dispatch_assignment<false>(x);
694 return *this;
695 }
696 // Set to std::complex.
697 #if defined(MPPP_HAVE_CONCEPTS)
698 template <cpp_complex T>
699 #else
700 template <typename T, detail::enable_if_t<is_cpp_complex<T>::value, int> = 0>
701 #endif
set(const T & c)702 real &set(const T &c)
703 {
704 if (mppp_unlikely(c.imag() != 0)) {
705 throw std::domain_error("Cannot set a real to a complex C++ value with a non-zero imaginary part of "
706 + detail::to_string(c.imag()));
707 }
708 return set(c.real());
709 }
710
711 private:
712 // Implementation of string setters.
713 MPPP_DLL_LOCAL void string_assignment_impl(const char *, int);
714 real &set_impl(const char *, int);
715 real &set_impl(const std::string &, int);
716 #if defined(MPPP_HAVE_STRING_VIEW)
717 real &set_impl(const std::string_view &, int);
718 #endif
719
720 public:
721 // Setter to string.
722 #if defined(MPPP_HAVE_CONCEPTS)
723 template <string_type T>
724 #else
725 template <typename T, detail::enable_if_t<is_string_type<T>::value, int> = 0>
726 #endif
set(const T & s,int base=10)727 real &set(const T &s, int base = 10)
728 {
729 return set_impl(s, base);
730 }
731 // Set to character range.
732 real &set(const char *begin, const char *end, int base = 10);
733
734 // Set to an mpfr_t.
735 real &set(const ::mpfr_t);
736
737 // Set to NaN.
738 real &set_nan();
739 // Set to infinity.
740 real &set_inf(int sign = 0);
741 // Set to zero.
742 real &set_zero(int sign = 0);
743
744 // Const reference to the internal mpfr_t.
get_mpfr_t() const745 MPPP_NODISCARD const mpfr_struct_t *get_mpfr_t() const
746 {
747 return &m_mpfr;
748 }
749 // Mutable reference to the internal mpfr_t.
_get_mpfr_t()750 mpfr_struct_t *_get_mpfr_t()
751 {
752 return &m_mpfr;
753 }
754
755 // Detect NaN.
nan_p() const756 MPPP_NODISCARD bool nan_p() const
757 {
758 return mpfr_nan_p(&m_mpfr) != 0;
759 }
760 // Detect infinity.
inf_p() const761 MPPP_NODISCARD bool inf_p() const
762 {
763 return mpfr_inf_p(&m_mpfr) != 0;
764 }
765 // Detect finite number.
number_p() const766 MPPP_NODISCARD bool number_p() const
767 {
768 return mpfr_number_p(&m_mpfr) != 0;
769 }
770 // Detect zero.
zero_p() const771 MPPP_NODISCARD bool zero_p() const
772 {
773 return mpfr_zero_p(&m_mpfr) != 0;
774 }
775 // Detect regular number.
regular_p() const776 MPPP_NODISCARD bool regular_p() const
777 {
778 return mpfr_regular_p(&m_mpfr) != 0;
779 }
780 // Detect integer.
781 MPPP_NODISCARD bool integer_p() const;
782 // Detect one.
783 MPPP_NODISCARD bool is_one() const;
784
785 // Detect sign.
sgn() const786 MPPP_NODISCARD int sgn() const
787 {
788 if (mppp_unlikely(nan_p())) {
789 // NOTE: mpfr_sgn() in this case would set an error flag, and we generally
790 // handle error flags as exceptions.
791 throw std::domain_error("Cannot determine the sign of a real NaN");
792 }
793 return mpfr_sgn(&m_mpfr);
794 }
795 // Get the sign bit.
signbit() const796 MPPP_NODISCARD bool signbit() const
797 {
798 return mpfr_signbit(&m_mpfr) != 0;
799 }
800
801 // Get the precision of this.
get_prec() const802 MPPP_NODISCARD ::mpfr_prec_t get_prec() const
803 {
804 return mpfr_get_prec(&m_mpfr);
805 }
806
807 private:
808 // Utility function to check precision in set_prec().
check_set_prec(::mpfr_prec_t p)809 static ::mpfr_prec_t check_set_prec(::mpfr_prec_t p)
810 {
811 if (mppp_unlikely(!detail::real_prec_check(p))) {
812 throw std::invalid_argument("Cannot set the precision of a real to the value " + detail::to_string(p)
813 + ": the maximum allowed precision is " + detail::to_string(real_prec_max())
814 + ", the minimum allowed precision is " + detail::to_string(real_prec_min()));
815 }
816 return p;
817 }
818 // mpfr_set_prec() wrapper, with or without prec checking.
819 template <bool Check>
set_prec_impl(::mpfr_prec_t p)820 void set_prec_impl(::mpfr_prec_t p)
821 {
822 ::mpfr_set_prec(&m_mpfr, Check ? check_set_prec(p) : p);
823 }
824 // mpfr_prec_round() wrapper, with or without prec checking.
825 template <bool Check>
prec_round_impl(::mpfr_prec_t p)826 void prec_round_impl(::mpfr_prec_t p)
827 {
828 ::mpfr_prec_round(&m_mpfr, Check ? check_set_prec(p) : p, MPFR_RNDN);
829 }
830
831 public:
832 // Destructively set the precision.
833 real &set_prec(::mpfr_prec_t);
834 // Set the precision maintaining the current value.
835 real &prec_round(::mpfr_prec_t);
836
837 private:
838 // Generic conversion.
839 // integer.
840 template <typename T, detail::enable_if_t<detail::is_integer<T>::value, int> = 0>
dispatch_conversion() const841 MPPP_NODISCARD T dispatch_conversion() const
842 {
843 if (mppp_unlikely(!number_p())) {
844 throw std::domain_error("Cannot convert a non-finite real to an integer");
845 }
846 MPPP_MAYBE_TLS detail::mpz_raii mpz;
847 // Truncate the value when converting to integer.
848 ::mpfr_get_z(&mpz.m_mpz, &m_mpfr, MPFR_RNDZ);
849 return T{&mpz.m_mpz};
850 }
851 // rational.
852 template <std::size_t SSize>
rational_conversion(rational<SSize> & rop) const853 bool rational_conversion(rational<SSize> &rop) const
854 {
855 #if defined(MPPP_MPFR_HAVE_MPFR_GET_Q)
856 MPPP_MAYBE_TLS detail::mpq_raii mpq;
857 // NOTE: we already checked outside
858 // that rop is a finite number, hence
859 // this function cannot fail.
860 ::mpfr_get_q(&mpq.m_mpq, &m_mpfr);
861 rop = &mpq.m_mpq;
862 return true;
863 #else
864 // Clear the range error flag before attempting the conversion.
865 ::mpfr_clear_erangeflag();
866 // NOTE: this call can fail if the exponent of this is very close to the upper/lower limits of the exponent
867 // type. If the call fails (signalled by a range flag being set), we will return error.
868 MPPP_MAYBE_TLS detail::mpz_raii mpz;
869 const ::mpfr_exp_t exp2 = ::mpfr_get_z_2exp(&mpz.m_mpz, &m_mpfr);
870 // NOTE: not sure at the moment how to trigger this, let's leave it for now.
871 // LCOV_EXCL_START
872 if (mppp_unlikely(::mpfr_erangeflag_p())) {
873 // Let's first reset the error flag.
874 ::mpfr_clear_erangeflag();
875 return false;
876 }
877 // LCOV_EXCL_STOP
878 // The conversion to n * 2**exp succeeded. We will build a rational
879 // from n and exp.
880 rop._get_num() = &mpz.m_mpz;
881 rop._get_den().set_one();
882 if (exp2 >= ::mpfr_exp_t(0)) {
883 // The output value will be an integer.
884 rop._get_num() <<= detail::make_unsigned(exp2);
885 } else {
886 // The output value will be a rational. Canonicalisation will be needed.
887 rop._get_den() <<= detail::nint_abs(exp2);
888 canonicalise(rop);
889 }
890 return true;
891 #endif
892 }
893 template <typename T, detail::enable_if_t<detail::is_rational<T>::value, int> = 0>
dispatch_conversion() const894 MPPP_NODISCARD T dispatch_conversion() const
895 {
896 if (mppp_unlikely(!number_p())) {
897 throw std::domain_error("Cannot convert a non-finite real to a rational");
898 }
899 T rop;
900 // LCOV_EXCL_START
901 if (mppp_unlikely(!rational_conversion(rop))) {
902 throw std::overflow_error("The exponent of a real is too large for conversion to rational");
903 }
904 // LCOV_EXCL_STOP
905 return rop;
906 }
907 // C++ floating-point.
908 template <typename T, detail::enable_if_t<std::is_floating_point<T>::value, int> = 0>
dispatch_conversion() const909 MPPP_NODISCARD T dispatch_conversion() const
910 {
911 if (std::is_same<T, float>::value) {
912 return static_cast<T>(::mpfr_get_flt(&m_mpfr, MPFR_RNDN));
913 }
914 if (std::is_same<T, double>::value) {
915 return static_cast<T>(::mpfr_get_d(&m_mpfr, MPFR_RNDN));
916 }
917 assert((std::is_same<T, long double>::value));
918 return static_cast<T>(::mpfr_get_ld(&m_mpfr, MPFR_RNDN));
919 }
920 // Small helper to raise an exception when converting to C++ integrals.
921 template <typename T>
raise_overflow_error() const922 [[noreturn]] void raise_overflow_error() const
923 {
924 throw std::overflow_error("Conversion of the real " + to_string() + " to the type '" + type_name<T>()
925 + "' results in overflow");
926 }
927 // Unsigned integrals, excluding bool.
928 template <typename T>
uint_conversion(T & rop) const929 bool uint_conversion(T &rop) const
930 {
931 // Clear the range error flag before attempting the conversion.
932 ::mpfr_clear_erangeflag();
933 // NOTE: this will handle correctly the case in which this is negative but greater than -1.
934 const unsigned long candidate = ::mpfr_get_ui(&m_mpfr, MPFR_RNDZ);
935 if (::mpfr_erangeflag_p()) {
936 // If the range error flag is set, it means the conversion failed because this is outside
937 // the range of unsigned long. Let's clear the error flag first.
938 ::mpfr_clear_erangeflag();
939 // If the value is positive, and the target type has a range greater than unsigned long,
940 // we will attempt the conversion again via integer.
941 if (detail::nl_max<T>() > detail::nl_max<unsigned long>() && sgn() > 0) {
942 return mppp::get(rop, static_cast<integer<2>>(*this));
943 }
944 return false;
945 }
946 if (candidate <= detail::nl_max<T>()) {
947 rop = static_cast<T>(candidate);
948 return true;
949 }
950 return false;
951 }
952 template <typename T,
953 detail::enable_if_t<detail::conjunction<detail::negation<std::is_same<bool, T>>, detail::is_integral<T>,
954 detail::is_unsigned<T>>::value,
955 int> = 0>
dispatch_conversion() const956 MPPP_NODISCARD T dispatch_conversion() const
957 {
958 if (mppp_unlikely(!number_p())) {
959 throw std::domain_error("Cannot convert a non-finite real to a C++ unsigned integral type");
960 }
961 T rop;
962 if (mppp_unlikely(!uint_conversion(rop))) {
963 raise_overflow_error<T>();
964 }
965 return rop;
966 }
967 // bool.
968 template <typename T, detail::enable_if_t<std::is_same<bool, T>::value, int> = 0>
dispatch_conversion() const969 MPPP_NODISCARD T dispatch_conversion() const
970 {
971 // NOTE: in C/C++ the conversion of NaN to bool gives true:
972 // https://stackoverflow.com/questions/9158567/nan-to-bool-conversion-true-or-false
973 return !zero_p();
974 }
975 // Signed integrals.
976 template <typename T>
sint_conversion(T & rop) const977 bool sint_conversion(T &rop) const
978 {
979 ::mpfr_clear_erangeflag();
980 const long candidate = ::mpfr_get_si(&m_mpfr, MPFR_RNDZ);
981 if (::mpfr_erangeflag_p()) {
982 // If the range error flag is set, it means the conversion failed because this is outside
983 // the range of long. Let's clear the error flag first.
984 ::mpfr_clear_erangeflag();
985 // If the target type has a range greater than long,
986 // we will attempt the conversion again via integer.
987 if (detail::nl_min<T>() < detail::nl_min<long>() && detail::nl_max<T>() > detail::nl_max<long>()) {
988 return mppp::get(rop, static_cast<integer<2>>(*this));
989 }
990 return false;
991 }
992 if (candidate >= detail::nl_min<T>() && candidate <= detail::nl_max<T>()) {
993 rop = static_cast<T>(candidate);
994 return true;
995 }
996 return false;
997 }
998 template <typename T,
999 detail::enable_if_t<detail::conjunction<detail::is_integral<T>, detail::is_signed<T>>::value, int> = 0>
dispatch_conversion() const1000 MPPP_NODISCARD T dispatch_conversion() const
1001 {
1002 if (mppp_unlikely(!number_p())) {
1003 throw std::domain_error("Cannot convert a non-finite real to a C++ signed integral type");
1004 }
1005 T rop;
1006 if (mppp_unlikely(!sint_conversion(rop))) {
1007 raise_overflow_error<T>();
1008 }
1009 return rop;
1010 }
1011 #if defined(MPPP_WITH_QUADMATH)
1012 template <typename T, detail::enable_if_t<std::is_same<real128, T>::value, int> = 0>
dispatch_conversion() const1013 MPPP_NODISCARD T dispatch_conversion() const
1014 {
1015 return convert_to_real128();
1016 }
1017 MPPP_NODISCARD real128 convert_to_real128() const;
1018 #endif
1019
1020 public:
1021 // Generic conversion operator.
1022 #if defined(MPPP_HAVE_CONCEPTS)
1023 template <real_interoperable T>
1024 #else
1025 template <typename T, detail::enable_if_t<is_real_interoperable<T>::value, int> = 0>
1026 #endif
operator T() const1027 explicit operator T() const
1028 {
1029 return dispatch_conversion<T>();
1030 }
1031 // Conversion to std::complex.
1032 #if defined(MPPP_HAVE_CONCEPTS)
1033 template <cpp_complex T>
1034 #else
1035 template <typename T, detail::enable_if_t<is_cpp_complex<T>::value, int> = 0>
1036 #endif
operator T() const1037 explicit operator T() const
1038 {
1039 using value_type = typename T::value_type;
1040
1041 return T{static_cast<value_type>(*this), value_type(0)};
1042 }
1043
1044 private:
1045 template <std::size_t SSize>
dispatch_get(integer<SSize> & rop) const1046 bool dispatch_get(integer<SSize> &rop) const
1047 {
1048 if (!number_p()) {
1049 return false;
1050 }
1051 MPPP_MAYBE_TLS detail::mpz_raii mpz;
1052 // Truncate the value when converting to integer.
1053 ::mpfr_get_z(&mpz.m_mpz, &m_mpfr, MPFR_RNDZ);
1054 rop = &mpz.m_mpz;
1055 return true;
1056 }
1057 template <std::size_t SSize>
dispatch_get(rational<SSize> & rop) const1058 bool dispatch_get(rational<SSize> &rop) const
1059 {
1060 if (!number_p()) {
1061 return false;
1062 }
1063 return rational_conversion(rop);
1064 }
dispatch_get(bool & b) const1065 bool dispatch_get(bool &b) const
1066 {
1067 b = !zero_p();
1068 return true;
1069 }
1070 template <typename T,
1071 detail::enable_if_t<detail::conjunction<detail::is_integral<T>, detail::is_unsigned<T>>::value, int> = 0>
dispatch_get(T & n) const1072 bool dispatch_get(T &n) const
1073 {
1074 if (!number_p()) {
1075 return false;
1076 }
1077 return uint_conversion(n);
1078 }
1079 template <typename T,
1080 detail::enable_if_t<detail::conjunction<detail::is_integral<T>, detail::is_signed<T>>::value, int> = 0>
dispatch_get(T & n) const1081 bool dispatch_get(T &n) const
1082 {
1083 if (!number_p()) {
1084 return false;
1085 }
1086 return sint_conversion(n);
1087 }
1088 template <typename T, detail::enable_if_t<std::is_floating_point<T>::value, int> = 0>
dispatch_get(T & x) const1089 bool dispatch_get(T &x) const
1090 {
1091 x = static_cast<T>(*this);
1092 return true;
1093 }
1094 #if defined(MPPP_WITH_QUADMATH)
1095 bool dispatch_get(real128 &) const;
1096 #endif
1097
1098 public:
1099 // Generic conversion function.
1100 #if defined(MPPP_HAVE_CONCEPTS)
1101 template <real_interoperable T>
1102 #else
1103 template <typename T, detail::enable_if_t<is_real_interoperable<T>::value, int> = 0>
1104 #endif
get(T & rop) const1105 bool get(T &rop) const
1106 {
1107 return dispatch_get(rop);
1108 }
1109 // Conversion function to std::complex.
1110 #if defined(MPPP_HAVE_CONCEPTS)
1111 template <cpp_complex T>
1112 #else
1113 template <typename T, detail::enable_if_t<is_cpp_complex<T>::value, int> = 0>
1114 #endif
get(T & rop) const1115 bool get(T &rop) const
1116 {
1117 rop = static_cast<T>(*this);
1118 return true;
1119 }
1120
1121 // Convert to string.
1122 MPPP_NODISCARD std::string to_string(int base = 10) const;
1123
1124 private:
1125 template <typename T>
1126 MPPP_DLL_LOCAL real &self_mpfr_unary(T &&);
1127 // Wrapper to apply the input unary MPFR function to this.
1128 // f must not need a rounding mode. Returns a reference to this.
1129 template <typename T>
self_mpfr_unary_nornd(T && f)1130 MPPP_DLL_LOCAL real &self_mpfr_unary_nornd(T &&f)
1131 {
1132 std::forward<T>(f)(&m_mpfr, &m_mpfr);
1133 return *this;
1134 }
1135
1136 public:
1137 // Negate in-place.
1138 real &neg();
1139 // In-place absolute value.
1140 real &abs();
1141
1142 // In-place square root.
1143 real &sqrt();
1144 // In-place reciprocal square root.
1145 real &rec_sqrt();
1146 #if defined(MPPP_WITH_ARB)
1147 // In-place sqrt1pm1.
1148 real &sqrt1pm1();
1149 #endif
1150 // In-place cubic root.
1151 real &cbrt();
1152
1153 // In-place squaring.
1154 real &sqr();
1155
1156 // In-place sine.
1157 real &sin();
1158 // In-place cosine.
1159 real &cos();
1160 // In-place tangent.
1161 real &tan();
1162 // In-place secant.
1163 real &sec();
1164 // In-place cosecant.
1165 real &csc();
1166 // In-place cotangent.
1167 real &cot();
1168 #if defined(MPPP_WITH_ARB)
1169 // Trig functions from Arb.
1170 real &sin_pi();
1171 real &cos_pi();
1172 real &tan_pi();
1173 real &cot_pi();
1174 real &sinc();
1175 real &sinc_pi();
1176 #endif
1177
1178 // In-place arccosine.
1179 real &acos();
1180 // In-place arcsine.
1181 real &asin();
1182 // In-place arctangent.
1183 real &atan();
1184
1185 // In-place hyperbolic cosine.
1186 real &cosh();
1187 // In-place hyperbolic sine.
1188 real &sinh();
1189 // In-place hyperbolic tangent.
1190 real &tanh();
1191 // In-place hyperbolic secant.
1192 real &sech();
1193 // In-place hyperbolic cosecant.
1194 real &csch();
1195 // In-place hyperbolic cotangent.
1196 real &coth();
1197
1198 // In-place inverse hyperbolic cosine.
1199 real &acosh();
1200 // In-place inverse hyperbolic sine.
1201 real &asinh();
1202 // In-place inverse hyperbolic tangent.
1203 real &atanh();
1204
1205 // In-place exponential.
1206 real &exp();
1207 // In-place base-2 exponential.
1208 real &exp2();
1209 // In-place base-10 exponential.
1210 real &exp10();
1211 // In-place exponential minus 1.
1212 real &expm1();
1213 // In-place logarithm.
1214 real &log();
1215 // In-place base-2 logarithm.
1216 real &log2();
1217 // In-place base-10 logarithm.
1218 real &log10();
1219 // In-place augmented logarithm.
1220 real &log1p();
1221
1222 // In-place Gamma function.
1223 real &gamma();
1224 // In-place logarithm of the Gamma function.
1225 real &lngamma();
1226 // In-place logarithm of the absolute value of the Gamma function.
1227 real &lgamma();
1228 // In-place Digamma function.
1229 real &digamma();
1230
1231 // In-place Bessel function of the first kind of order 0.
1232 real &j0();
1233 // In-place Bessel function of the first kind of order 1.
1234 real &j1();
1235 // In-place Bessel function of the second kind of order 0.
1236 real &y0();
1237 // In-place Bessel function of the second kind of order 1.
1238 real &y1();
1239
1240 // In-place exponential integral.
1241 real &eint();
1242 // In-place dilogarithm.
1243 real &li2();
1244 // In-place Riemann Zeta function.
1245 real &zeta();
1246 // In-place error function.
1247 real &erf();
1248 // In-place complementary error function.
1249 real &erfc();
1250 // In-place Airy function.
1251 real &ai();
1252 #if defined(MPPP_WITH_ARB)
1253 // In-place Lambert W function.
1254 real &lambert_w0();
1255 real &lambert_wm1();
1256 #endif
1257
1258 // In-place integer and remainder-related functions.
1259 real &ceil();
1260 real &floor();
1261 real &round();
1262 #if defined(MPPP_MPFR_HAVE_MPFR_ROUNDEVEN)
1263 real &roundeven();
1264 #endif
1265 real &trunc();
1266 real &frac();
1267
1268 #if defined(MPPP_MPFR_HAVE_MPFR_GET_STR_NDIGITS)
1269 // Get the number of significant digits required for a round-tripping representation.
1270 MPPP_NODISCARD std::size_t get_str_ndigits(int = 10) const;
1271 #endif
1272
1273 // Size of the serialised binary representation.
1274 MPPP_NODISCARD std::size_t binary_size() const;
1275
1276 private:
1277 void binary_save_impl(char *, std::size_t) const;
1278
1279 MPPP_DLL_LOCAL std::size_t binary_load_impl(const char *);
1280 std::size_t binary_load_impl(const char *, std::size_t, const char *);
1281
1282 public:
1283 std::size_t binary_save(char *) const;
1284 std::size_t binary_save(std::vector<char> &) const;
1285 template <std::size_t S>
binary_save(std::array<char,S> & dest) const1286 std::size_t binary_save(std::array<char, S> &dest) const
1287 {
1288 const auto bs = binary_size();
1289 if (bs > S) {
1290 return 0;
1291 }
1292 binary_save_impl(dest.data(), bs);
1293 return bs;
1294 }
1295 std::size_t binary_save(std::ostream &) const;
1296
1297 std::size_t binary_load(const char *);
1298 std::size_t binary_load(const std::vector<char> &);
1299 template <std::size_t S>
binary_load(const std::array<char,S> & src)1300 std::size_t binary_load(const std::array<char, S> &src)
1301 {
1302 return binary_load_impl(src.data(), detail::safe_cast<std::size_t>(src.size()), "std::array");
1303 }
1304 std::size_t binary_load(std::istream &);
1305
1306 private:
1307 mpfr_struct_t m_mpfr;
1308 };
1309
1310 template <typename T, typename U>
1311 using are_real_op_types
1312 = detail::disjunction<detail::conjunction<is_cvr_real<T>, is_cvr_real<U>>,
1313 detail::conjunction<is_cvr_real<T>, is_real_interoperable<detail::uncvref_t<U>>>,
1314 detail::conjunction<is_cvr_real<U>, is_real_interoperable<detail::uncvref_t<T>>>>;
1315
1316 template <typename T, typename U>
1317 using are_real_in_place_op_types
1318 = detail::conjunction<detail::negation<std::is_const<detail::unref_t<T>>>, are_real_op_types<T, U>>;
1319
1320 #if defined(MPPP_HAVE_CONCEPTS)
1321
1322 template <typename T, typename U>
1323 MPPP_CONCEPT_DECL real_op_types = are_real_op_types<T, U>::value;
1324
1325 template <typename T, typename U>
1326 MPPP_CONCEPT_DECL real_in_place_op_types = are_real_in_place_op_types<T, U>::value;
1327
1328 #endif
1329
1330 // Destructively set the precision.
set_prec(real & r,::mpfr_prec_t p)1331 inline void set_prec(real &r, ::mpfr_prec_t p)
1332 {
1333 r.set_prec(p);
1334 }
1335
1336 // Set the precision.
prec_round(real & r,::mpfr_prec_t p)1337 inline void prec_round(real &r, ::mpfr_prec_t p)
1338 {
1339 r.prec_round(p);
1340 }
1341
1342 // Get the precision.
get_prec(const real & r)1343 inline mpfr_prec_t get_prec(const real &r)
1344 {
1345 return r.get_prec();
1346 }
1347
1348 namespace detail
1349 {
1350
1351 template <typename... Args>
1352 using real_set_t = decltype(std::declval<real &>().set(std::declval<const Args &>()...));
1353
1354 } // namespace detail
1355
1356 #if defined(MPPP_HAVE_CONCEPTS)
1357
1358 template <typename... Args>
1359 MPPP_CONCEPT_DECL real_set_args = detail::is_detected<detail::real_set_t, Args...>::value;
1360
1361 #endif
1362
1363 // Generic setter.
1364 #if defined(MPPP_HAVE_CONCEPTS)
1365 template <real_set_args... Args>
1366 #else
1367 template <typename... Args, detail::enable_if_t<detail::is_detected<detail::real_set_t, Args...>::value, int> = 0>
1368 #endif
set(real & r,const Args &...args)1369 inline real &set(real &r, const Args &...args)
1370 {
1371 return r.set(args...);
1372 }
1373
1374 // Set to n*2**e.
1375 template <std::size_t SSize>
set_z_2exp(real & r,const integer<SSize> & n,::mpfr_exp_t e)1376 inline real &set_z_2exp(real &r, const integer<SSize> &n, ::mpfr_exp_t e)
1377 {
1378 ::mpfr_set_z_2exp(r._get_mpfr_t(), n.get_mpz_view(), e, MPFR_RNDN);
1379 return r;
1380 }
1381
1382 MPPP_DLL_PUBLIC real &set_ui_2exp(real &, unsigned long, ::mpfr_exp_t);
1383 MPPP_DLL_PUBLIC real &set_si_2exp(real &, long, ::mpfr_exp_t);
1384
1385 // Implementation of the constructor from n*2**e, integer overload.
1386 // Place it here so that set_z_2exp() is visible.
1387 template <std::size_t SSize>
1388 // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init, hicpp-member-init)
real(const integer<SSize> & n,::mpfr_exp_t e,::mpfr_prec_t p)1389 inline real::real(const integer<SSize> &n, ::mpfr_exp_t e, ::mpfr_prec_t p)
1390 {
1391 ::mpfr_init2(&m_mpfr, check_init_prec(p));
1392 set_z_2exp(*this, n, e);
1393 }
1394
1395 // Set to NaN.
set_nan(real & r)1396 inline real &set_nan(real &r)
1397 {
1398 return r.set_nan();
1399 }
1400
1401 // Set to infinity.
set_inf(real & r,int sign=0)1402 inline real &set_inf(real &r, int sign = 0)
1403 {
1404 return r.set_inf(sign);
1405 }
1406
1407 // Set to zero.
set_zero(real & r,int sign=0)1408 inline real &set_zero(real &r, int sign = 0)
1409 {
1410 return r.set_zero(sign);
1411 }
1412
1413 // Swap.
swap(real & a,real & b)1414 inline void swap(real &a, real &b) noexcept
1415 {
1416 ::mpfr_swap(a._get_mpfr_t(), b._get_mpfr_t());
1417 }
1418
1419 // Generic conversion functions.
1420 #if defined(MPPP_HAVE_CONCEPTS)
1421 template <real_interoperable T>
1422 #else
1423 template <typename T, detail::enable_if_t<is_real_interoperable<T>::value, int> = 0>
1424 #endif
get(T & rop,const real & x)1425 inline bool get(T &rop, const real &x)
1426 {
1427 return x.get(rop);
1428 }
1429
1430 #if defined(MPPP_HAVE_CONCEPTS)
1431 template <cpp_complex T>
1432 #else
1433 template <typename T, detail::enable_if_t<is_cpp_complex<T>::value, int> = 0>
1434 #endif
get(T & rop,const real & x)1435 inline bool get(T &rop, const real &x)
1436 {
1437 return x.get(rop);
1438 }
1439
1440 // Extract significand and exponent.
1441 template <std::size_t SSize>
get_z_2exp(integer<SSize> & n,const real & r)1442 inline mpfr_exp_t get_z_2exp(integer<SSize> &n, const real &r)
1443 {
1444 if (mppp_unlikely(!r.number_p())) {
1445 throw std::domain_error("Cannot extract the significand and the exponent of a non-finite real");
1446 }
1447 MPPP_MAYBE_TLS detail::mpz_raii m;
1448 ::mpfr_clear_erangeflag();
1449 auto retval = ::mpfr_get_z_2exp(&m.m_mpz, r.get_mpfr_t());
1450 // LCOV_EXCL_START
1451 if (mppp_unlikely(::mpfr_erangeflag_p())) {
1452 ::mpfr_clear_erangeflag();
1453 throw std::overflow_error("Cannot extract the exponent of the real value " + r.to_string()
1454 + ": the exponent's magnitude is too large");
1455 }
1456 // LCOV_EXCL_STOP
1457 n = &m.m_mpz;
1458 return retval;
1459 }
1460
1461 namespace detail
1462 {
1463
1464 // A small helper to init the pairs in the functions below. We need this because
1465 // we cannot take the address of a const real as a real *.
1466 template <typename Arg, enable_if_t<!is_ncrvr<Arg &&>::value, int> = 0>
mpfr_nary_op_init_pair(::mpfr_prec_t min_prec,Arg && arg)1467 inline std::pair<real *, ::mpfr_prec_t> mpfr_nary_op_init_pair(::mpfr_prec_t min_prec, Arg &&arg)
1468 {
1469 // arg is not a non-const rvalue ref, we cannot steal from it. Init with nullptr.
1470 return std::make_pair(static_cast<real *>(nullptr), c_max(arg.get_prec(), min_prec));
1471 }
1472
1473 template <typename Arg, enable_if_t<is_ncrvr<Arg &&>::value, int> = 0>
mpfr_nary_op_init_pair(::mpfr_prec_t min_prec,Arg && arg)1474 inline std::pair<real *, ::mpfr_prec_t> mpfr_nary_op_init_pair(::mpfr_prec_t min_prec, Arg &&arg)
1475 {
1476 // arg is a non-const rvalue ref, and a candidate for stealing resources.
1477 return std::make_pair(&arg, c_max(arg.get_prec(), min_prec));
1478 }
1479
1480 // A recursive function to determine, in an MPFR function call,
1481 // the largest argument we can steal resources from, and the max precision among
1482 // all the arguments.
mpfr_nary_op_check_steal(std::pair<real *,::mpfr_prec_t> &)1483 inline void mpfr_nary_op_check_steal(std::pair<real *, ::mpfr_prec_t> &) {}
1484
1485 // NOTE: we need 2 overloads for this, as we cannot extract a non-const pointer from
1486 // arg0 if arg0 is a const ref.
1487 template <typename Arg0, typename... Args, enable_if_t<!is_ncrvr<Arg0 &&>::value, int> = 0>
1488 void mpfr_nary_op_check_steal(std::pair<real *, ::mpfr_prec_t> &, Arg0 &&, Args &&...);
1489
1490 template <typename Arg0, typename... Args, enable_if_t<is_ncrvr<Arg0 &&>::value, int> = 0>
1491 void mpfr_nary_op_check_steal(std::pair<real *, ::mpfr_prec_t> &, Arg0 &&, Args &&...);
1492
1493 template <typename Arg0, typename... Args, enable_if_t<!is_ncrvr<Arg0 &&>::value, int>>
mpfr_nary_op_check_steal(std::pair<real *,::mpfr_prec_t> & p,Arg0 && arg0,Args &&...args)1494 inline void mpfr_nary_op_check_steal(std::pair<real *, ::mpfr_prec_t> &p, Arg0 &&arg0, Args &&...args)
1495 {
1496 // arg0 is not a non-const rvalue ref, we won't be able to steal from it regardless. Just
1497 // update the max prec.
1498 p.second = c_max(arg0.get_prec(), p.second);
1499 mpfr_nary_op_check_steal(p, std::forward<Args>(args)...);
1500 }
1501
1502 template <typename Arg0, typename... Args, enable_if_t<is_ncrvr<Arg0 &&>::value, int>>
mpfr_nary_op_check_steal(std::pair<real *,::mpfr_prec_t> & p,Arg0 && arg0,Args &&...args)1503 inline void mpfr_nary_op_check_steal(std::pair<real *, ::mpfr_prec_t> &p, Arg0 &&arg0, Args &&...args)
1504 {
1505 const auto prec0 = arg0.get_prec();
1506 if (!p.first || prec0 > p.first->get_prec()) {
1507 // The current argument arg0 is a non-const rvalue reference, and either it's
1508 // the first argument we encounter we can steal from, or it has a precision
1509 // larger than the current candidate for stealing resources from. This means that
1510 // arg0 is the new candidate.
1511 p.first = &arg0;
1512 }
1513 // Update the max precision among the arguments, if necessary.
1514 p.second = c_max(prec0, p.second);
1515 mpfr_nary_op_check_steal(p, std::forward<Args>(args)...);
1516 }
1517
1518 // A small wrapper to call an MPFR function f with arguments args. If the first param is true_type,
1519 // the rounding mode MPFR_RNDN will be appended at the end of the function arguments list.
1520 template <typename F, typename... Args>
mpfr_nary_func_wrapper(const std::true_type &,const F & f,Args &&...args)1521 inline void mpfr_nary_func_wrapper(const std::true_type &, const F &f, Args &&...args)
1522 {
1523 f(std::forward<Args>(args)..., MPFR_RNDN);
1524 }
1525
1526 template <typename F, typename... Args>
mpfr_nary_func_wrapper(const std::false_type &,const F & f,Args &&...args)1527 inline void mpfr_nary_func_wrapper(const std::false_type &, const F &f, Args &&...args)
1528 {
1529 f(std::forward<Args>(args)...);
1530 }
1531
1532 // The goal of this helper is to invoke the MPFR-like function object f with signature
1533 //
1534 // void f(mpfr_t out, const mpfr_t x0, const mpfr_t x1, ...)
1535 //
1536 // on the mpfr_t instances contained in the input real objects,
1537 //
1538 // f(rop._get_mpfr_t(), arg0.get_mpfr_t(), arg1.get_mpfr_t(), ...)
1539 //
1540 // The helper will ensure that, before the invocation, the precision
1541 // of rop is set to max(min_prec, arg0.get_prec(), arg1.get_prec(), ...).
1542 //
1543 // One of the input arguments may be used as return value in the invocation
1544 // instead of rop if it provides enough precision and it is passed as a non-const
1545 // rvalue reference. In such a case, the selected input argument will be swapped
1546 // into rop after the invocation and before the function returns.
1547 //
1548 // The Rnd flag controls whether to add the rounding mode (MPFR_RNDN) at the end
1549 // of the MPFR-like function object arguments list or not.
1550 //
1551 // This function requires that the MPFR-like function object being called supports
1552 // overlapping arguments (both input and output).
1553 template <bool Rnd, typename F, typename Arg0, typename... Args>
mpfr_nary_op_impl(::mpfr_prec_t min_prec,const F & f,real & rop,Arg0 && arg0,Args &&...args)1554 inline real &mpfr_nary_op_impl(::mpfr_prec_t min_prec, const F &f, real &rop, Arg0 &&arg0, Args &&...args)
1555 {
1556 // Make sure min_prec is valid.
1557 // NOTE: min_prec == 0 is ok, it just means
1558 // p below will be inited with arg0's precision
1559 // rather than min_prec.
1560 assert(min_prec == 0 || real_prec_check(min_prec));
1561
1562 // This pair will contain:
1563 //
1564 // - a pointer to the largest-precision arg from which we can steal resources (may be nullptr),
1565 // - the largest precision among all args and min_prec (i.e., the target precision
1566 // for rop).
1567 //
1568 // It is inited with arg0's precision (but no less than min_prec), and a pointer to arg0, if arg0 is a nonconst
1569 // rvalue ref (a nullptr otherwise).
1570 auto p = mpfr_nary_op_init_pair(min_prec, std::forward<Arg0>(arg0));
1571 // Finish setting up p by examining the remaining arguments.
1572 mpfr_nary_op_check_steal(p, std::forward<Args>(args)...);
1573
1574 // Cache for convenience.
1575 const auto r_prec = rop.get_prec();
1576
1577 if (p.second == r_prec) {
1578 // The target precision and the precision of the return value
1579 // match. No need to steal, just execute the function.
1580 mpfr_nary_func_wrapper(std::integral_constant<bool, Rnd>{}, f, rop._get_mpfr_t(), arg0.get_mpfr_t(),
1581 args.get_mpfr_t()...);
1582 } else {
1583 if (r_prec > p.second) {
1584 // The precision of the return value is larger than the target precision.
1585 // We can reset its precision destructively
1586 // because we know it does not overlap with any operand.
1587 rop.set_prec_impl<false>(p.second);
1588 mpfr_nary_func_wrapper(std::integral_constant<bool, Rnd>{}, f, rop._get_mpfr_t(), arg0.get_mpfr_t(),
1589 args.get_mpfr_t()...);
1590 } else if (p.first && p.first->get_prec() == p.second) {
1591 // The precision of the return value is smaller than the target precision,
1592 // and we have a candidate for stealing with enough precision: we will use it as return
1593 // value and then swap out the result to rop.
1594 mpfr_nary_func_wrapper(std::integral_constant<bool, Rnd>{}, f, p.first->_get_mpfr_t(), arg0.get_mpfr_t(),
1595 args.get_mpfr_t()...);
1596 swap(*p.first, rop);
1597 } else {
1598 // The precision of the return value is smaller than the target precision,
1599 // and either:
1600 //
1601 // - we cannot steal from any argument, or
1602 // - we can steal from an argument but the selected argument
1603 // does not have enough precision.
1604 //
1605 // In these cases, we will just set the precision of rop and call the function.
1606 //
1607 // NOTE: we need to set the precision without destroying the rop, as rop might
1608 // overlap with one of the arguments. Since this will be an increase in precision,
1609 // it should not entail a rounding operation.
1610 //
1611 // NOTE: we assume all the precs in the operands and min_prec are valid, so
1612 // we will not need to check them.
1613 rop.prec_round_impl<false>(p.second);
1614 mpfr_nary_func_wrapper(std::integral_constant<bool, Rnd>{}, f, rop._get_mpfr_t(), arg0.get_mpfr_t(),
1615 args.get_mpfr_t()...);
1616 }
1617 }
1618
1619 return rop;
1620 }
1621
1622 // The goal of this helper is to invoke the MPFR-like function object f with signature
1623 //
1624 // void f(mpfr_t out, const mpfr_t x0, const mpfr_t x1, ...)
1625 //
1626 // on the mpfr_t instances contained in the input real objects,
1627 //
1628 // f(rop._get_mpfr_t(), arg0.get_mpfr_t(), arg1.get_mpfr_t(), ...)
1629 //
1630 // and then return rop.
1631 //
1632 // The rop object will either be created within the helper with a precision
1633 // set to max(min_prec, arg0.get_prec(), arg1.get_prec(), ...),
1634 // or it will be one of the input arguments if it provides enough precision and
1635 // it is passed as a non-const rvalue reference.
1636 //
1637 // The Rnd flag controls whether to add the rounding mode (MPFR_RNDN) at the end
1638 // of the MPFR-like function object arguments list or not.
1639 //
1640 // This function requires that the MPFR-like function object being called supports
1641 // overlapping arguments (both input and output).
1642 template <bool Rnd, typename F, typename Arg0, typename... Args>
mpfr_nary_op_return_impl(::mpfr_prec_t min_prec,const F & f,Arg0 && arg0,Args &&...args)1643 inline real mpfr_nary_op_return_impl(::mpfr_prec_t min_prec, const F &f, Arg0 &&arg0, Args &&...args)
1644 {
1645 // Make sure min_prec is valid.
1646 // NOTE: min_prec == 0 is ok, it just means
1647 // p below will be inited with arg0's precision
1648 // rather than min_prec.
1649 assert(min_prec == 0 || real_prec_check(min_prec));
1650
1651 // This pair will contain:
1652 //
1653 // - a pointer to the largest-precision arg from which we can steal resources (may be nullptr),
1654 // - the largest precision among all args and min_prec (i.e., the target precision
1655 // for the return value).
1656 //
1657 // It is inited with arg0's precision (but no less than min_prec), and a pointer to arg0, if arg0 is a nonconst
1658 // rvalue ref (a nullptr otherwise).
1659 auto p = mpfr_nary_op_init_pair(min_prec, std::forward<Arg0>(arg0));
1660 // Finish setting up p by examining the remaining arguments.
1661 mpfr_nary_op_check_steal(p, std::forward<Args>(args)...);
1662
1663 if (p.first && p.first->get_prec() == p.second) {
1664 // We can steal from one or more args, and the precision of
1665 // the largest-precision arg we can steal from matches
1666 // the target precision. Use it.
1667 mpfr_nary_func_wrapper(std::integral_constant<bool, Rnd>{}, f, p.first->_get_mpfr_t(), arg0.get_mpfr_t(),
1668 args.get_mpfr_t()...);
1669 return std::move(*p.first);
1670 } else {
1671 // Either we cannot steal from any arg, or the candidate does not have
1672 // enough precision. Init a new value and use it instead.
1673 real retval{real::ptag{}, p.second, true};
1674 mpfr_nary_func_wrapper(std::integral_constant<bool, Rnd>{}, f, retval._get_mpfr_t(), arg0.get_mpfr_t(),
1675 args.get_mpfr_t()...);
1676 return retval;
1677 }
1678 }
1679
1680 } // namespace detail
1681
1682 // Ternary addition.
1683 #if defined(MPPP_HAVE_CONCEPTS)
1684 template <cvr_real T, cvr_real U>
1685 #else
1686 template <typename T, typename U, cvr_real_enabler<T, U> = 0>
1687 #endif
add(real & rop,T && a,U && b)1688 inline real &add(real &rop, T &&a, U &&b)
1689 {
1690 return detail::mpfr_nary_op_impl<true>(0, ::mpfr_add, rop, std::forward<T>(a), std::forward<U>(b));
1691 }
1692
1693 // Ternary subtraction.
1694 #if defined(MPPP_HAVE_CONCEPTS)
1695 template <cvr_real T, cvr_real U>
1696 #else
1697 template <typename T, typename U, cvr_real_enabler<T, U> = 0>
1698 #endif
sub(real & rop,T && a,U && b)1699 inline real &sub(real &rop, T &&a, U &&b)
1700 {
1701 return detail::mpfr_nary_op_impl<true>(0, ::mpfr_sub, rop, std::forward<T>(a), std::forward<U>(b));
1702 }
1703
1704 // Ternary multiplication.
1705 #if defined(MPPP_HAVE_CONCEPTS)
1706 template <cvr_real T, cvr_real U>
1707 #else
1708 template <typename T, typename U, cvr_real_enabler<T, U> = 0>
1709 #endif
mul(real & rop,T && a,U && b)1710 inline real &mul(real &rop, T &&a, U &&b)
1711 {
1712 return detail::mpfr_nary_op_impl<true>(0, ::mpfr_mul, rop, std::forward<T>(a), std::forward<U>(b));
1713 }
1714
1715 // Ternary division.
1716 #if defined(MPPP_HAVE_CONCEPTS)
1717 template <cvr_real T, cvr_real U>
1718 #else
1719 template <typename T, typename U, cvr_real_enabler<T, U> = 0>
1720 #endif
div(real & rop,T && a,U && b)1721 inline real &div(real &rop, T &&a, U &&b)
1722 {
1723 return detail::mpfr_nary_op_impl<true>(0, ::mpfr_div, rop, std::forward<T>(a), std::forward<U>(b));
1724 }
1725
1726 // Quaternary fused multiply–add.
1727 #if defined(MPPP_HAVE_CONCEPTS)
1728 template <cvr_real T, cvr_real U, cvr_real V>
1729 #else
1730 template <typename T, typename U, typename V, cvr_real_enabler<T, U, V> = 0>
1731 #endif
fma(real & rop,T && a,U && b,V && c)1732 inline real &fma(real &rop, T &&a, U &&b, V &&c)
1733 {
1734 return detail::mpfr_nary_op_impl<true>(0, ::mpfr_fma, rop, std::forward<T>(a), std::forward<U>(b),
1735 std::forward<V>(c));
1736 }
1737
1738 // Quaternary fused multiply–sub.
1739 #if defined(MPPP_HAVE_CONCEPTS)
1740 template <cvr_real T, cvr_real U, cvr_real V>
1741 #else
1742 template <typename T, typename U, typename V, cvr_real_enabler<T, U, V> = 0>
1743 #endif
fms(real & rop,T && a,U && b,V && c)1744 inline real &fms(real &rop, T &&a, U &&b, V &&c)
1745 {
1746 return detail::mpfr_nary_op_impl<true>(0, ::mpfr_fms, rop, std::forward<T>(a), std::forward<U>(b),
1747 std::forward<V>(c));
1748 }
1749
1750 // Ternary fused multiply–add.
1751 #if defined(MPPP_HAVE_CONCEPTS)
1752 template <cvr_real T, cvr_real U, cvr_real V>
1753 #else
1754 template <typename T, typename U, typename V, cvr_real_enabler<T, U, V> = 0>
1755 #endif
fma(T && a,U && b,V && c)1756 inline real fma(T &&a, U &&b, V &&c)
1757 {
1758 return detail::mpfr_nary_op_return_impl<true>(0, ::mpfr_fma, std::forward<T>(a), std::forward<U>(b),
1759 std::forward<V>(c));
1760 }
1761
1762 // Ternary fused multiply–sub.
1763 #if defined(MPPP_HAVE_CONCEPTS)
1764 template <cvr_real T, cvr_real U, cvr_real V>
1765 #else
1766 template <typename T, typename U, typename V, cvr_real_enabler<T, U, V> = 0>
1767 #endif
fms(T && a,U && b,V && c)1768 inline real fms(T &&a, U &&b, V &&c)
1769 {
1770 return detail::mpfr_nary_op_return_impl<true>(0, ::mpfr_fms, std::forward<T>(a), std::forward<U>(b),
1771 std::forward<V>(c));
1772 }
1773
1774 // mul2/div2 primitives.
1775 #if defined(MPPP_HAVE_CONCEPTS)
1776 template <cvr_real T>
1777 #else
1778 template <typename T, cvr_real_enabler<T> = 0>
1779 #endif
mul_2ui(T && x,unsigned long n)1780 inline real mul_2ui(T &&x, unsigned long n)
1781 {
1782 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2ui(r, o, n, MPFR_RNDN); };
1783 return detail::mpfr_nary_op_return_impl<false>(0, wrapper, std::forward<T>(x));
1784 }
1785
1786 #if defined(MPPP_HAVE_CONCEPTS)
1787 template <cvr_real T>
1788 #else
1789 template <typename T, cvr_real_enabler<T> = 0>
1790 #endif
mul_2ui(real & rop,T && x,unsigned long n)1791 inline real &mul_2ui(real &rop, T &&x, unsigned long n)
1792 {
1793 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2ui(r, o, n, MPFR_RNDN); };
1794 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(x));
1795 }
1796
1797 #if defined(MPPP_HAVE_CONCEPTS)
1798 template <cvr_real T>
1799 #else
1800 template <typename T, cvr_real_enabler<T> = 0>
1801 #endif
mul_2si(T && x,long n)1802 inline real mul_2si(T &&x, long n)
1803 {
1804 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2si(r, o, n, MPFR_RNDN); };
1805 return detail::mpfr_nary_op_return_impl<false>(0, wrapper, std::forward<T>(x));
1806 }
1807
1808 #if defined(MPPP_HAVE_CONCEPTS)
1809 template <cvr_real T>
1810 #else
1811 template <typename T, cvr_real_enabler<T> = 0>
1812 #endif
mul_2si(real & rop,T && x,long n)1813 inline real &mul_2si(real &rop, T &&x, long n)
1814 {
1815 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2si(r, o, n, MPFR_RNDN); };
1816 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(x));
1817 }
1818
1819 #if defined(MPPP_HAVE_CONCEPTS)
1820 template <cvr_real T>
1821 #else
1822 template <typename T, cvr_real_enabler<T> = 0>
1823 #endif
div_2ui(T && x,unsigned long n)1824 inline real div_2ui(T &&x, unsigned long n)
1825 {
1826 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2ui(r, o, n, MPFR_RNDN); };
1827 return detail::mpfr_nary_op_return_impl<false>(0, wrapper, std::forward<T>(x));
1828 }
1829
1830 #if defined(MPPP_HAVE_CONCEPTS)
1831 template <cvr_real T>
1832 #else
1833 template <typename T, cvr_real_enabler<T> = 0>
1834 #endif
div_2ui(real & rop,T && x,unsigned long n)1835 inline real &div_2ui(real &rop, T &&x, unsigned long n)
1836 {
1837 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2ui(r, o, n, MPFR_RNDN); };
1838 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(x));
1839 }
1840
1841 #if defined(MPPP_HAVE_CONCEPTS)
1842 template <cvr_real T>
1843 #else
1844 template <typename T, cvr_real_enabler<T> = 0>
1845 #endif
div_2si(T && x,long n)1846 inline real div_2si(T &&x, long n)
1847 {
1848 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2si(r, o, n, MPFR_RNDN); };
1849 return detail::mpfr_nary_op_return_impl<false>(0, wrapper, std::forward<T>(x));
1850 }
1851
1852 #if defined(MPPP_HAVE_CONCEPTS)
1853 template <cvr_real T>
1854 #else
1855 template <typename T, cvr_real_enabler<T> = 0>
1856 #endif
div_2si(real & rop,T && x,long n)1857 inline real &div_2si(real &rop, T &&x, long n)
1858 {
1859 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2si(r, o, n, MPFR_RNDN); };
1860 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(x));
1861 }
1862
1863 // Detect NaN.
nan_p(const real & r)1864 inline bool nan_p(const real &r)
1865 {
1866 return r.nan_p();
1867 }
1868
1869 // Detect inf.
inf_p(const real & r)1870 inline bool inf_p(const real &r)
1871 {
1872 return r.inf_p();
1873 }
1874
1875 // Detect finite.
number_p(const real & r)1876 inline bool number_p(const real &r)
1877 {
1878 return r.number_p();
1879 }
1880
1881 // Detect zero.
zero_p(const real & r)1882 inline bool zero_p(const real &r)
1883 {
1884 return r.zero_p();
1885 }
1886
1887 // Detect regular number.
regular_p(const real & r)1888 inline bool regular_p(const real &r)
1889 {
1890 return r.regular_p();
1891 }
1892
1893 // Detect integral value.
integer_p(const real & r)1894 inline bool integer_p(const real &r)
1895 {
1896 return r.integer_p();
1897 }
1898
1899 // Detect one.
is_one(const real & r)1900 inline bool is_one(const real &r)
1901 {
1902 return r.is_one();
1903 }
1904
1905 // Detect the sign.
sgn(const real & r)1906 inline int sgn(const real &r)
1907 {
1908 return r.sgn();
1909 }
1910
1911 // Get the sign bit.
signbit(const real & r)1912 inline bool signbit(const real &r)
1913 {
1914 return r.signbit();
1915 }
1916
1917 // Comparison.
1918 MPPP_DLL_PUBLIC int cmp(const real &, const real &);
1919
1920 // Comparison of absolute values.
1921 MPPP_DLL_PUBLIC int cmpabs(const real &, const real &);
1922
1923 // Comparison with integral multiples of powers of 2.
1924 MPPP_DLL_PUBLIC int cmp_ui_2exp(const real &, unsigned long, ::mpfr_exp_t);
1925 MPPP_DLL_PUBLIC int cmp_si_2exp(const real &, long, ::mpfr_exp_t);
1926
1927 // Equality predicate with special NaN handling.
1928 MPPP_DLL_PUBLIC bool real_equal_to(const real &, const real &);
1929
1930 // Less-than predicate with special NaN and moved-from handling.
1931 MPPP_DLL_PUBLIC bool real_lt(const real &, const real &);
1932
1933 // Greater-than predicate with special NaN and moved-from handling.
1934 MPPP_DLL_PUBLIC bool real_gt(const real &, const real &);
1935
1936 // These are helper macros to reduce typing when dealing with the common case
1937 // of exposing MPFR-like functions with a single argument (both variants with retval
1938 // and with return). "name" will be the name of the mppp function, "fname" is
1939 // the name of the MPFR-like function and "rnd" is a boolean flag that signals whether
1940 // fname requires a rounding mode argument or not.
1941 // The fname function must accept only mpfr_t arguments in input (plus the rounding mode if
1942 // rnd is true).
1943
1944 // These are the headers of the overloads that will be produced. They are different depending
1945 // on whether concepts are available or not.
1946 #if defined(MPPP_HAVE_CONCEPTS)
1947 #define MPPP_REAL_MPFR_UNARY_HEADER template <cvr_real T>
1948 #else
1949 #define MPPP_REAL_MPFR_UNARY_HEADER template <typename T, cvr_real_enabler<T> = 0>
1950 #endif
1951
1952 #define MPPP_REAL_MPFR_UNARY_IMPL(name, fname, rnd) \
1953 MPPP_REAL_MPFR_UNARY_HEADER inline real &name(real &rop, T &&op) \
1954 { \
1955 return detail::mpfr_nary_op_impl<rnd>(0, fname, rop, std::forward<T>(op)); \
1956 } \
1957 MPPP_REAL_MPFR_UNARY_HEADER inline real name(T &&r) \
1958 { \
1959 return detail::mpfr_nary_op_return_impl<rnd>(0, fname, std::forward<T>(r)); \
1960 }
1961
1962 // Machinery to expose the binary MPFR-like function fname as an mppp function called "name".
1963 //
1964 // Two overloads of "name" will be provided:
1965 //
1966 // - an overload in which the return value is passed by
1967 // reference as the first argument,
1968 // - an overload which returns the result.
1969 //
1970 // The first overload accepts only real arguments.
1971 //
1972 // The second overload is generic and it accepts 2 input arguments, at least one of which must be a
1973 // real. The other argument, if not real, will be converted to real in the usual way.
1974 //
1975 // The rnd param (a boolean) indicates if fname requires a rounding mode as last argument or not.
1976 //
1977 // The fname function must accept only mpfr_t arguments in input (plus the rounding mode if
1978 // rnd is true).
1979
1980 // These are the headers of the overloads that will be produced. They are different depending
1981 // on whether concepts are available or not.
1982 #if defined(MPPP_HAVE_CONCEPTS)
1983 #define MPPP_REAL_MPFR_BINARY_HEADER1 template <cvr_real T, cvr_real U>
1984 #define MPPP_REAL_MPFR_BINARY_HEADER2 template <typename T, real_op_types<T> U>
1985 #else
1986 #define MPPP_REAL_MPFR_BINARY_HEADER1 template <typename T, typename U, cvr_real_enabler<T, U> = 0>
1987 #define MPPP_REAL_MPFR_BINARY_HEADER2 \
1988 template <typename T, typename U, detail::enable_if_t<are_real_op_types<U, T>::value, int> = 0>
1989 #endif
1990
1991 // The actual macro.
1992 #define MPPP_REAL_MPFR_BINARY_IMPL(name, fname, rnd) \
1993 /* The overload which accepts the return value in input. */ \
1994 MPPP_REAL_MPFR_BINARY_HEADER1 inline real &name(real &rop, T &&y, U &&x) \
1995 { \
1996 return detail::mpfr_nary_op_impl<rnd>(0, fname, rop, std::forward<T>(y), std::forward<U>(x)); \
1997 } \
1998 /* Implementation details of the other overload. */ \
1999 namespace detail \
2000 { \
2001 /* Both arguments are real. */ \
2002 template <typename T, typename U, cvr_real_enabler<T, U> = 0> \
2003 inline real dispatch_##name(T &&y, U &&x) \
2004 { \
2005 return mpfr_nary_op_return_impl<rnd>(0, fname, std::forward<T>(y), std::forward<U>(x)); \
2006 } \
2007 /* Only the first argument is real. */ \
2008 template <typename T, typename U, \
2009 enable_if_t<conjunction<is_cvr_real<T>, is_real_interoperable<U>>::value, int> = 0> \
2010 inline real dispatch_##name(T &&a, const U &x) \
2011 { \
2012 MPPP_MAYBE_TLS real tmp; \
2013 tmp.set_prec(c_max(a.get_prec(), real_deduce_precision(x))); \
2014 tmp.set(x); \
2015 return dispatch_##name(std::forward<T>(a), tmp); \
2016 } \
2017 /* Only the second argument is real. */ \
2018 template <typename T, typename U, \
2019 enable_if_t<conjunction<is_real_interoperable<T>, is_cvr_real<U>>::value, int> = 0> \
2020 inline real dispatch_##name(const T &x, U &&a) \
2021 { \
2022 MPPP_MAYBE_TLS real tmp; \
2023 tmp.set_prec(c_max(a.get_prec(), real_deduce_precision(x))); \
2024 tmp.set(x); \
2025 return dispatch_##name(tmp, std::forward<U>(a)); \
2026 } \
2027 } \
2028 /* The overload which returns the result. */ \
2029 MPPP_REAL_MPFR_BINARY_HEADER2 inline real name(T &&y, U &&x) \
2030 { \
2031 return detail::dispatch_##name(std::forward<T>(y), std::forward<U>(x)); \
2032 }
2033
2034 // Neg and abs.
MPPP_REAL_MPFR_UNARY_IMPL(neg,::mpfr_neg,true)2035 MPPP_REAL_MPFR_UNARY_IMPL(neg, ::mpfr_neg, true)
2036 MPPP_REAL_MPFR_UNARY_IMPL(abs, ::mpfr_abs, true)
2037
2038 // Positive difference.
2039 MPPP_REAL_MPFR_BINARY_IMPL(dim, ::mpfr_dim, true)
2040
2041 // Square root.
2042 MPPP_REAL_MPFR_UNARY_IMPL(sqrt, ::mpfr_sqrt, true)
2043
2044 #if defined(MPPP_WITH_ARB)
2045
2046 // sqrt1pm1.
2047 MPPP_REAL_MPFR_UNARY_IMPL(sqrt1pm1, detail::arb_sqrt1pm1, false)
2048
2049 #endif
2050
2051 // Reciprocal square root.
2052 MPPP_REAL_MPFR_UNARY_IMPL(rec_sqrt, ::mpfr_rec_sqrt, true)
2053
2054 // Cubic root.
2055 MPPP_REAL_MPFR_UNARY_IMPL(cbrt, ::mpfr_cbrt, true)
2056
2057 #if defined(MPPP_MPFR_HAVE_MPFR_ROOTN_UI)
2058
2059 // K-th root.
2060 #if defined(MPPP_HAVE_CONCEPTS)
2061 template <cvr_real T>
2062 #else
2063 template <typename T, cvr_real_enabler<T> = 0>
2064 #endif
2065 inline real &rootn_ui(real &rop, T &&op, unsigned long k)
2066 {
2067 auto wrapper = [k](::mpfr_t r, const ::mpfr_t o) { ::mpfr_rootn_ui(r, o, k, MPFR_RNDN); };
2068 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(op));
2069 }
2070
2071 #if defined(MPPP_HAVE_CONCEPTS)
2072 template <cvr_real T>
2073 #else
2074 template <typename T, cvr_real_enabler<T> = 0>
2075 #endif
rootn_ui(T && r,unsigned long k)2076 inline real rootn_ui(T &&r, unsigned long k)
2077 {
2078 auto wrapper = [k](::mpfr_t rop, const ::mpfr_t op) { ::mpfr_rootn_ui(rop, op, k, MPFR_RNDN); };
2079 return detail::mpfr_nary_op_return_impl<false>(0, wrapper, std::forward<T>(r));
2080 }
2081
2082 #endif
2083
2084 // Ternary exponentiation.
2085 #if defined(MPPP_HAVE_CONCEPTS)
2086 template <cvr_real T, cvr_real U>
2087 #else
2088 template <typename T, typename U, cvr_real_enabler<T, U> = 0>
2089 #endif
pow(real & rop,T && op1,U && op2)2090 inline real &pow(real &rop, T &&op1, U &&op2)
2091 {
2092 return detail::mpfr_nary_op_impl<true>(0, ::mpfr_pow, rop, std::forward<T>(op1), std::forward<U>(op2));
2093 }
2094
2095 namespace detail
2096 {
2097
2098 // real-real.
2099 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cvr_real<U>>::value, int> = 0>
dispatch_real_pow(T && op1,U && op2)2100 inline real dispatch_real_pow(T &&op1, U &&op2)
2101 {
2102 return mpfr_nary_op_return_impl<true>(0, ::mpfr_pow, std::forward<T>(op1), std::forward<U>(op2));
2103 }
2104
2105 // real-integer.
2106 template <typename T, std::size_t SSize, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_pow(T && a,const integer<SSize> & n)2107 inline real dispatch_real_pow(T &&a, const integer<SSize> &n)
2108 {
2109 auto wrapper = [&n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_pow_z(r, o, n.get_mpz_view(), MPFR_RNDN); };
2110
2111 // NOTE: in these mpfr_nary_op_return_impl() invocations, we are passing a min_prec
2112 // which is by definition valid because it is produced by an invocation of
2113 // real_deduce_precision() (which does clamping).
2114 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2115 }
2116
2117 // real-unsigned integral.
2118 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_unsigned_integral<U>>::value, int> = 0>
dispatch_real_pow(T && a,const U & n)2119 inline real dispatch_real_pow(T &&a, const U &n)
2120 {
2121 if (n <= nl_max<unsigned long>()) {
2122 auto wrapper
2123 = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_pow_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
2124
2125 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2126 } else {
2127 return dispatch_real_pow(std::forward<T>(a), integer<2>{n});
2128 }
2129 }
2130
2131 // Special casing for bool
2132 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_pow(T && a,const bool & n)2133 inline real dispatch_real_pow(T &&a, const bool &n)
2134 {
2135 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_pow_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
2136
2137 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2138 }
2139
2140 // real-signed integral.
2141 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_signed_integral<U>>::value, int> = 0>
dispatch_real_pow(T && a,const U & n)2142 inline real dispatch_real_pow(T &&a, const U &n)
2143 {
2144 if (n <= nl_max<long>() && n >= nl_min<long>()) {
2145 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_pow_si(r, o, static_cast<long>(n), MPFR_RNDN); };
2146
2147 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2148 } else {
2149 return dispatch_real_pow(std::forward<T>(a), integer<2>{n});
2150 }
2151 }
2152
2153 // real-(floating point, rational, real128).
2154 template <typename T, typename U,
2155 enable_if_t<conjunction<is_cvr_real<T>, disjunction<is_cpp_floating_point<U>, is_rational<U>
2156 #if defined(MPPP_WITH_QUADMATH)
2157 ,
2158 std::is_same<U, real128>
2159 #endif
2160 >>::value,
2161 int> = 0>
dispatch_real_pow(T && a,const U & x)2162 inline real dispatch_real_pow(T &&a, const U &x)
2163 {
2164 MPPP_MAYBE_TLS real tmp;
2165 tmp.set_prec(c_max(a.get_prec(), real_deduce_precision(x)));
2166 tmp.set(x);
2167 return dispatch_real_pow(std::forward<T>(a), tmp);
2168 }
2169
2170 // (everything but unsigned integral)-real.
2171 template <
2172 typename T, typename U,
2173 enable_if_t<conjunction<is_real_interoperable<T>, negation<is_cpp_unsigned_integral<T>>, is_cvr_real<U>>::value,
2174 int> = 0>
dispatch_real_pow(const T & x,U && a)2175 inline real dispatch_real_pow(const T &x, U &&a)
2176 {
2177 MPPP_MAYBE_TLS real tmp;
2178 tmp.set_prec(c_max(a.get_prec(), real_deduce_precision(x)));
2179 tmp.set(x);
2180 return dispatch_real_pow(tmp, std::forward<U>(a));
2181 }
2182
2183 // unsigned integral-real.
2184 template <
2185 typename T, typename U,
2186 enable_if_t<conjunction<is_real_interoperable<T>, is_cpp_unsigned_integral<T>, is_cvr_real<U>>::value, int> = 0>
dispatch_real_pow(const T & n,U && a)2187 inline real dispatch_real_pow(const T &n, U &&a)
2188 {
2189 if (n <= nl_max<unsigned long>()) {
2190 auto wrapper
2191 = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_ui_pow(r, static_cast<unsigned long>(n), o, MPFR_RNDN); };
2192
2193 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<U>(a));
2194 } else {
2195 return dispatch_real_pow(integer<2>{n}, std::forward<U>(a));
2196 }
2197 }
2198
2199 // Special casing for bool.
2200 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_pow(const bool & n,T && a)2201 inline real dispatch_real_pow(const bool &n, T &&a)
2202 {
2203 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_ui_pow(r, static_cast<unsigned long>(n), o, MPFR_RNDN); };
2204
2205 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2206 }
2207
2208 } // namespace detail
2209
2210 // Binary exponentiation.
2211 #if defined(MPPP_HAVE_CONCEPTS)
2212 template <typename T, typename U>
2213 requires real_op_types<T, U>
2214 #else
2215 template <typename T, typename U, detail::enable_if_t<are_real_op_types<T, U>::value, int> = 0>
2216 #endif
pow(T && op1,U && op2)2217 inline real pow(T &&op1, U &&op2)
2218 {
2219 return detail::dispatch_real_pow(std::forward<T>(op1), std::forward<U>(op2));
2220 }
2221
2222 // Squaring.
MPPP_REAL_MPFR_UNARY_IMPL(sqr,::mpfr_sqr,true)2223 MPPP_REAL_MPFR_UNARY_IMPL(sqr, ::mpfr_sqr, true)
2224
2225 // Trigonometric functions.
2226 MPPP_REAL_MPFR_UNARY_IMPL(sin, ::mpfr_sin, true)
2227 MPPP_REAL_MPFR_UNARY_IMPL(cos, ::mpfr_cos, true)
2228 MPPP_REAL_MPFR_UNARY_IMPL(tan, ::mpfr_tan, true)
2229 MPPP_REAL_MPFR_UNARY_IMPL(sec, ::mpfr_sec, true)
2230 MPPP_REAL_MPFR_UNARY_IMPL(csc, ::mpfr_csc, true)
2231 MPPP_REAL_MPFR_UNARY_IMPL(cot, ::mpfr_cot, true)
2232
2233 #if defined(MPPP_WITH_ARB)
2234
2235 MPPP_REAL_MPFR_UNARY_IMPL(sin_pi, detail::arb_sin_pi, false)
2236 MPPP_REAL_MPFR_UNARY_IMPL(cos_pi, detail::arb_cos_pi, false)
2237 MPPP_REAL_MPFR_UNARY_IMPL(tan_pi, detail::arb_tan_pi, false)
2238 MPPP_REAL_MPFR_UNARY_IMPL(cot_pi, detail::arb_cot_pi, false)
2239 MPPP_REAL_MPFR_UNARY_IMPL(sinc, detail::arb_sinc, false)
2240 MPPP_REAL_MPFR_UNARY_IMPL(sinc_pi, detail::arb_sinc_pi, false)
2241
2242 #endif
2243
2244 MPPP_REAL_MPFR_UNARY_IMPL(asin, ::mpfr_asin, true)
2245 MPPP_REAL_MPFR_UNARY_IMPL(acos, ::mpfr_acos, true)
2246 MPPP_REAL_MPFR_UNARY_IMPL(atan, ::mpfr_atan, true)
2247
2248 // sin and cos at the same time.
2249 // NOTE: we don't have the machinery to steal resources
2250 // for multiple retvals, thus we do a manual implementation
2251 // of this function. We keep the signature with cvr_real
2252 // for consistency with the other functions.
2253 #if defined(MPPP_HAVE_CONCEPTS)
2254 template <cvr_real T>
2255 #else
2256 template <typename T, cvr_real_enabler<T> = 0>
2257 #endif
2258 inline void sin_cos(real &sop, real &cop, T &&op)
2259 {
2260 if (mppp_unlikely(&sop == &cop)) {
2261 throw std::invalid_argument(
2262 "In the real sin_cos() function, the return values 'sop' and 'cop' must be distinct objects");
2263 }
2264
2265 // Set the precision of sop and cop to the
2266 // precision of op.
2267 const auto op_prec = op.get_prec();
2268 // NOTE: use prec_round() to avoid issues in case
2269 // sop/cop overlap with op.
2270 sop.prec_round(op_prec);
2271 cop.prec_round(op_prec);
2272
2273 // Run the mpfr function.
2274 ::mpfr_sin_cos(sop._get_mpfr_t(), cop._get_mpfr_t(), op.get_mpfr_t(), MPFR_RNDN);
2275 }
2276
MPPP_REAL_MPFR_BINARY_IMPL(atan2,::mpfr_atan2,true)2277 MPPP_REAL_MPFR_BINARY_IMPL(atan2, ::mpfr_atan2, true)
2278
2279 // Hyperbolic functions.
2280 MPPP_REAL_MPFR_UNARY_IMPL(sinh, ::mpfr_sinh, true)
2281 MPPP_REAL_MPFR_UNARY_IMPL(cosh, ::mpfr_cosh, true)
2282 MPPP_REAL_MPFR_UNARY_IMPL(tanh, ::mpfr_tanh, true)
2283 MPPP_REAL_MPFR_UNARY_IMPL(sech, ::mpfr_sech, true)
2284 MPPP_REAL_MPFR_UNARY_IMPL(csch, ::mpfr_csch, true)
2285 MPPP_REAL_MPFR_UNARY_IMPL(coth, ::mpfr_coth, true)
2286 MPPP_REAL_MPFR_UNARY_IMPL(asinh, ::mpfr_asinh, true)
2287 MPPP_REAL_MPFR_UNARY_IMPL(acosh, ::mpfr_acosh, true)
2288 MPPP_REAL_MPFR_UNARY_IMPL(atanh, ::mpfr_atanh, true)
2289
2290 // sinh and cosh at the same time.
2291 // NOTE: we don't have the machinery to steal resources
2292 // for multiple retvals, thus we do a manual implementation
2293 // of this function. We keep the signature with cvr_real
2294 // for consistency with the other functions.
2295 #if defined(MPPP_HAVE_CONCEPTS)
2296 template <cvr_real T>
2297 #else
2298 template <typename T, cvr_real_enabler<T> = 0>
2299 #endif
2300 inline void sinh_cosh(real &sop, real &cop, T &&op)
2301 {
2302 if (mppp_unlikely(&sop == &cop)) {
2303 throw std::invalid_argument(
2304 "In the real sinh_cosh() function, the return values 'sop' and 'cop' must be distinct objects");
2305 }
2306
2307 // Set the precision of sop and cop to the
2308 // precision of op.
2309 const auto op_prec = op.get_prec();
2310 // NOTE: use prec_round() to avoid issues in case
2311 // sop/cop overlap with op.
2312 sop.prec_round(op_prec);
2313 cop.prec_round(op_prec);
2314
2315 // Run the mpfr function.
2316 ::mpfr_sinh_cosh(sop._get_mpfr_t(), cop._get_mpfr_t(), op.get_mpfr_t(), MPFR_RNDN);
2317 }
2318
2319 // Exponentials and logarithms.
MPPP_REAL_MPFR_UNARY_IMPL(exp,::mpfr_exp,true)2320 MPPP_REAL_MPFR_UNARY_IMPL(exp, ::mpfr_exp, true)
2321 MPPP_REAL_MPFR_UNARY_IMPL(exp2, ::mpfr_exp2, true)
2322 MPPP_REAL_MPFR_UNARY_IMPL(exp10, ::mpfr_exp10, true)
2323 MPPP_REAL_MPFR_UNARY_IMPL(expm1, ::mpfr_expm1, true)
2324 MPPP_REAL_MPFR_UNARY_IMPL(log, ::mpfr_log, true)
2325 MPPP_REAL_MPFR_UNARY_IMPL(log2, ::mpfr_log2, true)
2326 MPPP_REAL_MPFR_UNARY_IMPL(log10, ::mpfr_log10, true)
2327 MPPP_REAL_MPFR_UNARY_IMPL(log1p, ::mpfr_log1p, true)
2328
2329 #if defined(MPPP_WITH_ARB)
2330
2331 // log_hypot.
2332 MPPP_REAL_MPFR_BINARY_IMPL(log_hypot, detail::arb_log_hypot, false)
2333
2334 // log_base_ui.
2335 #if defined(MPPP_HAVE_CONCEPTS)
2336 template <cvr_real T>
2337 #else
2338 template <typename T, cvr_real_enabler<T> = 0>
2339 #endif
2340 inline real &log_base_ui(real &rop, T &&op, unsigned long b)
2341 {
2342 auto wrapper = [b](::mpfr_t r, const ::mpfr_t o) { detail::arb_log_base_ui(r, o, b); };
2343 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(op));
2344 }
2345
2346 #if defined(MPPP_HAVE_CONCEPTS)
2347 template <cvr_real T>
2348 #else
2349 template <typename T, cvr_real_enabler<T> = 0>
2350 #endif
log_base_ui(T && r,unsigned long b)2351 inline real log_base_ui(T &&r, unsigned long b)
2352 {
2353 auto wrapper = [b](::mpfr_t rop, const ::mpfr_t op) { detail::arb_log_base_ui(rop, op, b); };
2354 return detail::mpfr_nary_op_return_impl<false>(0, wrapper, std::forward<T>(r));
2355 }
2356
2357 #endif
2358
2359 // Gamma functions.
MPPP_REAL_MPFR_UNARY_IMPL(gamma,::mpfr_gamma,true)2360 MPPP_REAL_MPFR_UNARY_IMPL(gamma, ::mpfr_gamma, true)
2361 MPPP_REAL_MPFR_UNARY_IMPL(lngamma, ::mpfr_lngamma, true)
2362 MPPP_REAL_MPFR_UNARY_IMPL(lgamma, detail::real_lgamma_wrapper, true)
2363 MPPP_REAL_MPFR_UNARY_IMPL(digamma, ::mpfr_digamma, true)
2364
2365 #if defined(MPPP_MPFR_HAVE_MPFR_GAMMA_INC)
2366
2367 // gamma_inc.
2368 MPPP_REAL_MPFR_BINARY_IMPL(gamma_inc, ::mpfr_gamma_inc, true)
2369
2370 #endif
2371
2372 // Bessel functions.
2373 MPPP_REAL_MPFR_UNARY_IMPL(j0, ::mpfr_j0, true)
2374 MPPP_REAL_MPFR_UNARY_IMPL(j1, ::mpfr_j1, true)
2375
2376 // Bessel function of the first kind of order n.
2377 #if defined(MPPP_HAVE_CONCEPTS)
2378 template <cvr_real T>
2379 #else
2380 template <typename T, cvr_real_enabler<T> = 0>
2381 #endif
2382 inline real &jn(real &rop, long n, T &&op)
2383 {
2384 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_jn(r, n, o, MPFR_RNDN); };
2385 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(op));
2386 }
2387
2388 #if defined(MPPP_HAVE_CONCEPTS)
2389 template <cvr_real T>
2390 #else
2391 template <typename T, cvr_real_enabler<T> = 0>
2392 #endif
jn(long n,T && r)2393 inline real jn(long n, T &&r)
2394 {
2395 auto wrapper = [n](::mpfr_t rop, const ::mpfr_t op) { ::mpfr_jn(rop, n, op, MPFR_RNDN); };
2396 return detail::mpfr_nary_op_return_impl<false>(0, wrapper, std::forward<T>(r));
2397 }
2398
MPPP_REAL_MPFR_UNARY_IMPL(y0,::mpfr_y0,true)2399 MPPP_REAL_MPFR_UNARY_IMPL(y0, ::mpfr_y0, true)
2400 MPPP_REAL_MPFR_UNARY_IMPL(y1, ::mpfr_y1, true)
2401
2402 // Bessel function of the second kind of order n.
2403 #if defined(MPPP_HAVE_CONCEPTS)
2404 template <cvr_real T>
2405 #else
2406 template <typename T, cvr_real_enabler<T> = 0>
2407 #endif
2408 inline real &yn(real &rop, long n, T &&op)
2409 {
2410 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_yn(r, n, o, MPFR_RNDN); };
2411 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(op));
2412 }
2413
2414 #if defined(MPPP_HAVE_CONCEPTS)
2415 template <cvr_real T>
2416 #else
2417 template <typename T, cvr_real_enabler<T> = 0>
2418 #endif
yn(long n,T && r)2419 inline real yn(long n, T &&r)
2420 {
2421 auto wrapper = [n](::mpfr_t rop, const ::mpfr_t op) { ::mpfr_yn(rop, n, op, MPFR_RNDN); };
2422 return detail::mpfr_nary_op_return_impl<false>(0, wrapper, std::forward<T>(r));
2423 }
2424
2425 #if defined(MPPP_WITH_ARB)
2426
MPPP_REAL_MPFR_BINARY_IMPL(jx,detail::arb_hypgeom_bessel_j,false)2427 MPPP_REAL_MPFR_BINARY_IMPL(jx, detail::arb_hypgeom_bessel_j, false)
2428 MPPP_REAL_MPFR_BINARY_IMPL(yx, detail::arb_hypgeom_bessel_y, false)
2429
2430 #endif
2431
2432 // Polylogarithms.
2433 MPPP_REAL_MPFR_UNARY_IMPL(li2, detail::real_li2_wrapper, true)
2434
2435 #if defined(MPPP_WITH_ARB)
2436
2437 // Polylogarithm, integer order.
2438 #if defined(MPPP_HAVE_CONCEPTS)
2439 template <cvr_real T>
2440 #else
2441 template <typename T, cvr_real_enabler<T> = 0>
2442 #endif
2443 inline real &polylog_si(real &rop, long n, T &&op)
2444 {
2445 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { detail::arb_polylog_si(r, n, o); };
2446 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(op));
2447 }
2448
2449 #if defined(MPPP_HAVE_CONCEPTS)
2450 template <cvr_real T>
2451 #else
2452 template <typename T, cvr_real_enabler<T> = 0>
2453 #endif
polylog_si(long n,T && r)2454 inline real polylog_si(long n, T &&r)
2455 {
2456 auto wrapper = [n](::mpfr_t rop, const ::mpfr_t op) { detail::arb_polylog_si(rop, n, op); };
2457 return detail::mpfr_nary_op_return_impl<false>(0, wrapper, std::forward<T>(r));
2458 }
2459
2460 // Polylogarithm, real order.
MPPP_REAL_MPFR_BINARY_IMPL(polylog,detail::arb_polylog,false)2461 MPPP_REAL_MPFR_BINARY_IMPL(polylog, detail::arb_polylog, false)
2462
2463 #endif
2464
2465 // Other special functions.
2466 MPPP_REAL_MPFR_UNARY_IMPL(eint, ::mpfr_eint, true)
2467 MPPP_REAL_MPFR_UNARY_IMPL(zeta, ::mpfr_zeta, true)
2468 MPPP_REAL_MPFR_UNARY_IMPL(erf, ::mpfr_erf, true)
2469 MPPP_REAL_MPFR_UNARY_IMPL(erfc, ::mpfr_erfc, true)
2470 MPPP_REAL_MPFR_UNARY_IMPL(ai, ::mpfr_ai, true)
2471
2472 #if defined(MPPP_WITH_ARB)
2473
2474 MPPP_REAL_MPFR_UNARY_IMPL(lambert_w0, detail::arb_lambert_w0, false)
2475 MPPP_REAL_MPFR_UNARY_IMPL(lambert_wm1, detail::arb_lambert_wm1, false)
2476
2477 #endif
2478
2479 #if defined(MPPP_MPFR_HAVE_MPFR_BETA)
2480
2481 // beta.
2482 MPPP_REAL_MPFR_BINARY_IMPL(beta, ::mpfr_beta, true)
2483
2484 #endif
2485
2486 // hypot.
2487 MPPP_REAL_MPFR_BINARY_IMPL(hypot, ::mpfr_hypot, true)
2488
2489 // agm.
2490 MPPP_REAL_MPFR_BINARY_IMPL(agm, ::mpfr_agm, true)
2491
2492 // Integer and remainder-related functions.
2493 MPPP_REAL_MPFR_UNARY_IMPL(ceil, detail::real_ceil_wrapper, false)
2494 MPPP_REAL_MPFR_UNARY_IMPL(floor, detail::real_floor_wrapper, false)
2495 MPPP_REAL_MPFR_UNARY_IMPL(round, detail::real_round_wrapper, false)
2496 #if defined(MPPP_MPFR_HAVE_MPFR_ROUNDEVEN)
2497 MPPP_REAL_MPFR_UNARY_IMPL(roundeven, detail::real_roundeven_wrapper, false)
2498 #endif
2499 MPPP_REAL_MPFR_UNARY_IMPL(trunc, detail::real_trunc_wrapper, false)
2500 MPPP_REAL_MPFR_UNARY_IMPL(frac, detail::real_frac_wrapper, false)
2501
2502 // modf.
2503 // NOTE: we don't have the machinery to steal resources
2504 // for multiple retvals, thus we do a manual implementation
2505 // of this function. We keep the signature with cvr_real
2506 // for consistency with the other functions.
2507 #if defined(MPPP_HAVE_CONCEPTS)
2508 template <cvr_real T>
2509 #else
2510 template <typename T, cvr_real_enabler<T> = 0>
2511 #endif
2512 inline void modf(real &iop, real &fop, T &&op)
2513 {
2514 if (mppp_unlikely(&iop == &fop)) {
2515 throw std::invalid_argument(
2516 "In the real modf() function, the return values 'iop' and 'fop' must be distinct objects");
2517 }
2518 if (mppp_unlikely(op.nan_p())) {
2519 throw std::domain_error("In the real modf() function, the input argument cannot be NaN");
2520 }
2521
2522 // Set the precision of iop and fop to the
2523 // precision of op.
2524 const auto op_prec = op.get_prec();
2525 // NOTE: use prec_round() to avoid issues in case
2526 // iop/fop overlap with op.
2527 iop.prec_round(op_prec);
2528 fop.prec_round(op_prec);
2529
2530 // Run the mpfr function.
2531 ::mpfr_modf(iop._get_mpfr_t(), fop._get_mpfr_t(), op.get_mpfr_t(), MPFR_RNDN);
2532 }
2533
MPPP_REAL_MPFR_BINARY_IMPL(fmod,::mpfr_fmod,true)2534 MPPP_REAL_MPFR_BINARY_IMPL(fmod, ::mpfr_fmod, true)
2535 MPPP_REAL_MPFR_BINARY_IMPL(remainder, ::mpfr_remainder, true)
2536
2537 #if defined(MPPP_HAVE_CONCEPTS)
2538 template <cvr_real T, cvr_real U>
2539 #else
2540 template <typename T, typename U, cvr_real_enabler<T, U> = 0>
2541 #endif
2542 inline real &remquo(real &rop, long *q, T &&x, U &&y)
2543 {
2544 auto wrapper = [q](::mpfr_t r, const ::mpfr_t o1, const ::mpfr_t o2) { ::mpfr_remquo(r, q, o1, o2, MPFR_RNDN); };
2545 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(x), std::forward<U>(y));
2546 }
2547
2548 #if defined(MPPP_MPFR_HAVE_MPFR_FMODQUO)
2549
2550 #if defined(MPPP_HAVE_CONCEPTS)
2551 template <cvr_real T, cvr_real U>
2552 #else
2553 template <typename T, typename U, cvr_real_enabler<T, U> = 0>
2554 #endif
fmodquo(real & rop,long * q,T && x,U && y)2555 inline real &fmodquo(real &rop, long *q, T &&x, U &&y)
2556 {
2557 auto wrapper = [q](::mpfr_t r, const ::mpfr_t o1, const ::mpfr_t o2) { ::mpfr_fmodquo(r, q, o1, o2, MPFR_RNDN); };
2558 return detail::mpfr_nary_op_impl<false>(0, wrapper, rop, std::forward<T>(x), std::forward<U>(y));
2559 }
2560
2561 #endif
2562
2563 #undef MPPP_REAL_MPFR_UNARY_HEADER
2564 #undef MPPP_REAL_MPFR_UNARY_IMPL
2565 #undef MPPP_REAL_MPFR_BINARY_HEADER1
2566 #undef MPPP_REAL_MPFR_BINARY_HEADER2
2567 #undef MPPP_REAL_MPFR_BINARY_IMPL
2568
2569 // Output stream operator.
2570 MPPP_DLL_PUBLIC std::ostream &operator<<(std::ostream &, const real &);
2571
2572 #if defined(MPPP_MPFR_HAVE_MPFR_GET_STR_NDIGITS)
2573
2574 // Get the number of significant digits required for a round-tripping representation.
2575 MPPP_DLL_PUBLIC std::size_t get_str_ndigits(const real &, int = 10);
2576
2577 #endif
2578
2579 // Binary serialization.
2580 MPPP_DLL_PUBLIC std::size_t binary_size(const real &);
2581
2582 // Save in binary format.
2583 template <typename T>
binary_save(const real & x,T && dest)2584 inline auto binary_save(const real &x, T &&dest) -> decltype(x.binary_save(std::forward<T>(dest)))
2585 {
2586 return x.binary_save(std::forward<T>(dest));
2587 }
2588
2589 // Load in binary format.
2590 template <typename T>
binary_load(real & x,T && src)2591 inline auto binary_load(real &x, T &&src) -> decltype(x.binary_load(std::forward<T>(src)))
2592 {
2593 return x.binary_load(std::forward<T>(src));
2594 }
2595
2596 // Constants.
2597 MPPP_DLL_PUBLIC real real_pi(::mpfr_prec_t);
2598 MPPP_DLL_PUBLIC real &real_pi(real &);
2599 MPPP_DLL_PUBLIC real real_log2(::mpfr_prec_t);
2600 MPPP_DLL_PUBLIC real &real_log2(real &);
2601 MPPP_DLL_PUBLIC real real_euler(::mpfr_prec_t);
2602 MPPP_DLL_PUBLIC real &real_euler(real &);
2603 MPPP_DLL_PUBLIC real real_catalan(::mpfr_prec_t);
2604 MPPP_DLL_PUBLIC real &real_catalan(real &);
2605
2606 // Identity operator.
2607 #if defined(MPPP_HAVE_CONCEPTS)
2608 template <cvr_real T>
2609 #else
2610 template <typename T, cvr_real_enabler<T> = 0>
2611 #endif
operator +(T && r)2612 inline real operator+(T &&r)
2613 {
2614 return std::forward<T>(r);
2615 }
2616
2617 namespace detail
2618 {
2619
2620 // real-real.
2621 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cvr_real<U>>::value, int> = 0>
dispatch_real_binary_add(T && a,U && b)2622 inline real dispatch_real_binary_add(T &&a, U &&b)
2623 {
2624 return mpfr_nary_op_return_impl<true>(0, ::mpfr_add, std::forward<T>(a), std::forward<U>(b));
2625 }
2626
2627 // real-integer.
2628 template <typename T, std::size_t SSize>
dispatch_real_binary_add(T && a,const integer<SSize> & n)2629 inline real dispatch_real_binary_add(T &&a, const integer<SSize> &n)
2630 {
2631 auto wrapper = [&n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_add_z(r, o, n.get_mpz_view(), MPFR_RNDN); };
2632
2633 // NOTE: in these mpfr_nary_op_return_impl() invocations, we are passing a min_prec
2634 // which is by definition valid because it is produced by an invocation of
2635 // real_deduce_precision() (which does clamping).
2636 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2637 }
2638
2639 // integer-real.
2640 template <typename T, std::size_t SSize>
dispatch_real_binary_add(const integer<SSize> & n,T && a)2641 inline real dispatch_real_binary_add(const integer<SSize> &n, T &&a)
2642 {
2643 return dispatch_real_binary_add(std::forward<T>(a), n);
2644 }
2645
2646 // real-unsigned integral.
2647 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_unsigned_integral<U>>::value, int> = 0>
dispatch_real_binary_add(T && a,const U & n)2648 inline real dispatch_real_binary_add(T &&a, const U &n)
2649 {
2650 if (n <= nl_max<unsigned long>()) {
2651 auto wrapper
2652 = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_add_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
2653
2654 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2655 } else {
2656 return dispatch_real_binary_add(std::forward<T>(a), integer<2>{n});
2657 }
2658 }
2659
2660 // real-bool.
2661 // NOTE: make this explicit (rather than letting bool fold into
2662 // the unsigned integrals overload) in order to avoid MSVC warnings.
2663 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_binary_add(T && a,const bool & n)2664 inline real dispatch_real_binary_add(T &&a, const bool &n)
2665 {
2666 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_add_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
2667 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2668 }
2669
2670 // unsigned integral-real.
2671 template <typename T, typename U, enable_if_t<conjunction<is_cpp_unsigned_integral<T>, is_cvr_real<U>>::value, int> = 0>
dispatch_real_binary_add(const T & n,U && a)2672 inline real dispatch_real_binary_add(const T &n, U &&a)
2673 {
2674 return dispatch_real_binary_add(std::forward<U>(a), n);
2675 }
2676
2677 // real-signed integral.
2678 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_signed_integral<U>>::value, int> = 0>
dispatch_real_binary_add(T && a,const U & n)2679 inline real dispatch_real_binary_add(T &&a, const U &n)
2680 {
2681 if (n <= nl_max<long>() && n >= nl_min<long>()) {
2682 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_add_si(r, o, static_cast<long>(n), MPFR_RNDN); };
2683
2684 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2685 } else {
2686 return dispatch_real_binary_add(std::forward<T>(a), integer<2>{n});
2687 }
2688 }
2689
2690 // signed integral-real.
2691 template <typename T, typename U, enable_if_t<conjunction<is_cpp_signed_integral<T>, is_cvr_real<U>>::value, int> = 0>
dispatch_real_binary_add(const T & n,U && a)2692 inline real dispatch_real_binary_add(const T &n, U &&a)
2693 {
2694 return dispatch_real_binary_add(std::forward<U>(a), n);
2695 }
2696
2697 // real-rational.
2698 template <typename T, std::size_t SSize>
dispatch_real_binary_add(T && a,const rational<SSize> & q)2699 inline real dispatch_real_binary_add(T &&a, const rational<SSize> &q)
2700 {
2701 const auto qv = detail::get_mpq_view(q);
2702
2703 auto wrapper = [&qv](::mpfr_t r, const ::mpfr_t o) { ::mpfr_add_q(r, o, &qv, MPFR_RNDN); };
2704
2705 return mpfr_nary_op_return_impl<false>(real_deduce_precision(q), wrapper, std::forward<T>(a));
2706 }
2707
2708 // rational-real.
2709 template <typename T, std::size_t SSize>
dispatch_real_binary_add(const rational<SSize> & q,T && a)2710 inline real dispatch_real_binary_add(const rational<SSize> &q, T &&a)
2711 {
2712 return dispatch_real_binary_add(std::forward<T>(a), q);
2713 }
2714
2715 // real-(float, double).
2716 template <typename T, typename U,
2717 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, float>, std::is_same<U, double>>>::value,
2718 int> = 0>
dispatch_real_binary_add(T && a,const U & x)2719 inline real dispatch_real_binary_add(T &&a, const U &x)
2720 {
2721 // NOTE: the MPFR docs state that mpfr_add_d() assumes that
2722 // the radix of double is a power of 2. If we ever run into platforms
2723 // for which this is not true, we can add a compile-time dispatch
2724 // that uses the long double implementation instead.
2725 constexpr auto dradix = static_cast<unsigned>(std::numeric_limits<double>::radix);
2726 static_assert(!(dradix & (dradix - 1)), "mpfr_add_d() requires the radix of the 'double' type to be a power of 2.");
2727
2728 auto wrapper = [x](::mpfr_t r, const ::mpfr_t o) { ::mpfr_add_d(r, o, static_cast<double>(x), MPFR_RNDN); };
2729
2730 return mpfr_nary_op_return_impl<false>(real_deduce_precision(x), wrapper, std::forward<T>(a));
2731 }
2732
2733 // (float, double)-real.
2734 template <typename T, typename U,
2735 enable_if_t<conjunction<is_cvr_real<U>, disjunction<std::is_same<T, float>, std::is_same<T, double>>>::value,
2736 int> = 0>
dispatch_real_binary_add(const T & x,U && a)2737 inline real dispatch_real_binary_add(const T &x, U &&a)
2738 {
2739 return dispatch_real_binary_add(std::forward<U>(a), x);
2740 }
2741
2742 // real-(long double, real128).
2743 template <typename T, typename U,
2744 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, long double>
2745 #if defined(MPPP_WITH_QUADMATH)
2746 ,
2747 std::is_same<U, real128>
2748 #endif
2749 >>::value,
2750 int> = 0>
dispatch_real_binary_add(T && a,const U & x)2751 inline real dispatch_real_binary_add(T &&a, const U &x)
2752 {
2753 MPPP_MAYBE_TLS real tmp;
2754 tmp.set_prec(c_max(a.get_prec(), real_deduce_precision(x)));
2755 tmp.set(x);
2756 return dispatch_real_binary_add(std::forward<T>(a), tmp);
2757 }
2758
2759 // (long double, real128)-real.
2760 template <typename T, typename U,
2761 enable_if_t<conjunction<disjunction<std::is_same<T, long double>
2762 #if defined(MPPP_WITH_QUADMATH)
2763 ,
2764 std::is_same<T, real128>
2765 #endif
2766 >,
2767 is_cvr_real<U>>::value,
2768 int> = 0>
dispatch_real_binary_add(const T & x,U && a)2769 inline real dispatch_real_binary_add(const T &x, U &&a)
2770 {
2771 return dispatch_real_binary_add(std::forward<U>(a), x);
2772 }
2773
2774 } // namespace detail
2775
2776 // Binary addition.
2777 #if defined(MPPP_HAVE_CONCEPTS)
2778 template <typename T, typename U>
2779 requires real_op_types<T, U>
2780 #else
2781 template <typename T, typename U, detail::enable_if_t<are_real_op_types<T, U>::value, int> = 0>
2782 #endif
operator +(T && a,U && b)2783 inline real operator+(T &&a, U &&b)
2784 {
2785 return detail::dispatch_real_binary_add(std::forward<T>(a), std::forward<U>(b));
2786 }
2787
2788 namespace detail
2789 {
2790
2791 // real-real.
2792 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_in_place_add(real & a,T && b)2793 inline void dispatch_real_in_place_add(real &a, T &&b)
2794 {
2795 add(a, a, std::forward<T>(b));
2796 }
2797
2798 MPPP_DLL_PUBLIC void dispatch_real_in_place_add_integer_impl(real &, const ::mpz_t, ::mpfr_prec_t);
2799
2800 // real-integer.
2801 template <std::size_t SSize>
dispatch_real_in_place_add(real & a,const integer<SSize> & n)2802 inline void dispatch_real_in_place_add(real &a, const integer<SSize> &n)
2803 {
2804 dispatch_real_in_place_add_integer_impl(a, n.get_mpz_view(), real_deduce_precision(n));
2805 }
2806
2807 // real-unsigned C++ integral.
2808 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_in_place_add(real & a,const T & n)2809 inline void dispatch_real_in_place_add(real &a, const T &n)
2810 {
2811 if (n <= nl_max<unsigned long>()) {
2812 auto wrapper
2813 = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_add_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
2814
2815 mpfr_nary_op_impl<false>(real_deduce_precision(n), wrapper, a, a);
2816 } else {
2817 dispatch_real_in_place_add(a, integer<2>{n});
2818 }
2819 }
2820
2821 // real-bool.
2822 // NOTE: make this explicit (rather than letting bool fold into
2823 // the unsigned integrals overload) in order to avoid MSVC warnings.
2824 MPPP_DLL_PUBLIC void dispatch_real_in_place_add(real &, bool);
2825
2826 // real-signed C++ integral.
2827 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_in_place_add(real & a,const T & n)2828 inline void dispatch_real_in_place_add(real &a, const T &n)
2829 {
2830 if (n <= nl_max<long>() && n >= nl_min<long>()) {
2831 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_add_si(r, o, static_cast<long>(n), MPFR_RNDN); };
2832
2833 mpfr_nary_op_impl<false>(real_deduce_precision(n), wrapper, a, a);
2834 } else {
2835 dispatch_real_in_place_add(a, integer<2>{n});
2836 }
2837 }
2838
2839 MPPP_DLL_PUBLIC void dispatch_real_in_place_add_rational_impl(real &, const ::mpq_t, ::mpfr_prec_t);
2840
2841 // real-rational.
2842 template <std::size_t SSize>
dispatch_real_in_place_add(real & a,const rational<SSize> & q)2843 inline void dispatch_real_in_place_add(real &a, const rational<SSize> &q)
2844 {
2845 const auto qv = detail::get_mpq_view(q);
2846 dispatch_real_in_place_add_rational_impl(a, &qv, real_deduce_precision(q));
2847 }
2848
2849 // real-(float, double).
2850 MPPP_DLL_PUBLIC void dispatch_real_in_place_add(real &, const float &);
2851 MPPP_DLL_PUBLIC void dispatch_real_in_place_add(real &, const double &);
2852
2853 // real-(long double, real128).
2854 MPPP_DLL_PUBLIC void dispatch_real_in_place_add(real &, const long double &);
2855 #if defined(MPPP_WITH_QUADMATH)
2856 MPPP_DLL_PUBLIC void dispatch_real_in_place_add(real &, const real128 &);
2857 #endif
2858
2859 // (real interoperable)-real.
2860 template <typename T, typename U, enable_if_t<is_real_interoperable<T>::value, int> = 0>
dispatch_real_in_place_add(T & x,U && a)2861 inline void dispatch_real_in_place_add(T &x, U &&a)
2862 {
2863 x = static_cast<T>(x + std::forward<U>(a));
2864 }
2865
2866 } // namespace detail
2867
2868 // In-place addition.
2869 #if defined(MPPP_HAVE_CONCEPTS)
2870 template <typename T, typename U>
2871 requires real_in_place_op_types<T, U>
2872 #else
2873 template <typename T, typename U, detail::enable_if_t<are_real_in_place_op_types<T, U>::value, int> = 0>
2874 #endif
operator +=(T & a,U && b)2875 inline T &operator+=(T &a, U &&b)
2876 {
2877 detail::dispatch_real_in_place_add(a, std::forward<U>(b));
2878 return a;
2879 }
2880
2881 // Prefix increment.
2882 MPPP_DLL_PUBLIC real &operator++(real &);
2883
2884 // Suffix increment.
2885 MPPP_DLL_PUBLIC real operator++(real &, int);
2886
2887 // Negated copy.
2888 #if defined(MPPP_HAVE_CONCEPTS)
2889 template <cvr_real T>
2890 #else
2891 template <typename T, cvr_real_enabler<T> = 0>
2892 #endif
operator -(T && r)2893 inline real operator-(T &&r)
2894 {
2895 real retval{std::forward<T>(r)};
2896 retval.neg();
2897 return retval;
2898 }
2899
2900 namespace detail
2901 {
2902
2903 // real-real.
2904 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cvr_real<U>>::value, int> = 0>
dispatch_real_binary_sub(T && a,U && b)2905 inline real dispatch_real_binary_sub(T &&a, U &&b)
2906 {
2907 return mpfr_nary_op_return_impl<true>(0, ::mpfr_sub, std::forward<T>(a), std::forward<U>(b));
2908 }
2909
2910 // real-integer.
2911 template <typename T, std::size_t SSize>
dispatch_real_binary_sub(T && a,const integer<SSize> & n)2912 inline real dispatch_real_binary_sub(T &&a, const integer<SSize> &n)
2913 {
2914 auto wrapper = [&n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_sub_z(r, o, n.get_mpz_view(), MPFR_RNDN); };
2915
2916 // NOTE: in these mpfr_nary_op_return_impl() invocations, we are passing a min_prec
2917 // which is by definition valid because it is produced by an invocation of
2918 // real_deduce_precision() (which does clamping).
2919 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2920 }
2921
2922 // integer-real.
2923 template <typename T, std::size_t SSize>
dispatch_real_binary_sub(const integer<SSize> & n,T && a)2924 inline real dispatch_real_binary_sub(const integer<SSize> &n, T &&a)
2925 {
2926 auto wrapper = [&n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_z_sub(r, n.get_mpz_view(), o, MPFR_RNDN); };
2927
2928 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2929 }
2930
2931 // real-unsigned integral.
2932 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_unsigned_integral<U>>::value, int> = 0>
dispatch_real_binary_sub(T && a,const U & n)2933 inline real dispatch_real_binary_sub(T &&a, const U &n)
2934 {
2935 if (n <= nl_max<unsigned long>()) {
2936 auto wrapper
2937 = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_sub_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
2938
2939 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2940 } else {
2941 return dispatch_real_binary_sub(std::forward<T>(a), integer<2>{n});
2942 }
2943 }
2944
2945 // unsigned integral-real.
2946 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_unsigned_integral<U>>::value, int> = 0>
dispatch_real_binary_sub(const U & n,T && a)2947 inline real dispatch_real_binary_sub(const U &n, T &&a)
2948 {
2949 if (n <= nl_max<unsigned long>()) {
2950 auto wrapper
2951 = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_ui_sub(r, static_cast<unsigned long>(n), o, MPFR_RNDN); };
2952
2953 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2954 } else {
2955 return dispatch_real_binary_sub(integer<2>{n}, std::forward<T>(a));
2956 }
2957 }
2958
2959 // real-bool.
2960 // NOTE: make this explicit (rather than letting bool fold into
2961 // the unsigned integrals overload) in order to avoid MSVC warnings.
2962 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_binary_sub(T && a,const bool & n)2963 inline real dispatch_real_binary_sub(T &&a, const bool &n)
2964 {
2965 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_sub_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
2966 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2967 }
2968
2969 // bool-real.
2970 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_binary_sub(const bool & n,T && a)2971 inline real dispatch_real_binary_sub(const bool &n, T &&a)
2972 {
2973 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_ui_sub(r, static_cast<unsigned long>(n), o, MPFR_RNDN); };
2974 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2975 }
2976
2977 // real-signed integral.
2978 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_signed_integral<U>>::value, int> = 0>
dispatch_real_binary_sub(T && a,const U & n)2979 inline real dispatch_real_binary_sub(T &&a, const U &n)
2980 {
2981 if (n <= nl_max<long>() && n >= nl_min<long>()) {
2982 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_sub_si(r, o, static_cast<long>(n), MPFR_RNDN); };
2983
2984 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2985 } else {
2986 return dispatch_real_binary_sub(std::forward<T>(a), integer<2>{n});
2987 }
2988 }
2989
2990 // signed integral-real.
2991 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_signed_integral<U>>::value, int> = 0>
dispatch_real_binary_sub(const U & n,T && a)2992 inline real dispatch_real_binary_sub(const U &n, T &&a)
2993 {
2994 if (n <= nl_max<long>() && n >= nl_min<long>()) {
2995 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_si_sub(r, static_cast<long>(n), o, MPFR_RNDN); };
2996
2997 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
2998 } else {
2999 return dispatch_real_binary_sub(integer<2>{n}, std::forward<T>(a));
3000 }
3001 }
3002
3003 // real-rational.
3004 template <typename T, std::size_t SSize>
dispatch_real_binary_sub(T && a,const rational<SSize> & q)3005 inline real dispatch_real_binary_sub(T &&a, const rational<SSize> &q)
3006 {
3007 const auto qv = detail::get_mpq_view(q);
3008
3009 auto wrapper = [&qv](::mpfr_t r, const ::mpfr_t o) { ::mpfr_sub_q(r, o, &qv, MPFR_RNDN); };
3010
3011 return mpfr_nary_op_return_impl<false>(real_deduce_precision(q), wrapper, std::forward<T>(a));
3012 }
3013
3014 // rational-real.
3015 template <typename T, std::size_t SSize>
dispatch_real_binary_sub(const rational<SSize> & q,T && a)3016 inline real dispatch_real_binary_sub(const rational<SSize> &q, T &&a)
3017 {
3018 // NOTE: apparently there's no mpfr_q_sub() primitive.
3019 return -dispatch_real_binary_sub(std::forward<T>(a), q);
3020 }
3021
3022 // real-(float, double).
3023 template <typename T, typename U,
3024 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, float>, std::is_same<U, double>>>::value,
3025 int> = 0>
dispatch_real_binary_sub(T && a,const U & x)3026 inline real dispatch_real_binary_sub(T &&a, const U &x)
3027 {
3028 // NOTE: the MPFR docs state that mpfr_sub_d() assumes that
3029 // the radix of double is a power of 2. If we ever run into platforms
3030 // for which this is not true, we can add a compile-time dispatch
3031 // that uses the long double implementation instead.
3032 constexpr auto dradix = static_cast<unsigned>(std::numeric_limits<double>::radix);
3033 static_assert(!(dradix & (dradix - 1)), "mpfr_sub_d() requires the radix of the 'double' type to be a power of 2.");
3034
3035 auto wrapper = [x](::mpfr_t r, const ::mpfr_t o) { ::mpfr_sub_d(r, o, static_cast<double>(x), MPFR_RNDN); };
3036
3037 return mpfr_nary_op_return_impl<false>(real_deduce_precision(x), wrapper, std::forward<T>(a));
3038 }
3039
3040 // (float, double)-real.
3041 template <typename T, typename U,
3042 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, float>, std::is_same<U, double>>>::value,
3043 int> = 0>
dispatch_real_binary_sub(const U & x,T && a)3044 inline real dispatch_real_binary_sub(const U &x, T &&a)
3045 {
3046 // NOTE: the MPFR docs state that mpfr_d_sub() assumes that
3047 // the radix of double is a power of 2. If we ever run into platforms
3048 // for which this is not true, we can add a compile-time dispatch
3049 // that uses the long double implementation instead.
3050 constexpr auto dradix = static_cast<unsigned>(std::numeric_limits<double>::radix);
3051 static_assert(!(dradix & (dradix - 1)), "mpfr_d_sub() requires the radix of the 'double' type to be a power of 2.");
3052
3053 auto wrapper = [x](::mpfr_t r, const ::mpfr_t o) { ::mpfr_d_sub(r, static_cast<double>(x), o, MPFR_RNDN); };
3054
3055 return mpfr_nary_op_return_impl<false>(real_deduce_precision(x), wrapper, std::forward<T>(a));
3056 }
3057
3058 // real-(long double, real128).
3059 template <typename T, typename U,
3060 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, long double>
3061 #if defined(MPPP_WITH_QUADMATH)
3062 ,
3063 std::is_same<U, real128>
3064 #endif
3065 >>::value,
3066 int> = 0>
dispatch_real_binary_sub(T && a,const U & x)3067 inline real dispatch_real_binary_sub(T &&a, const U &x)
3068 {
3069 MPPP_MAYBE_TLS real tmp;
3070 tmp.set_prec(c_max(a.get_prec(), real_deduce_precision(x)));
3071 tmp.set(x);
3072 return dispatch_real_binary_sub(std::forward<T>(a), tmp);
3073 }
3074
3075 // (long double, real128)-real.
3076 template <typename T, typename U,
3077 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, long double>
3078 #if defined(MPPP_WITH_QUADMATH)
3079 ,
3080 std::is_same<U, real128>
3081 #endif
3082 >>::value,
3083 int> = 0>
dispatch_real_binary_sub(const U & x,T && a)3084 inline real dispatch_real_binary_sub(const U &x, T &&a)
3085 {
3086 MPPP_MAYBE_TLS real tmp;
3087 tmp.set_prec(c_max(a.get_prec(), real_deduce_precision(x)));
3088 tmp.set(x);
3089 return dispatch_real_binary_sub(tmp, std::forward<T>(a));
3090 }
3091
3092 } // namespace detail
3093
3094 // Binary subtraction.
3095 #if defined(MPPP_HAVE_CONCEPTS)
3096 template <typename T, typename U>
3097 requires real_op_types<T, U>
3098 #else
3099 template <typename T, typename U, detail::enable_if_t<are_real_op_types<T, U>::value, int> = 0>
3100 #endif
operator -(T && a,U && b)3101 inline real operator-(T &&a, U &&b)
3102 {
3103 return detail::dispatch_real_binary_sub(std::forward<T>(a), std::forward<U>(b));
3104 }
3105
3106 namespace detail
3107 {
3108
3109 // real-real.
3110 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_in_place_sub(real & a,T && b)3111 inline void dispatch_real_in_place_sub(real &a, T &&b)
3112 {
3113 sub(a, a, std::forward<T>(b));
3114 }
3115
3116 MPPP_DLL_PUBLIC void dispatch_real_in_place_sub_integer_impl(real &, const ::mpz_t, ::mpfr_prec_t);
3117
3118 // real-integer.
3119 template <std::size_t SSize>
dispatch_real_in_place_sub(real & a,const integer<SSize> & n)3120 inline void dispatch_real_in_place_sub(real &a, const integer<SSize> &n)
3121 {
3122 dispatch_real_in_place_sub_integer_impl(a, n.get_mpz_view(), real_deduce_precision(n));
3123 }
3124
3125 // real-unsigned C++ integral.
3126 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_in_place_sub(real & a,const T & n)3127 inline void dispatch_real_in_place_sub(real &a, const T &n)
3128 {
3129 if (n <= nl_max<unsigned long>()) {
3130 auto wrapper
3131 = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_sub_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
3132
3133 mpfr_nary_op_impl<false>(real_deduce_precision(n), wrapper, a, a);
3134 } else {
3135 dispatch_real_in_place_sub(a, integer<2>{n});
3136 }
3137 }
3138
3139 // real-bool.
3140 // NOTE: make this explicit (rather than letting bool fold into
3141 // the unsigned integrals overload) in order to avoid MSVC warnings.
3142 MPPP_DLL_PUBLIC void dispatch_real_in_place_sub(real &, bool);
3143
3144 // real-signed C++ integral.
3145 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_in_place_sub(real & a,const T & n)3146 inline void dispatch_real_in_place_sub(real &a, const T &n)
3147 {
3148 if (n <= nl_max<long>() && n >= nl_min<long>()) {
3149 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_sub_si(r, o, static_cast<long>(n), MPFR_RNDN); };
3150
3151 mpfr_nary_op_impl<false>(real_deduce_precision(n), wrapper, a, a);
3152 } else {
3153 dispatch_real_in_place_sub(a, integer<2>{n});
3154 }
3155 }
3156
3157 MPPP_DLL_PUBLIC void dispatch_real_in_place_sub_rational_impl(real &, const ::mpq_t, ::mpfr_prec_t);
3158
3159 // real-rational.
3160 template <std::size_t SSize>
dispatch_real_in_place_sub(real & a,const rational<SSize> & q)3161 inline void dispatch_real_in_place_sub(real &a, const rational<SSize> &q)
3162 {
3163 const auto qv = detail::get_mpq_view(q);
3164 dispatch_real_in_place_sub_rational_impl(a, &qv, real_deduce_precision(q));
3165 }
3166
3167 // real-(float, double).
3168 MPPP_DLL_PUBLIC void dispatch_real_in_place_sub(real &, const float &);
3169 MPPP_DLL_PUBLIC void dispatch_real_in_place_sub(real &, const double &);
3170
3171 // real-(long double, real128).
3172 MPPP_DLL_PUBLIC void dispatch_real_in_place_sub(real &, const long double &);
3173 #if defined(MPPP_WITH_QUADMATH)
3174 MPPP_DLL_PUBLIC void dispatch_real_in_place_sub(real &, const real128 &);
3175 #endif
3176
3177 // (real interoperable)-real.
3178 template <typename T, typename U, enable_if_t<is_real_interoperable<T>::value, int> = 0>
dispatch_real_in_place_sub(T & x,U && a)3179 inline void dispatch_real_in_place_sub(T &x, U &&a)
3180 {
3181 x = static_cast<T>(x - std::forward<U>(a));
3182 }
3183
3184 } // namespace detail
3185
3186 // In-place subtraction.
3187 #if defined(MPPP_HAVE_CONCEPTS)
3188 template <typename T, typename U>
3189 requires real_in_place_op_types<T, U>
3190 #else
3191 template <typename T, typename U, detail::enable_if_t<are_real_in_place_op_types<T, U>::value, int> = 0>
3192 #endif
operator -=(T & a,U && b)3193 inline T &operator-=(T &a, U &&b)
3194 {
3195 detail::dispatch_real_in_place_sub(a, std::forward<U>(b));
3196 return a;
3197 }
3198
3199 // Prefix decrement.
3200 MPPP_DLL_PUBLIC real &operator--(real &);
3201
3202 // Suffix decrement.
3203 MPPP_DLL_PUBLIC real operator--(real &, int);
3204
3205 namespace detail
3206 {
3207
3208 // real-real.
3209 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cvr_real<U>>::value, int> = 0>
dispatch_real_binary_mul(T && a,U && b)3210 inline real dispatch_real_binary_mul(T &&a, U &&b)
3211 {
3212 return mpfr_nary_op_return_impl<true>(0, ::mpfr_mul, std::forward<T>(a), std::forward<U>(b));
3213 }
3214
3215 // real-integer.
3216 template <typename T, std::size_t SSize>
dispatch_real_binary_mul(T && a,const integer<SSize> & n)3217 inline real dispatch_real_binary_mul(T &&a, const integer<SSize> &n)
3218 {
3219 auto wrapper = [&n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_z(r, o, n.get_mpz_view(), MPFR_RNDN); };
3220
3221 // NOTE: in these mpfr_nary_op_return_impl() invocations, we are passing a min_prec
3222 // which is by definition valid because it is produced by an invocation of
3223 // real_deduce_precision() (which does clamping).
3224 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3225 }
3226
3227 // integer-real.
3228 template <typename T, std::size_t SSize>
dispatch_real_binary_mul(const integer<SSize> & n,T && a)3229 inline real dispatch_real_binary_mul(const integer<SSize> &n, T &&a)
3230 {
3231 return dispatch_real_binary_mul(std::forward<T>(a), n);
3232 }
3233
3234 // real-unsigned integral.
3235 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_unsigned_integral<U>>::value, int> = 0>
dispatch_real_binary_mul(T && a,const U & n)3236 inline real dispatch_real_binary_mul(T &&a, const U &n)
3237 {
3238 if (n <= nl_max<unsigned long>()) {
3239 auto wrapper
3240 = [n](::mpfr_t r, const ::mpfr_t o) { mpfr_mul_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
3241
3242 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3243 } else {
3244 return dispatch_real_binary_mul(std::forward<T>(a), integer<2>{n});
3245 }
3246 }
3247
3248 // real-bool.
3249 // NOTE: make this explicit (rather than letting bool fold into
3250 // the unsigned integrals overload) in order to avoid MSVC warnings.
3251 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_binary_mul(T && a,const bool & n)3252 inline real dispatch_real_binary_mul(T &&a, const bool &n)
3253 {
3254 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { mpfr_mul_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
3255 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3256 }
3257
3258 // unsigned integral-real.
3259 template <typename T, typename U, enable_if_t<conjunction<is_cpp_unsigned_integral<T>, is_cvr_real<U>>::value, int> = 0>
dispatch_real_binary_mul(const T & n,U && a)3260 inline real dispatch_real_binary_mul(const T &n, U &&a)
3261 {
3262 return dispatch_real_binary_mul(std::forward<U>(a), n);
3263 }
3264
3265 // real-signed integral.
3266 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_signed_integral<U>>::value, int> = 0>
dispatch_real_binary_mul(T && a,const U & n)3267 inline real dispatch_real_binary_mul(T &&a, const U &n)
3268 {
3269 if (n <= nl_max<long>() && n >= nl_min<long>()) {
3270 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { mpfr_mul_si(r, o, static_cast<long>(n), MPFR_RNDN); };
3271
3272 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3273 } else {
3274 return dispatch_real_binary_mul(std::forward<T>(a), integer<2>{n});
3275 }
3276 }
3277
3278 // signed integral-real.
3279 template <typename T, typename U, enable_if_t<conjunction<is_cpp_signed_integral<T>, is_cvr_real<U>>::value, int> = 0>
dispatch_real_binary_mul(const T & n,U && a)3280 inline real dispatch_real_binary_mul(const T &n, U &&a)
3281 {
3282 return dispatch_real_binary_mul(std::forward<U>(a), n);
3283 }
3284
3285 // real-rational.
3286 template <typename T, std::size_t SSize>
dispatch_real_binary_mul(T && a,const rational<SSize> & q)3287 inline real dispatch_real_binary_mul(T &&a, const rational<SSize> &q)
3288 {
3289 const auto qv = detail::get_mpq_view(q);
3290
3291 auto wrapper = [&qv](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_q(r, o, &qv, MPFR_RNDN); };
3292
3293 return mpfr_nary_op_return_impl<false>(real_deduce_precision(q), wrapper, std::forward<T>(a));
3294 }
3295
3296 // rational-real.
3297 template <typename T, std::size_t SSize>
dispatch_real_binary_mul(const rational<SSize> & q,T && a)3298 inline real dispatch_real_binary_mul(const rational<SSize> &q, T &&a)
3299 {
3300 return dispatch_real_binary_mul(std::forward<T>(a), q);
3301 }
3302
3303 // real-(float, double).
3304 template <typename T, typename U,
3305 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, float>, std::is_same<U, double>>>::value,
3306 int> = 0>
dispatch_real_binary_mul(T && a,const U & x)3307 inline real dispatch_real_binary_mul(T &&a, const U &x)
3308 {
3309 // NOTE: the MPFR docs state that mpfr_mul_d() assumes that
3310 // the radix of double is a power of 2. If we ever run into platforms
3311 // for which this is not true, we can add a compile-time dispatch
3312 // that uses the long double implementation instead.
3313 constexpr auto dradix = static_cast<unsigned>(std::numeric_limits<double>::radix);
3314 static_assert(!(dradix & (dradix - 1)), "mpfr_mul_d() requires the radix of the 'double' type to be a power of 2.");
3315
3316 auto wrapper = [x](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_d(r, o, static_cast<double>(x), MPFR_RNDN); };
3317
3318 return mpfr_nary_op_return_impl<false>(real_deduce_precision(x), wrapper, std::forward<T>(a));
3319 }
3320
3321 // (float, double)-real.
3322 template <typename T, typename U,
3323 enable_if_t<conjunction<is_cvr_real<U>, disjunction<std::is_same<T, float>, std::is_same<T, double>>>::value,
3324 int> = 0>
dispatch_real_binary_mul(const T & x,U && a)3325 inline real dispatch_real_binary_mul(const T &x, U &&a)
3326 {
3327 return dispatch_real_binary_mul(std::forward<U>(a), x);
3328 }
3329
3330 // real-(long double, real128).
3331 template <typename T, typename U,
3332 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, long double>
3333 #if defined(MPPP_WITH_QUADMATH)
3334 ,
3335 std::is_same<U, real128>
3336 #endif
3337 >>::value,
3338 int> = 0>
dispatch_real_binary_mul(T && a,const U & x)3339 inline real dispatch_real_binary_mul(T &&a, const U &x)
3340 {
3341 MPPP_MAYBE_TLS real tmp;
3342 tmp.set_prec(c_max(a.get_prec(), real_deduce_precision(x)));
3343 tmp.set(x);
3344 return dispatch_real_binary_mul(std::forward<T>(a), tmp);
3345 }
3346
3347 // (long double, real128)-real.
3348 template <typename T, typename U,
3349 enable_if_t<conjunction<disjunction<std::is_same<T, long double>
3350 #if defined(MPPP_WITH_QUADMATH)
3351 ,
3352 std::is_same<T, real128>
3353 #endif
3354 >,
3355 is_cvr_real<U>>::value,
3356 int> = 0>
dispatch_real_binary_mul(const T & x,U && a)3357 inline real dispatch_real_binary_mul(const T &x, U &&a)
3358 {
3359 return dispatch_real_binary_mul(std::forward<U>(a), x);
3360 }
3361
3362 } // namespace detail
3363
3364 // Binary multiplication.
3365 #if defined(MPPP_HAVE_CONCEPTS)
3366 template <typename T, typename U>
3367 requires real_op_types<T, U>
3368 #else
3369 template <typename T, typename U, detail::enable_if_t<are_real_op_types<T, U>::value, int> = 0>
3370 #endif
operator *(T && a,U && b)3371 inline real operator*(T &&a, U &&b)
3372 {
3373 return detail::dispatch_real_binary_mul(std::forward<T>(a), std::forward<U>(b));
3374 }
3375
3376 namespace detail
3377 {
3378
3379 // real-real.
3380 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_in_place_mul(real & a,T && b)3381 inline void dispatch_real_in_place_mul(real &a, T &&b)
3382 {
3383 mul(a, a, std::forward<T>(b));
3384 }
3385
3386 MPPP_DLL_PUBLIC void dispatch_real_in_place_mul_integer_impl(real &, const ::mpz_t, ::mpfr_prec_t);
3387
3388 // real-integer.
3389 template <std::size_t SSize>
dispatch_real_in_place_mul(real & a,const integer<SSize> & n)3390 inline void dispatch_real_in_place_mul(real &a, const integer<SSize> &n)
3391 {
3392 dispatch_real_in_place_mul_integer_impl(a, n.get_mpz_view(), real_deduce_precision(n));
3393 }
3394
3395 // real-unsigned C++ integral.
3396 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_in_place_mul(real & a,const T & n)3397 inline void dispatch_real_in_place_mul(real &a, const T &n)
3398 {
3399 if (n <= nl_max<unsigned long>()) {
3400 auto wrapper
3401 = [n](::mpfr_t r, const ::mpfr_t o) { mpfr_mul_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
3402
3403 mpfr_nary_op_impl<false>(real_deduce_precision(n), wrapper, a, a);
3404 } else {
3405 dispatch_real_in_place_mul(a, integer<2>{n});
3406 }
3407 }
3408
3409 // real-bool.
3410 // NOTE: make this explicit (rather than letting bool fold into
3411 // the unsigned integrals overload) in order to avoid MSVC warnings.
3412 MPPP_DLL_PUBLIC void dispatch_real_in_place_mul(real &, bool);
3413
3414 // real-signed C++ integral.
3415 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_in_place_mul(real & a,const T & n)3416 inline void dispatch_real_in_place_mul(real &a, const T &n)
3417 {
3418 if (n <= nl_max<long>() && n >= nl_min<long>()) {
3419 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { mpfr_mul_si(r, o, static_cast<long>(n), MPFR_RNDN); };
3420
3421 mpfr_nary_op_impl<false>(real_deduce_precision(n), wrapper, a, a);
3422 } else {
3423 dispatch_real_in_place_mul(a, integer<2>{n});
3424 }
3425 }
3426
3427 MPPP_DLL_PUBLIC void dispatch_real_in_place_mul_rational_impl(real &, const ::mpq_t, ::mpfr_prec_t);
3428
3429 // real-rational.
3430 template <std::size_t SSize>
dispatch_real_in_place_mul(real & a,const rational<SSize> & q)3431 inline void dispatch_real_in_place_mul(real &a, const rational<SSize> &q)
3432 {
3433 const auto qv = detail::get_mpq_view(q);
3434 dispatch_real_in_place_mul_rational_impl(a, &qv, real_deduce_precision(q));
3435 }
3436
3437 // real-(float, double).
3438 MPPP_DLL_PUBLIC void dispatch_real_in_place_mul(real &, const float &);
3439 MPPP_DLL_PUBLIC void dispatch_real_in_place_mul(real &, const double &);
3440
3441 // real-(long double, real128).
3442 MPPP_DLL_PUBLIC void dispatch_real_in_place_mul(real &, const long double &);
3443 #if defined(MPPP_WITH_QUADMATH)
3444 MPPP_DLL_PUBLIC void dispatch_real_in_place_mul(real &, const real128 &);
3445 #endif
3446
3447 // (real interoperable)-real.
3448 template <typename T, typename U, enable_if_t<is_real_interoperable<T>::value, int> = 0>
dispatch_real_in_place_mul(T & x,U && a)3449 inline void dispatch_real_in_place_mul(T &x, U &&a)
3450 {
3451 x = static_cast<T>(x * std::forward<U>(a));
3452 }
3453
3454 } // namespace detail
3455
3456 // In-place multiplication.
3457 #if defined(MPPP_HAVE_CONCEPTS)
3458 template <typename T, typename U>
3459 requires real_in_place_op_types<T, U>
3460 #else
3461 template <typename T, typename U, detail::enable_if_t<are_real_in_place_op_types<T, U>::value, int> = 0>
3462 #endif
operator *=(T & a,U && b)3463 inline T &operator*=(T &a, U &&b)
3464 {
3465 detail::dispatch_real_in_place_mul(a, std::forward<U>(b));
3466 return a;
3467 }
3468
3469 namespace detail
3470 {
3471
3472 // real-real.
3473 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cvr_real<U>>::value, int> = 0>
dispatch_real_binary_div(T && a,U && b)3474 inline real dispatch_real_binary_div(T &&a, U &&b)
3475 {
3476 return mpfr_nary_op_return_impl<true>(0, ::mpfr_div, std::forward<T>(a), std::forward<U>(b));
3477 }
3478
3479 // (long double, real128, integer, rational)-real.
3480 // NOTE: place it here because it is used in the
3481 // implementations below.
3482 template <
3483 typename T, typename U,
3484 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, long double>, is_integer<U>, is_rational<U>
3485 #if defined(MPPP_WITH_QUADMATH)
3486 ,
3487 std::is_same<U, real128>
3488 #endif
3489 >>::value,
3490 int> = 0>
dispatch_real_binary_div(const U & x,T && a)3491 inline real dispatch_real_binary_div(const U &x, T &&a)
3492 {
3493 MPPP_MAYBE_TLS real tmp;
3494 tmp.set_prec(c_max(a.get_prec(), real_deduce_precision(x)));
3495 tmp.set(x);
3496 return dispatch_real_binary_div(tmp, std::forward<T>(a));
3497 }
3498
3499 // real-integer.
3500 template <typename T, std::size_t SSize>
dispatch_real_binary_div(T && a,const integer<SSize> & n)3501 inline real dispatch_real_binary_div(T &&a, const integer<SSize> &n)
3502 {
3503 auto wrapper = [&n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_z(r, o, n.get_mpz_view(), MPFR_RNDN); };
3504
3505 // NOTE: in these mpfr_nary_op_return_impl() invocations, we are passing a min_prec
3506 // which is by definition valid because it is produced by an invocation of
3507 // real_deduce_precision() (which does clamping).
3508 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3509 }
3510
3511 // real-unsigned integral.
3512 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_unsigned_integral<U>>::value, int> = 0>
dispatch_real_binary_div(T && a,const U & n)3513 inline real dispatch_real_binary_div(T &&a, const U &n)
3514 {
3515 if (n <= nl_max<unsigned long>()) {
3516 auto wrapper
3517 = [n](::mpfr_t r, const ::mpfr_t o) { mpfr_div_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
3518
3519 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3520 } else {
3521 return dispatch_real_binary_div(std::forward<T>(a), integer<2>{n});
3522 }
3523 }
3524
3525 // unsigned integral-real.
3526 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_unsigned_integral<U>>::value, int> = 0>
dispatch_real_binary_div(const U & n,T && a)3527 inline real dispatch_real_binary_div(const U &n, T &&a)
3528 {
3529 if (n <= nl_max<unsigned long>()) {
3530 auto wrapper
3531 = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_ui_div(r, static_cast<unsigned long>(n), o, MPFR_RNDN); };
3532
3533 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3534 } else {
3535 return dispatch_real_binary_div(integer<2>{n}, std::forward<T>(a));
3536 }
3537 }
3538
3539 // real-bool.
3540 // NOTE: make this explicit (rather than letting bool fold into
3541 // the unsigned integrals overload) in order to avoid MSVC warnings.
3542 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_binary_div(T && a,const bool & n)3543 inline real dispatch_real_binary_div(T &&a, const bool &n)
3544 {
3545 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { mpfr_div_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
3546 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3547 }
3548
3549 // bool-real.
3550 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_binary_div(const bool & n,T && a)3551 inline real dispatch_real_binary_div(const bool &n, T &&a)
3552 {
3553 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_ui_div(r, static_cast<unsigned long>(n), o, MPFR_RNDN); };
3554 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3555 }
3556
3557 // real-signed integral.
3558 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_signed_integral<U>>::value, int> = 0>
dispatch_real_binary_div(T && a,const U & n)3559 inline real dispatch_real_binary_div(T &&a, const U &n)
3560 {
3561 if (n <= nl_max<long>() && n >= nl_min<long>()) {
3562 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { mpfr_div_si(r, o, static_cast<long>(n), MPFR_RNDN); };
3563
3564 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3565 } else {
3566 return dispatch_real_binary_div(std::forward<T>(a), integer<2>{n});
3567 }
3568 }
3569
3570 // signed integral-real.
3571 template <typename T, typename U, enable_if_t<conjunction<is_cvr_real<T>, is_cpp_signed_integral<U>>::value, int> = 0>
dispatch_real_binary_div(const U & n,T && a)3572 inline real dispatch_real_binary_div(const U &n, T &&a)
3573 {
3574 if (n <= nl_max<long>() && n >= nl_min<long>()) {
3575 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_si_div(r, static_cast<long>(n), o, MPFR_RNDN); };
3576
3577 return mpfr_nary_op_return_impl<false>(real_deduce_precision(n), wrapper, std::forward<T>(a));
3578 } else {
3579 return dispatch_real_binary_div(integer<2>{n}, std::forward<T>(a));
3580 }
3581 }
3582
3583 // real-rational.
3584 template <typename T, std::size_t SSize>
dispatch_real_binary_div(T && a,const rational<SSize> & q)3585 inline real dispatch_real_binary_div(T &&a, const rational<SSize> &q)
3586 {
3587 const auto qv = detail::get_mpq_view(q);
3588
3589 auto wrapper = [&qv](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_q(r, o, &qv, MPFR_RNDN); };
3590
3591 return mpfr_nary_op_return_impl<false>(real_deduce_precision(q), wrapper, std::forward<T>(a));
3592 }
3593
3594 // real-(float, double).
3595 template <typename T, typename U,
3596 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, float>, std::is_same<U, double>>>::value,
3597 int> = 0>
dispatch_real_binary_div(T && a,const U & x)3598 inline real dispatch_real_binary_div(T &&a, const U &x)
3599 {
3600 // NOTE: the MPFR docs state that mpfr_div_d() assumes that
3601 // the radix of double is a power of 2. If we ever run into platforms
3602 // for which this is not true, we can add a compile-time dispatch
3603 // that uses the long double implementation instead.
3604 constexpr auto dradix = static_cast<unsigned>(std::numeric_limits<double>::radix);
3605 static_assert(!(dradix & (dradix - 1)), "mpfr_div_d() requires the radix of the 'double' type to be a power of 2.");
3606
3607 auto wrapper = [x](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_d(r, o, static_cast<double>(x), MPFR_RNDN); };
3608
3609 return mpfr_nary_op_return_impl<false>(real_deduce_precision(x), wrapper, std::forward<T>(a));
3610 }
3611
3612 // (float, double)-real.
3613 template <typename T, typename U,
3614 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, float>, std::is_same<U, double>>>::value,
3615 int> = 0>
dispatch_real_binary_div(const U & x,T && a)3616 inline real dispatch_real_binary_div(const U &x, T &&a)
3617 {
3618 // NOTE: the MPFR docs state that mpfr_d_div() assumes that
3619 // the radix of double is a power of 2. If we ever run into platforms
3620 // for which this is not true, we can add a compile-time dispatch
3621 // that uses the long double implementation instead.
3622 constexpr auto dradix = static_cast<unsigned>(std::numeric_limits<double>::radix);
3623 static_assert(!(dradix & (dradix - 1)), "mpfr_d_div() requires the radix of the 'double' type to be a power of 2.");
3624
3625 auto wrapper = [x](::mpfr_t r, const ::mpfr_t o) { ::mpfr_d_div(r, static_cast<double>(x), o, MPFR_RNDN); };
3626
3627 return mpfr_nary_op_return_impl<false>(real_deduce_precision(x), wrapper, std::forward<T>(a));
3628 }
3629
3630 // real-(long double, real128).
3631 template <typename T, typename U,
3632 enable_if_t<conjunction<is_cvr_real<T>, disjunction<std::is_same<U, long double>
3633 #if defined(MPPP_WITH_QUADMATH)
3634 ,
3635 std::is_same<U, real128>
3636 #endif
3637 >>::value,
3638 int> = 0>
dispatch_real_binary_div(T && a,const U & x)3639 inline real dispatch_real_binary_div(T &&a, const U &x)
3640 {
3641 MPPP_MAYBE_TLS real tmp;
3642 tmp.set_prec(c_max(a.get_prec(), real_deduce_precision(x)));
3643 tmp.set(x);
3644 return dispatch_real_binary_div(std::forward<T>(a), tmp);
3645 }
3646
3647 } // namespace detail
3648
3649 // Binary division.
3650 #if defined(MPPP_HAVE_CONCEPTS)
3651 template <typename U, real_op_types<U> T>
3652 #else
3653 template <typename T, typename U, detail::enable_if_t<are_real_op_types<T, U>::value, int> = 0>
3654 #endif
operator /(T && a,U && b)3655 inline real operator/(T &&a, U &&b)
3656 {
3657 return detail::dispatch_real_binary_div(std::forward<T>(a), std::forward<U>(b));
3658 }
3659
3660 namespace detail
3661 {
3662
3663 // real-real.
3664 template <typename T, enable_if_t<is_cvr_real<T>::value, int> = 0>
dispatch_real_in_place_div(real & a,T && b)3665 inline void dispatch_real_in_place_div(real &a, T &&b)
3666 {
3667 div(a, a, std::forward<T>(b));
3668 }
3669
3670 MPPP_DLL_PUBLIC void dispatch_real_in_place_div_integer_impl(real &, const ::mpz_t, ::mpfr_prec_t);
3671
3672 // real-integer.
3673 template <std::size_t SSize>
dispatch_real_in_place_div(real & a,const integer<SSize> & n)3674 inline void dispatch_real_in_place_div(real &a, const integer<SSize> &n)
3675 {
3676 dispatch_real_in_place_div_integer_impl(a, n.get_mpz_view(), real_deduce_precision(n));
3677 }
3678
3679 // real-unsigned C++ integral.
3680 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_in_place_div(real & a,const T & n)3681 inline void dispatch_real_in_place_div(real &a, const T &n)
3682 {
3683 if (n <= nl_max<unsigned long>()) {
3684 auto wrapper
3685 = [n](::mpfr_t r, const ::mpfr_t o) { mpfr_div_ui(r, o, static_cast<unsigned long>(n), MPFR_RNDN); };
3686
3687 mpfr_nary_op_impl<false>(real_deduce_precision(n), wrapper, a, a);
3688 } else {
3689 dispatch_real_in_place_div(a, integer<2>{n});
3690 }
3691 }
3692
3693 // real-bool.
3694 // NOTE: make this explicit (rather than letting bool fold into
3695 // the unsigned integrals overload) in order to avoid MSVC warnings.
3696 MPPP_DLL_PUBLIC void dispatch_real_in_place_div(real &, bool);
3697
3698 // real-signed C++ integral.
3699 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_in_place_div(real & a,const T & n)3700 inline void dispatch_real_in_place_div(real &a, const T &n)
3701 {
3702 if (n <= nl_max<long>() && n >= nl_min<long>()) {
3703 auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { mpfr_div_si(r, o, static_cast<long>(n), MPFR_RNDN); };
3704
3705 mpfr_nary_op_impl<false>(real_deduce_precision(n), wrapper, a, a);
3706 } else {
3707 dispatch_real_in_place_div(a, integer<2>{n});
3708 }
3709 }
3710
3711 MPPP_DLL_PUBLIC void dispatch_real_in_place_div_rational_impl(real &, const ::mpq_t, ::mpfr_prec_t);
3712
3713 // real-rational.
3714 template <std::size_t SSize>
dispatch_real_in_place_div(real & a,const rational<SSize> & q)3715 inline void dispatch_real_in_place_div(real &a, const rational<SSize> &q)
3716 {
3717 const auto qv = detail::get_mpq_view(q);
3718 dispatch_real_in_place_div_rational_impl(a, &qv, real_deduce_precision(q));
3719 }
3720
3721 // real-(float, double).
3722 MPPP_DLL_PUBLIC void dispatch_real_in_place_div(real &, const float &);
3723 MPPP_DLL_PUBLIC void dispatch_real_in_place_div(real &, const double &);
3724
3725 // real-(long double, real128).
3726 MPPP_DLL_PUBLIC void dispatch_real_in_place_div(real &, const long double &);
3727 #if defined(MPPP_WITH_QUADMATH)
3728 MPPP_DLL_PUBLIC void dispatch_real_in_place_div(real &, const real128 &);
3729 #endif
3730
3731 // (real interoperable)-real.
3732 template <typename T, typename U, enable_if_t<is_real_interoperable<T>::value, int> = 0>
dispatch_real_in_place_div(T & x,U && a)3733 inline void dispatch_real_in_place_div(T &x, U &&a)
3734 {
3735 x = static_cast<T>(x / std::forward<U>(a));
3736 }
3737
3738 } // namespace detail
3739
3740 // In-place division.
3741 #if defined(MPPP_HAVE_CONCEPTS)
3742 template <typename T, typename U>
3743 requires real_in_place_op_types<T, U>
3744 #else
3745 template <typename T, typename U, detail::enable_if_t<are_real_in_place_op_types<T, U>::value, int> = 0>
3746 #endif
operator /=(T & a,U && b)3747 inline T &operator/=(T &a, U &&b)
3748 {
3749 detail::dispatch_real_in_place_div(a, std::forward<U>(b));
3750 return a;
3751 }
3752
3753 template <typename T, typename U>
3754 using are_real_eq_op_types
3755 = detail::disjunction<are_real_op_types<T, U>, detail::conjunction<std::is_same<T, real>, is_cpp_complex<U>>,
3756 detail::conjunction<std::is_same<U, real>, is_cpp_complex<T>>>;
3757
3758 #if defined(MPPP_HAVE_CONCEPTS)
3759
3760 template <typename T, typename U>
3761 MPPP_CONCEPT_DECL real_eq_op_types = are_real_eq_op_types<T, U>::value;
3762
3763 #endif
3764
3765 namespace detail
3766 {
3767
3768 // real-real.
3769 MPPP_DLL_PUBLIC bool dispatch_real_equality(const real &, const real &);
3770
3771 MPPP_DLL_PUBLIC bool dispatch_real_equality_integer_impl(const real &, const ::mpz_t);
3772
3773 // real-integer.
3774 template <std::size_t SSize>
dispatch_real_equality(const real & r,const integer<SSize> & n)3775 inline bool dispatch_real_equality(const real &r, const integer<SSize> &n)
3776 {
3777 return dispatch_real_equality_integer_impl(r, n.get_mpz_view());
3778 }
3779
3780 // real-unsigned c++ integral.
3781 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_equality(const real & r,const T & n)3782 inline bool dispatch_real_equality(const real &r, const T &n)
3783 {
3784 if (n <= nl_max<unsigned long>()) {
3785 if (r.nan_p()) {
3786 return false;
3787 } else {
3788 return mpfr_cmp_ui(r.get_mpfr_t(), static_cast<unsigned long>(n)) == 0;
3789 }
3790 } else {
3791 return dispatch_real_equality(r, integer<2>{n});
3792 }
3793 }
3794
3795 // NOTE: treat bool explicitly in order to avoid MSVC warnings.
3796 MPPP_DLL_PUBLIC bool dispatch_real_equality(const real &, bool);
3797
3798 // real-signed c++ integral.
3799 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_equality(const real & r,const T & n)3800 inline bool dispatch_real_equality(const real &r, const T &n)
3801 {
3802 if (n <= nl_max<long>() && n >= nl_min<long>()) {
3803 if (r.nan_p()) {
3804 return false;
3805 } else {
3806 return mpfr_cmp_si(r.get_mpfr_t(), static_cast<long>(n)) == 0;
3807 }
3808 } else {
3809 return dispatch_real_equality(r, integer<2>{n});
3810 }
3811 }
3812
3813 MPPP_DLL_PUBLIC bool dispatch_real_equality_rational_impl(const real &, const ::mpq_t);
3814
3815 // real-rational.
3816 template <std::size_t SSize>
dispatch_real_equality(const real & r,const rational<SSize> & q)3817 inline bool dispatch_real_equality(const real &r, const rational<SSize> &q)
3818 {
3819 const auto qv = detail::get_mpq_view(q);
3820
3821 return dispatch_real_equality_rational_impl(r, &qv);
3822 }
3823
3824 // real-C++ floating point.
3825 MPPP_DLL_PUBLIC bool dispatch_real_equality(const real &, const float &);
3826 MPPP_DLL_PUBLIC bool dispatch_real_equality(const real &, const double &);
3827 MPPP_DLL_PUBLIC bool dispatch_real_equality(const real &, const long double &);
3828
3829 #if defined(MPPP_WITH_QUADMATH)
3830
3831 // real-real128.
3832 MPPP_DLL_PUBLIC bool dispatch_real_equality(const real &, const real128 &);
3833
3834 #endif
3835
3836 // real-std::complex.
3837 template <typename T>
dispatch_real_equality(const real & r,const std::complex<T> & c)3838 inline bool dispatch_real_equality(const real &r, const std::complex<T> &c)
3839 {
3840 return c.imag() == T(0) && dispatch_real_equality(r, c.real());
3841 }
3842
3843 // (anything)-real.
3844 template <typename T>
dispatch_real_equality(const T & x,const real & r)3845 inline bool dispatch_real_equality(const T &x, const real &r)
3846 {
3847 return dispatch_real_equality(r, x);
3848 }
3849
3850 } // namespace detail
3851
3852 // Equality operator.
3853 #if defined(MPPP_HAVE_CONCEPTS)
3854 template <typename T, typename U>
3855 requires real_eq_op_types<T, U>
3856 #else
3857 template <typename T, typename U, detail::enable_if_t<are_real_eq_op_types<T, U>::value, int> = 0>
3858 #endif
operator ==(const T & a,const U & b)3859 inline bool operator==(const T &a, const U &b)
3860 {
3861 return detail::dispatch_real_equality(a, b);
3862 }
3863
3864 // Inequality operator.
3865 #if defined(MPPP_HAVE_CONCEPTS)
3866 template <typename T, typename U>
3867 requires real_eq_op_types<T, U>
3868 #else
3869 template <typename T, typename U, detail::enable_if_t<are_real_eq_op_types<T, U>::value, int> = 0>
3870 #endif
operator !=(const T & a,const U & b)3871 inline bool operator!=(const T &a, const U &b)
3872 {
3873 return !(a == b);
3874 }
3875
3876 namespace detail
3877 {
3878
3879 // real-real.
3880 MPPP_DLL_PUBLIC bool dispatch_real_gt(const real &, const real &);
3881
3882 MPPP_DLL_PUBLIC bool dispatch_real_gt_integer_impl(const real &, const ::mpz_t);
3883
3884 // real-integer.
3885 template <std::size_t SSize>
dispatch_real_gt(const real & r,const integer<SSize> & n)3886 inline bool dispatch_real_gt(const real &r, const integer<SSize> &n)
3887 {
3888 return dispatch_real_gt_integer_impl(r, n.get_mpz_view());
3889 }
3890
3891 MPPP_DLL_PUBLIC bool dispatch_real_gt_integer_impl(const ::mpz_t, const real &);
3892
3893 // integer-real.
3894 template <std::size_t SSize>
dispatch_real_gt(const integer<SSize> & n,const real & r)3895 inline bool dispatch_real_gt(const integer<SSize> &n, const real &r)
3896 {
3897 return dispatch_real_gt_integer_impl(n.get_mpz_view(), r);
3898 }
3899
3900 // real-unsigned c++ integral.
3901 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_gt(const real & r,const T & n)3902 inline bool dispatch_real_gt(const real &r, const T &n)
3903 {
3904 if (n <= nl_max<unsigned long>()) {
3905 if (r.nan_p()) {
3906 return false;
3907 } else {
3908 return mpfr_cmp_ui(r.get_mpfr_t(), static_cast<unsigned long>(n)) > 0;
3909 }
3910 } else {
3911 return dispatch_real_gt(r, integer<2>{n});
3912 }
3913 }
3914
3915 // unsigned c++ integral-real.
3916 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_gt(const T & n,const real & r)3917 inline bool dispatch_real_gt(const T &n, const real &r)
3918 {
3919 if (n <= nl_max<unsigned long>()) {
3920 if (r.nan_p()) {
3921 return false;
3922 } else {
3923 return mpfr_cmp_ui(r.get_mpfr_t(), static_cast<unsigned long>(n)) < 0;
3924 }
3925 } else {
3926 return dispatch_real_gt(integer<2>{n}, r);
3927 }
3928 }
3929
3930 // NOTE: treat bool explicitly in order to avoid MSVC warnings.
3931 MPPP_DLL_PUBLIC bool dispatch_real_gt(const real &, bool);
3932 MPPP_DLL_PUBLIC bool dispatch_real_gt(bool, const real &);
3933
3934 // real-signed c++ integral.
3935 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_gt(const real & r,const T & n)3936 inline bool dispatch_real_gt(const real &r, const T &n)
3937 {
3938 if (n <= nl_max<long>() && n >= nl_min<long>()) {
3939 if (r.nan_p()) {
3940 return false;
3941 } else {
3942 return mpfr_cmp_si(r.get_mpfr_t(), static_cast<long>(n)) > 0;
3943 }
3944 } else {
3945 return dispatch_real_gt(r, integer<2>{n});
3946 }
3947 }
3948
3949 // signed c++ integral-real.
3950 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_gt(const T & n,const real & r)3951 inline bool dispatch_real_gt(const T &n, const real &r)
3952 {
3953 if (n <= nl_max<long>() && n >= nl_min<long>()) {
3954 if (r.nan_p()) {
3955 return false;
3956 } else {
3957 return mpfr_cmp_si(r.get_mpfr_t(), static_cast<long>(n)) < 0;
3958 }
3959 } else {
3960 return dispatch_real_gt(integer<2>{n}, r);
3961 }
3962 }
3963
3964 MPPP_DLL_PUBLIC bool dispatch_real_gt_rational_impl(const real &, const ::mpq_t);
3965
3966 // real-rational.
3967 template <std::size_t SSize>
dispatch_real_gt(const real & r,const rational<SSize> & q)3968 inline bool dispatch_real_gt(const real &r, const rational<SSize> &q)
3969 {
3970 const auto qv = detail::get_mpq_view(q);
3971
3972 return dispatch_real_gt_rational_impl(r, &qv);
3973 }
3974
3975 MPPP_DLL_PUBLIC bool dispatch_real_gt_rational_impl(const ::mpq_t, const real &);
3976
3977 // rational-real.
3978 template <std::size_t SSize>
dispatch_real_gt(const rational<SSize> & q,const real & r)3979 inline bool dispatch_real_gt(const rational<SSize> &q, const real &r)
3980 {
3981 const auto qv = detail::get_mpq_view(q);
3982
3983 return dispatch_real_gt_rational_impl(&qv, r);
3984 }
3985
3986 // real-C++ floating point.
3987 MPPP_DLL_PUBLIC bool dispatch_real_gt(const real &, const float &);
3988 MPPP_DLL_PUBLIC bool dispatch_real_gt(const real &, const double &);
3989 MPPP_DLL_PUBLIC bool dispatch_real_gt(const real &, const long double &);
3990
3991 // C++ floating point-real.
3992 MPPP_DLL_PUBLIC bool dispatch_real_gt(const float &, const real &);
3993 MPPP_DLL_PUBLIC bool dispatch_real_gt(const double &, const real &);
3994 MPPP_DLL_PUBLIC bool dispatch_real_gt(const long double &, const real &);
3995
3996 #if defined(MPPP_WITH_QUADMATH)
3997
3998 // real-real128.
3999 MPPP_DLL_PUBLIC bool dispatch_real_gt(const real &, const real128 &);
4000
4001 // real128-real.
4002 MPPP_DLL_PUBLIC bool dispatch_real_gt(const real128 &, const real &);
4003
4004 #endif
4005
4006 } // namespace detail
4007
4008 // Greater-than operator.
4009 #if defined(MPPP_HAVE_CONCEPTS)
4010 template <typename T, typename U>
4011 requires real_op_types<T, U>
4012 #else
4013 template <typename T, typename U, detail::enable_if_t<are_real_op_types<T, U>::value, int> = 0>
4014 #endif
operator >(const T & a,const U & b)4015 inline bool operator>(const T &a, const U &b)
4016 {
4017 return detail::dispatch_real_gt(a, b);
4018 }
4019
4020 namespace detail
4021 {
4022
4023 // real-real.
4024 MPPP_DLL_PUBLIC bool dispatch_real_gte(const real &, const real &);
4025
4026 MPPP_DLL_PUBLIC bool dispatch_real_gte_integer_impl(const real &, const ::mpz_t);
4027
4028 // real-integer.
4029 template <std::size_t SSize>
dispatch_real_gte(const real & r,const integer<SSize> & n)4030 inline bool dispatch_real_gte(const real &r, const integer<SSize> &n)
4031 {
4032 return dispatch_real_gte_integer_impl(r, n.get_mpz_view());
4033 }
4034
4035 MPPP_DLL_PUBLIC bool dispatch_real_gte_integer_impl(const ::mpz_t, const real &);
4036
4037 // integer-real.
4038 template <std::size_t SSize>
dispatch_real_gte(const integer<SSize> & n,const real & r)4039 inline bool dispatch_real_gte(const integer<SSize> &n, const real &r)
4040 {
4041 return dispatch_real_gte_integer_impl(n.get_mpz_view(), r);
4042 }
4043
4044 // real-unsigned c++ integral.
4045 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_gte(const real & r,const T & n)4046 inline bool dispatch_real_gte(const real &r, const T &n)
4047 {
4048 if (n <= nl_max<unsigned long>()) {
4049 if (r.nan_p()) {
4050 return false;
4051 } else {
4052 return mpfr_cmp_ui(r.get_mpfr_t(), static_cast<unsigned long>(n)) >= 0;
4053 }
4054 } else {
4055 return dispatch_real_gte(r, integer<2>{n});
4056 }
4057 }
4058
4059 // unsigned c++ integral-real.
4060 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_gte(const T & n,const real & r)4061 inline bool dispatch_real_gte(const T &n, const real &r)
4062 {
4063 if (n <= nl_max<unsigned long>()) {
4064 if (r.nan_p()) {
4065 return false;
4066 } else {
4067 return mpfr_cmp_ui(r.get_mpfr_t(), static_cast<unsigned long>(n)) <= 0;
4068 }
4069 } else {
4070 return dispatch_real_gte(integer<2>{n}, r);
4071 }
4072 }
4073
4074 // NOTE: treat bool explicitly in order to avoid MSVC warnings.
4075 MPPP_DLL_PUBLIC bool dispatch_real_gte(const real &, bool);
4076 MPPP_DLL_PUBLIC bool dispatch_real_gte(bool, const real &);
4077
4078 // real-signed c++ integral.
4079 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_gte(const real & r,const T & n)4080 inline bool dispatch_real_gte(const real &r, const T &n)
4081 {
4082 if (n <= nl_max<long>() && n >= nl_min<long>()) {
4083 if (r.nan_p()) {
4084 return false;
4085 } else {
4086 return mpfr_cmp_si(r.get_mpfr_t(), static_cast<long>(n)) >= 0;
4087 }
4088 } else {
4089 return dispatch_real_gte(r, integer<2>{n});
4090 }
4091 }
4092
4093 // signed c++ integral-real.
4094 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_gte(const T & n,const real & r)4095 inline bool dispatch_real_gte(const T &n, const real &r)
4096 {
4097 if (n <= nl_max<long>() && n >= nl_min<long>()) {
4098 if (r.nan_p()) {
4099 return false;
4100 } else {
4101 return mpfr_cmp_si(r.get_mpfr_t(), static_cast<long>(n)) <= 0;
4102 }
4103 } else {
4104 return dispatch_real_gte(integer<2>{n}, r);
4105 }
4106 }
4107
4108 MPPP_DLL_PUBLIC bool dispatch_real_gte_rational_impl(const real &, const ::mpq_t);
4109
4110 // real-rational.
4111 template <std::size_t SSize>
dispatch_real_gte(const real & r,const rational<SSize> & q)4112 inline bool dispatch_real_gte(const real &r, const rational<SSize> &q)
4113 {
4114 const auto qv = detail::get_mpq_view(q);
4115
4116 return dispatch_real_gte_rational_impl(r, &qv);
4117 }
4118
4119 MPPP_DLL_PUBLIC bool dispatch_real_gte_rational_impl(const ::mpq_t, const real &);
4120
4121 // rational-real.
4122 template <std::size_t SSize>
dispatch_real_gte(const rational<SSize> & q,const real & r)4123 inline bool dispatch_real_gte(const rational<SSize> &q, const real &r)
4124 {
4125 const auto qv = detail::get_mpq_view(q);
4126
4127 return dispatch_real_gte_rational_impl(&qv, r);
4128 }
4129
4130 // real-C++ floating point.
4131 MPPP_DLL_PUBLIC bool dispatch_real_gte(const real &, const float &);
4132 MPPP_DLL_PUBLIC bool dispatch_real_gte(const real &, const double &);
4133 MPPP_DLL_PUBLIC bool dispatch_real_gte(const real &, const long double &);
4134
4135 // C++ floating point-real.
4136 MPPP_DLL_PUBLIC bool dispatch_real_gte(const float &, const real &);
4137 MPPP_DLL_PUBLIC bool dispatch_real_gte(const double &, const real &);
4138 MPPP_DLL_PUBLIC bool dispatch_real_gte(const long double &, const real &);
4139
4140 #if defined(MPPP_WITH_QUADMATH)
4141
4142 // real-real128.
4143 MPPP_DLL_PUBLIC bool dispatch_real_gte(const real &, const real128 &);
4144
4145 // real128-real.
4146 MPPP_DLL_PUBLIC bool dispatch_real_gte(const real128 &, const real &);
4147
4148 #endif
4149
4150 } // namespace detail
4151
4152 // Greater-than or equal operator.
4153 #if defined(MPPP_HAVE_CONCEPTS)
4154 template <typename T, typename U>
4155 requires real_op_types<T, U>
4156 #else
4157 template <typename T, typename U, detail::enable_if_t<are_real_op_types<T, U>::value, int> = 0>
4158 #endif
operator >=(const T & a,const U & b)4159 inline bool operator>=(const T &a, const U &b)
4160 {
4161 return detail::dispatch_real_gte(a, b);
4162 }
4163
4164 namespace detail
4165 {
4166
4167 // real-real.
4168 MPPP_DLL_PUBLIC bool dispatch_real_lt(const real &, const real &);
4169
4170 MPPP_DLL_PUBLIC bool dispatch_real_lt_integer_impl(const real &, const ::mpz_t);
4171
4172 // real-integer.
4173 template <std::size_t SSize>
dispatch_real_lt(const real & r,const integer<SSize> & n)4174 inline bool dispatch_real_lt(const real &r, const integer<SSize> &n)
4175 {
4176 return dispatch_real_lt_integer_impl(r, n.get_mpz_view());
4177 }
4178
4179 MPPP_DLL_PUBLIC bool dispatch_real_lt_integer_impl(const ::mpz_t, const real &);
4180
4181 // integer-real.
4182 template <std::size_t SSize>
dispatch_real_lt(const integer<SSize> & n,const real & r)4183 inline bool dispatch_real_lt(const integer<SSize> &n, const real &r)
4184 {
4185 return dispatch_real_lt_integer_impl(n.get_mpz_view(), r);
4186 }
4187
4188 // real-unsigned c++ integral.
4189 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_lt(const real & r,const T & n)4190 inline bool dispatch_real_lt(const real &r, const T &n)
4191 {
4192 if (n <= nl_max<unsigned long>()) {
4193 if (r.nan_p()) {
4194 return false;
4195 } else {
4196 return mpfr_cmp_ui(r.get_mpfr_t(), static_cast<unsigned long>(n)) < 0;
4197 }
4198 } else {
4199 return dispatch_real_lt(r, integer<2>{n});
4200 }
4201 }
4202
4203 // unsigned c++ integral-real.
4204 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_lt(const T & n,const real & r)4205 inline bool dispatch_real_lt(const T &n, const real &r)
4206 {
4207 if (n <= nl_max<unsigned long>()) {
4208 if (r.nan_p()) {
4209 return false;
4210 } else {
4211 return mpfr_cmp_ui(r.get_mpfr_t(), static_cast<unsigned long>(n)) > 0;
4212 }
4213 } else {
4214 return dispatch_real_lt(integer<2>{n}, r);
4215 }
4216 }
4217
4218 // NOTE: treat bool explicitly in order to avoid MSVC warnings.
4219 MPPP_DLL_PUBLIC bool dispatch_real_lt(const real &, bool);
4220 MPPP_DLL_PUBLIC bool dispatch_real_lt(bool, const real &);
4221
4222 // real-signed c++ integral.
4223 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_lt(const real & r,const T & n)4224 inline bool dispatch_real_lt(const real &r, const T &n)
4225 {
4226 if (n <= nl_max<long>() && n >= nl_min<long>()) {
4227 if (r.nan_p()) {
4228 return false;
4229 } else {
4230 return mpfr_cmp_si(r.get_mpfr_t(), static_cast<long>(n)) < 0;
4231 }
4232 } else {
4233 return dispatch_real_lt(r, integer<2>{n});
4234 }
4235 }
4236
4237 // signed c++ integral-real.
4238 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_lt(const T & n,const real & r)4239 inline bool dispatch_real_lt(const T &n, const real &r)
4240 {
4241 if (n <= nl_max<long>() && n >= nl_min<long>()) {
4242 if (r.nan_p()) {
4243 return false;
4244 } else {
4245 return mpfr_cmp_si(r.get_mpfr_t(), static_cast<long>(n)) > 0;
4246 }
4247 } else {
4248 return dispatch_real_lt(integer<2>{n}, r);
4249 }
4250 }
4251
4252 MPPP_DLL_PUBLIC bool dispatch_real_lt_rational_impl(const real &, const ::mpq_t);
4253
4254 // real-rational.
4255 template <std::size_t SSize>
dispatch_real_lt(const real & r,const rational<SSize> & q)4256 inline bool dispatch_real_lt(const real &r, const rational<SSize> &q)
4257 {
4258 const auto qv = detail::get_mpq_view(q);
4259
4260 return dispatch_real_lt_rational_impl(r, &qv);
4261 }
4262
4263 MPPP_DLL_PUBLIC bool dispatch_real_lt_rational_impl(const ::mpq_t, const real &);
4264
4265 // rational-real.
4266 template <std::size_t SSize>
dispatch_real_lt(const rational<SSize> & q,const real & r)4267 inline bool dispatch_real_lt(const rational<SSize> &q, const real &r)
4268 {
4269 const auto qv = detail::get_mpq_view(q);
4270
4271 return dispatch_real_lt_rational_impl(&qv, r);
4272 }
4273
4274 // real-C++ floating point.
4275 MPPP_DLL_PUBLIC bool dispatch_real_lt(const real &, const float &);
4276 MPPP_DLL_PUBLIC bool dispatch_real_lt(const real &, const double &);
4277 MPPP_DLL_PUBLIC bool dispatch_real_lt(const real &, const long double &);
4278
4279 // C++ floating point-real.
4280 MPPP_DLL_PUBLIC bool dispatch_real_lt(const float &, const real &);
4281 MPPP_DLL_PUBLIC bool dispatch_real_lt(const double &, const real &);
4282 MPPP_DLL_PUBLIC bool dispatch_real_lt(const long double &, const real &);
4283
4284 #if defined(MPPP_WITH_QUADMATH)
4285
4286 // real-real128.
4287 MPPP_DLL_PUBLIC bool dispatch_real_lt(const real &, const real128 &);
4288
4289 // real128-real.
4290 MPPP_DLL_PUBLIC bool dispatch_real_lt(const real128 &, const real &);
4291
4292 #endif
4293
4294 } // namespace detail
4295
4296 // Less-than operator.
4297 #if defined(MPPP_HAVE_CONCEPTS)
4298 template <typename T, typename U>
4299 requires real_op_types<T, U>
4300 #else
4301 template <typename T, typename U, detail::enable_if_t<are_real_op_types<T, U>::value, int> = 0>
4302 #endif
operator <(const T & a,const U & b)4303 inline bool operator<(const T &a, const U &b)
4304 {
4305 return detail::dispatch_real_lt(a, b);
4306 }
4307
4308 namespace detail
4309 {
4310
4311 // real-real.
4312 MPPP_DLL_PUBLIC bool dispatch_real_lte(const real &, const real &);
4313
4314 MPPP_DLL_PUBLIC bool dispatch_real_lte_integer_impl(const real &, const ::mpz_t);
4315
4316 // real-integer.
4317 template <std::size_t SSize>
dispatch_real_lte(const real & r,const integer<SSize> & n)4318 inline bool dispatch_real_lte(const real &r, const integer<SSize> &n)
4319 {
4320 return dispatch_real_lte_integer_impl(r, n.get_mpz_view());
4321 }
4322
4323 MPPP_DLL_PUBLIC bool dispatch_real_lte_integer_impl(const ::mpz_t, const real &);
4324
4325 // integer-real.
4326 template <std::size_t SSize>
dispatch_real_lte(const integer<SSize> & n,const real & r)4327 inline bool dispatch_real_lte(const integer<SSize> &n, const real &r)
4328 {
4329 return dispatch_real_lte_integer_impl(n.get_mpz_view(), r);
4330 }
4331
4332 // real-unsigned c++ integral.
4333 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_lte(const real & r,const T & n)4334 inline bool dispatch_real_lte(const real &r, const T &n)
4335 {
4336 if (n <= nl_max<unsigned long>()) {
4337 if (r.nan_p()) {
4338 return false;
4339 } else {
4340 return mpfr_cmp_ui(r.get_mpfr_t(), static_cast<unsigned long>(n)) <= 0;
4341 }
4342 } else {
4343 return dispatch_real_lte(r, integer<2>{n});
4344 }
4345 }
4346
4347 // unsigned c++ integral-real.
4348 template <typename T, enable_if_t<is_cpp_unsigned_integral<T>::value, int> = 0>
dispatch_real_lte(const T & n,const real & r)4349 inline bool dispatch_real_lte(const T &n, const real &r)
4350 {
4351 if (n <= nl_max<unsigned long>()) {
4352 if (r.nan_p()) {
4353 return false;
4354 } else {
4355 return mpfr_cmp_ui(r.get_mpfr_t(), static_cast<unsigned long>(n)) >= 0;
4356 }
4357 } else {
4358 return dispatch_real_lte(integer<2>{n}, r);
4359 }
4360 }
4361
4362 // NOTE: treat bool explicitly in order to avoid MSVC warnings.
4363 MPPP_DLL_PUBLIC bool dispatch_real_lte(const real &, bool);
4364 MPPP_DLL_PUBLIC bool dispatch_real_lte(bool, const real &);
4365
4366 // real-signed c++ integral.
4367 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_lte(const real & r,const T & n)4368 inline bool dispatch_real_lte(const real &r, const T &n)
4369 {
4370 if (n <= nl_max<long>() && n >= nl_min<long>()) {
4371 if (r.nan_p()) {
4372 return false;
4373 } else {
4374 return mpfr_cmp_si(r.get_mpfr_t(), static_cast<long>(n)) <= 0;
4375 }
4376 } else {
4377 return dispatch_real_lte(r, integer<2>{n});
4378 }
4379 }
4380
4381 // signed c++ integral-real.
4382 template <typename T, enable_if_t<is_cpp_signed_integral<T>::value, int> = 0>
dispatch_real_lte(const T & n,const real & r)4383 inline bool dispatch_real_lte(const T &n, const real &r)
4384 {
4385 if (n <= nl_max<long>() && n >= nl_min<long>()) {
4386 if (r.nan_p()) {
4387 return false;
4388 } else {
4389 return mpfr_cmp_si(r.get_mpfr_t(), static_cast<long>(n)) >= 0;
4390 }
4391 } else {
4392 return dispatch_real_lte(integer<2>{n}, r);
4393 }
4394 }
4395
4396 MPPP_DLL_PUBLIC bool dispatch_real_lte_rational_impl(const real &, const ::mpq_t);
4397
4398 // real-rational.
4399 template <std::size_t SSize>
dispatch_real_lte(const real & r,const rational<SSize> & q)4400 inline bool dispatch_real_lte(const real &r, const rational<SSize> &q)
4401 {
4402 const auto qv = detail::get_mpq_view(q);
4403
4404 return dispatch_real_lte_rational_impl(r, &qv);
4405 }
4406
4407 MPPP_DLL_PUBLIC bool dispatch_real_lte_rational_impl(const ::mpq_t, const real &);
4408
4409 // rational-real.
4410 template <std::size_t SSize>
dispatch_real_lte(const rational<SSize> & q,const real & r)4411 inline bool dispatch_real_lte(const rational<SSize> &q, const real &r)
4412 {
4413 const auto qv = detail::get_mpq_view(q);
4414
4415 return dispatch_real_lte_rational_impl(&qv, r);
4416 }
4417
4418 // real-C++ floating point.
4419 MPPP_DLL_PUBLIC bool dispatch_real_lte(const real &, const float &);
4420 MPPP_DLL_PUBLIC bool dispatch_real_lte(const real &, const double &);
4421 MPPP_DLL_PUBLIC bool dispatch_real_lte(const real &, const long double &);
4422
4423 // C++ floating point-real.
4424 MPPP_DLL_PUBLIC bool dispatch_real_lte(const float &, const real &);
4425 MPPP_DLL_PUBLIC bool dispatch_real_lte(const double &, const real &);
4426 MPPP_DLL_PUBLIC bool dispatch_real_lte(const long double &, const real &);
4427
4428 #if defined(MPPP_WITH_QUADMATH)
4429
4430 // real-real128.
4431 MPPP_DLL_PUBLIC bool dispatch_real_lte(const real &, const real128 &);
4432
4433 // real128-real.
4434 MPPP_DLL_PUBLIC bool dispatch_real_lte(const real128 &, const real &);
4435
4436 #endif
4437
4438 } // namespace detail
4439
4440 // Less-than or equal operator.
4441 #if defined(MPPP_HAVE_CONCEPTS)
4442 template <typename T, typename U>
4443 requires real_op_types<T, U>
4444 #else
4445 template <typename T, typename U, detail::enable_if_t<are_real_op_types<T, U>::value, int> = 0>
4446 #endif
operator <=(const T & a,const U & b)4447 inline bool operator<=(const T &a, const U &b)
4448 {
4449 return detail::dispatch_real_lte(a, b);
4450 }
4451
4452 // Implementation of integer's assignment
4453 // from real.
4454 template <std::size_t SSize>
operator =(const real & x)4455 inline integer<SSize> &integer<SSize>::operator=(const real &x)
4456 {
4457 // NOLINTNEXTLINE(cppcoreguidelines-c-copy-assignment-signature, misc-unconventional-assign-operator)
4458 return *this = static_cast<integer<SSize>>(x);
4459 }
4460
4461 // Implementation of rational's assignment
4462 // from real.
4463 template <std::size_t SSize>
operator =(const real & x)4464 inline rational<SSize> &rational<SSize>::operator=(const real &x)
4465 {
4466 // NOLINTNEXTLINE(cppcoreguidelines-c-copy-assignment-signature, misc-unconventional-assign-operator)
4467 return *this = static_cast<rational<SSize>>(x);
4468 }
4469
4470 } // namespace mppp
4471
4472 #if defined(MPPP_WITH_BOOST_S11N)
4473
4474 // Never track the address of real objects
4475 // during serialization.
4476 BOOST_CLASS_TRACKING(mppp::real, boost::serialization::track_never)
4477
4478 #endif
4479
4480 #include <mp++/detail/real_literals.hpp>
4481
4482 // Support for pretty printing in xeus-cling.
4483 #if defined(__CLING__)
4484
4485 #if __has_include(<nlohmann/json.hpp>)
4486
4487 #include <nlohmann/json.hpp>
4488
4489 namespace mppp
4490 {
4491
mime_bundle_repr(const real & x)4492 inline nlohmann::json mime_bundle_repr(const real &x)
4493 {
4494 auto bundle = nlohmann::json::object();
4495
4496 bundle["text/plain"] = x.to_string();
4497
4498 return bundle;
4499 }
4500
4501 } // namespace mppp
4502
4503 #endif
4504
4505 #endif
4506
4507 #else
4508
4509 #error The real.hpp header was included but mp++ was not configured with the MPPP_WITH_MPFR option.
4510
4511 #endif
4512
4513 #endif
4514