Warning, /include/eigen3/unsupported/Eigen/MPRealSupport is written in an unsupported language. File is not indexed.
0001 // This file is part of a joint effort between Eigen, a lightweight C++ template library
0002 // for linear algebra, and MPFR C++, a C++ interface to MPFR library (http://www.holoborodko.com/pavel/)
0003 //
0004 // Copyright (C) 2010-2012 Pavel Holoborodko <pavel@holoborodko.com>
0005 // Copyright (C) 2010 Konstantin Holoborodko <konstantin@holoborodko.com>
0006 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
0007 //
0008 // This Source Code Form is subject to the terms of the Mozilla
0009 // Public License v. 2.0. If a copy of the MPL was not distributed
0010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0011
0012 #ifndef EIGEN_MPREALSUPPORT_MODULE_H
0013 #define EIGEN_MPREALSUPPORT_MODULE_H
0014
0015 #include "../../Eigen/Core"
0016 #include <mpreal.h>
0017
0018 namespace Eigen {
0019
0020 /**
0021 * \defgroup MPRealSupport_Module MPFRC++ Support module
0022 * \code
0023 * #include <Eigen/MPRealSupport>
0024 * \endcode
0025 *
0026 * This module provides support for multi precision floating point numbers
0027 * via the <a href="http://www.holoborodko.com/pavel/mpfr">MPFR C++</a>
0028 * library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>.
0029 *
0030 * \warning MPFR C++ is licensed under the GPL.
0031 *
0032 * You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder.
0033 *
0034 * Here is an example:
0035 *
0036 \code
0037 #include <iostream>
0038 #include <Eigen/MPRealSupport>
0039 #include <Eigen/LU>
0040 using namespace mpfr;
0041 using namespace Eigen;
0042 int main()
0043 {
0044 // set precision to 256 bits (double has only 53 bits)
0045 mpreal::set_default_prec(256);
0046 // Declare matrix and vector types with multi-precision scalar type
0047 typedef Matrix<mpreal,Dynamic,Dynamic> MatrixXmp;
0048 typedef Matrix<mpreal,Dynamic,1> VectorXmp;
0049
0050 MatrixXmp A = MatrixXmp::Random(100,100);
0051 VectorXmp b = VectorXmp::Random(100);
0052
0053 // Solve Ax=b using LU
0054 VectorXmp x = A.lu().solve(b);
0055 std::cout << "relative error: " << (A*x - b).norm() / b.norm() << std::endl;
0056 return 0;
0057 }
0058 \endcode
0059 *
0060 */
0061
0062 template<> struct NumTraits<mpfr::mpreal>
0063 : GenericNumTraits<mpfr::mpreal>
0064 {
0065 enum {
0066 IsInteger = 0,
0067 IsSigned = 1,
0068 IsComplex = 0,
0069 RequireInitialization = 1,
0070 ReadCost = HugeCost,
0071 AddCost = HugeCost,
0072 MulCost = HugeCost
0073 };
0074
0075 typedef mpfr::mpreal Real;
0076 typedef mpfr::mpreal NonInteger;
0077
0078 static inline Real highest (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(Precision); }
0079 static inline Real lowest (long Precision = mpfr::mpreal::get_default_prec()) { return -mpfr::maxval(Precision); }
0080
0081 // Constants
0082 static inline Real Pi (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_pi(Precision); }
0083 static inline Real Euler (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_euler(Precision); }
0084 static inline Real Log2 (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_log2(Precision); }
0085 static inline Real Catalan (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_catalan(Precision); }
0086
0087 static inline Real epsilon (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(Precision); }
0088 static inline Real epsilon (const Real& x) { return mpfr::machine_epsilon(x); }
0089
0090 #ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS
0091 static inline int digits10 (long Precision = mpfr::mpreal::get_default_prec()) { return std::numeric_limits<Real>::digits10(Precision); }
0092 static inline int digits10 (const Real& x) { return std::numeric_limits<Real>::digits10(x); }
0093
0094 static inline int digits () { return std::numeric_limits<Real>::digits(); }
0095 static inline int digits (const Real& x) { return std::numeric_limits<Real>::digits(x); }
0096 #endif
0097
0098 static inline Real dummy_precision()
0099 {
0100 mpfr_prec_t weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100;
0101 return mpfr::machine_epsilon(weak_prec);
0102 }
0103 };
0104
0105 namespace internal {
0106
0107 template<> inline mpfr::mpreal random<mpfr::mpreal>()
0108 {
0109 return mpfr::random();
0110 }
0111
0112 template<> inline mpfr::mpreal random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b)
0113 {
0114 return a + (b-a) * random<mpfr::mpreal>();
0115 }
0116
0117 inline bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
0118 {
0119 return mpfr::abs(a) <= mpfr::abs(b) * eps;
0120 }
0121
0122 inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
0123 {
0124 return mpfr::isEqualFuzzy(a,b,eps);
0125 }
0126
0127 inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
0128 {
0129 return a <= b || mpfr::isEqualFuzzy(a,b,eps);
0130 }
0131
0132 template<> inline long double cast<mpfr::mpreal,long double>(const mpfr::mpreal& x)
0133 { return x.toLDouble(); }
0134
0135 template<> inline double cast<mpfr::mpreal,double>(const mpfr::mpreal& x)
0136 { return x.toDouble(); }
0137
0138 template<> inline long cast<mpfr::mpreal,long>(const mpfr::mpreal& x)
0139 { return x.toLong(); }
0140
0141 template<> inline int cast<mpfr::mpreal,int>(const mpfr::mpreal& x)
0142 { return int(x.toLong()); }
0143
0144 // Specialize GEBP kernel and traits for mpreal (no need for peeling, nor complicated stuff)
0145 // This also permits to directly call mpfr's routines and avoid many temporaries produced by mpreal
0146 template<>
0147 class gebp_traits<mpfr::mpreal, mpfr::mpreal, false, false>
0148 {
0149 public:
0150 typedef mpfr::mpreal ResScalar;
0151 enum {
0152 Vectorizable = false,
0153 LhsPacketSize = 1,
0154 RhsPacketSize = 1,
0155 ResPacketSize = 1,
0156 NumberOfRegisters = 1,
0157 nr = 1,
0158 mr = 1,
0159 LhsProgress = 1,
0160 RhsProgress = 1
0161 };
0162 typedef ResScalar LhsPacket;
0163 typedef ResScalar RhsPacket;
0164 typedef ResScalar ResPacket;
0165 typedef LhsPacket LhsPacket4Packing;
0166
0167 };
0168
0169
0170
0171 template<typename Index, typename DataMapper, bool ConjugateLhs, bool ConjugateRhs>
0172 struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,DataMapper,1,1,ConjugateLhs,ConjugateRhs>
0173 {
0174 typedef mpfr::mpreal mpreal;
0175
0176 EIGEN_DONT_INLINE
0177 void operator()(const DataMapper& res, const mpreal* blockA, const mpreal* blockB,
0178 Index rows, Index depth, Index cols, const mpreal& alpha,
0179 Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0)
0180 {
0181 if(rows==0 || cols==0 || depth==0)
0182 return;
0183
0184 mpreal acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())),
0185 tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr()));
0186
0187 if(strideA==-1) strideA = depth;
0188 if(strideB==-1) strideB = depth;
0189
0190 for(Index i=0; i<rows; ++i)
0191 {
0192 for(Index j=0; j<cols; ++j)
0193 {
0194 const mpreal *A = blockA + i*strideA + offsetA;
0195 const mpreal *B = blockB + j*strideB + offsetB;
0196
0197 acc1 = 0;
0198 for(Index k=0; k<depth; k++)
0199 {
0200 mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_srcptr(), B[k].mpfr_srcptr(), mpreal::get_default_rnd());
0201 mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd());
0202 }
0203
0204 mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_srcptr(), alpha.mpfr_srcptr(), mpreal::get_default_rnd());
0205 mpfr_add(res(i,j).mpfr_ptr(), res(i,j).mpfr_srcptr(), acc1.mpfr_srcptr(), mpreal::get_default_rnd());
0206 }
0207 }
0208 }
0209 };
0210 } // end namespace internal
0211 }
0212
0213 #endif // EIGEN_MPREALSUPPORT_MODULE_H
