Using boost multiprecision with RcppEigen

I have written a library using RcppArmadillo. The problem I have is that for some parameters, the arma::solve function does not give me the exact solution, since the "A" matrix has an rcond close to 0. If arma::solve could solve the linear equation exactly, that would not be a problem. But it gives me an approximate solution, which is not good enough for me.

Then, I have thought in using RcppEigen, and use boost multiprecision variables. If I understand correctly Eigen, the Eigen solver will give me a solution in multiprecision, and it is quite likely this solution will be good enough for me (even with a float128).

But when I try to implement this plan, I have an error, and I do not know what to do. For example, the following works:

#include <RcppArmadillo.h>
#include <RcppEigen.h>
#include <cmath>
#include <cstdint>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/eigen.hpp>

// Correctly setup the build environment 
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(BH)]]
using namespace Rcpp;
using namespace arma;
using namespace Eigen;
using namespace boost;
using namespace boost::multiprecision;
using Eigen::Map;
using Eigen::MatrixXd;                  // variable size matrix, double precision
using Eigen::VectorXd;
using Eigen::Matrix;
using Eigen::Dynamic;
namespace mp = boost::multiprecision;


    // [[Rcpp::export]]
    Eigen::MatrixXd onesfgh_LPPLS_RcppEigen(int t2, int t1, double tc, double m, double w) {
      //Definitions & pre-computations:
      unsigned int tdim;
      double ttc, powttc, wlog, sinwlog, coswlog;
      tdim = t2 - t1 + 1;
      ttc = tc - t1 + 1;
      Eigen::MatrixXd output(tdim, 4); //first dimension is time, second is onesfgh
      //Main calculations:
      for (unsigned int i = 0; i < tdim; i++) {
        output(i, 0) = 1.0;
        ttc--;
        powttc = pow(ttc, m);
        wlog = w * log(ttc);
        sinwlog = sin(wlog);
        coswlog = cos(wlog);
        //sincos(wlog, &sinwlog, &coswlog);
        output(i, 1) = powttc;
        output(i, 2) = coswlog * powttc;
        output(i, 3) = sinwlog * powttc;
      }
      return output;
    }

But:

// [[Rcpp::export]]
Eigen::Matrix<boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off>, Eigen::Dynamic, 4> onesfgh_LPPLS_RcppEigen(int t2, int t1, double tc, double m, double w) {
  //Definitions & pre-computations:
  unsigned int tdim;
  boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off> tc128(tc), m128(m), w128(w);
  boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off> ttc, powttc, wlog, sinwlog, coswlog;
  boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off> t1128(double(t1));
  tdim = t2 - t1 + 1;
  ttc = tc128 - t1128 + 1.0;
  Eigen::Matrix<boost::multiprecision::number<boost::multiprecision::float128, boost::multiprecision::et_off>, Eigen::Dynamic, 4> output(tdim, 4); //first dimension is time, second is onesfgh
  //Main calculations:
  for (unsigned int i = 0; i < tdim; i++) {
    output(i, 0) = 1.0;
    ttc--;
    powttc = boost::multiprecision::pow(ttc, m128);
    wlog = w128 * boost::multiprecision::log(ttc);
    sinwlog = boost::multiprecision::sin(wlog);
    coswlog = boost::multiprecision::cos(wlog);
    //sincos(wlog, &sinwlog, &coswlog);
    output(i, 1) = powttc;
    output(i, 2) = coswlog * powttc;
    output(i, 3) = sinwlog * powttc;
  }
  return output;
}

it does not work. I get the error message:

RcppExports.cpp:11:15: error: 'boost' was not declared in this scope

In particular, I have tried the advice given in https://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/math_toolkit/high_precision/using_test.html but without using a typedef, since in another post, Dirk Eddelbuettel suggests typedefs are not easy to deal with (and they should be place in another .h file).

Is there any suggestion how to proceed?

EDIT:

I have modified the function, by not exporting. Now the code is:

  #include <RcppArmadillo.h>
    #include <RcppEigen.h>
    #include <cmath>
    #include <cstdint>
    #include <boost/multiprecision/float128.hpp>
    #include <boost/multiprecision/eigen.hpp>

    // Correctly setup the build environment 
    // [[Rcpp::depends(RcppArmadillo)]]
    // [[Rcpp::depends(RcppEigen)]]
    // [[Rcpp::depends(BH)]]
    using namespace Rcpp;
    using namespace arma;
    using namespace Eigen;
    using namespace boost;
    using namespace boost::multiprecision;
    using Eigen::Map;
    using Eigen::MatrixXd;                  
    using Eigen::VectorXd;
    using Eigen::Matrix;
    using Eigen::Dynamic;
    namespace mp = boost::multiprecision;

    Eigen::Matrix<mp::float128, Eigen::Dynamic, 4> onesfgh_LPPLS_RcppEigen(int `t2, int t1, double tc, double m, double w) {`
      //Definitions & pre-computations:
      unsigned int tdim;
      mp::float128 tc128(tc), m128(m), w128(w);
      mp::float128 ttc, powttc, wlog, sinwlog, coswlog;
      mp::float128 t1128(double(t1));
      tdim = t2 - t1 + 1;
      ttc = tc128 - t1128 + 1.0;
      Eigen::Matrix<mp::float128, Eigen::Dynamic, 4> output(tdim, 4); 
      //Main calculations:
      for (unsigned int i = 0; i < tdim; i++) {
        output(i, 0) = 1.0;
        ttc--;
        powttc = mp::pow(ttc, m128);
        wlog = w128 * mp::log(ttc);
        sinwlog = mp::sin(wlog);
        coswlog = mp::cos(wlog);
        //sincos(wlog, &sinwlog, &coswlog);
        output(i, 1) = powttc;
        output(i, 2) = coswlog * powttc;
        output(i, 3) = sinwlog * powttc;
      }
      return output;
    }

And now I get the following error message:

babel_RcppArmadillo.cpp:6:42: fatal error: boost/multiprecision/eigen.hpp: 
No such file or directory
#include <boost/multiprecision/eigen.hpp>                                          
compilation terminated.

Which is odd, since I assume that BH includes the subfolders of boost.

Solution

Update: This should be unnecessary by now, since the BH package includes the required header since January 2019.

The boost/multiprecision/eigen.hpp header was added in version 1.68, while the BH package currently provides boost 1.66. You have to install an updated boost separately. The following should work but is untested:

Download boost 1.68 and unpack it into some suitable directory. On Linux and other Unix-like systems I would probably use /usr/local/include. Otherwise I would use any path without spaces in the name.
Do not depend on the BH package, i.e. remove // [[Rcpp::depends(BH)]] for sourceCpp et al. or Imports: BH for package use.
If boost was installed in a non-standard location in step 1., you have to tell the compiler to look there. This can be done with
```
PKG_CPPFLAGS=-I<path-to-boost>
```
either via Sys.setenv for sourceCpp et al. or within src/Makevars(.win) for package use.