1 // This file is part of a joint effort between Eigen, a lightweight C++ template library 2 // for linear algebra, and MPFR C++, a C++ interface to MPFR library (http://www.holoborodko.com/pavel/) 3 // 4 // Copyright (C) 2010-2012 Pavel Holoborodko <pavel (a] holoborodko.com> 5 // Copyright (C) 2010 Konstantin Holoborodko <konstantin (a] holoborodko.com> 6 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud (a] inria.fr> 7 // 8 // This Source Code Form is subject to the terms of the Mozilla 9 // Public License v. 2.0. If a copy of the MPL was not distributed 10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 11 12 #ifndef EIGEN_MPREALSUPPORT_MODULE_H 13 #define EIGEN_MPREALSUPPORT_MODULE_H 14 15 #include <Eigen/Core> 16 #include <mpreal.h> 17 18 namespace Eigen { 19 20 /** 21 * \defgroup MPRealSupport_Module MPFRC++ Support module 22 * \code 23 * #include <Eigen/MPRealSupport> 24 * \endcode 25 * 26 * This module provides support for multi precision floating point numbers 27 * via the <a href="http://www.holoborodko.com/pavel/mpfr">MPFR C++</a> 28 * library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>. 29 * 30 * You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder. 31 * 32 * Here is an example: 33 * 34 \code 35 #include <iostream> 36 #include <Eigen/MPRealSupport> 37 #include <Eigen/LU> 38 using namespace mpfr; 39 using namespace Eigen; 40 int main() 41 { 42 // set precision to 256 bits (double has only 53 bits) 43 mpreal::set_default_prec(256); 44 // Declare matrix and vector types with multi-precision scalar type 45 typedef Matrix<mpreal,Dynamic,Dynamic> MatrixXmp; 46 typedef Matrix<mpreal,Dynamic,1> VectorXmp; 47 48 MatrixXmp A = MatrixXmp::Random(100,100); 49 VectorXmp b = VectorXmp::Random(100); 50 51 // Solve Ax=b using LU 52 VectorXmp x = A.lu().solve(b); 53 std::cout << "relative error: " << (A*x - b).norm() / b.norm() << std::endl; 54 return 0; 55 } 56 \endcode 57 * 58 */ 59 60 template<> struct NumTraits<mpfr::mpreal> 61 : GenericNumTraits<mpfr::mpreal> 62 { 63 enum { 64 IsInteger = 0, 65 IsSigned = 1, 66 IsComplex = 0, 67 RequireInitialization = 1, 68 ReadCost = 10, 69 AddCost = 10, 70 MulCost = 40 71 }; 72 73 typedef mpfr::mpreal Real; 74 typedef mpfr::mpreal NonInteger; 75 76 inline static Real highest (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(Precision); } 77 inline static Real lowest (long Precision = mpfr::mpreal::get_default_prec()) { return -mpfr::maxval(Precision); } 78 79 // Constants 80 inline static Real Pi (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_pi(Precision); } 81 inline static Real Euler (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_euler(Precision); } 82 inline static Real Log2 (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_log2(Precision); } 83 inline static Real Catalan (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_catalan(Precision); } 84 85 inline static Real epsilon (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(Precision); } 86 inline static Real epsilon (const Real& x) { return mpfr::machine_epsilon(x); } 87 88 inline static Real dummy_precision() 89 { 90 unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100; 91 return mpfr::machine_epsilon(weak_prec); 92 } 93 }; 94 95 namespace internal { 96 97 template<> inline mpfr::mpreal random<mpfr::mpreal>() 98 { 99 return mpfr::random(); 100 } 101 102 template<> inline mpfr::mpreal random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b) 103 { 104 return a + (b-a) * random<mpfr::mpreal>(); 105 } 106 107 inline bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) 108 { 109 return mpfr::abs(a) <= mpfr::abs(b) * eps; 110 } 111 112 inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) 113 { 114 return mpfr::isEqualFuzzy(a,b,eps); 115 } 116 117 inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) 118 { 119 return a <= b || mpfr::isEqualFuzzy(a,b,eps); 120 } 121 122 template<> inline long double cast<mpfr::mpreal,long double>(const mpfr::mpreal& x) 123 { return x.toLDouble(); } 124 125 template<> inline double cast<mpfr::mpreal,double>(const mpfr::mpreal& x) 126 { return x.toDouble(); } 127 128 template<> inline long cast<mpfr::mpreal,long>(const mpfr::mpreal& x) 129 { return x.toLong(); } 130 131 template<> inline int cast<mpfr::mpreal,int>(const mpfr::mpreal& x) 132 { return int(x.toLong()); } 133 134 // Specialize GEBP kernel and traits for mpreal (no need for peeling, nor complicated stuff) 135 // This also permits to directly call mpfr's routines and avoid many temporaries produced by mpreal 136 template<> 137 class gebp_traits<mpfr::mpreal, mpfr::mpreal, false, false> 138 { 139 public: 140 typedef mpfr::mpreal ResScalar; 141 enum { 142 nr = 2, // must be 2 for proper packing... 143 mr = 1, 144 WorkSpaceFactor = nr, 145 LhsProgress = 1, 146 RhsProgress = 1 147 }; 148 }; 149 150 template<typename Index, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> 151 struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,mr,nr,ConjugateLhs,ConjugateRhs> 152 { 153 typedef mpfr::mpreal mpreal; 154 155 EIGEN_DONT_INLINE 156 void operator()(mpreal* res, Index resStride, const mpreal* blockA, const mpreal* blockB, Index rows, Index depth, Index cols, mpreal alpha, 157 Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0, mpreal* /*unpackedB*/ = 0) 158 { 159 mpreal acc1, acc2, tmp; 160 161 if(strideA==-1) strideA = depth; 162 if(strideB==-1) strideB = depth; 163 164 for(Index j=0; j<cols; j+=nr) 165 { 166 Index actual_nr = (std::min<Index>)(nr,cols-j); 167 mpreal *C1 = res + j*resStride; 168 mpreal *C2 = res + (j+1)*resStride; 169 for(Index i=0; i<rows; i++) 170 { 171 mpreal *B = const_cast<mpreal*>(blockB) + j*strideB + offsetB*actual_nr; 172 mpreal *A = const_cast<mpreal*>(blockA) + i*strideA + offsetA; 173 acc1 = 0; 174 acc2 = 0; 175 for(Index k=0; k<depth; k++) 176 { 177 mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[0].mpfr_ptr(), mpreal::get_default_rnd()); 178 mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd()); 179 180 if(actual_nr==2) { 181 mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[1].mpfr_ptr(), mpreal::get_default_rnd()); 182 mpfr_add(acc2.mpfr_ptr(), acc2.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd()); 183 } 184 185 B+=actual_nr; 186 } 187 188 mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd()); 189 mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_ptr(), acc1.mpfr_ptr(), mpreal::get_default_rnd()); 190 191 if(actual_nr==2) { 192 mpfr_mul(acc2.mpfr_ptr(), acc2.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd()); 193 mpfr_add(C2[i].mpfr_ptr(), C2[i].mpfr_ptr(), acc2.mpfr_ptr(), mpreal::get_default_rnd()); 194 } 195 } 196 } 197 } 198 }; 199 200 } // end namespace internal 201 } 202 203 #endif // EIGEN_MPREALSUPPORT_MODULE_H 204
