diff --git a/Eigen/src/Geometry/arch/Geometry_SSE.h b/Eigen/src/Geometry/arch/Geometry_SSE.h index 080b87d4b..0078e4aab 100644 --- a/Eigen/src/Geometry/arch/Geometry_SSE.h +++ b/Eigen/src/Geometry/arch/Geometry_SSE.h @@ -64,4 +64,58 @@ struct ei_cross3_impl } }; + + + +template +struct ei_quat_product +{ + inline static Quaternion run(const QuaternionBase& _a, const QuaternionBase& _b) + { + const Packet2d mask = _mm_castsi128_pd(_mm_set_epi32(0x0,0x0,0x80000000,0x0)); + + Quaternion res; + + const double* a = _a.coeffs().data(); + Packet2d b_xy = _b.coeffs().template packet(0); + Packet2d b_zw = _b.coeffs().template packet(2); + Packet2d a_xx = ei_pset1(a[0]); + Packet2d a_yy = ei_pset1(a[1]); + Packet2d a_zz = ei_pset1(a[2]); + Packet2d a_ww = ei_pset1(a[3]); + + // two temporaries: + Packet2d t1, t2; + + /* + * t1 = ww*xy + yy*zw + * t2 = zz*xy - xx*zw + * res.xy = t1 +/- swap(t2) + */ + t1 = ei_padd(ei_pmul(a_ww, b_xy), ei_pmul(a_yy, b_zw)); + t2 = ei_psub(ei_pmul(a_zz, b_xy), ei_pmul(a_xx, b_zw)); +#ifdef __SSE3__ + ei_pstore(&res.x(), _mm_addsub_pd(t1, ei_preverse(t2))); +#else + ei_pstore(&res.x(), ei_padd(t1, ei_pxor(mask,ei_preverse(t2)))); +#endif + + /* + * t1 = ww*zw - yy*xy + * t2 = zz*zw + xx*xy + * res.zw = t1 -/+ swap(t2) = swap( swap(t1) +/- t2) + */ + t1 = ei_psub(ei_pmul(a_ww, b_zw), ei_pmul(a_yy, b_xy)); + t2 = ei_padd(ei_pmul(a_zz, b_zw), ei_pmul(a_xx, b_xy)); +#ifdef __SSE3__ + ei_pstore(&res.z(), ei_preverse(_mm_addsub_pd(ei_preverse(t1), t2))); +#else + ei_pstore(&res.z(), ei_psub(t1, ei_pxor(mask,ei_preverse(t2)))); +#endif + + return res; +} +}; + + #endif // EIGEN_GEOMETRY_SSE_H diff --git a/bench/quatmul.cpp b/bench/quatmul.cpp new file mode 100644 index 000000000..d91a4b01b --- /dev/null +++ b/bench/quatmul.cpp @@ -0,0 +1,47 @@ +#include +#include +#include +#include + +using namespace Eigen; + +template +EIGEN_DONT_INLINE void quatmul_default(const Quat& a, const Quat& b, Quat& c) +{ + c = a * b; +} + +template +EIGEN_DONT_INLINE void quatmul_novec(const Quat& a, const Quat& b, Quat& c) +{ + c = ei_quat_product<0, Quat, Quat, typename Quat::Scalar, Aligned>::run(a,b); +} + +template void bench(const std::string& label) +{ + int tries = 10; + int rep = 1000000; + BenchTimer t; + + Quat a(4, 1, 2, 3); + Quat b(2, 3, 4, 5); + Quat c; + + std::cout.precision(3); + + BENCH(t, tries, rep, quatmul_default(a,b,c)); + std::cout << label << " default " << 1e3*t.best(CPU_TIMER) << "ms \t" << 1e-6*double(rep)/(t.best(CPU_TIMER)) << " M mul/s\n"; + + BENCH(t, tries, rep, quatmul_novec(a,b,c)); + std::cout << label << " novec " << 1e3*t.best(CPU_TIMER) << "ms \t" << 1e-6*double(rep)/(t.best(CPU_TIMER)) << " M mul/s\n"; +} + +int main() +{ + bench("float "); + bench("double"); + + return 0; + +} + diff --git a/cmake/FindGMP.cmake b/cmake/FindGMP.cmake new file mode 100644 index 000000000..1f0273960 --- /dev/null +++ b/cmake/FindGMP.cmake @@ -0,0 +1,21 @@ +# Try to find the GNU Multiple Precision Arithmetic Library (GMP) +# See http://gmplib.org/ + +if (GMP_INCLUDES AND GMP_LIBRARIES) + set(GMP_FIND_QUIETLY TRUE) +endif (GMP_INCLUDES AND GMP_LIBRARIES) + +find_path(GMP_INCLUDES + NAMES + gmp.h + PATHS + $ENV{GMPDIR} + ${INCLUDE_INSTALL_DIR} +) + +find_library(GMP_LIBRARIES gmp PATHS $ENV{GMPDIR} ${LIB_INSTALL_DIR}) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(GMP DEFAULT_MSG + GMP_INCLUDES GMP_LIBRARIES) +mark_as_advanced(GMP_INCLUDES GMP_LIBRARIES) diff --git a/cmake/FindMPFR.cmake b/cmake/FindMPFR.cmake new file mode 100644 index 000000000..3aa97f0b4 --- /dev/null +++ b/cmake/FindMPFR.cmake @@ -0,0 +1,21 @@ +# Try to find the MPFR library +# See http://www.mpfr.org/ + +if (MPFR_INCLUDES AND MPFR_LIBRARIES) + set(MPFR_FIND_QUIETLY TRUE) +endif (MPFR_INCLUDES AND MPFR_LIBRARIES) + +find_path(MPFR_INCLUDES + NAMES + mpfr.h + PATHS + $ENV{GMPDIR} + ${INCLUDE_INSTALL_DIR} +) + +find_library(MPFR_LIBRARIES mpfr PATHS $ENV{GMPDIR} ${LIB_INSTALL_DIR}) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(MPFR DEFAULT_MSG + MPFR_INCLUDES MPFR_LIBRARIES) +mark_as_advanced(MPFR_INCLUDES MPFR_LIBRARIES) diff --git a/doc/unsupported_modules.dox b/doc/unsupported_modules.dox index dff74ca1b..9f0517c4c 100644 --- a/doc/unsupported_modules.dox +++ b/doc/unsupported_modules.dox @@ -61,4 +61,7 @@ namespace Eigen { /** \ingroup Unsupported_modules * \defgroup SparseExtra_Module */ +/** \ingroup Unsupported_modules + * \defgroup MpfrcxxSupport_Module */ + } // namespace Eigen diff --git a/unsupported/Eigen/MPRealSupport b/unsupported/Eigen/MPRealSupport new file mode 100644 index 000000000..abc719ba7 --- /dev/null +++ b/unsupported/Eigen/MPRealSupport @@ -0,0 +1,160 @@ +// This file is part of a joint effort between Eigen, a lightweight C++ template library +// for linear algebra, and MPFR C++, a C++ interface to MPFR library (http://www.holoborodko.com/pavel/) +// +// Copyright (C) 2010 Pavel Holoborodko +// Copyright (C) 2010 Konstantin Holoborodko +// Copyright (C) 2010 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// this library. If not, see . +// +// Contributors: +// Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen, Pere Constans + +#ifndef EIGEN_MPREALSUPPORT_MODULE_H +#define EIGEN_MPREALSUPPORT_MODULE_H + +#include +#include + +namespace Eigen { + + /** \defgroup MPRealSupport_Module MPFRC++ Support module + * + * \code + * #include + * \endcode + * + * This module provides support for multi precision floating point numbers + * via the MPFR C++ + * library which itself is built upon MPFR/GMP. + * + * Here is an example: + * +\code +#include +#include +#include +using namespace mpfr; +using namespace Eigen; +int main() +{ + // set precision to 256 bits (double has only 53 bits) + mpreal::set_default_prec(256); + // Declare matrix and vector types with multi-precision scalar type + typedef Matrix MatrixXmp; + typedef Matrix VectorXmp; + + MatrixXmp A = MatrixXmp::Random(100,100); + VectorXmp b = VectorXmp::Random(100); + + // Solve Ax=b using LU + VectorXmp x = A.lu().solve(b); + std::cout << "relative error: " << (A*x - b).norm() / b.norm() << std::endl; + return 0; +} +\endcode + * + */ + + template<> struct NumTraits + : GenericNumTraits + { + enum { + IsInteger = 0, + IsSigned = 1, + IsComplex = 0, + ReadCost = 10, + AddCost = 10, + MulCost = 40 + }; + + typedef mpfr::mpreal Real; + typedef mpfr::mpreal NonInteger; + + inline static mpfr::mpreal highest() { return mpfr::mpreal_max(mpfr::mpreal::get_default_prec()); } + inline static mpfr::mpreal lowest() { return -mpfr::mpreal_max(mpfr::mpreal::get_default_prec()); } + + inline static Real epsilon() + { + return mpfr::machine_epsilon(mpfr::mpreal::get_default_prec()); + } + inline static Real dummy_precision() + { + unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1)*90)/100; + return mpfr::machine_epsilon(weak_prec); + } + }; + + template<> mpfr::mpreal ei_random() + { +#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) + static gmp_randstate_t state; + static bool isFirstTime = true; + + if(isFirstTime) + { + gmp_randinit_default(state); + gmp_randseed_ui(state,(unsigned)time(NULL)); + isFirstTime = false; + } + + return mpfr::urandom(state)*2-1; +#else + return mpfr::mpreal(ei_random()); +#endif + } + + template<> mpfr::mpreal ei_random(const mpfr::mpreal& a, const mpfr::mpreal& b) + { + return a + (b-a) * ei_random(); + } +} + +namespace mpfr { + + inline const mpreal& ei_conj(const mpreal& x) { return x; } + inline const mpreal& ei_real(const mpreal& x) { return x; } + inline mpreal ei_imag(const mpreal&) { return 0.0; } + inline mpreal ei_abs(const mpreal& x) { return fabs(x); } + inline mpreal ei_abs2(const mpreal& x) { return x*x; } + inline mpreal ei_sqrt(const mpreal& x) { return sqrt(x); } + inline mpreal ei_exp(const mpreal& x) { return exp(x); } + inline mpreal ei_log(const mpreal& x) { return log(x); } + inline mpreal ei_sin(const mpreal& x) { return sin(x); } + inline mpreal ei_cos(const mpreal& x) { return cos(x); } + inline mpreal ei_pow(const mpreal& x, mpreal& y) { return pow(x, y); } + + bool ei_isMuchSmallerThan(const mpreal& a, const mpreal& b, const mpreal& prec) + { + return ei_abs(a) <= abs(b) * prec; + } + + inline bool ei_isApprox(const mpreal& a, const mpreal& b, const mpreal& prec) + { + return ei_abs(a - b) <= min(abs(a), abs(b)) * prec; + } + + inline bool ei_isApproxOrLessThan(const mpreal& a, const mpreal& b, const mpreal& prec) + { + return a <= b || ei_isApprox(a, b, prec); + } +} + +#endif // EIGEN_MPREALSUPPORT_MODULE_H diff --git a/unsupported/doc/examples/CMakeLists.txt b/unsupported/doc/examples/CMakeLists.txt index 61af3e5a5..978f9afd8 100644 --- a/unsupported/doc/examples/CMakeLists.txt +++ b/unsupported/doc/examples/CMakeLists.txt @@ -2,6 +2,8 @@ FILE(GLOB examples_SRCS "*.cpp") ADD_CUSTOM_TARGET(unsupported_examples) +INCLUDE_DIRECTORIES(../../../unsupported ../../../unsupported/test) + FOREACH(example_src ${examples_SRCS}) GET_FILENAME_COMPONENT(example ${example_src} NAME_WE) ADD_EXECUTABLE(example_${example} ${example_src}) diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index c7cfc8881..db051caec 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -55,10 +55,10 @@ include_directories(../../test ../../unsupported ../../Eigen) if(ADOLC_FOUND) include_directories(${ADOLC_INCLUDES}) - ei_add_property(EIGEN_TESTED_BACKENDS "Adolc") + ei_add_property(EIGEN_TESTED_BACKENDS "Adolc, ") ei_add_test(forward_adolc " " ${ADOLC_LIBRARIES}) else(ADOLC_FOUND) - ei_add_property(EIGEN_MISSING_BACKENDS "Adolc") + ei_add_property(EIGEN_MISSING_BACKENDS "Adolc, ") endif(ADOLC_FOUND) ei_add_test(NonLinearOptimization) @@ -70,6 +70,15 @@ ei_add_test(matrix_function) ei_add_test(alignedvector3) ei_add_test(FFT) +find_package(MPFR) +if(MPFR_FOUND) + include_directories(${MPFR_INCLUDES}) + ei_add_property(EIGEN_TESTED_BACKENDS "MPFR C++, ") + ei_add_test(mpreal_support " " ${MPFR_LIBRARIES} ) +else() + ei_add_property(EIGEN_MISSING_BACKENDS "MPFR C++, ") +endif() + ei_add_test(sparse_llt " " "${SPARSE_LIBS}") ei_add_test(sparse_ldlt " " "${SPARSE_LIBS}") ei_add_test(sparse_lu " " "${SPARSE_LIBS}") diff --git a/unsupported/test/mpreal.cpp b/unsupported/test/mpreal.cpp new file mode 100644 index 000000000..bb34c246d --- /dev/null +++ b/unsupported/test/mpreal.cpp @@ -0,0 +1,408 @@ +/* + Multi-precision real number class. C++ wrapper fo MPFR library. + Project homepage: http://www.holoborodko.com/pavel/ + Contact e-mail: pavel@holoborodko.com + + Copyright (c) 2008-2010 Pavel Holoborodko + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + + Contributors: + Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, + Heinz van Saanen, Pere Constans, Dmitriy Gubanov +*/ + +#include +#include "mpreal.h" + +using std::ws; +using std::cerr; +using std::endl; +using std::string; +using std::ostream; +using std::istream; + +namespace mpfr{ + +mp_rnd_t mpreal::default_rnd = mpfr_get_default_rounding_mode(); +mp_prec_t mpreal::default_prec = mpfr_get_default_prec(); +int mpreal::default_base = 10; +int mpreal::double_bits = -1; + +// Default constructor: creates mp number and initializes it to 0. +mpreal::mpreal() +{ + mpfr_init2(mp,default_prec); + mpfr_set_ui(mp,0,default_rnd); +} + +mpreal::mpreal(const mpreal& u) +{ + mpfr_init2(mp,mpfr_get_prec(u.mp)); + mpfr_set(mp,u.mp,default_rnd); +} + +mpreal::mpreal(const mpfr_t u) +{ + mpfr_init2(mp,mpfr_get_prec(u)); + mpfr_set(mp,u,default_rnd); +} + +mpreal::mpreal(const mpf_t u) +{ + mpfr_init2(mp,mpf_get_prec(u)); + mpfr_set_f(mp,u,default_rnd); +} + +mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode) +{ + mpfr_init2(mp,prec); + mpfr_set_z(mp,u,mode); +} + +mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode) +{ + mpfr_init2(mp,prec); + mpfr_set_q(mp,u,mode); +} + +mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) +{ + if(double_bits == -1 || fits_in_bits(u, double_bits)) + { + mpfr_init2(mp,prec); + mpfr_set_d(mp,u,mode); + } + else + throw conversion_overflow(); +} + +mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode) +{ + mpfr_init2(mp,prec); + mpfr_set_ld(mp,u,mode); +} + +mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) +{ + mpfr_init2(mp,prec); + mpfr_set_ui(mp,u,mode); +} + +mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) +{ + mpfr_init2(mp,prec); + mpfr_set_ui(mp,u,mode); +} + +mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) +{ + mpfr_init2(mp,prec); + mpfr_set_si(mp,u,mode); +} + +mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) +{ + mpfr_init2(mp,prec); + mpfr_set_si(mp,u,mode); +} + +mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) +{ + mpfr_init2(mp,prec); + mpfr_set_str(mp, s, base, mode); +} + +mpreal::~mpreal() +{ + mpfr_clear(mp); +} + +// Operators - Assignment +mpreal& mpreal::operator=(const char* s) +{ + mpfr_t t; + if(0==mpfr_init_set_str(t,s,default_base,default_rnd)) + { + mpfr_set(mp,t,mpreal::default_rnd); + mpfr_clear(t); + }else{ + mpfr_clear(t); + // cerr<<"fail to convert string"<p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); + + mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode); + return a; +} + +const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode) +{ + mpreal a; + mp_prec_t p1, p2, p3; + + p1 = v1.get_prec(); + p2 = v2.get_prec(); + p3 = v3.get_prec(); + + a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); + + mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode); + return a; +} + +const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode) +{ + mpreal a; + mp_prec_t p1, p2; + + p1 = v1.get_prec(); + p2 = v2.get_prec(); + + a.set_prec(p1>p2?p1:p2); + + mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode); + + return a; +} + +const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) +{ + mpreal a; + mp_prec_t yp, xp; + + yp = y.get_prec(); + xp = x.get_prec(); + + a.set_prec(yp>xp?yp:xp); + + mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode); + + return a; +} + +const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode) +{ + mpreal x; + mpfr_ptr* t; + unsigned long int i; + + t = new mpfr_ptr[n]; + for (i=0;ixp?yp:xp); + + mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode); + + return a; +} + +const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) +{ + mpreal a; + mp_prec_t yp, xp; + + yp = y.get_prec(); + xp = x.get_prec(); + + a.set_prec(yp>xp?yp:xp); + + mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode); + + return a; +} + +template +std::string to_string(T t, std::ios_base & (*f)(std::ios_base&)) +{ + std::ostringstream oss; + oss << f << t; + return oss.str(); +} + +mpreal::operator std::string() const +{ + return to_string(); +} + +string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const +{ + char *s, *ns = NULL; + size_t slen, nslen; + mp_exp_t exp; + string out; + + if(mpfr_inf_p(mp)) + { + if(mpfr_sgn(mp)>0) return "+@Inf@"; + else return "-@Inf@"; + } + + if(mpfr_zero_p(mp)) return "0"; + if(mpfr_nan_p(mp)) return "@NaN@"; + + + s = mpfr_get_str(NULL,&exp,b,0,mp,mode); + ns = mpfr_get_str(NULL,&exp,b,n,mp,mode); + + if(s!=NULL && ns!=NULL) + { + slen = strlen(s); + nslen = strlen(ns); + if(nslen<=slen) + { + mpfr_free_str(s); + s = ns; + slen = nslen; + } + else { + mpfr_free_str(ns); + } + + // Make human eye-friendly formatting if possible + if (exp>0 && static_cast(exp)s+exp) ptr--; + + if(ptr==s+exp) out = string(s,exp+1); + else out = string(s,exp+1)+'.'+string(s+exp+1,ptr-(s+exp+1)+1); + + //out = string(s,exp+1)+'.'+string(s+exp+1); + } + else + { + // Remove zeros starting from right end + char* ptr = s+slen-1; + while (*ptr=='0' && ptr>s+exp-1) ptr--; + + if(ptr==s+exp-1) out = string(s,exp); + else out = string(s,exp)+'.'+string(s+exp,ptr-(s+exp)+1); + + //out = string(s,exp)+'.'+string(s+exp); + } + + }else{ // exp<0 || exp>slen + if(s[0]=='-') + { + // Remove zeros starting from right end + char* ptr = s+slen-1; + while (*ptr=='0' && ptr>s+1) ptr--; + + if(ptr==s+1) out = string(s,2); + else out = string(s,2)+'.'+string(s+2,ptr-(s+2)+1); + + //out = string(s,2)+'.'+string(s+2); + } + else + { + // Remove zeros starting from right end + char* ptr = s+slen-1; + while (*ptr=='0' && ptr>s) ptr--; + + if(ptr==s) out = string(s,1); + else out = string(s,1)+'.'+string(s+1,ptr-(s+1)+1); + + //out = string(s,1)+'.'+string(s+1); + } + + // Make final string + if(--exp) + { + if(exp>0) out += "e+"+mpfr::to_string(exp,std::dec); + else out += "e"+mpfr::to_string(exp,std::dec); + } + } + + mpfr_free_str(s); + return out; + }else{ + return "conversion error!"; + } +} + +////////////////////////////////////////////////////////////////////////// +// I/O +ostream& operator<<(ostream& os, const mpreal& v) +{ + return os<>(istream &is, mpreal& v) +{ + char c; + string s = ""; + mpfr_t t; + + if(is.good()) + { + is>>ws; + while ((c = is.get())!=EOF) + { + if(c ==' ' || c == '\t' || c == '\n' || c == '\r') + { + is.putback(c); + break; + } + s += c; + } + + if(s.size() != 0) + { + // Protect current value from alternation in case of input error + // so some error handling(roll back) procedure can be used + if(0==mpfr_init_set_str(t,s.c_str(),mpreal::default_base,mpreal::default_rnd)) + { + mpfr_set(v.mp,t,mpreal::default_rnd); + mpfr_clear(t); + + }else{ + mpfr_clear(t); + cerr<<"error reading from istream"< +#include +#include +#include +#include +#include + +#include + +// Detect compiler using signatures from http://predef.sourceforge.net/ +#if defined(__GNUC__) && defined(__INTEL_COMPILER) + #define IsInf(x) isinf(x) // GNU C/C++ + Intel ICC compiler + +#elif defined(__GNUC__) + #define IsInf(x) std::isinf(x) // GNU C/C++ + +#elif defined(_MSC_VER) + #define IsInf(x) (!_finite(x)) // Microsoft Visual C++ + +#else + #define IsInf(x) std::isinf(x) // C99 conformance +#endif + +namespace mpfr { + +class mpreal { +private: + mpfr_t mp; + +public: + static mp_rnd_t default_rnd; + static mp_prec_t default_prec; + static int default_base; + static int double_bits; + +public: + // Constructors && type conversion + mpreal(); + mpreal(const mpreal& u); + + mpreal(const mpfr_t u); + mpreal(const mpf_t u); + + mpreal(const mpz_t u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const mpq_t u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const double u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const long double u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const unsigned long int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const unsigned int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const long int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const char* s, mp_prec_t prec = default_prec, int base = default_base, mp_rnd_t mode = default_rnd); + + ~mpreal(); + + // Operations + // = + // +, -, *, /, ++, --, <<, >> + // *=, +=, -=, /=, + // <, >, ==, <=, >= + + // = + mpreal& operator=(const mpreal& v); + mpreal& operator=(const mpf_t v); + mpreal& operator=(const mpz_t v); + mpreal& operator=(const mpq_t v); + mpreal& operator=(const long double v); + mpreal& operator=(const double v); + mpreal& operator=(const unsigned long int v); + mpreal& operator=(const unsigned int v); + mpreal& operator=(const long int v); + mpreal& operator=(const int v); + mpreal& operator=(const char* s); + + // + + mpreal& operator+=(const mpreal& v); + mpreal& operator+=(const mpf_t v); + mpreal& operator+=(const mpz_t v); + mpreal& operator+=(const mpq_t v); + mpreal& operator+=(const long double u); + mpreal& operator+=(const double u); + mpreal& operator+=(const unsigned long int u); + mpreal& operator+=(const unsigned int u); + mpreal& operator+=(const long int u); + mpreal& operator+=(const int u); + const mpreal operator+() const; + mpreal& operator++ (); + const mpreal operator++ (int); + + // - + mpreal& operator-=(const mpreal& v); + mpreal& operator-=(const mpz_t v); + mpreal& operator-=(const mpq_t v); + mpreal& operator-=(const long double u); + mpreal& operator-=(const double u); + mpreal& operator-=(const unsigned long int u); + mpreal& operator-=(const unsigned int u); + mpreal& operator-=(const long int u); + mpreal& operator-=(const int u); + const mpreal operator-() const; + friend const mpreal operator-(const unsigned long int b, const mpreal& a); + friend const mpreal operator-(const unsigned int b, const mpreal& a); + friend const mpreal operator-(const long int b, const mpreal& a); + friend const mpreal operator-(const int b, const mpreal& a); + friend const mpreal operator-(const double b, const mpreal& a); + mpreal& operator-- (); + const mpreal operator-- (int); + + // * + mpreal& operator*=(const mpreal& v); + mpreal& operator*=(const mpz_t v); + mpreal& operator*=(const mpq_t v); + mpreal& operator*=(const long double v); + mpreal& operator*=(const double v); + mpreal& operator*=(const unsigned long int v); + mpreal& operator*=(const unsigned int v); + mpreal& operator*=(const long int v); + mpreal& operator*=(const int v); + + // / + mpreal& operator/=(const mpreal& v); + mpreal& operator/=(const mpz_t v); + mpreal& operator/=(const mpq_t v); + mpreal& operator/=(const long double v); + mpreal& operator/=(const double v); + mpreal& operator/=(const unsigned long int v); + mpreal& operator/=(const unsigned int v); + mpreal& operator/=(const long int v); + mpreal& operator/=(const int v); + friend const mpreal operator/(const unsigned long int b, const mpreal& a); + friend const mpreal operator/(const unsigned int b, const mpreal& a); + friend const mpreal operator/(const long int b, const mpreal& a); + friend const mpreal operator/(const int b, const mpreal& a); + friend const mpreal operator/(const double b, const mpreal& a); + + //<<= Fast Multiplication by 2^u + mpreal& operator<<=(const unsigned long int u); + mpreal& operator<<=(const unsigned int u); + mpreal& operator<<=(const long int u); + mpreal& operator<<=(const int u); + + //>>= Fast Division by 2^u + mpreal& operator>>=(const unsigned long int u); + mpreal& operator>>=(const unsigned int u); + mpreal& operator>>=(const long int u); + mpreal& operator>>=(const int u); + + // Boolean Operators + friend bool operator > (const mpreal& a, const mpreal& b); + friend bool operator >= (const mpreal& a, const mpreal& b); + friend bool operator < (const mpreal& a, const mpreal& b); + friend bool operator <= (const mpreal& a, const mpreal& b); + friend bool operator == (const mpreal& a, const mpreal& b); + friend bool operator != (const mpreal& a, const mpreal& b); + + // Type Conversion operators + inline operator long double() const; + inline operator double() const; + inline operator float() const; + inline operator unsigned long() const; + inline operator unsigned int() const; + inline operator long() const; + operator std::string() const; + inline operator mpfr_ptr(); + + // Math Functions + friend const mpreal sqr(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend int cmpabs(const mpreal& a,const mpreal& b); + + friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + + friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + + friend const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + + friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal fac_ui (unsigned long int v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal log1p (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal expm1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + + friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _j0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _j1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _jn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _y0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _y1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _yn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend int sgn(const mpreal& v); // -1 or +1 + +// MPFR 2.4.0 Specifics +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + friend int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal li2(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); +#endif + +// MPFR 3.0.0 Specifics +#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) + friend const mpreal digamma(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal urandom (gmp_randstate_t& state,mp_rnd_t rnd_mode = mpreal::default_rnd); // use gmp_randinit_default() to init state, gmp_randclear() to clear + friend bool _isregular(const mpreal& v); +#endif + + // Exponent and mantissa manipulation + friend const mpreal frexp(const mpreal& v, mp_exp_t* exp); + friend const mpreal ldexp(const mpreal& v, mp_exp_t exp); + + // Splits mpreal value into fractional and integer parts. + // Returns fractional part and stores integer part in n. + friend const mpreal modf(const mpreal& v, mpreal& n); + + // Constants + // don't forget to call mpfr_free_cache() for every thread where you are using const-functions + friend const mpreal const_log2 (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal const_pi (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal const_euler (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal const_catalan (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); + // returns +inf iff sign>=0 otherwise -inf + friend const mpreal const_infinity(int sign = 1, mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); + + // Output/ Input + friend std::ostream& operator<<(std::ostream& os, const mpreal& v); + friend std::istream& operator>>(std::istream& is, mpreal& v); + + // Integer Related Functions + friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal ceil (const mpreal& v); + friend const mpreal floor(const mpreal& v); + friend const mpreal round(const mpreal& v); + friend const mpreal trunc(const mpreal& v); + friend const mpreal rint_ceil (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal rint_floor(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal rint_round(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal rint_trunc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); + + // Miscellaneous Functions + friend const mpreal nexttoward (const mpreal& x, const mpreal& y); + friend const mpreal nextabove (const mpreal& x); + friend const mpreal nextbelow (const mpreal& x); + + // use gmp_randinit_default() to init state, gmp_randclear() to clear + friend const mpreal urandomb (gmp_randstate_t& state); + +// MPFR < 2.4.2 Specifics +#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2)) + friend const mpreal random2 (mp_size_t size, mp_exp_t exp); +#endif + + // Instance Checkers + friend bool _isnan(const mpreal& v); + friend bool _isinf(const mpreal& v); + friend bool _isnum(const mpreal& v); + friend bool _iszero(const mpreal& v); + friend bool _isint(const mpreal& v); + + // Set/Get instance properties + inline mp_prec_t get_prec() const; + inline void set_prec(mp_prec_t prec, mp_rnd_t rnd_mode = default_rnd); // Change precision with rounding mode + + // Set mpreal to +-inf, NaN + void set_inf(int sign = +1); + void set_nan(); + + // sign = -1 or +1 + void set_sign(int sign, mp_rnd_t rnd_mode = default_rnd); + + //Exponent + mp_exp_t get_exp(); + int set_exp(mp_exp_t e); + int check_range (int t, mp_rnd_t rnd_mode = default_rnd); + int subnormalize (int t,mp_rnd_t rnd_mode = default_rnd); + + // Inexact conversion from float + inline bool fits_in_bits(double x, int n); + + // Set/Get global properties + static void set_default_prec(mp_prec_t prec); + static mp_prec_t get_default_prec(); + static void set_default_base(int base); + static int get_default_base(); + static void set_double_bits(int dbits); + static int get_double_bits(); + static void set_default_rnd(mp_rnd_t rnd_mode); + static mp_rnd_t get_default_rnd(); + static mp_exp_t get_emin (void); + static mp_exp_t get_emax (void); + static mp_exp_t get_emin_min (void); + static mp_exp_t get_emin_max (void); + static mp_exp_t get_emax_min (void); + static mp_exp_t get_emax_max (void); + static int set_emin (mp_exp_t exp); + static int set_emax (mp_exp_t exp); + + // Get/Set conversions + // Convert mpreal to string with n significant digits in base b + // n = 0 -> convert with the maximum available digits + std::string to_string(size_t n = 0, int b = default_base, mp_rnd_t mode = default_rnd) const; + + // Efficient swapping of two mpreal values + friend void swap(mpreal& x, mpreal& y); + + //Min Max - macros is evil. Needed for systems which defines max and min globally as macros (e.g. Windows) + //Hope that globally defined macros use > < operations only + #ifndef max + friend const mpreal max(const mpreal& x, const mpreal& y); + #endif + + #ifndef min + friend const mpreal min(const mpreal& x, const mpreal& y); + #endif +}; + +////////////////////////////////////////////////////////////////////////// +// Exceptions +class conversion_overflow : public std::exception { +public: + std::string why() { return "inexact conversion from floating point"; } +}; + +////////////////////////////////////////////////////////////////////////// +// + Addition +const mpreal operator+(const mpreal& a, const mpreal& b); + +// + Fast specialized addition - implemented through fast += operations +const mpreal operator+(const mpreal& a, const mpz_t b); +const mpreal operator+(const mpreal& a, const mpq_t b); +const mpreal operator+(const mpreal& a, const long double b); +const mpreal operator+(const mpreal& a, const double b); +const mpreal operator+(const mpreal& a, const unsigned long int b); +const mpreal operator+(const mpreal& a, const unsigned int b); +const mpreal operator+(const mpreal& a, const long int b); +const mpreal operator+(const mpreal& a, const int b); +const mpreal operator+(const mpreal& a, const char* b); +const mpreal operator+(const char* a, const mpreal& b); +const std::string operator+(const mpreal& a, const std::string b); +const std::string operator+(const std::string a, const mpreal& b); + +const mpreal operator+(const mpz_t b, const mpreal& a); +const mpreal operator+(const mpq_t b, const mpreal& a); +const mpreal operator+(const long double b, const mpreal& a); +const mpreal operator+(const double b, const mpreal& a); +const mpreal operator+(const unsigned long int b, const mpreal& a); +const mpreal operator+(const unsigned int b, const mpreal& a); +const mpreal operator+(const long int b, const mpreal& a); +const mpreal operator+(const int b, const mpreal& a); + +////////////////////////////////////////////////////////////////////////// +// - Subtraction +const mpreal operator-(const mpreal& a, const mpreal& b); + +// - Fast specialized subtraction - implemented through fast -= operations +const mpreal operator-(const mpreal& a, const mpz_t b); +const mpreal operator-(const mpreal& a, const mpq_t b); +const mpreal operator-(const mpreal& a, const long double b); +const mpreal operator-(const mpreal& a, const double b); +const mpreal operator-(const mpreal& a, const unsigned long int b); +const mpreal operator-(const mpreal& a, const unsigned int b); +const mpreal operator-(const mpreal& a, const long int b); +const mpreal operator-(const mpreal& a, const int b); +const mpreal operator-(const mpreal& a, const char* b); +const mpreal operator-(const char* a, const mpreal& b); + +const mpreal operator-(const mpz_t b, const mpreal& a); +const mpreal operator-(const mpq_t b, const mpreal& a); +const mpreal operator-(const long double b, const mpreal& a); +//const mpreal operator-(const double b, const mpreal& a); + +////////////////////////////////////////////////////////////////////////// +// * Multiplication +const mpreal operator*(const mpreal& a, const mpreal& b); + +// * Fast specialized multiplication - implemented through fast *= operations +const mpreal operator*(const mpreal& a, const mpz_t b); +const mpreal operator*(const mpreal& a, const mpq_t b); +const mpreal operator*(const mpreal& a, const long double b); +const mpreal operator*(const mpreal& a, const double b); +const mpreal operator*(const mpreal& a, const unsigned long int b); +const mpreal operator*(const mpreal& a, const unsigned int b); +const mpreal operator*(const mpreal& a, const long int b); +const mpreal operator*(const mpreal& a, const int b); + +const mpreal operator*(const mpz_t b, const mpreal& a); +const mpreal operator*(const mpq_t b, const mpreal& a); +const mpreal operator*(const long double b, const mpreal& a); +const mpreal operator*(const double b, const mpreal& a); +const mpreal operator*(const unsigned long int b, const mpreal& a); +const mpreal operator*(const unsigned int b, const mpreal& a); +const mpreal operator*(const long int b, const mpreal& a); +const mpreal operator*(const int b, const mpreal& a); + +////////////////////////////////////////////////////////////////////////// +// / Division +const mpreal operator/(const mpreal& a, const mpreal& b); + +// / Fast specialized division - implemented through fast /= operations +const mpreal operator/(const mpreal& a, const mpz_t b); +const mpreal operator/(const mpreal& a, const mpq_t b); +const mpreal operator/(const mpreal& a, const long double b); +const mpreal operator/(const mpreal& a, const double b); +const mpreal operator/(const mpreal& a, const unsigned long int b); +const mpreal operator/(const mpreal& a, const unsigned int b); +const mpreal operator/(const mpreal& a, const long int b); +const mpreal operator/(const mpreal& a, const int b); + +const mpreal operator/(const long double b, const mpreal& a); + +////////////////////////////////////////////////////////////////////////// +// Shifts operators - Multiplication/Division by a power of 2 +const mpreal operator<<(const mpreal& v, const unsigned long int k); +const mpreal operator<<(const mpreal& v, const unsigned int k); +const mpreal operator<<(const mpreal& v, const long int k); +const mpreal operator<<(const mpreal& v, const int k); + +const mpreal operator>>(const mpreal& v, const unsigned long int k); +const mpreal operator>>(const mpreal& v, const unsigned int k); +const mpreal operator>>(const mpreal& v, const long int k); +const mpreal operator>>(const mpreal& v, const int k); + +////////////////////////////////////////////////////////////////////////// +// Boolean operators +bool operator < (const mpreal& a, const unsigned long int b); +bool operator < (const mpreal& a, const unsigned int b); +bool operator < (const mpreal& a, const long int b); +bool operator < (const mpreal& a, const int b); +bool operator < (const mpreal& a, const long double b); +bool operator < (const mpreal& a, const double b); + +bool operator < (const unsigned long int a,const mpreal& b); +bool operator < (const unsigned int a, const mpreal& b); +bool operator < (const long int a, const mpreal& b); +bool operator < (const int a, const mpreal& b); +bool operator < (const long double a, const mpreal& b); +bool operator < (const double a, const mpreal& b); + +bool operator > (const mpreal& a, const unsigned long int b); +bool operator > (const mpreal& a, const unsigned int b); +bool operator > (const mpreal& a, const long int b); +bool operator > (const mpreal& a, const int b); +bool operator > (const mpreal& a, const long double b); +bool operator > (const mpreal& a, const double b); + +bool operator > (const unsigned long int a,const mpreal& b); +bool operator > (const unsigned int a, const mpreal& b); +bool operator > (const long int a, const mpreal& b); +bool operator > (const int a, const mpreal& b); +bool operator > (const long double a, const mpreal& b); +bool operator > (const double a, const mpreal& b); + +bool operator >= (const mpreal& a, const unsigned long int b); +bool operator >= (const mpreal& a, const unsigned int b); +bool operator >= (const mpreal& a, const long int b); +bool operator >= (const mpreal& a, const int b); +bool operator >= (const mpreal& a, const long double b); +bool operator >= (const mpreal& a, const double b); + +bool operator >= (const unsigned long int a,const mpreal& b); +bool operator >= (const unsigned int a, const mpreal& b); +bool operator >= (const long int a, const mpreal& b); +bool operator >= (const int a, const mpreal& b); +bool operator >= (const long double a, const mpreal& b); +bool operator >= (const double a, const mpreal& b); + +bool operator <= (const mpreal& a, const unsigned long int b); +bool operator <= (const mpreal& a, const unsigned int b); +bool operator <= (const mpreal& a, const long int b); +bool operator <= (const mpreal& a, const int b); +bool operator <= (const mpreal& a, const long double b); +bool operator <= (const mpreal& a, const double b); + +bool operator <= (const unsigned long int a,const mpreal& b); +bool operator <= (const unsigned int a, const mpreal& b); +bool operator <= (const long int a, const mpreal& b); +bool operator <= (const int a, const mpreal& b); +bool operator <= (const long double a, const mpreal& b); +bool operator <= (const double a, const mpreal& b); + +bool operator == (const mpreal& a, const unsigned long int b); +bool operator == (const mpreal& a, const unsigned int b); +bool operator == (const mpreal& a, const long int b); +bool operator == (const mpreal& a, const int b); +bool operator == (const mpreal& a, const long double b); +bool operator == (const mpreal& a, const double b); + +bool operator == (const unsigned long int a,const mpreal& b); +bool operator == (const unsigned int a, const mpreal& b); +bool operator == (const long int a, const mpreal& b); +bool operator == (const int a, const mpreal& b); +bool operator == (const long double a, const mpreal& b); +bool operator == (const double a, const mpreal& b); + +bool operator != (const mpreal& a, const unsigned long int b); +bool operator != (const mpreal& a, const unsigned int b); +bool operator != (const mpreal& a, const long int b); +bool operator != (const mpreal& a, const int b); +bool operator != (const mpreal& a, const long double b); +bool operator != (const mpreal& a, const double b); + +bool operator != (const unsigned long int a,const mpreal& b); +bool operator != (const unsigned int a, const mpreal& b); +bool operator != (const long int a, const mpreal& b); +bool operator != (const int a, const mpreal& b); +bool operator != (const long double a, const mpreal& b); +bool operator != (const double a, const mpreal& b); + +////////////////////////////////////////////////////////////////////////// +// sqrt +const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal sqrt(const long int v, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::default_rnd); + +////////////////////////////////////////////////////////////////////////// +// pow +const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +////////////////////////////////////////////////////////////////////////// +// Estimate machine epsilon for the given precision +inline const mpreal machine_epsilon(mp_prec_t prec); +inline const mpreal mpreal_min(mp_prec_t prec); +inline const mpreal mpreal_max(mp_prec_t prec); + +////////////////////////////////////////////////////////////////////////// +// Implementation of inline functions +////////////////////////////////////////////////////////////////////////// + +////////////////////////////////////////////////////////////////////////// +// Operators - Assignment +inline mpreal& mpreal::operator=(const mpreal& v) +{ + if (this!= &v) mpfr_set(mp,v.mp,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const mpf_t v) +{ + mpfr_set_f(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const mpz_t v) +{ + mpfr_set_z(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const mpq_t v) +{ + mpfr_set_q(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const long double v) +{ + mpfr_set_ld(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const double v) +{ + if(double_bits == -1 || fits_in_bits(v, double_bits)) + { + mpfr_set_d(mp,v,default_rnd); + } + else + throw conversion_overflow(); + + return *this; +} + +inline mpreal& mpreal::operator=(const unsigned long int v) +{ + mpfr_set_ui(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const unsigned int v) +{ + mpfr_set_ui(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const long int v) +{ + mpfr_set_si(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const int v) +{ + mpfr_set_si(mp,v,default_rnd); + return *this; +} + +////////////////////////////////////////////////////////////////////////// +// + Addition +inline mpreal& mpreal::operator+=(const mpreal& v) +{ + mpfr_add(mp,mp,v.mp,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+=(const mpf_t u) +{ + *this += mpreal(u); + return *this; +} + +inline mpreal& mpreal::operator+=(const mpz_t u) +{ + mpfr_add_z(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+=(const mpq_t u) +{ + mpfr_add_q(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+= (const long double u) +{ + return *this += mpreal(u); +} + +inline mpreal& mpreal::operator+= (const double u) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpfr_add_d(mp,mp,u,default_rnd); + return *this; +#else + return *this += mpreal(u); +#endif +} + +inline mpreal& mpreal::operator+=(const unsigned long int u) +{ + mpfr_add_ui(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+=(const unsigned int u) +{ + mpfr_add_ui(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+=(const long int u) +{ + mpfr_add_si(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+=(const int u) +{ + mpfr_add_si(mp,mp,u,default_rnd); + return *this; +} + +inline const mpreal mpreal::operator+()const +{ + return mpreal(*this); +} + +inline const mpreal operator+(const mpreal& a, const mpreal& b) +{ + // prec(a+b) = max(prec(a),prec(b)) + if(a.get_prec()>b.get_prec()) return mpreal(a) += b; + else return mpreal(b) += a; +} + +inline const std::string operator+(const mpreal& a, const std::string b) +{ + return (std::string)a+b; +} + +inline const std::string operator+(const std::string a, const mpreal& b) +{ + return a+(std::string)b; +} + +inline const mpreal operator+(const mpreal& a, const mpz_t b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const char* b) +{ + return a+mpreal(b); +} + +inline const mpreal operator+(const char* a, const mpreal& b) +{ + return mpreal(a)+b; + +} + +inline const mpreal operator+(const mpreal& a, const mpq_t b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const long double b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const double b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const unsigned long int b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const unsigned int b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const long int b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const int b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpz_t b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpq_t b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const long double b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const double b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const unsigned long int b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const unsigned int b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const long int b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const int b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline mpreal& mpreal::operator++() +{ + *this += 1; + return *this; +} + +inline const mpreal mpreal::operator++ (int) +{ + mpreal x(*this); + *this += 1; + return x; +} + +inline mpreal& mpreal::operator--() +{ + *this -= 1; + return *this; +} + +inline const mpreal mpreal::operator-- (int) +{ + mpreal x(*this); + *this -= 1; + return x; +} + +////////////////////////////////////////////////////////////////////////// +// - Subtraction +inline mpreal& mpreal::operator-= (const mpreal& v) +{ + mpfr_sub(mp,mp,v.mp,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const mpz_t v) +{ + mpfr_sub_z(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const mpq_t v) +{ + mpfr_sub_q(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const long double v) +{ + return *this -= mpreal(v); +} + +inline mpreal& mpreal::operator-=(const double v) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpfr_sub_d(mp,mp,v,default_rnd); + return *this; +#else + return *this -= mpreal(v); +#endif +} + +inline mpreal& mpreal::operator-=(const unsigned long int v) +{ + mpfr_sub_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const unsigned int v) +{ + mpfr_sub_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const long int v) +{ + mpfr_sub_si(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const int v) +{ + mpfr_sub_si(mp,mp,v,default_rnd); + return *this; +} + +inline const mpreal mpreal::operator-()const +{ + mpreal u(*this); + mpfr_neg(u.mp,u.mp,default_rnd); + return u; +} + +inline const mpreal operator-(const mpreal& a, const mpreal& b) +{ + // prec(a-b) = max(prec(a),prec(b)) + if(a.get_prec()>b.get_prec()) return mpreal(a) -= b; + else return -(mpreal(b) -= a); +} + +inline const mpreal operator-(const mpreal& a, const mpz_t b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const mpq_t b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const long double b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const double b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const unsigned long int b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const unsigned int b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const long int b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const int b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpz_t b, const mpreal& a) +{ + return -(mpreal(a) -= b); +} + +inline const mpreal operator-(const mpq_t b, const mpreal& a) +{ + return -(mpreal(a) -= b); +} + +inline const mpreal operator-(const long double b, const mpreal& a) +{ + return -(mpreal(a) -= b); +} + +inline const mpreal operator-(const double b, const mpreal& a) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpreal x(a); + mpfr_d_sub(x.mp,b,a.mp,mpreal::default_rnd); + return x; +#else + return -(mpreal(a) -= b); +#endif +} + +inline const mpreal operator-(const unsigned long int b, const mpreal& a) +{ + mpreal x(a); + mpfr_ui_sub(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator-(const unsigned int b, const mpreal& a) +{ + mpreal x(a); + mpfr_ui_sub(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator-(const long int b, const mpreal& a) +{ + mpreal x(a); + mpfr_si_sub(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator-(const int b, const mpreal& a) +{ + mpreal x(a); + mpfr_si_sub(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator-(const mpreal& a, const char* b) +{ + return a-mpreal(b); +} + +inline const mpreal operator-(const char* a, const mpreal& b) +{ + return mpreal(a)-b; +} + +////////////////////////////////////////////////////////////////////////// +// * Multiplication +inline mpreal& mpreal::operator*= (const mpreal& v) +{ + mpfr_mul(mp,mp,v.mp,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const mpz_t v) +{ + mpfr_mul_z(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const mpq_t v) +{ + mpfr_mul_q(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const long double v) +{ + return *this *= mpreal(v); +} + +inline mpreal& mpreal::operator*=(const double v) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpfr_mul_d(mp,mp,v,default_rnd); + return *this; +#else + return *this *= mpreal(v); +#endif +} + +inline mpreal& mpreal::operator*=(const unsigned long int v) +{ + mpfr_mul_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const unsigned int v) +{ + mpfr_mul_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const long int v) +{ + mpfr_mul_si(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const int v) +{ + mpfr_mul_si(mp,mp,v,default_rnd); + return *this; +} + +inline const mpreal operator*(const mpreal& a, const mpreal& b) +{ + // prec(a*b) = max(prec(a),prec(b)) + if(a.get_prec()>b.get_prec()) return mpreal(a) *= b; + else return mpreal(b) *= a; +} + +inline const mpreal operator*(const mpreal& a, const mpz_t b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const mpq_t b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const long double b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const double b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const unsigned long int b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const unsigned int b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const long int b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const int b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpz_t b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpq_t b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const long double b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const double b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const unsigned long int b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const unsigned int b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const long int b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const int b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +////////////////////////////////////////////////////////////////////////// +// / Division +inline mpreal& mpreal::operator/=(const mpreal& v) +{ + mpfr_div(mp,mp,v.mp,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const mpz_t v) +{ + mpfr_div_z(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const mpq_t v) +{ + mpfr_div_q(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const long double v) +{ + return *this /= mpreal(v); +} + +inline mpreal& mpreal::operator/=(const double v) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpfr_div_d(mp,mp,v,default_rnd); + return *this; +#else + return *this /= mpreal(v); +#endif +} + +inline mpreal& mpreal::operator/=(const unsigned long int v) +{ + mpfr_div_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const unsigned int v) +{ + mpfr_div_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const long int v) +{ + mpfr_div_si(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const int v) +{ + mpfr_div_si(mp,mp,v,default_rnd); + return *this; +} + +inline const mpreal operator/(const mpreal& a, const mpreal& b) +{ + mpreal x(a); + mp_prec_t pb; + mp_prec_t pa; + + // prec(a/b) = max(prec(a),prec(b)) + pa = a.get_prec(); + pb = b.get_prec(); + if(pb>pa) x.set_prec(pb); + + return x /= b; +} + +inline const mpreal operator/(const mpreal& a, const mpz_t b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const mpq_t b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const long double b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const double b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const unsigned long int b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const unsigned int b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const long int b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const int b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const unsigned long int b, const mpreal& a) +{ + mpreal x(a); + mpfr_ui_div(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator/(const unsigned int b, const mpreal& a) +{ + mpreal x(a); + mpfr_ui_div(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator/(const long int b, const mpreal& a) +{ + mpreal x(a); + mpfr_si_div(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator/(const int b, const mpreal& a) +{ + mpreal x(a); + mpfr_si_div(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator/(const long double b, const mpreal& a) +{ + mpreal x(b); + return x/a; +} + +inline const mpreal operator/(const double b, const mpreal& a) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpreal x(a); + mpfr_d_div(x.mp,b,a.mp,mpreal::default_rnd); + return x; +#else + mpreal x(b); + return x/a; +#endif +} + +////////////////////////////////////////////////////////////////////////// +// Shifts operators - Multiplication/Division by power of 2 +inline mpreal& mpreal::operator<<=(const unsigned long int u) +{ + mpfr_mul_2ui(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator<<=(const unsigned int u) +{ + mpfr_mul_2ui(mp,mp,static_cast(u),default_rnd); + return *this; +} + +inline mpreal& mpreal::operator<<=(const long int u) +{ + mpfr_mul_2si(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator<<=(const int u) +{ + mpfr_mul_2si(mp,mp,static_cast(u),default_rnd); + return *this; +} + +inline mpreal& mpreal::operator>>=(const unsigned long int u) +{ + mpfr_div_2ui(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator>>=(const unsigned int u) +{ + mpfr_div_2ui(mp,mp,static_cast(u),default_rnd); + return *this; +} + +inline mpreal& mpreal::operator>>=(const long int u) +{ + mpfr_div_2si(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator>>=(const int u) +{ + mpfr_div_2si(mp,mp,static_cast(u),default_rnd); + return *this; +} + +inline const mpreal operator<<(const mpreal& v, const unsigned long int k) +{ + return mul_2ui(v,k); +} + +inline const mpreal operator<<(const mpreal& v, const unsigned int k) +{ + return mul_2ui(v,static_cast(k)); +} + +inline const mpreal operator<<(const mpreal& v, const long int k) +{ + return mul_2si(v,k); +} + +inline const mpreal operator<<(const mpreal& v, const int k) +{ + return mul_2si(v,static_cast(k)); +} + +inline const mpreal operator>>(const mpreal& v, const unsigned long int k) +{ + return div_2ui(v,k); +} + +inline const mpreal operator>>(const mpreal& v, const long int k) +{ + return div_2si(v,k); +} + +inline const mpreal operator>>(const mpreal& v, const unsigned int k) +{ + return div_2ui(v,static_cast(k)); +} + +inline const mpreal operator>>(const mpreal& v, const int k) +{ + return div_2si(v,static_cast(k)); +} + +// mul_2ui +inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_mul_2ui(x.mp,v.mp,k,rnd_mode); + return x; +} + +// mul_2si +inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_mul_2si(x.mp,v.mp,k,rnd_mode); + return x; +} + +inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_div_2ui(x.mp,v.mp,k,rnd_mode); + return x; +} + +inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_div_2si(x.mp,v.mp,k,rnd_mode); + return x; +} + +////////////////////////////////////////////////////////////////////////// +//Boolean operators +inline bool operator > (const mpreal& a, const mpreal& b) +{ + return (mpfr_greater_p(a.mp,b.mp)!=0); +} + +inline bool operator > (const mpreal& a, const unsigned long int b) +{ + return a>mpreal(b); +} + +inline bool operator > (const mpreal& a, const unsigned int b) +{ + return a>mpreal(b); +} + +inline bool operator > (const mpreal& a, const long int b) +{ + return a>mpreal(b); +} + +inline bool operator > (const mpreal& a, const int b) +{ + return a>mpreal(b); +} + +inline bool operator > (const mpreal& a, const long double b) +{ + return a>mpreal(b); +} + +inline bool operator > (const mpreal& a, const double b) +{ + return a>mpreal(b); +} + +inline bool operator > (const unsigned long int a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator > (const unsigned int a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator > (const long int a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator > (const int a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator > (const long double a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator > (const double a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator >= (const mpreal& a, const mpreal& b) +{ + return (mpfr_greaterequal_p(a.mp,b.mp)!=0); +} + +inline bool operator >= (const mpreal& a, const unsigned long int b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const mpreal& a, const unsigned int b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const mpreal& a, const long int b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const mpreal& a, const int b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const mpreal& a, const long double b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const mpreal& a, const double b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const unsigned long int a,const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator >= (const unsigned int a, const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator >= (const long int a, const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator >= (const int a, const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator >= (const long double a, const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator >= (const double a, const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator < (const mpreal& a, const mpreal& b) +{ + return (mpfr_less_p(a.mp,b.mp)!=0); +} + +inline bool operator < (const mpreal& a, const unsigned long int b) +{ + return a= MPFR_VERSION_NUM(3,0,0)) +inline bool _isregular(const mpreal& v) +{ + return (mpfr_regular_p(v.mp)); +} +#endif // MPFR 3.0.0 Specifics + +////////////////////////////////////////////////////////////////////////// +// Type Converters +inline mpreal::operator double() const +{ + return mpfr_get_d(mp,default_rnd); +} + +inline mpreal::operator float() const +{ + return (float)mpfr_get_d(mp,default_rnd); +} + +inline mpreal::operator long double() const +{ + return mpfr_get_ld(mp,default_rnd); +} + +inline mpreal::operator unsigned long() const +{ + return mpfr_get_ui(mp,default_rnd); +} + +inline mpreal::operator unsigned int() const +{ + return static_cast(mpfr_get_ui(mp,default_rnd)); +} + +inline mpreal::operator long() const +{ + return mpfr_get_si(mp,default_rnd); +} + +inline mpreal::operator mpfr_ptr() +{ + return mp; +} + +////////////////////////////////////////////////////////////////////////// +// Set/Get number properties +inline int sgn(const mpreal& v) +{ + int r = mpfr_signbit(v.mp); + return (r>0?-1:1); +} + +inline void mpreal::set_sign(int sign, mp_rnd_t rnd_mode) +{ + mpfr_setsign(mp,mp,(sign<0?1:0),rnd_mode); +} + +inline mp_prec_t mpreal::get_prec() const +{ + return mpfr_get_prec(mp); +} + +inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpfr_prec_round(mp,prec,rnd_mode); +} + +inline void mpreal::set_inf(int sign) +{ + mpfr_set_inf(mp,sign); +} + +inline void mpreal::set_nan() +{ + mpfr_set_nan(mp); +} + +inline mp_exp_t mpreal::get_exp () +{ + return mpfr_get_exp(mp); +} + +inline int mpreal::set_exp (mp_exp_t e) +{ + return mpfr_set_exp(mp,e); +} + +inline const mpreal frexp(const mpreal& v, mp_exp_t* exp) +{ + mpreal x(v); + *exp = x.get_exp(); + x.set_exp(0); + return x; +} + +inline const mpreal ldexp(const mpreal& v, mp_exp_t exp) +{ + mpreal x(v); + + // rounding is not important since we just increasing the exponent + mpfr_mul_2si(x.mp,x.mp,exp,mpreal::default_rnd); + return x; +} + +inline const mpreal machine_epsilon(mp_prec_t prec) +{ + // smallest eps such that 1.0+eps != 1.0 + // depends (of cause) on the precision + mpreal x(1,prec); + return nextabove(x)-x; +} + +inline const mpreal mpreal_min(mp_prec_t prec) +{ + // min = 1/2*2^emin = 2^(emin-1) + + mpreal x(1,prec); + return x <<= mpreal::get_emin()-1; +} + +inline const mpreal mpreal_max(mp_prec_t prec) +{ + // max = (1-eps)*2^emax, assume eps = 0?, + // and use emax-1 to prevent value to be +inf + // max = 2^(emax-1) + + mpreal x(1,prec); + return x <<= mpreal::get_emax()-1; +} + +inline const mpreal modf(const mpreal& v, mpreal& n) +{ + mpreal frac(v); + + // rounding is not important since we are using the same number + mpfr_frac(frac.mp,frac.mp,mpreal::default_rnd); + mpfr_trunc(n.mp,v.mp); + return frac; +} + +inline int mpreal::check_range (int t, mp_rnd_t rnd_mode) +{ + return mpfr_check_range(mp,t,rnd_mode); +} + +inline int mpreal::subnormalize (int t,mp_rnd_t rnd_mode) +{ + return mpfr_subnormalize(mp,t,rnd_mode); +} + +inline mp_exp_t mpreal::get_emin (void) +{ + return mpfr_get_emin(); +} + +inline int mpreal::set_emin (mp_exp_t exp) +{ + return mpfr_set_emin(exp); +} + +inline mp_exp_t mpreal::get_emax (void) +{ + return mpfr_get_emax(); +} + +inline int mpreal::set_emax (mp_exp_t exp) +{ + return mpfr_set_emax(exp); +} + +inline mp_exp_t mpreal::get_emin_min (void) +{ + return mpfr_get_emin_min(); +} + +inline mp_exp_t mpreal::get_emin_max (void) +{ + return mpfr_get_emin_max(); +} + +inline mp_exp_t mpreal::get_emax_min (void) +{ + return mpfr_get_emax_min(); +} + +inline mp_exp_t mpreal::get_emax_max (void) +{ + return mpfr_get_emax_max(); +} + +////////////////////////////////////////////////////////////////////////// +// Mathematical Functions +////////////////////////////////////////////////////////////////////////// +inline const mpreal sqr(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sqr(x.mp,x.mp,rnd_mode); + return x; +} + +inline const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sqrt(x.mp,x.mp,rnd_mode); + return x; +} + +inline const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode) +{ + mpreal x; + mpfr_sqrt_ui(x.mp,v,rnd_mode); + return x; +} + +inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode) +{ + return sqrt(static_cast(v),rnd_mode); +} + +inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode) +{ + if (v>=0) return sqrt(static_cast(v),rnd_mode); + else return mpreal(); // NaN +} + +inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) +{ + if (v>=0) return sqrt(static_cast(v),rnd_mode); + else return mpreal(); // NaN +} + +inline const mpreal sqrt(const long double v, mp_rnd_t rnd_mode) +{ + return sqrt(mpreal(v),rnd_mode); +} + +inline const mpreal sqrt(const double v, mp_rnd_t rnd_mode) +{ + return sqrt(mpreal(v),rnd_mode); +} + +inline const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_cbrt(x.mp,x.mp,rnd_mode); + return x; +} + +inline const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_root(x.mp,x.mp,k,rnd_mode); + return x; +} + +inline const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_abs(x.mp,x.mp,rnd_mode); + return x; +} + +inline const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_abs(x.mp,x.mp,rnd_mode); + return x; +} + +inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode) +{ + mpreal x(a); + mpfr_dim(x.mp,a.mp,b.mp,rnd_mode); + return x; +} + +inline int cmpabs(const mpreal& a,const mpreal& b) +{ + return mpfr_cmpabs(a.mp,b.mp); +} + +inline const mpreal log (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_log(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal log2(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_log2(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_log10(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal exp(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_exp(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal exp2(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_exp2(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_exp10(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_cos(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sin(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_tan(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sec(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_csc(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_cot(x.mp,v.mp,rnd_mode); + return x; +} + +inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode) +{ + return mpfr_sin_cos(s.mp,c.mp,v.mp,rnd_mode); +} + +inline const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_acos(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_asin(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_atan(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode) +{ + mpreal a; + mp_prec_t yp, xp; + + yp = y.get_prec(); + xp = x.get_prec(); + + a.set_prec(yp>xp?yp:xp); + + mpfr_atan2(a.mp, y.mp, x.mp, rnd_mode); + + return a; +} + +inline const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_cosh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sinh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_tanh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sech(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_csch(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_coth(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_acosh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_asinh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_atanh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal fac_ui (unsigned long int v, mp_rnd_t rnd_mode) +{ + mpreal x; + mpfr_fac_ui(x.mp,v,rnd_mode); + return x; +} + +inline const mpreal log1p (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_log1p(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal expm1 (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_expm1(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_eint(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_gamma(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_lngamma(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_lgamma(x.mp,signp,v.mp,rnd_mode); + return x; +} + +inline const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_zeta(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_erf(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_erfc(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal _j0 (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_j0(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal _j1 (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_j1(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal _jn (long n, const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_jn(x.mp,n,v.mp,rnd_mode); + return x; +} + +inline const mpreal _y0 (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_y0(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal _y1 (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_y1(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal _yn (long n, const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_yn(x.mp,n,v.mp,rnd_mode); + return x; +} + +////////////////////////////////////////////////////////////////////////// +// MPFR 2.4.0 Specifics +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + +inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode) +{ + return mpfr_sinh_cosh(s.mp,c.mp,v.mp,rnd_mode); +} + +inline const mpreal li2(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_li2(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) +{ + mpreal a; + mp_prec_t yp, xp; + + yp = y.get_prec(); + xp = x.get_prec(); + + a.set_prec(yp>xp?yp:xp); + + mpfr_fmod(a.mp, x.mp, y.mp, rnd_mode); + + return a; +} + +inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rec_sqrt(x.mp,v.mp,rnd_mode); + return x; +} +#endif // MPFR 2.4.0 Specifics + +////////////////////////////////////////////////////////////////////////// +// MPFR 3.0.0 Specifics +#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) +inline const mpreal digamma(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_digamma(x.mp,v.mp,rnd_mode); + return x; +} +#endif // MPFR 3.0.0 Specifics + +////////////////////////////////////////////////////////////////////////// +// Constants +inline const mpreal const_log2 (mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpreal x; + x.set_prec(prec); + mpfr_const_log2(x.mp,rnd_mode); + return x; +} + +inline const mpreal const_pi (mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpreal x; + x.set_prec(prec); + mpfr_const_pi(x.mp,rnd_mode); + return x; +} + +inline const mpreal const_euler (mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpreal x; + x.set_prec(prec); + mpfr_const_euler(x.mp,rnd_mode); + return x; +} + +inline const mpreal const_catalan (mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpreal x; + x.set_prec(prec); + mpfr_const_catalan(x.mp,rnd_mode); + return x; +} + +inline const mpreal const_infinity (int sign, mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpreal x; + x.set_prec(prec,rnd_mode); + mpfr_set_inf(x.mp, sign); + return x; +} + +////////////////////////////////////////////////////////////////////////// +// Integer Related Functions +inline const mpreal rint(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rint(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal ceil(const mpreal& v) +{ + mpreal x(v); + mpfr_ceil(x.mp,v.mp); + return x; + +} + +inline const mpreal floor(const mpreal& v) +{ + mpreal x(v); + mpfr_floor(x.mp,v.mp); + return x; +} + +inline const mpreal round(const mpreal& v) +{ + mpreal x(v); + mpfr_round(x.mp,v.mp); + return x; +} + +inline const mpreal trunc(const mpreal& v) +{ + mpreal x(v); + mpfr_trunc(x.mp,v.mp); + return x; +} + +inline const mpreal rint_ceil (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rint_ceil(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal rint_floor(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rint_floor(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal rint_round(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rint_round(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal rint_trunc(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rint_trunc(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_frac(x.mp,v.mp,rnd_mode); + return x; +} + +////////////////////////////////////////////////////////////////////////// +// Miscellaneous Functions +inline void swap(mpreal& a, mpreal& b) +{ + mpfr_swap(a.mp,b.mp); +} + +#ifndef max +inline const mpreal max(const mpreal& x, const mpreal& y) +{ + return (x>y?x:y); +} +#endif + +#ifndef min +inline const mpreal min(const mpreal& x, const mpreal& y) +{ + return (x= MPFR_VERSION_NUM(3,0,0)) +// use gmp_randinit_default() to init state, gmp_randclear() to clear +inline const mpreal urandom (gmp_randstate_t& state,mp_rnd_t rnd_mode) +{ + mpreal x; + mpfr_urandom(x.mp,state,rnd_mode); + return x; +} +#endif + +#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2)) +inline const mpreal random2 (mp_size_t size, mp_exp_t exp) +{ + mpreal x; + mpfr_random2(x.mp,size,exp); + return x; +} +#endif + +////////////////////////////////////////////////////////////////////////// +// Set/Get global properties +inline void mpreal::set_default_prec(mp_prec_t prec) +{ + default_prec = prec; + mpfr_set_default_prec(prec); +} + +inline mp_prec_t mpreal::get_default_prec() +{ + return mpfr_get_default_prec(); +} + +inline void mpreal::set_default_base(int base) +{ + default_base = base; +} + +inline int mpreal::get_default_base() +{ + return default_base; +} + +inline void mpreal::set_default_rnd(mp_rnd_t rnd_mode) +{ + default_rnd = rnd_mode; + mpfr_set_default_rounding_mode(rnd_mode); +} + +inline mp_rnd_t mpreal::get_default_rnd() +{ + return mpfr_get_default_rounding_mode(); +} + +inline void mpreal::set_double_bits(int dbits) +{ + double_bits = dbits; +} + +inline int mpreal::get_double_bits() +{ + return double_bits; +} + +inline bool mpreal::fits_in_bits(double x, int n) +{ + int i; + double t; + return IsInf(x) || (std::modf ( std::ldexp ( std::frexp ( x, &i ), n ), &t ) == 0.0); +} + +inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode) +{ + mpreal x(a); + mpfr_pow(x.mp,x.mp,b.mp,rnd_mode); + return x; +} + +inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode) +{ + mpreal x(a); + mpfr_pow_z(x.mp,x.mp,b,rnd_mode); + return x; +} + +inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode) +{ + mpreal x(a); + mpfr_pow_ui(x.mp,x.mp,b,rnd_mode); + return x; +} + +inline const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode) +{ + return pow(a,static_cast(b),rnd_mode); +} + +inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode) +{ + mpreal x(a); + mpfr_pow_si(x.mp,x.mp,b,rnd_mode); + return x; +} + +inline const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode) +{ + return pow(a,static_cast(b),rnd_mode); +} + +inline const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode) +{ + return pow(a,mpreal(b),rnd_mode); +} + +inline const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode) +{ + return pow(a,mpreal(b),rnd_mode); +} + +inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode) +{ + mpreal x(a); + mpfr_ui_pow(x.mp,a,b.mp,rnd_mode); + return x; +} + +inline const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode) +{ + return pow(static_cast(a),b,rnd_mode); +} + +inline const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode) +{ + if (a>=0) return pow(static_cast(a),b,rnd_mode); + else return pow(mpreal(a),b,rnd_mode); +} + +inline const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode) +{ + if (a>=0) return pow(static_cast(a),b,rnd_mode); + else return pow(mpreal(a),b,rnd_mode); +} + +inline const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),b,rnd_mode); +} + +inline const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),b,rnd_mode); +} + +// pow unsigned long int +inline const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode) +{ + mpreal x(a); + mpfr_ui_pow_ui(x.mp,a,b,rnd_mode); + return x; +} + +inline const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode) +{ + return pow(a,static_cast(b),rnd_mode); //mpfr_ui_pow_ui +} + +inline const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode) +{ + if(b>0) return pow(a,static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow +} + +inline const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode) +{ + if(b>0) return pow(a,static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow +} + +inline const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode) +{ + return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow +} + +inline const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode) +{ + return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow +} + +// pow unsigned int +inline const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode) +{ + return pow(static_cast(a),b,rnd_mode); //mpfr_ui_pow_ui +} + +inline const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode) +{ + return pow(static_cast(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui +} + +inline const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode) +{ + if(b>0) return pow(static_cast(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow +} + +inline const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode) +{ + if(b>0) return pow(static_cast(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow +} + +inline const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode) +{ + return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow +} + +inline const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode) +{ + return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow +} + +// pow long int +inline const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode) +{ + if (a>0) return pow(static_cast(a),b,rnd_mode); //mpfr_ui_pow_ui + else return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui +} + +inline const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode) +{ + if (a>0) return pow(static_cast(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(mpreal(a),static_cast(b),rnd_mode); //mpfr_pow_ui +} + +inline const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode) +{ + if (a>0) + { + if(b>0) return pow(static_cast(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow + }else{ + return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si + } +} + +inline const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode) +{ + if (a>0) + { + if(b>0) return pow(static_cast(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow + }else{ + return pow(mpreal(a),static_cast(b),rnd_mode); // mpfr_pow_si + } +} + +inline const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode) +{ + if (a>=0) return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow + else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow +} + +inline const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode) +{ + if (a>=0) return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow + else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow +} + +// pow int +inline const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode) +{ + if (a>0) return pow(static_cast(a),b,rnd_mode); //mpfr_ui_pow_ui + else return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui +} + +inline const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode) +{ + if (a>0) return pow(static_cast(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(mpreal(a),static_cast(b),rnd_mode); //mpfr_pow_ui +} + +inline const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode) +{ + if (a>0) + { + if(b>0) return pow(static_cast(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow + }else{ + return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si + } +} + +inline const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode) +{ + if (a>0) + { + if(b>0) return pow(static_cast(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow + }else{ + return pow(mpreal(a),static_cast(b),rnd_mode); // mpfr_pow_si + } +} + +inline const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode) +{ + if (a>=0) return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow + else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow +} + +inline const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode) +{ + if (a>=0) return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow + else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow +} + +// pow long double +inline const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),mpreal(b),rnd_mode); +} + +inline const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui +} + +inline const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),static_cast(b),rnd_mode); //mpfr_pow_ui +} + +inline const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si +} + +inline const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),static_cast(b),rnd_mode); // mpfr_pow_si +} + +inline const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),mpreal(b),rnd_mode); +} + +inline const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),b,rnd_mode); // mpfr_pow_ui +} + +inline const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),static_cast(b),rnd_mode); // mpfr_pow_ui +} + +inline const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si +} + +inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode) +{ + return pow(mpreal(a),static_cast(b),rnd_mode); // mpfr_pow_si +} +} + +// Explicit specialization of std::swap for mpreal numbers +// Thus standard algorithms will use efficient version of swap (due to Koenig lookup) +// Non-throwing swap C++ idiom: http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Non-throwing_swap +namespace std +{ + template <> + inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) + { + return mpfr::swap(x, y); + } +} + +#endif /* __MP_REAL_H__ */ diff --git a/unsupported/test/mpreal_support.cpp b/unsupported/test/mpreal_support.cpp new file mode 100644 index 000000000..32570f092 --- /dev/null +++ b/unsupported/test/mpreal_support.cpp @@ -0,0 +1,42 @@ +#include "main.h" +#include + +using namespace mpfr; +using namespace std; +using namespace Eigen; + +void test_mpreal_support() +{ + // set precision to 256 bits (double has only 53 bits) + mpreal::set_default_prec(256); + typedef Matrix MatrixXmp; + + std::cerr << "epsilon = " << NumTraits::epsilon() << "\n"; + std::cerr << "dummy_precision = " << NumTraits::dummy_precision() << "\n"; + std::cerr << "highest = " << NumTraits::highest() << "\n"; + std::cerr << "lowest = " << NumTraits::lowest() << "\n"; + + for(int i = 0; i < g_repeat; i++) { + int s = ei_random(1,100); + MatrixXmp A = MatrixXmp::Random(s,s); + MatrixXmp B = MatrixXmp::Random(s,s); + MatrixXmp S = A.adjoint() * A; + MatrixXmp X; + + // Cholesky + X = S.selfadjointView().llt().solve(B); + VERIFY_IS_APPROX((S.selfadjointView()*X).eval(),B); + + // partial LU + X = A.lu().solve(B); + VERIFY_IS_APPROX((A*X).eval(),B); + + // symmetric eigenvalues + SelfAdjointEigenSolver eig(S); + VERIFY_IS_EQUAL(eig.info(), Success); + VERIFY_IS_APPROX((S.selfadjointView() * eig.eigenvectors()), + eig.eigenvectors() * eig.eigenvalues().asDiagonal()); + } +} + +#include "mpreal.cpp"