How to put Rmpfr values into a function in R?

Question

I am calculating the inverse of a Vandermonde Matrix. I have written the codes to calculate the inverse explicitly by its formula as below:

library(gtools)

#input is the generation vector of terms of Vandermonde matrix.
FMinv <- function(base){
  n=length(base)
  inv=matrix(nrow=n,ncol=n)
  for (i in 1:n){
    for (j in 1:n){
      if(j<n){
        a=as.matrix(combinations(n,n-j,repeats.allowed = F))
        arow.tmp=nrow(a) #this is in fact a[,1]
        b=which(a==i)%%length(a[,1])
        nrowdel=length(b)
        b=replace(b,b==0,length(a[,1]))
        a=a[-b,]
        if(arow.tmp-nrowdel>1){
          a=as.matrix(a)
          nrowa=nrow(a)
          prod=vector()
          for(k in 1:nrowa){
            prod[k]=prod(base[a[k,]])
          }
          num=sum(prod)
        }
        if(arow.tmp-nrowdel==1){
          num=prod(base[a])
        }
        den=base[i]*prod(base[-i]-base[i])
        inv[i,j]=(-1)^(j-1)*num/den
      }
      if(j==n){
        inv[i,j]=1/(base[i]*prod(base[i]-base[-i]))
      }
    }
  }
  return(inv)
}

And I define a base as follow:

> library(Rmpfr)
> a=mpfr(c(10:1),1000)/Rmpfr::mpfr(sum(1:10),1000)
> a
10 'mpfr' numbers of precision  1000   bits 
 [1]  0.18181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181819
 [2]  0.16363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363637
 [3]  0.14545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545456
 [4]  0.12727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727274
 [5]  0.10909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909091
 [6] 0.090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909094
 [7] 0.072727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727278
 [8] 0.054545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545455
[9] 0.036363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363639
[10] 0.018181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181819

However, when I was attempting to put the "a" into the function, I got:

> FMinv(a)
Error in sum(prod) : invalid 'type' (list) of argument

By checking its type,

> typeof(a)
[1] "list"

The only thing that I know to transform it to values is asNumeric() in Rmpfr. However,

> asNumeric(a)
 [1] 0.18181818 0.16363636 0.14545455 0.12727273 0.10909091 0.09090909 0.07272727 0.05454545 0.03636364 0.01818182

it lost the remaining digits.

Is there anyway to put the "a" into my function without losing decimals?

Thanks!

Answer 1

The trick is to use S3 methods.
Define a generic, a default method to be called with your "normal" numbers, meaning, objects of class "numeric" and the function the question is asking for.

That is the problem one. It took a while but I believe the code below is right.

library(OBsMD)

FMinv <- function(...) UseMethod("FMinv")

FMinv.default <- function(base) {
    # Your function
    # unchanged
}



FMinv.mpfr <- function(base, precBits = getPrec(base)) {
    n <- length(base)
    inv <- mpfr(rep(0, n*n), precBits = precBits)
    inv <- matrix(inv, nrow = n, ncol = n)
    for (i in 1:n) {
        for (j in 1:n) {
            if (j < n) {
                a <- combinations(n, n - j, repeats.allowed = F)
                a <- as.matrix(a)
                arow.tmp <- nrow(a)  # this is in fact a[, 1]
                b <- which(a == i) %% length(a[, 1])
                nrowdel <- length(b)
                b <- replace(b, b == 0, length(a[, 1]))
                a <- a[-b, ]
                num <- mpfr(0, precBits[1])
                if (arow.tmp - nrowdel > 1) {
                  a <- as.matrix(a)
                  nrowa <- nrow(a)
                  for (k in 1:nrowa) {
                    num <- num + prod(base[a[k, ]])
                  }
                }
                if (arow.tmp - nrowdel == 1) {
                  num <- num + prod(base[a])
                }
                den <- base[i] * prod(base[-i] - base[i])
                inv[i, j] <- (-1)^(j - 1) * num/den
            }
            if (j == n) {
                inv[i, j] <- 1/(base[i] * prod(base[i] - base[-i]))
            }
        }
    }
    return(inv)
}

Now test both methods and compare some of the results' values.

library(Rmpfr)

a <- mpfr(c(10:1),1000)/Rmpfr::mpfr(sum(1:10),1000)

inv1 <- FMinv(asNumeric(a))
inv2 <- FMinv(a)

inv1[10, 10]
#[1] -6.98014e+11

inv2[10, 10]
#1 'mpfr' number of precision  1000   bits 
#[1] -698013564040.84166942239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906474