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