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  * 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>
38using namespace mpfr;
39using namespace Eigen;
40int 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