blob: c4ea4ec5f6d80f7c97670371e3b3cb64fee2106b [file] [log] [blame]
 // 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-2012 Pavel Holoborodko // Copyright (C) 2010 Konstantin Holoborodko // Copyright (C) 2010 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_MPREALSUPPORT_MODULE_H #define EIGEN_MPREALSUPPORT_MODULE_H #include "../../Eigen/Core" #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. * * \warning MPFR C++ is licensed under the GPL. * * You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder. * * 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, RequireInitialization = 1, ReadCost = HugeCost, AddCost = HugeCost, MulCost = HugeCost }; typedef mpfr::mpreal Real; typedef mpfr::mpreal NonInteger; static inline Real highest (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(Precision); } static inline Real lowest (long Precision = mpfr::mpreal::get_default_prec()) { return -mpfr::maxval(Precision); } // Constants static inline Real Pi (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_pi(Precision); } static inline Real Euler (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_euler(Precision); } static inline Real Log2 (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_log2(Precision); } static inline Real Catalan (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_catalan(Precision); } static inline Real epsilon (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(Precision); } static inline Real epsilon (const Real& x) { return mpfr::machine_epsilon(x); } #ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS static inline int digits10 (long Precision = mpfr::mpreal::get_default_prec()) { return std::numeric_limits::digits10(Precision); } static inline int digits10 (const Real& x) { return std::numeric_limits::digits10(x); } static inline int digits () { return std::numeric_limits::digits(); } static inline int digits (const Real& x) { return std::numeric_limits::digits(x); } #endif static inline Real dummy_precision() { mpfr_prec_t weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100; return mpfr::machine_epsilon(weak_prec); } }; namespace internal { template<> inline mpfr::mpreal random() { return mpfr::random(); } template<> inline mpfr::mpreal random(const mpfr::mpreal& a, const mpfr::mpreal& b) { return a + (b-a) * random(); } inline bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) { return mpfr::abs(a) <= mpfr::abs(b) * eps; } inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) { return mpfr::isEqualFuzzy(a,b,eps); } inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) { return a <= b || mpfr::isEqualFuzzy(a,b,eps); } template<> inline long double cast(const mpfr::mpreal& x) { return x.toLDouble(); } template<> inline double cast(const mpfr::mpreal& x) { return x.toDouble(); } template<> inline long cast(const mpfr::mpreal& x) { return x.toLong(); } template<> inline int cast(const mpfr::mpreal& x) { return int(x.toLong()); } // Specialize GEBP kernel and traits for mpreal (no need for peeling, nor complicated stuff) // This also permits to directly call mpfr's routines and avoid many temporaries produced by mpreal template<> class gebp_traits { public: typedef mpfr::mpreal ResScalar; enum { Vectorizable = false, LhsPacketSize = 1, RhsPacketSize = 1, ResPacketSize = 1, NumberOfRegisters = 1, nr = 1, mr = 1, LhsProgress = 1, RhsProgress = 1 }; typedef ResScalar LhsPacket; typedef ResScalar RhsPacket; typedef ResScalar ResPacket; typedef LhsPacket LhsPacket4Packing; }; template struct gebp_kernel { typedef mpfr::mpreal mpreal; EIGEN_DONT_INLINE void operator()(const DataMapper& res, const mpreal* blockA, const mpreal* blockB, Index rows, Index depth, Index cols, const mpreal& alpha, Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0) { if(rows==0 || cols==0 || depth==0) return; mpreal acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())), tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr())); if(strideA==-1) strideA = depth; if(strideB==-1) strideB = depth; for(Index i=0; i