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.