1

I am using R, and I have a fraction that is a complex number:

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

Notice, I am calling erfz ... the complex error function. This is the source from the pracma library for that function:

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


You will notice that the N and D for each is getting large, but the fraction N/D is staying relatively manageable.

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

At some point, R breaks:

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

Now it breaks:

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

So I believe I could plug into some C++ libraries that could perform this calculation. The numerator is likely easier that the denominator because of the complex error function.

I have included standalone C++ code within R before, using this link as a guide: fast large matrix multiplication in R

Something like:

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

I am hoping for some help on the Numerator, which just seems like a large number, so a library that allows that.

For the Denominator, also handling of large numbers, but with access to the erfz function in a C++ library.

Then I can source it using Rcpp

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

How to multiply (and take exponent) of large complex numbers in the numerator? How to multiple and apply the complex error function in the denominator?

Maybe a function for each?


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

I am familiar with gmp https://gmplib.org/ but would prefer not having to implement that entire library. Nor do I know if it will meet my needs.

Dharman
  • 30,962
  • 25
  • 85
  • 135
mshaffer
  • 959
  • 1
  • 9
  • 19
  • Are you aware that [GMP](https://gmplib.org/) has R packages such as [Rmpfr](https://cran.r-project.org/web/packages/Rmpfr/index.html) and [gmp](https://cran.r-project.org/web/packages/gmp/index.html)? You could look at those and their reverse depends. – Dirk Eddelbuettel Dec 03 '20 at 13:38
  • ```> si = mpfr(2+3i, precBits=80); Error in mpfr.default(2 + (0+3i), precBits = 80) : invalid 'x'. Must be numeric (logical, raw) or character``` This is library(Rmpfr) when I try and load a complex number. `si = mpfr(2, precBits=80);` [no complexity] works fine. – mshaffer Dec 04 '20 at 16:50
  • Sorry, was worth a try. Next up, maybe check Boost via package BH? – Dirk Eddelbuettel Dec 04 '20 at 16:54

0 Answers0