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#endif
94
95    static inline Real dummy_precision()
96    {
97      mpfr_prec_t weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100;
98      return mpfr::machine_epsilon(weak_prec);
99    }
100  };
101
102  namespace internal {
103
104  template<> inline mpfr::mpreal random<mpfr::mpreal>()
105  {
106    return mpfr::random();
107  }
108
109  template<> inline mpfr::mpreal random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b)
110  {
111    return a + (b-a) * random<mpfr::mpreal>();
112  }
113
114  inline bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
115  {
116    return mpfr::abs(a) <= mpfr::abs(b) * eps;
117  }
118
119  inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
120  {
121    return mpfr::isEqualFuzzy(a,b,eps);
122  }
123
124  inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
125  {
126    return a <= b || mpfr::isEqualFuzzy(a,b,eps);
127  }
128
129  template<> inline long double cast<mpfr::mpreal,long double>(const mpfr::mpreal& x)
130  { return x.toLDouble(); }
131
132  template<> inline double cast<mpfr::mpreal,double>(const mpfr::mpreal& x)
133  { return x.toDouble(); }
134
135  template<> inline long cast<mpfr::mpreal,long>(const mpfr::mpreal& x)
136  { return x.toLong(); }
137
138  template<> inline int cast<mpfr::mpreal,int>(const mpfr::mpreal& x)
139  { return int(x.toLong()); }
140
141  // Specialize GEBP kernel and traits for mpreal (no need for peeling, nor complicated stuff)
142  // This also permits to directly call mpfr's routines and avoid many temporaries produced by mpreal
143    template<>
144    class gebp_traits<mpfr::mpreal, mpfr::mpreal, false, false>
145    {
146    public:
147      typedef mpfr::mpreal ResScalar;
148      enum {
149        Vectorizable = false,
150        LhsPacketSize = 1,
151        RhsPacketSize = 1,
152        ResPacketSize = 1,
153        NumberOfRegisters = 1,
154        nr = 1,
155        mr = 1,
156        LhsProgress = 1,
157        RhsProgress = 1
158      };
159      typedef ResScalar LhsPacket;
160      typedef ResScalar RhsPacket;
161      typedef ResScalar ResPacket;
162
163    };
164
165
166
167    template<typename Index, typename DataMapper, bool ConjugateLhs, bool ConjugateRhs>
168    struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,DataMapper,1,1,ConjugateLhs,ConjugateRhs>
169    {
170      typedef mpfr::mpreal mpreal;
171
172      EIGEN_DONT_INLINE
173      void operator()(const DataMapper& res, const mpreal* blockA, const mpreal* blockB,
174                      Index rows, Index depth, Index cols, const mpreal& alpha,
175                      Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0)
176      {
177        if(rows==0 || cols==0 || depth==0)
178          return;
179
180        mpreal  acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())),
181                tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr()));
182
183        if(strideA==-1) strideA = depth;
184        if(strideB==-1) strideB = depth;
185
186        for(Index i=0; i<rows; ++i)
187        {
188          for(Index j=0; j<cols; ++j)
189          {
190            const mpreal *A = blockA + i*strideA + offsetA;
191            const mpreal *B = blockB + j*strideB + offsetB;
192
193            acc1 = 0;
194            for(Index k=0; k<depth; k++)
195            {
196              mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_srcptr(), B[k].mpfr_srcptr(), mpreal::get_default_rnd());
197              mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(),  mpreal::get_default_rnd());
198            }
199
200            mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_srcptr(), alpha.mpfr_srcptr(), mpreal::get_default_rnd());
201            mpfr_add(res(i,j).mpfr_ptr(), res(i,j).mpfr_srcptr(), acc1.mpfr_srcptr(),  mpreal::get_default_rnd());
202          }
203        }
204      }
205    };
206  } // end namespace internal
207}
208
209#endif // EIGEN_MPREALSUPPORT_MODULE_H
210