Summations zu.mpfr and zstar based on a probability sample
First suggestion would be to write a function for the high precision conversion because you keep using mpfr(..., precBits = precBits) over and over:

Z &lt;- function(x) mpfr(x, precBits = 1000000)


It won't make your code faster but it will save you a lot of typing and make it easier on the eye.

Next, if you do not know why your code is slow, you need to learn to use the profiler:

Rprof(tmp &lt;- tempfile())
[INSERT_YOUR_CODE_HERE]
Rprof()
summaryRprof(tmp)
unlink(tmp)


It will show you that 90% of the time is spent doing the division in your prod = prod * ... line of code. So a solution would be to compute the numerator (num) and denominator (den) by iteration, then do a single division at the end:

zu.mpfr &lt;- function(freq, u) {
  tot &lt;- Z(0)
  n &lt;- sum(freq)
  for (i in seq_along(freq)) {
    num &lt;- Z(1)
    den &lt;- Z(1)
    for (j in (0:(u-1))) {
      num &lt;- num * Z(freq[i] - j)
      den &lt;- den * Z(n - j)
    }
    prod &lt;- num / den
    tot &lt;- tot + prod
  }
  return(tot)
}


This alone makes your code three to four times faster.

Next, I would recommend looking for places where you can vectorize. I understand that your code is dealing with high precision numbers so maybe there is no other way than using for loops if only the simple operators (+, *, /) are available... However, it appears that if x is a vector of mpfr numbers, prod(x) also returns a mpfr number so I am ready to trust it is doing the right thing. If so, then you could do:

zu.mpfr &lt;- function(freq, u) {
  tot &lt;- Z(0)
  n &lt;- sum(freq)
  J &lt;- 0:(u-1)
  for (i in seq_along(freq)) {
    tot &lt;- tot + prod(Z(freq[i] - J)) / prod(Z(n - J))
  }
  return(tot)
}


Have you noticed that we are computing the same denominator over and over? Move it outside the loop to save a little more:

zu.mpfr &lt;- function(freq, u) {
  tot &lt;- Z(0)
  n &lt;- sum(freq)
  J &lt;- 0:(u-1)
  den &lt;- prod(Z(n - J))
  for (i in seq_along(freq)) {
    tot &lt;- tot + prod(Z(freq[i] - J)) / den
  }
  return(tot)
}


Last, if memory size is not an issue (i.e. K remains small), you could save even more time by working with matrices. It's a bit harder to digest (including the computation of row-wise cumulative products), but it's elegant how short your code becomes (notice we don't need zu.mpfr anymore). This version is 5 to 6 times faster than your original one:

  sample &lt;- c(rmultinom(1, n, asNumeric(p)))
  x &lt;- Z(outer(sample, 0:(K-1), FUN = "-"))
  y &lt;- Z(sum(sample) - 0:(K-1))
  z &lt;- apply(x, 1, cumprod) / cumprod(y)
  zstar &lt;- tail(rowSums(z), -1)