High precision sum with R involving binomial coefficients and logs

Question

I am trying to get the sum of terms obtained from the products of binomial coefficients (very large integers) and logarithms (small reals), each terms having alternating signs.

For example:

library(Rmpfr)
binom <- function(n,i) {factorial(n)/(factorial(n-i)*factorial(i))}

i <- 30
n <-60
Ui <- rep(0,i)
for (k in (0:(i-1))) {
    Ui[k+1] <-  (-1)^(i-1-k) * binom(i-1,k)/(n-k) * log(n-k)
}
U <- sum(mpfr(Ui, 1024))

returns 7.2395....e-10, which is very far from the actual response returned by Mathematica, that is, -5.11...e-20.

How can I make the sum be accurate? I checked manually the Ui and all individually seem accurate to many digits.

Edit

Here is the Mathematica code that computes the same sum. It works on integers and only convert to reals once the sum is over. I increased the number of reported decimals.

Reason for this?

In the end, I need to get the ratio of two numbers obtained with similar computations. When the two numbers are off a couple of order of magnitude, the ratio obtained is just simply unpredictable.

Answer 1

You need to work with the mpfr objects during the whole calculation, rather than just at the summation:

library(Rmpfr)

i <- 30
n <- 60
k <- 0:(i - 1)

nk <- mpfr(n - k, 128)
(U <- sum((-1)^(i-1-k)*choose(i-1,k)/(nk)*log(nk)))
#> 1 'mpfr' number of precision  128   bits 
#> [1] -5.110333215290518581300810256453669394729e-20

nk <- mpfr(n - k, 256)
(U <- sum((-1)^(i-1-k)*choose(i-1,k)/(nk)*log(nk)))
#> 1 'mpfr' number of precision  256   bits 
#> [1] -5.110333215285320173235309727002720346864555872897902728222060861935229197560667e-20

nk <- mpfr(n - k, 512)
(U <- sum((-1)^(i-1-k)*choose(i-1,k)/(nk)*log(nk)))
#> 1 'mpfr' number of precision  512   bits 
#> [1] -5.1103332152853201732353097270027203468645558728979134452318939958128833820370490135678222208577277855238767473116630391351888405531035522832949562601913591e-20