kfr

Fast, modern C++ DSP framework, FFT, Sample Rate Conversion, FIR/IIR/Biquad Filters (SSE, AVX, AVX-512, ARM NEON)
Log | Files | Refs | README

mpfrplus.hpp (29735B)


      1 #pragma once
      2 
      3 #define MPFR_STRING(str) #str
      4 
      5 #if defined __clang__
      6 #define MPFR_DIAG_PRAGMA(pragma) _Pragma(MPFR_STRING(clang diagnostic pragma))
      7 #else
      8 #define MPFR_DIAG_PRAGMA(pragma)
      9 #endif
     10 
     11 MPFR_DIAG_PRAGMA(push)
     12 MPFR_DIAG_PRAGMA(ignored "-Wreserved-id-macro")
     13 MPFR_DIAG_PRAGMA(ignored "-Wpadded")
     14 MPFR_DIAG_PRAGMA(ignored "-Wdeprecated")
     15 MPFR_DIAG_PRAGMA(ignored "-Wshorten-64-to-32")
     16 MPFR_DIAG_PRAGMA(ignored "-Wsign-conversion")
     17 #include "mpfr.h"
     18 MPFR_DIAG_PRAGMA(pop)
     19 #include <cmath>
     20 #include <limits>
     21 #include <string>
     22 #include <type_traits>
     23 
     24 namespace mpfr
     25 {
     26 
     27 /// Rounding mode values
     28 enum rounding_mode
     29 {
     30     nearest        = MPFR_RNDN, // round to nearest, with ties to even
     31     toward_zero    = MPFR_RNDZ, // round toward zero
     32     toward_pinf    = MPFR_RNDU, // round toward +Inf
     33     toward_ninf    = MPFR_RNDD, // round toward -Inf
     34     away_from_zero = MPFR_RNDA // round away from zero
     35 };
     36 
     37 constexpr int extra_precision = 5;
     38 
     39 namespace internal
     40 {
     41 
     42 struct with_precision_t
     43 {
     44 };
     45 
     46 constexpr with_precision_t with_precision{};
     47 } // namespace internal
     48 
     49 namespace internal
     50 {
     51 inline mpfr_prec_t& precision()
     52 {
     53     static mpfr_prec_t prec = mpfr_get_default_prec();
     54     return prec;
     55 }
     56 inline mpfr_rnd_t& rounding_mode()
     57 {
     58     static mpfr_rnd_t rnd = mpfr_get_default_rounding_mode();
     59     return rnd;
     60 }
     61 } // namespace internal
     62 
     63 /// Temporarily sets the precision
     64 struct scoped_precision
     65 {
     66     inline scoped_precision(mpfr_prec_t prec)
     67     {
     68         saved_precision       = internal::precision();
     69         internal::precision() = prec;
     70     }
     71     inline ~scoped_precision() { internal::precision() = saved_precision; }
     72     mpfr_prec_t saved_precision;
     73 };
     74 
     75 /// Temporarily sets the rounding mode
     76 struct scoped_rounding_mode
     77 {
     78     inline scoped_rounding_mode(mpfr_rnd_t rnd)
     79     {
     80         saved_rounding_mode       = internal::rounding_mode();
     81         internal::rounding_mode() = rnd;
     82     }
     83     inline scoped_rounding_mode(rounding_mode rnd)
     84     {
     85         saved_rounding_mode       = internal::rounding_mode();
     86         internal::rounding_mode() = static_cast<mpfr_rnd_t>(rnd);
     87     }
     88     inline ~scoped_rounding_mode() { internal::rounding_mode() = saved_rounding_mode; }
     89     mpfr_rnd_t saved_rounding_mode;
     90 };
     91 
     92 #define MPFR_MACRO_CONCAT2(x, y) x##y
     93 #define MPFR_MACRO_CONCAT(x, y) MPFR_MACRO_CONCAT2(x, y)
     94 
     95 #define MPFR_CXX_ROUNDING_MODE(rnd)                                                                          \
     96     MPFR_MACRO_CONCAT(::mpfr::scoped_rounding_mode rounding_mode, __COUNTER__)(rnd)
     97 #define MPFR_CXX_PRECISION(prec) MPFR_MACRO_CONCAT(::mpfr::scoped_precision precision, __COUNTER__)(prec)
     98 
     99 #define MPFR_FN(fn)                                                                                          \
    100     struct fn_##fn                                                                                           \
    101     {                                                                                                        \
    102         template <typename... Args>                                                                          \
    103         inline decltype(fn(std::declval<Args>()...)) operator()(Args&&... args) const                        \
    104         {                                                                                                    \
    105             return fn(std::forward<Args>(args)...);                                                          \
    106         }                                                                                                    \
    107     };
    108 
    109 #define MPFR_CXX_CMP(op, mpfr_fn)                                                                            \
    110     inline bool operator op(const number& x, const number& y)                                                \
    111     {                                                                                                        \
    112         return mpfr_fn(x.mpfr_val(), y.mpfr_val()) op 0;                                                     \
    113     }
    114 #define MPFR_CXX_CMP_RT(op, mpfr_fn, type)                                                                   \
    115     inline bool operator op(const number& x, type y) { return mpfr_fn(x.mpfr_val(), y) op 0; }
    116 #define MPFR_CXX_CMP_TR(op, mpfr_fn, type)                                                                   \
    117     inline bool operator op(type x, const number& y) { return 0 op mpfr_fn(y.mpfr_val(), x); }
    118 #define MPFR_CXX_UNARY(fn, mpfr_fn)                                                                          \
    119     inline number fn(const number& x)                                                                        \
    120     {                                                                                                        \
    121         number result;                                                                                       \
    122         mpfr_fn(result.mpfr_val(), x.mpfr_val(), internal::rounding_mode());                                 \
    123         return result;                                                                                       \
    124     }
    125 #define MPFR_CXX_UNARY_RND(fn, mpfr_fn, rnd)                                                                 \
    126     inline number fn(const number& x)                                                                        \
    127     {                                                                                                        \
    128         number result;                                                                                       \
    129         mpfr_fn(result.mpfr_val(), x.mpfr_val() rnd);                                                        \
    130         return result;                                                                                       \
    131     }
    132 #define MPFR_CXX_BINARY(fn, mpfr_fn)                                                                         \
    133     inline number fn(const number& x, const number& y)                                                       \
    134     {                                                                                                        \
    135         number result;                                                                                       \
    136         mpfr_fn(result.mpfr_val(), x.mpfr_val(), y.mpfr_val(), internal::rounding_mode());                   \
    137         return result;                                                                                       \
    138     }
    139 #define MPFR_CXX_TERNARY(fn, mpfr_fn)                                                                        \
    140     inline number fn(const number& x, const number& y, const number& z)                                      \
    141     {                                                                                                        \
    142         number result;                                                                                       \
    143         mpfr_fn(result.mpfr_val(), x.mpfr_val(), y.mpfr_val(), z.mpfr_val(), internal::rounding_mode());     \
    144         return result;                                                                                       \
    145     }
    146 #define MPFR_CXX_BINARY_RT(fn, mpfr_fn, type)                                                                \
    147     inline number fn(const number& x, type y)                                                                \
    148     {                                                                                                        \
    149         number result;                                                                                       \
    150         mpfr_fn(result.mpfr_val(), x.mpfr_val(), y, internal::rounding_mode());                              \
    151         return result;                                                                                       \
    152     }
    153 #define MPFR_CXX_BINARY_TR(fn, mpfr_fn, type)                                                                \
    154     inline number fn(type x, const number& y)                                                                \
    155     {                                                                                                        \
    156         number result;                                                                                       \
    157         mpfr_fn(result.mpfr_val(), x, y.mpfr_val(), internal::rounding_mode());                              \
    158         return result;                                                                                       \
    159     }
    160 #define MPFR_CXX_INPLACE_BINARY(fn, mpfr_fn)                                                                 \
    161     inline number& fn(number& x, const number& y)                                                            \
    162     {                                                                                                        \
    163         mpfr_fn(x.mpfr_val(), x.mpfr_val(), y.mpfr_val(), internal::rounding_mode());                        \
    164         return x;                                                                                            \
    165     }
    166 #define MPFR_CXX_INPLACE_BINARY_T(fn, mpfr_fn, type)                                                         \
    167     inline number& fn(number& x, type y)                                                                     \
    168     {                                                                                                        \
    169         mpfr_fn(x.mpfr_val(), x.mpfr_val(), y, internal::rounding_mode());                                   \
    170         return x;                                                                                            \
    171     }
    172 #define MPFR_CXX_CONST(fn, mpfr_fn)                                                                          \
    173     inline number fn()                                                                                       \
    174     {                                                                                                        \
    175         number result;                                                                                       \
    176         mpfr_fn(result.mpfr_val(), internal::rounding_mode());                                               \
    177         return result;                                                                                       \
    178     }
    179 
    180 /// Floating-point value with arbitrary precision
    181 struct number
    182 {
    183 private:
    184     mpfr_t val;
    185     bool owns;
    186 
    187 public:
    188     inline number() noexcept : owns(true) { mpfr_init2(val, internal::precision()); }
    189     inline ~number()
    190     {
    191         if (owns)
    192             mpfr_clear(val);
    193     }
    194     inline number(internal::with_precision_t, mpfr_prec_t precision) : owns(true)
    195     {
    196         mpfr_init2(val, precision);
    197     }
    198     /// Copy-constructor
    199     inline number(const number& value) noexcept : owns(true)
    200     {
    201         mpfr_init2(val, internal::precision());
    202         mpfr_set(val, value.val, internal::rounding_mode());
    203     }
    204     /// Move-constructor
    205     inline number(number&& value) noexcept : owns(true)
    206     {
    207         val[0]     = value.val[0];
    208         value.owns = false;
    209     }
    210     /**
    211     Construct from string
    212     `str`  C-string containing string representation of the number
    213     `base` Radix (defaults to 10)
    214 
    215     See also operator""_mpfr
    216     */
    217     inline number(const char* str, int base = 10) noexcept : owns(true)
    218     {
    219         mpfr_init2(val, internal::precision());
    220         mpfr_set_str(val, str, base, internal::rounding_mode());
    221     }
    222     inline number(std::nullptr_t) = delete;
    223 #define MPFR_CXX_CTOR_T(mpfr_fn, type)                                                                       \
    224     inline number(type value) noexcept : owns(true)                                                          \
    225     {                                                                                                        \
    226         mpfr_init2(val, internal::precision());                                                              \
    227         mpfr_fn(val, value, internal::rounding_mode());                                                      \
    228     }
    229 #define MPFR_CXX_ASGN_T(mpfr_fn, type)                                                                       \
    230     inline number& operator=(type value) noexcept                                                            \
    231     {                                                                                                        \
    232         mpfr_fn(val, value, internal::rounding_mode());                                                      \
    233         return *this;                                                                                        \
    234     }
    235     MPFR_CXX_CTOR_T(mpfr_set_d, double)
    236     MPFR_CXX_CTOR_T(mpfr_set_ld, long double)
    237     MPFR_CXX_CTOR_T(mpfr_set_flt, float)
    238     MPFR_CXX_CTOR_T(mpfr_set_si, int)
    239     MPFR_CXX_CTOR_T(mpfr_set_ui, unsigned int)
    240     MPFR_CXX_CTOR_T(mpfr_set_si, long int)
    241     MPFR_CXX_CTOR_T(mpfr_set_ui, unsigned long int)
    242 #ifdef _MPFR_H_HAVE_INTMAX_T
    243     MPFR_CXX_CTOR_T(mpfr_set_sj, intmax_t)
    244     MPFR_CXX_CTOR_T(mpfr_set_uj, uintmax_t)
    245 #endif
    246 
    247     MPFR_CXX_ASGN_T(mpfr_set_d, double)
    248     MPFR_CXX_ASGN_T(mpfr_set_ld, long double)
    249     MPFR_CXX_ASGN_T(mpfr_set_flt, float)
    250     MPFR_CXX_ASGN_T(mpfr_set_si, int)
    251     MPFR_CXX_ASGN_T(mpfr_set_ui, unsigned int)
    252     MPFR_CXX_ASGN_T(mpfr_set_si, long int)
    253     MPFR_CXX_ASGN_T(mpfr_set_ui, unsigned long int)
    254 #ifdef _MPFR_H_HAVE_INTMAX_T
    255     MPFR_CXX_ASGN_T(mpfr_set_sj, intmax_t)
    256     MPFR_CXX_ASGN_T(mpfr_set_uj, uintmax_t)
    257 #endif
    258 
    259     inline number& operator=(const number& value) noexcept
    260     {
    261         mpfr_set(val, value.val, internal::rounding_mode());
    262         return *this;
    263     }
    264     /// Get internal mpfr value
    265     inline mpfr_ptr mpfr_val() { return &val[0]; }
    266     inline mpfr_srcptr mpfr_val() const { return &val[0]; }
    267 
    268     MPFR_DIAG_PRAGMA(push)
    269     MPFR_DIAG_PRAGMA(ignored "-Wold-style-cast")
    270     /// Current precision
    271     inline mpfr_prec_t prec() const { return mpfr_get_prec(val); }
    272     /// Is this number zero?
    273     inline bool iszero() const { return mpfr_zero_p(val); }
    274     /// Is this number NaN?
    275     inline bool isnan() const { return mpfr_nan_p(val); }
    276     /// Is this number Infinity?
    277     inline bool isinfinity() const { return mpfr_inf_p(val); }
    278     /// Is this number finite?
    279     inline bool isfinite() const { return !isinfinity() && !isnan(); }
    280     /// Returns sign of the number
    281     inline int sign() const { return MPFR_SIGN(val); }
    282 
    283     MPFR_DIAG_PRAGMA(pop)
    284     /// Converts precision
    285     inline number withprec(mpfr_prec_t prec) const
    286     {
    287         number r(internal::with_precision, prec);
    288         r = *this;
    289         return r;
    290     }
    291 
    292     /// Explicit conversion to float
    293     inline explicit operator float() const noexcept { return mpfr_get_flt(val, internal::rounding_mode()); }
    294     /// Explicit conversion to double
    295     inline explicit operator double() const noexcept { return mpfr_get_d(val, internal::rounding_mode()); }
    296     /// Explicit conversion to long double
    297     inline explicit operator long double() const noexcept
    298     {
    299         return mpfr_get_ld(val, internal::rounding_mode());
    300     }
    301 
    302     std::string to_string() const
    303     {
    304         char* str;
    305         mpfr_asprintf(&str, "%.*Rg", prec(), val);
    306         std::string result = str;
    307         mpfr_free_str(str);
    308         return result;
    309     }
    310 };
    311 
    312 #ifdef MPFR_USE_UDL
    313 /**
    314 Constructs number from string using User-Defined-Literal (C++11).
    315 #### Example:
    316 ```C++
    317 using namespace mpfr;
    318 number r = "1.1111111111111111111111111111111111111111"_mpfr;
    319 ```
    320 */
    321 inline number operator""_mpfr(const char* str, size_t) { return str; }
    322 /// Constructs number from long double using User-Defined-Literal (C++11)
    323 inline number operator""_mpfr(long double val) { return val; }
    324 #endif
    325 
    326 MPFR_CXX_CMP(==, mpfr_cmp)
    327 MPFR_CXX_CMP(!=, mpfr_cmp)
    328 MPFR_CXX_CMP(<, mpfr_cmp)
    329 MPFR_CXX_CMP(>, mpfr_cmp)
    330 MPFR_CXX_CMP(>=, mpfr_cmp)
    331 MPFR_CXX_CMP(<=, mpfr_cmp)
    332 MPFR_CXX_CMP_RT(==, mpfr_cmp_si, int)
    333 MPFR_CXX_CMP_RT(!=, mpfr_cmp_si, int)
    334 MPFR_CXX_CMP_RT(<, mpfr_cmp_si, int)
    335 MPFR_CXX_CMP_RT(>, mpfr_cmp_si, int)
    336 MPFR_CXX_CMP_RT(>=, mpfr_cmp_si, int)
    337 MPFR_CXX_CMP_RT(<=, mpfr_cmp_si, int)
    338 MPFR_CXX_CMP_RT(==, mpfr_cmp_ui, unsigned int)
    339 MPFR_CXX_CMP_RT(!=, mpfr_cmp_ui, unsigned int)
    340 MPFR_CXX_CMP_RT(<, mpfr_cmp_ui, unsigned int)
    341 MPFR_CXX_CMP_RT(>, mpfr_cmp_ui, unsigned int)
    342 MPFR_CXX_CMP_RT(>=, mpfr_cmp_ui, unsigned int)
    343 MPFR_CXX_CMP_RT(<=, mpfr_cmp_ui, unsigned int)
    344 MPFR_CXX_CMP_RT(==, mpfr_cmp_si, long int)
    345 MPFR_CXX_CMP_RT(!=, mpfr_cmp_si, long int)
    346 MPFR_CXX_CMP_RT(<, mpfr_cmp_si, long int)
    347 MPFR_CXX_CMP_RT(>, mpfr_cmp_si, long int)
    348 MPFR_CXX_CMP_RT(>=, mpfr_cmp_si, long int)
    349 MPFR_CXX_CMP_RT(<=, mpfr_cmp_si, long int)
    350 MPFR_CXX_CMP_RT(==, mpfr_cmp_ui, unsigned long int)
    351 MPFR_CXX_CMP_RT(!=, mpfr_cmp_ui, unsigned long int)
    352 MPFR_CXX_CMP_RT(<, mpfr_cmp_ui, unsigned long int)
    353 MPFR_CXX_CMP_RT(>, mpfr_cmp_ui, unsigned long int)
    354 MPFR_CXX_CMP_RT(>=, mpfr_cmp_ui, unsigned long int)
    355 MPFR_CXX_CMP_RT(<=, mpfr_cmp_ui, unsigned long int)
    356 MPFR_CXX_CMP_RT(==, mpfr_cmp_d, double)
    357 MPFR_CXX_CMP_RT(!=, mpfr_cmp_d, double)
    358 MPFR_CXX_CMP_RT(<, mpfr_cmp_d, double)
    359 MPFR_CXX_CMP_RT(>, mpfr_cmp_d, double)
    360 MPFR_CXX_CMP_RT(>=, mpfr_cmp_d, double)
    361 MPFR_CXX_CMP_RT(<=, mpfr_cmp_d, double)
    362 
    363 MPFR_CXX_CMP_TR(==, mpfr_cmp_si, int)
    364 MPFR_CXX_CMP_TR(!=, mpfr_cmp_si, int)
    365 MPFR_CXX_CMP_TR(<, mpfr_cmp_si, int)
    366 MPFR_CXX_CMP_TR(>, mpfr_cmp_si, int)
    367 MPFR_CXX_CMP_TR(>=, mpfr_cmp_si, int)
    368 MPFR_CXX_CMP_TR(<=, mpfr_cmp_si, int)
    369 MPFR_CXX_CMP_TR(==, mpfr_cmp_ui, unsigned int)
    370 MPFR_CXX_CMP_TR(!=, mpfr_cmp_ui, unsigned int)
    371 MPFR_CXX_CMP_TR(<, mpfr_cmp_ui, unsigned int)
    372 MPFR_CXX_CMP_TR(>, mpfr_cmp_ui, unsigned int)
    373 MPFR_CXX_CMP_TR(>=, mpfr_cmp_ui, unsigned int)
    374 MPFR_CXX_CMP_TR(<=, mpfr_cmp_ui, unsigned int)
    375 MPFR_CXX_CMP_TR(==, mpfr_cmp_si, long int)
    376 MPFR_CXX_CMP_TR(!=, mpfr_cmp_si, long int)
    377 MPFR_CXX_CMP_TR(<, mpfr_cmp_si, long int)
    378 MPFR_CXX_CMP_TR(>, mpfr_cmp_si, long int)
    379 MPFR_CXX_CMP_TR(>=, mpfr_cmp_si, long int)
    380 MPFR_CXX_CMP_TR(<=, mpfr_cmp_si, long int)
    381 MPFR_CXX_CMP_TR(==, mpfr_cmp_ui, unsigned long int)
    382 MPFR_CXX_CMP_TR(!=, mpfr_cmp_ui, unsigned long int)
    383 MPFR_CXX_CMP_TR(<, mpfr_cmp_ui, unsigned long int)
    384 MPFR_CXX_CMP_TR(>, mpfr_cmp_ui, unsigned long int)
    385 MPFR_CXX_CMP_TR(>=, mpfr_cmp_ui, unsigned long int)
    386 MPFR_CXX_CMP_TR(<=, mpfr_cmp_ui, unsigned long int)
    387 MPFR_CXX_CMP_TR(==, mpfr_cmp_d, double)
    388 MPFR_CXX_CMP_TR(!=, mpfr_cmp_d, double)
    389 MPFR_CXX_CMP_TR(<, mpfr_cmp_d, double)
    390 MPFR_CXX_CMP_TR(>, mpfr_cmp_d, double)
    391 MPFR_CXX_CMP_TR(>=, mpfr_cmp_d, double)
    392 MPFR_CXX_CMP_TR(<=, mpfr_cmp_d, double)
    393 
    394 MPFR_CXX_UNARY(operator-, mpfr_neg)
    395 MPFR_CXX_BINARY(operator+, mpfr_add)
    396 MPFR_CXX_BINARY(operator-, mpfr_sub)
    397 MPFR_CXX_BINARY(operator*, mpfr_mul)
    398 MPFR_CXX_BINARY(operator/, mpfr_div)
    399 MPFR_CXX_BINARY_RT(operator+, mpfr_add_si, int)
    400 MPFR_CXX_BINARY_RT(operator-, mpfr_sub_si, int)
    401 MPFR_CXX_BINARY_RT(operator*, mpfr_mul_si, int)
    402 MPFR_CXX_BINARY_RT(operator/, mpfr_div_si, int)
    403 MPFR_CXX_BINARY_RT(operator+, mpfr_add_si, long int)
    404 MPFR_CXX_BINARY_RT(operator-, mpfr_sub_si, long int)
    405 MPFR_CXX_BINARY_RT(operator*, mpfr_mul_si, long int)
    406 MPFR_CXX_BINARY_RT(operator/, mpfr_div_si, long int)
    407 MPFR_CXX_BINARY_RT(operator+, mpfr_add_ui, unsigned int)
    408 MPFR_CXX_BINARY_RT(operator-, mpfr_sub_ui, unsigned int)
    409 MPFR_CXX_BINARY_RT(operator*, mpfr_mul_ui, unsigned int)
    410 MPFR_CXX_BINARY_RT(operator/, mpfr_div_ui, unsigned int)
    411 MPFR_CXX_BINARY_RT(operator+, mpfr_add_ui, unsigned long int)
    412 MPFR_CXX_BINARY_RT(operator-, mpfr_sub_ui, unsigned long int)
    413 MPFR_CXX_BINARY_RT(operator*, mpfr_mul_ui, unsigned long int)
    414 MPFR_CXX_BINARY_RT(operator/, mpfr_div_ui, unsigned long int)
    415 MPFR_CXX_BINARY_RT(operator+, mpfr_add_d, double)
    416 MPFR_CXX_BINARY_RT(operator-, mpfr_sub_d, double)
    417 MPFR_CXX_BINARY_RT(operator*, mpfr_mul_d, double)
    418 MPFR_CXX_BINARY_RT(operator/, mpfr_div_d, double)
    419 
    420 MPFR_CXX_BINARY_TR(operator-, mpfr_si_sub, int)
    421 MPFR_CXX_BINARY_TR(operator/, mpfr_si_div, int)
    422 MPFR_CXX_BINARY_TR(operator-, mpfr_si_sub, long int)
    423 MPFR_CXX_BINARY_TR(operator/, mpfr_si_div, long int)
    424 MPFR_CXX_BINARY_TR(operator-, mpfr_ui_sub, unsigned int)
    425 MPFR_CXX_BINARY_TR(operator/, mpfr_ui_div, unsigned int)
    426 MPFR_CXX_BINARY_TR(operator-, mpfr_ui_sub, unsigned long int)
    427 MPFR_CXX_BINARY_TR(operator/, mpfr_ui_div, unsigned long int)
    428 MPFR_CXX_BINARY_TR(operator-, mpfr_d_sub, double)
    429 MPFR_CXX_BINARY_TR(operator/, mpfr_d_div, double)
    430 
    431 MPFR_CXX_INPLACE_BINARY(operator+=, mpfr_add)
    432 MPFR_CXX_INPLACE_BINARY(operator-=, mpfr_sub)
    433 MPFR_CXX_INPLACE_BINARY(operator*=, mpfr_mul)
    434 MPFR_CXX_INPLACE_BINARY(operator/=, mpfr_div)
    435 
    436 MPFR_CXX_INPLACE_BINARY_T(operator+=, mpfr_add_si, int)
    437 MPFR_CXX_INPLACE_BINARY_T(operator-=, mpfr_sub_si, int)
    438 MPFR_CXX_INPLACE_BINARY_T(operator*=, mpfr_mul_si, int)
    439 MPFR_CXX_INPLACE_BINARY_T(operator/=, mpfr_div_si, int)
    440 MPFR_CXX_INPLACE_BINARY_T(operator+=, mpfr_add_si, long int)
    441 MPFR_CXX_INPLACE_BINARY_T(operator-=, mpfr_sub_si, long int)
    442 MPFR_CXX_INPLACE_BINARY_T(operator*=, mpfr_mul_si, long int)
    443 MPFR_CXX_INPLACE_BINARY_T(operator/=, mpfr_div_si, long int)
    444 MPFR_CXX_INPLACE_BINARY_T(operator+=, mpfr_add_ui, unsigned int)
    445 MPFR_CXX_INPLACE_BINARY_T(operator-=, mpfr_sub_ui, unsigned int)
    446 MPFR_CXX_INPLACE_BINARY_T(operator*=, mpfr_mul_ui, unsigned int)
    447 MPFR_CXX_INPLACE_BINARY_T(operator/=, mpfr_div_ui, unsigned int)
    448 MPFR_CXX_INPLACE_BINARY_T(operator+=, mpfr_add_ui, unsigned long int)
    449 MPFR_CXX_INPLACE_BINARY_T(operator-=, mpfr_sub_ui, unsigned long int)
    450 MPFR_CXX_INPLACE_BINARY_T(operator*=, mpfr_mul_ui, unsigned long int)
    451 MPFR_CXX_INPLACE_BINARY_T(operator/=, mpfr_div_ui, unsigned long int)
    452 MPFR_CXX_INPLACE_BINARY_T(operator+=, mpfr_add_d, double)
    453 MPFR_CXX_INPLACE_BINARY_T(operator-=, mpfr_sub_d, double)
    454 MPFR_CXX_INPLACE_BINARY_T(operator*=, mpfr_mul_d, double)
    455 MPFR_CXX_INPLACE_BINARY_T(operator/=, mpfr_div_d, double)
    456 
    457 MPFR_CXX_UNARY(sqrt, mpfr_sqrt)
    458 MPFR_CXX_UNARY(recipsqrt, mpfr_rec_sqrt)
    459 MPFR_CXX_UNARY(cbrt, mpfr_cbrt)
    460 MPFR_CXX_UNARY(sin, mpfr_sin)
    461 MPFR_CXX_UNARY(cos, mpfr_cos)
    462 MPFR_CXX_UNARY(tan, mpfr_tan)
    463 MPFR_CXX_UNARY(cot, mpfr_cot)
    464 MPFR_CXX_UNARY(sec, mpfr_sec)
    465 MPFR_CXX_UNARY(csc, mpfr_csc)
    466 MPFR_CXX_UNARY(asin, mpfr_asin)
    467 MPFR_CXX_UNARY(acos, mpfr_acos)
    468 MPFR_CXX_UNARY(atan, mpfr_atan)
    469 MPFR_CXX_UNARY(sinh, mpfr_sinh)
    470 MPFR_CXX_UNARY(cosh, mpfr_cosh)
    471 MPFR_CXX_UNARY(tanh, mpfr_tanh)
    472 MPFR_CXX_UNARY(coth, mpfr_coth)
    473 MPFR_CXX_UNARY(asinh, mpfr_asinh)
    474 MPFR_CXX_UNARY(acosh, mpfr_acosh)
    475 MPFR_CXX_UNARY(atanh, mpfr_atanh)
    476 MPFR_CXX_UNARY(log, mpfr_log)
    477 MPFR_CXX_UNARY(log2, mpfr_log2)
    478 MPFR_CXX_UNARY(log10, mpfr_log10)
    479 MPFR_CXX_UNARY(log1p, mpfr_log1p)
    480 MPFR_CXX_UNARY(exp, mpfr_exp)
    481 MPFR_CXX_UNARY(exp2, mpfr_exp2)
    482 MPFR_CXX_UNARY(exp10, mpfr_exp10)
    483 MPFR_CXX_UNARY(expm1, mpfr_expm1)
    484 MPFR_CXX_UNARY_RND(floor, mpfr_floor, )
    485 MPFR_CXX_UNARY_RND(ceil, mpfr_ceil, )
    486 MPFR_CXX_UNARY_RND(trunc, mpfr_trunc, )
    487 MPFR_CXX_UNARY_RND(round, mpfr_round, )
    488 MPFR_CXX_UNARY(abs, mpfr_abs)
    489 MPFR_CXX_UNARY(erf, mpfr_erf)
    490 MPFR_CXX_UNARY(erfc, mpfr_erfc)
    491 MPFR_CXX_UNARY(gamma, mpfr_gamma)
    492 MPFR_CXX_UNARY(lngamma, mpfr_lngamma)
    493 MPFR_CXX_UNARY(digamma, mpfr_digamma)
    494 
    495 MPFR_CXX_BINARY(min, mpfr_min)
    496 MPFR_CXX_BINARY(max, mpfr_max)
    497 MPFR_CXX_BINARY(dim, mpfr_dim)
    498 MPFR_CXX_BINARY(hypot, mpfr_hypot)
    499 MPFR_CXX_BINARY(pow, mpfr_pow)
    500 MPFR_CXX_BINARY_RT(pow, mpfr_pow_si, int)
    501 MPFR_CXX_BINARY_RT(pow, mpfr_pow_ui, unsigned int)
    502 MPFR_CXX_BINARY_RT(pow, mpfr_pow_si, long int)
    503 MPFR_CXX_BINARY_RT(pow, mpfr_pow_ui, unsigned long int)
    504 MPFR_CXX_BINARY(atan2, mpfr_atan2)
    505 
    506 inline number sinc(const number& x)
    507 {
    508     number result;
    509     if (x.iszero())
    510     {
    511         return 1;
    512     }
    513     scoped_precision p(internal::precision() * 2);
    514     result = sin(x) / x;
    515     return result.withprec(internal::precision());
    516 }
    517 
    518 inline number logb(const number& x) { return floor(log2(x)); }
    519 
    520 inline number pow(const number& x, double y)
    521 {
    522     number result;
    523     number yy(y);
    524     mpfr_pow(result.mpfr_val(), x.mpfr_val(), yy.mpfr_val(), internal::rounding_mode());
    525     return result;
    526 }
    527 
    528 MPFR_CXX_TERNARY(fma, mpfr_fma)
    529 MPFR_CXX_TERNARY(fms, mpfr_fms)
    530 
    531 MPFR_CXX_CONST(const_pi, mpfr_const_pi)
    532 MPFR_CXX_CONST(pi, mpfr_const_pi)
    533 MPFR_CXX_CONST(const_euler, mpfr_const_euler)
    534 MPFR_CXX_CONST(const_log2, mpfr_const_log2)
    535 MPFR_CXX_CONST(const_catalan, mpfr_const_catalan)
    536 
    537 inline void sincos(number& sin, number& cos, const number& x)
    538 {
    539     mpfr_sin_cos(sin.mpfr_val(), cos.mpfr_val(), x.mpfr_val(), internal::rounding_mode());
    540 }
    541 inline void sinhcosh(number& sinh, number& cosh, const number& x)
    542 {
    543     mpfr_sinh_cosh(sinh.mpfr_val(), cosh.mpfr_val(), x.mpfr_val(), internal::rounding_mode());
    544 }
    545 
    546 /// Return Not-a-Number
    547 inline number nan()
    548 {
    549     number result;
    550     mpfr_set_nan(result.mpfr_val());
    551     return result;
    552 }
    553 /// Return Infinity (positive or negative)
    554 inline number infinity(int sign = 1)
    555 {
    556     number result;
    557     mpfr_set_inf(result.mpfr_val(), sign);
    558     return result;
    559 }
    560 /// Return zero (positive or negative)
    561 inline number zero(int sign = 1)
    562 {
    563     number result;
    564     mpfr_set_zero(result.mpfr_val(), sign);
    565     return result;
    566 }
    567 /// Return n / d
    568 inline number fraction(long n, long d)
    569 {
    570     number num(n);
    571     return num / d;
    572 }
    573 /// Return 1/x
    574 inline number reciprocal(const number& x) { return 1 / x; }
    575 
    576 inline number horner(number /*x*/) { return zero(); }
    577 inline number horner(number /*x*/, number c0) { return c0; }
    578 template <typename... Ts>
    579 inline number horner(number x, number c0, number c1, Ts... values)
    580 {
    581     return fma(horner(x, c1, values...), x, c0);
    582 }
    583 
    584 struct complex
    585 {
    586     inline complex() noexcept                          = default;
    587     inline complex(const complex&) noexcept            = default;
    588     inline complex(complex&&) noexcept                 = default;
    589     inline complex& operator=(const complex&) noexcept = default;
    590     inline complex& operator=(complex&&) noexcept      = default;
    591     inline complex(const number& real) noexcept : real(real), imag(zero()) {}
    592     inline complex(const number& real, const number& imag) noexcept : real(real), imag(imag) {}
    593     inline bool isreal() const { return imag.iszero(); }
    594     inline bool iszero() const { return real.iszero() && imag.iszero(); }
    595     inline bool isnan() const { return real.isnan() || imag.isnan(); }
    596     inline bool isinfinity() const { return real.isinfinity() || imag.isinfinity(); }
    597     inline bool isfinite() const { return real.isfinite() && imag.isfinite(); }
    598     /// Converts precision
    599     inline complex withprec(mpfr_prec_t prec) const { return { real.withprec(prec), imag.withprec(prec) }; }
    600     number real;
    601     number imag;
    602 };
    603 
    604 inline complex operator+(const complex& x, const complex& y) { return { x.real + y.real, x.imag + y.imag }; }
    605 inline complex operator+(const complex& x, const number& y) { return { x.real + y, x.imag }; }
    606 inline complex operator+(const number& x, const complex& y) { return y + x; }
    607 inline complex operator-(const complex& x, const complex& y) { return { x.real - y.real, x.imag - y.imag }; }
    608 inline complex operator-(const complex& x, const number& y) { return { x.real - y, x.imag }; }
    609 inline complex operator-(const number& x, const complex& y) { return { x - y.real, -y.imag }; }
    610 inline complex operator*(const complex& x, const complex& y)
    611 {
    612     return { x.real * y.real - x.imag * y.imag, x.real * y.imag + x.imag * y.real };
    613 }
    614 inline complex operator*(const complex& x, const number& y) { return { x.real * y, x.imag * y }; }
    615 inline complex operator*(const number& x, const complex& y) { return y * x; }
    616 
    617 inline number abs(const complex& x) { return hypot(x.real, x.imag); }
    618 inline number phase(const complex& x) { return atan2(x.imag, x.real); }
    619 inline complex conj(const complex& x) { return { x.real, -x.imag }; }
    620 inline complex reciprocal(const complex& x)
    621 {
    622     number a = abs(x);
    623     return { x.real / a, -x.imag / a };
    624 }
    625 inline complex operator/(const complex& x, const complex& y) { return x * reciprocal(y); }
    626 inline complex operator/(const complex& x, const number& y) { return x * reciprocal(y); }
    627 inline complex operator/(const number& x, const complex& y) { return x * reciprocal(y); }
    628 
    629 inline complex exp(const complex& x)
    630 {
    631     complex sc;
    632     sincos(sc.imag, sc.real, x.imag);
    633     number e = exp(x.real);
    634     return { e * sc.imag, e * sc.imag };
    635 }
    636 inline complex sin(const complex& x)
    637 {
    638     complex sc;
    639     complex sch;
    640     sincos(sc.imag, sc.real, x.real);
    641     sinhcosh(sch.imag, sch.real, x.imag);
    642     return { sc.imag * sch.real, sc.real * sch.imag };
    643 }
    644 inline complex cos(const complex& x)
    645 {
    646     complex sc;
    647     complex sch;
    648     sincos(sc.imag, sc.real, x.real);
    649     sinhcosh(sch.imag, sch.real, x.imag);
    650     return { sc.real * sch.real, -sc.imag * sch.imag };
    651 }
    652 inline complex tan(const complex& x) { return sin(x) / cos(x); }
    653 
    654 inline complex polar(const complex& x) { return { hypot(x.real, x.imag), atan2(x.imag, x.real) }; }
    655 inline complex cartesian(const complex& x)
    656 {
    657     number c, s;
    658     sincos(s, c, x.imag);
    659     return { x.real * c, x.real * s };
    660 }
    661 
    662 MPFR_FN(sqrt)
    663 MPFR_FN(recipsqrt)
    664 MPFR_FN(cbrt)
    665 MPFR_FN(sin)
    666 MPFR_FN(sinc)
    667 MPFR_FN(cos)
    668 MPFR_FN(tan)
    669 MPFR_FN(cot)
    670 MPFR_FN(sec)
    671 MPFR_FN(csc)
    672 MPFR_FN(asin)
    673 MPFR_FN(acos)
    674 MPFR_FN(atan)
    675 MPFR_FN(sinh)
    676 MPFR_FN(cosh)
    677 MPFR_FN(tanh)
    678 MPFR_FN(coth)
    679 MPFR_FN(asinh)
    680 MPFR_FN(acosh)
    681 MPFR_FN(atanh)
    682 MPFR_FN(log)
    683 MPFR_FN(log2)
    684 MPFR_FN(log10)
    685 MPFR_FN(log1p)
    686 MPFR_FN(logb)
    687 MPFR_FN(exp)
    688 MPFR_FN(exp2)
    689 MPFR_FN(exp10)
    690 MPFR_FN(expm1)
    691 MPFR_FN(floor)
    692 MPFR_FN(ceil)
    693 MPFR_FN(trunc)
    694 MPFR_FN(round)
    695 MPFR_FN(abs)
    696 MPFR_FN(erf)
    697 MPFR_FN(erfc)
    698 MPFR_FN(gamma)
    699 MPFR_FN(lngamma)
    700 MPFR_FN(digamma)
    701 MPFR_FN(min)
    702 MPFR_FN(max)
    703 MPFR_FN(dim)
    704 MPFR_FN(hypot)
    705 MPFR_FN(pow)
    706 MPFR_FN(atan2)
    707 MPFR_FN(fma)
    708 MPFR_FN(fms)
    709 MPFR_FN(horner)
    710 MPFR_FN(sincos)
    711 MPFR_FN(sinhcosh)
    712 MPFR_FN(nan)
    713 MPFR_FN(infinity)
    714 MPFR_FN(zero)
    715 MPFR_FN(fraction)
    716 MPFR_FN(reciprocal)
    717 } // namespace mpfr
    718 
    719 #undef MPFR_CXX_UNARY
    720 #undef MPFR_CXX_UNARY_RND
    721 #undef MPFR_CXX_BINARY
    722 #undef MPFR_CXX_TERNARY
    723 #undef MPFR_CXX_BINARY_RT
    724 #undef MPFR_CXX_BINARY_TR
    725 #undef MPFR_CXX_INPLACE_BINARY
    726 #undef MPFR_CXX_INPLACE_BINARY_T
    727 #undef MPFR_CXX_CTOR_T
    728 #undef MPFR_CXX_ASGN_T