Rcpp с вычислением четверной точности

Question

Каков наилучший способ для численного вычисления таких вещей, как следующие в Rcpp?

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

Во многих случаях вычисление может потребовать многократной точности (для опыта), но конечный результат может быть округлен до обычного двойного.

Через четверть? Через буст?

Если вы остаетесь в R (за пределами Rcpp), есть действительно удобные обертки для этой работы:

library(Rmpfr)

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

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

Но как получить доступ с помощью rcpp?

Answer 1

Быстрый способ начать работу - это установить BH упакуйте и используйте библиотеку Boost Multiprecision, которая предоставляет несколько типов с плавающей запятой повышенной точности. Например, этот код демонстрирует float128 а также mpf_float_100 типы:

// [[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>();
}

---

Последнее требует добавления флагов для ссылки на libmpfr а также libgmp перед компиляцией; первый не делает, так как это обертка вокруг встроенного в GCC __float128:

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

qexp()
# [1] "0.731058578630004879251159241821836351"

mpfr_exp()
# [1] "0.731058578630004879251159241821836274365144640165056519276365907919040453070204639387474532075981245292174466493140773"

По сравнению с 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

Answer 2

Боюсь, вы можете быть в основном сбиты с толку.

Rcpp - это мост между R и C++. Это не новый язык и не новая система.

А в C++ вы используете что-то вроде GNU gmplib, иначе говоря, арифметическая библиотека GNU Multiple Precision, которая, конечно, уже была упакована для R как пакет gmp.

Но если вы хотите написать функции C++, вызываемые из R, взаимодействующие с библиотекой GMP - вы, конечно. Возникает вопрос: "Как написать функцию C++, вызываемую из R, которая обращается к внешней библиотеке" - и этот вопрос также неоднократно освещался здесь.