具有四精度计算的 Rcpp [英] Rcpp with quad precision computation

问题描述

在 Rcpp 中数值计算如下内容的最佳方法是什么?

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

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

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

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

通过quadmath?通过升压?

Via quadmath? Via boost?

如果你留在 R(Rcpp 之外)，有非常舒适的包装器来完成这项工作:

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))

但是如何使用 rcpp 访问?

But how to access with rcpp?

推荐答案

一种快速启动和运行的方法是安装 BH 包并使用 Boost Multiprecision 库，它提供了 几种扩展精度浮点类型.例如，这段代码演示了 float128 和 mpf_float_100 类型:

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

<小时>

后者需要在编译前添加用于链接到libmpfr 和libgmp 的标志；前者没有，因为它是 GCC 内置 __float128 的包装器:

---

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"

与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的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！