具有四精度计算的 Rcpp

Question

在 Rcpp 中对以下内容进行数值计算的最佳方法是什么？

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

在许多情况下，计算可能需要多精度（对于 exp），但最终结果可以四舍五入为通常的双精度。

通过 quadmath？通过升压？

如果你留在 R 中（在 Rcpp 之外），确实有很舒服的包装器来完成这项工作：

library(Rmpfr)

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

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

但是如何用rcpp访问呢？

Answer 1

恐怕你会很困惑。

Rcpp 是 R 和 C++ 之间的桥梁。它既不是新语言，也不是新系统。

而在 C++ 中，您使用类似 GNU gmplib aka The GNU Multiple Precision Arithmetic Library which, of course, has already been wrapped for R as package gmp

的东西

但是，如果您想编写从 R 调用的 C++ 函数，连接 GMP 库——您当然可以。它回落到 'how do I write a C++ function callable from R which accesses an external library'——这个问题也在这里被多次提及。

Answer 2

快速起床和运行的方法是安装BH包并使用Boost Multiprecision库，它提供several extended precision floating point types. For example, this code demos the float128 and 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

具有四精度计算的 Rcpp

Rcpp with quad precision computation

r

rcpp