Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
 
Loading...
Searching...
No Matches
MPRealSupport
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@holoborodko.com>
5// Copyright (C) 2010 Konstantin Holoborodko <konstantin@holoborodko.com>
6// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@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
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