Rcpp具有四精度计算

Question

在Rcpp中以数字方式计算以下内容的最佳方法是什么？

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

在许多情况下,计算可能需要多精度(对于exp),但最终结果可以舍入到通常的double.

通过quadmath？通过提升？

如果你留在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