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.
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