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 <[email protected]> 5// Copyright (C) 2010 Konstantin Holoborodko <[email protected]> 6// Copyright (C) 2010 Gael Guennebaud <[email protected]> 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 18namespace 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 * \warning MPFR C++ is licensed under the GPL. 31 * 32 * You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder. 33 * 34 * Here is an example: 35 * 36\code 37#include <iostream> 38#include <Eigen/MPRealSupport> 39#include <Eigen/LU> 40using namespace mpfr; 41using namespace Eigen; 42int main() 43{ 44 // set precision to 256 bits (double has only 53 bits) 45 mpreal::set_default_prec(256); 46 // Declare matrix and vector types with multi-precision scalar type 47 typedef Matrix<mpreal,Dynamic,Dynamic> MatrixXmp; 48 typedef Matrix<mpreal,Dynamic,1> VectorXmp; 49 50 MatrixXmp A = MatrixXmp::Random(100,100); 51 VectorXmp b = VectorXmp::Random(100); 52 53 // Solve Ax=b using LU 54 VectorXmp x = A.lu().solve(b); 55 std::cout << "relative error: " << (A*x - b).norm() / b.norm() << std::endl; 56 return 0; 57} 58\endcode 59 * 60 */ 61 62 template<> struct NumTraits<mpfr::mpreal> 63 : GenericNumTraits<mpfr::mpreal> 64 { 65 enum { 66 IsInteger = 0, 67 IsSigned = 1, 68 IsComplex = 0, 69 RequireInitialization = 1, 70 ReadCost = HugeCost, 71 AddCost = HugeCost, 72 MulCost = HugeCost 73 }; 74 75 typedef mpfr::mpreal Real; 76 typedef mpfr::mpreal NonInteger; 77 78 static inline Real highest (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(Precision); } 79 static inline Real lowest (long Precision = mpfr::mpreal::get_default_prec()) { return -mpfr::maxval(Precision); } 80 81 // Constants 82 static inline Real Pi (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_pi(Precision); } 83 static inline Real Euler (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_euler(Precision); } 84 static inline Real Log2 (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_log2(Precision); } 85 static inline Real Catalan (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_catalan(Precision); } 86 87 static inline Real epsilon (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(Precision); } 88 static inline Real epsilon (const Real& x) { return mpfr::machine_epsilon(x); } 89 90#ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS 91 static inline int digits10 (long Precision = mpfr::mpreal::get_default_prec()) { return std::numeric_limits<Real>::digits10(Precision); } 92 static inline int digits10 (const Real& x) { return std::numeric_limits<Real>::digits10(x); } 93 94 static inline int digits () { return std::numeric_limits<Real>::digits(); } 95 static inline int digits (const Real& x) { return std::numeric_limits<Real>::digits(x); } 96#endif 97 98 static inline Real dummy_precision() 99 { 100 mpfr_prec_t weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100; 101 return mpfr::machine_epsilon(weak_prec); 102 } 103 }; 104 105 namespace internal { 106 107 template<> inline mpfr::mpreal random<mpfr::mpreal>() 108 { 109 return mpfr::random(); 110 } 111 112 template<> inline mpfr::mpreal random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b) 113 { 114 return a + (b-a) * random<mpfr::mpreal>(); 115 } 116 117 inline bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) 118 { 119 return mpfr::abs(a) <= mpfr::abs(b) * eps; 120 } 121 122 inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) 123 { 124 return mpfr::isEqualFuzzy(a,b,eps); 125 } 126 127 inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) 128 { 129 return a <= b || mpfr::isEqualFuzzy(a,b,eps); 130 } 131 132 template<> inline long double cast<mpfr::mpreal,long double>(const mpfr::mpreal& x) 133 { return x.toLDouble(); } 134 135 template<> inline double cast<mpfr::mpreal,double>(const mpfr::mpreal& x) 136 { return x.toDouble(); } 137 138 template<> inline long cast<mpfr::mpreal,long>(const mpfr::mpreal& x) 139 { return x.toLong(); } 140 141 template<> inline int cast<mpfr::mpreal,int>(const mpfr::mpreal& x) 142 { return int(x.toLong()); } 143 144 // Specialize GEBP kernel and traits for mpreal (no need for peeling, nor complicated stuff) 145 // This also permits to directly call mpfr's routines and avoid many temporaries produced by mpreal 146 template<> 147 class gebp_traits<mpfr::mpreal, mpfr::mpreal, false, false> 148 { 149 public: 150 typedef mpfr::mpreal ResScalar; 151 enum { 152 Vectorizable = false, 153 LhsPacketSize = 1, 154 RhsPacketSize = 1, 155 ResPacketSize = 1, 156 NumberOfRegisters = 1, 157 nr = 1, 158 mr = 1, 159 LhsProgress = 1, 160 RhsProgress = 1 161 }; 162 typedef ResScalar LhsPacket; 163 typedef ResScalar RhsPacket; 164 typedef ResScalar ResPacket; 165 typedef LhsPacket LhsPacket4Packing; 166 167 }; 168 169 170 171 template<typename Index, typename DataMapper, bool ConjugateLhs, bool ConjugateRhs> 172 struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,DataMapper,1,1,ConjugateLhs,ConjugateRhs> 173 { 174 typedef mpfr::mpreal mpreal; 175 176 EIGEN_DONT_INLINE 177 void operator()(const DataMapper& res, const mpreal* blockA, const mpreal* blockB, 178 Index rows, Index depth, Index cols, const mpreal& alpha, 179 Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0) 180 { 181 if(rows==0 || cols==0 || depth==0) 182 return; 183 184 mpreal acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())), 185 tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr())); 186 187 if(strideA==-1) strideA = depth; 188 if(strideB==-1) strideB = depth; 189 190 for(Index i=0; i<rows; ++i) 191 { 192 for(Index j=0; j<cols; ++j) 193 { 194 const mpreal *A = blockA + i*strideA + offsetA; 195 const mpreal *B = blockB + j*strideB + offsetB; 196 197 acc1 = 0; 198 for(Index k=0; k<depth; k++) 199 { 200 mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_srcptr(), B[k].mpfr_srcptr(), mpreal::get_default_rnd()); 201 mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd()); 202 } 203 204 mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_srcptr(), alpha.mpfr_srcptr(), mpreal::get_default_rnd()); 205 mpfr_add(res(i,j).mpfr_ptr(), res(i,j).mpfr_srcptr(), acc1.mpfr_srcptr(), mpreal::get_default_rnd()); 206 } 207 } 208 } 209 }; 210 } // end namespace internal 211} 212 213#endif // EIGEN_MPREALSUPPORT_MODULE_H 214