Rcpp with quad precision computation

Question

What is the best way to numerically compute things like the following within Rcpp?

exp(-1500)/(exp(-1500)+exp(-1501))

In many cases, the computation may require multi precision (for the exp), but the end result can be rounded to a usual double.

Via quadmath? Via boost?

If you stay in R (outside of Rcpp), there are really comfy wrappers to do the job:

library(Rmpfr)

a = mpfr(-1500,100)
b = mpfr(-1501,100)

exp(a)/(exp(a)+exp(b))

But how to access with rcpp?

Why not just compute the same thing in a way that *doesn't* require such an extreme exponent range? — user2357112, Feb 16 '17 at 01:03
For example, here, you could just add 1500 to all the arguments to `exp`. — user2357112, Feb 16 '17 at 01:06
this is just a minimal example ... I would like to apply in situations where I don't know what the numbers are. — Inferrator, Feb 16 '17 at 01:07
Then *compute* the adjustment! Take the maximum argument (-1500 here) and subtract it from all the arguments. — user2357112, Feb 16 '17 at 02:29
I can also do series expansion, but that wasn't the question — Inferrator, Feb 16 '17 at 02:30

score 8 · Accepted Answer · answered Feb 17 '17 at 20:56

A quick way to get up and running is to install the BH package and use the Boost Multiprecision library, which provides several extended precision floating point types. For example, this code demos the float128 and mpf_float_100 types:

// [[Rcpp::depends(BH)]]
#include <Rcpp.h>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/mpfr.hpp>

namespace mp = boost::multiprecision;

// [[Rcpp::export]]
std::string qexp(double da = -1500.0, double db = -1501.0)
{
    mp::float128 a(da), b(db);
    mp::float128 res = mp::exp(a) / (mp::exp(a) + mp::exp(b));
    return res.convert_to<std::string>();
}

// [[Rcpp::export]]
std::string mpfr_exp(double da = -1500.0, double db = -1501.0)
{
    mp::mpf_float_100 a(da), b(db);
    mp::mpf_float_100 res = mp::exp(a) / (mp::exp(a) + mp::exp(b));
    return res.convert_to<std::string>();
}

The latter requires adding flags for linking to the libmpfr and libgmp before compiling; the former does not, as it is a wrapper around GCC's built in __float128:

Sys.setenv("PKG_LIBS" = "-lmpfr -lgmp")
Rcpp::sourceCpp('/tmp/quadexp.cpp')

qexp()
# [1] "0.731058578630004879251159241821836351"

mpfr_exp()
# [1] "0.731058578630004879251159241821836274365144640165056519276365907919040453070204639387474532075981245292174466493140773"

Comparing to Rmpfr,

library(Rmpfr)

a <- mpfr(-1500, 100)
b <- mpfr(-1501, 100)

exp(a) / (exp(a) + exp(b))
# 1 'mpfr' number of precision  100   bits 
# [1] 7.3105857863000487925115924182206e-1

Dirk Eddelbuettel · Answer 2 · 2019-01-04T21:22:06.587

1

I fear you may be mostly confused.

Rcpp is a bridge between R and C++. It is neither a new language, nor a new system.

And in C++ you use something like GNU gmplib aka The GNU Multiple Precision Arithmetic Library which, of course, has already been wrapped for R as package gmp

But if you want to write C++ functions, called from R, interfacing the GMP library -- you certainly. It falls back to 'how do I write a C++ function callable from R which accesses an external library' -- and that question has been covered numerous times here too.

edited Jan 04 '19 at 21:22

answered Feb 16 '17 at 01:58

Dirk Eddelbuettel

360,940
56
644
725

1

I don't feel confused, I was just looking for a minimal example. I'll just take a look the gmp sources then, thanks. – Inferrator Feb 16 '17 at 02:39
It's just like anything else requiring an external library: entirely doable, but you need to ensure the library is present on all relevant OSs. – Dirk Eddelbuettel Feb 16 '17 at 02:57
For the record, the accepted answer states exactly what I said: _an external library is needed_. It also kind enough to show one. – Dirk Eddelbuettel Jan 04 '19 at 21:23

Rcpp with quad precision computation

2 Answers2

Linked