c++ - 使用 Rcpp，如何计算两个大复数的比率

Question

我正在使用 R，并且我有一个复数分数：

library(pracma);
sigma.jj = 1;
si = 2+3i;
N = 2*exp( 0.5*(- (sigma.jj^2)*(si^2)) );   # 23.39454+6.80796i
D = 1-pracma::erfz( (2^(-0.5))* ( si*sigma.jj )); # 2.225013-1.597002i
N/D;  # 5.489974+7.000163i

请注意，我正在调用erfz...复杂的错误函数。这是该pracma函数库的源代码：

# Complex error function
erfz    <- function(z)
{
    if (is.null(z)) return( NULL )
    else if (!is.numeric(z) && !is.complex(z))
        stop("Argument 'z' must be a numeric or complex scalar or vector.")

    a0 <- abs(z)
    c0 <- exp(-z * z)
    z1 <- ifelse (Re(z) < 0, -z, z) 

    i <- a0 <= 5.8
    work.i <- i
    cer <- rep(NA, length = length(z))
    if ( sum(work.i) > 0) {
        cs <- z1
        cr <- cs
        for (k in 1:120) {
            cr[work.i] <- cr[work.i] * z1[work.i] * z1[work.i]/(k + 0.5)
            cs[work.i] <- cs[work.i] + cr[work.i]
            work.i <- work.i & (abs(cr/cs) >= 1e-15)
        if (sum(work.i) == 0) break
        }
        cer[i] <- 2 * c0[i] * cs[i]/sqrt(pi)
    }
    work.i <- !i
    if( sum(work.i) > 0) {
        cl <- 1/z1
        cr <- cl
        for (k in 1:13) {
            cr[work.i] <- -cr[work.i] * (k - 0.5)/(z1[work.i] * z1[work.i])
            cl[work.i] <-  cl[work.i] + cr[work.i]
            work.i <- work.i & (abs(cr/cl) >= 1e-15)
        if (sum(work.i) == 0) break
        }
        cer[!i] <- 1 - c0[!i] * cl[!i]/sqrt(pi)
    }
    cer[ Re(z) < 0] <- -cer[ Re(z) < 0]
    return(cer)
}

---

您会注意到每个的 N 和 D 都在变大，但比例N/D保持相对可控。

sigma.jj = 1;
si = 2+13i;
N = 2*exp( 0.5*(- (sigma.jj^2)*(si^2)) );   # 8.73323e+35-1.029433e+36i
D = 1-pracma::erfz( (2^(-0.5))* ( si*sigma.jj )); # -2.692745e+34-3.114988e+34i
N/D;  # 5.04326+32.39578i

在某些时候，R 中断：

sigma.jj = 1;
si = 2+33i;
N = 2*exp( 0.5*(- (sigma.jj^2)*(si^2)) );   # -8.046987e+235+2.13732e+234i
D = 1-pracma::erfz( (2^(-0.5))* ( si*sigma.jj )); # -3.31367e+232+9.716936e+233i
N/D;  # 5.01787+82.64292i

现在它打破了：

sigma.jj = 1;
si = 2+38i;
N = 2*exp( 0.5*(- (sigma.jj^2)*(si^2)) );   # Inf-Infi
D = 1-pracma::erfz( (2^(-0.5))* ( si*sigma.jj )); # NaN-Infi
N/D;  # NaN+NaNi

所以我相信我可以插入一些可以执行此计算的 C++ 库。由于复杂的误差函数，分子可能比分母更容易。

我之前在 R 中包含了独立的 C++ 代码，使用此链接作为指南：fast large matrix multiplication in R

就像是：

// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]

#include <RcppArmadillo.h>
#include <RcppEigen.h>

// [[Rcpp::export]]
SEXP armaMatMult(arma::mat A, arma::mat B){
    arma::mat C = A * B;

    return Rcpp::wrap(C);
}

// [[Rcpp::export]]
SEXP eigenMatTrans(Eigen::MatrixXd A){
    Eigen::MatrixXd C = A.transpose();

    return Rcpp::wrap(C);
}

// [[Rcpp::export]]
SEXP eigenMatMult(Eigen::MatrixXd A, Eigen::MatrixXd B){
    Eigen::MatrixXd C = A * B;

    return Rcpp::wrap(C);
}

// [[Rcpp::export]]
SEXP eigenMapMatMult(const Eigen::Map<Eigen::MatrixXd> A, Eigen::Map<Eigen::MatrixXd> B){
    Eigen::MatrixXd C = A * B;

    return Rcpp::wrap(C);
}

我希望对 Numerator 有所帮助，这似乎是一个很大的数字，所以一个允许这样做的库。

对于分母，也可以处理大数，但可以访问erfzC++ 库中的函数。

然后我可以使用 Rcpp 获取它

library(Rcpp)
sourceCpp("multiply.cpp")

如何将分子中的大复数相乘（并取指数）？如何在分母中多重化和应用复杂的误差函数？

也许每个功能？

---

https://www.tutorialspoint.com/handling-large-numbers-in-cplusplus

https://www.tutorialspoint.com/cpp_standard_library/complex.htm

https://en.cppreference.com/w/cpp/numeric/math/erf

我熟悉 gmp https://gmplib.org/，但不希望不必实现整个库。我也不知道它是否能满足我的需要。

标签： c++rrcppcomplex-numberslargenumber