Computation in R

Question

I am trying to find the following expression for $n=50, m=40$ in the software R.

\begin{equation} m \binom{n}{m} \sum_{i=0}^{m-1}\frac{ \binom{m-1}{i}(-1)^i }{n-m+1+i} \end{equation} Ideally this expression should be evaluated to be 1 (in fact this expression gives value 1 for all $m <n$ and that can be proved theoretically). But instead it is giving 166254.90890021951. I understand it is precision problem. How do i get over it ?

Answer 1

My attempt:

n = 50;  m = 40;  k = m*choose(n,m);  i = 0:(m-1)
k*sum(choose(m-1, i)*(-1)^i/(n- m + 1 + i))
## -16608.58

There may be issues of computational stability because $k = 410,891,126,800$ and the sum computes to $-4.042088e-08.$

Trying the computation for smaller values of $n$ and $m$, all seems OK.

n = 5; m = 4; k =m*choose(n,m);  i = 0:(m-1)
k*sum(choose(m-1, i)*(-1)^i/(n-m+1 + i))
## 1

n = 10; m = 8; k =m*choose(n,m);  i = 0:(m-1)
k*sum(choose(m-1, i)*(-1)^i/(n-m+1 + i))
## 1

So you seem to be at (or over) the edge of overflow for the case with $n = 50, m = 40.$

Try using logs for the constant $k$ and for each term in the sum (exponentiating before summing). For example, you could use lchoose for parts of the computation:

choose(50, 40)
## 10272278170
lchoose(50, 40)
## 23.05271
exp(lchoose(50,40))
## 10272278170

Note: I do not share @Donald's reservations about R. However, in any computation in any computer language, you need to realize that you are not dealing with all the real numbers; you are dealing with a huge, but finite set of the rationals. Consequently you need to be aware of possibility of overflow and underflow in intermediate steps. This is especially true when dealing with factorials, binomial coefficients, Gamma functions, and exponentials.

Answer 2

If you want to calculate this directly as written, you'll need arbitrary precision via a library in R.

f <- function(n,m) {
  require(gmp)
  i <- 0:(m-1)
  m*as.bigq(chooseZ(n,m)) * sum(chooseZ(m-1, i)*(-1)^i / (n-m+1+i))
}

f(50,40)
# Big Rational ('bigq') :
# [1] 1