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