Here is my code:

library(gtools)
library(Rmpfr)
library(OBsMD)
precBits=1000000
K=9
n=1000
p=mpfr(c(K:1),precBits = precBits);p=mpfr(p,precBits = precBits)/mpfr(sum(p),precBits = precBits)

zu.mpfr <- function(freq, u, precBits){
  sum=mpfr(0,precBits = precBits)
  n=sum(freq)
  for (i in (1:length(freq))){
    prod=mpfr(1,precBits = precBits)
    for (j in (0:(u-1))){
      prod=prod*(mpfr(freq[i]-j,precBits = precBits))/mpfr(n-j,precBits=precBits)
    }
    sum=sum+prod
  }
  return(sum)
}

etastar<-rep(mpfr(0,precBits = precBits),K-1)
etastar=c(etastar)
for (i in 1:(K-1)){
  etastar[i]=sum(p^mpfr(i+1,precBits))
}

Following is the part that I need to make it faster:

start.time<-Sys.time()
sample=rmultinom(1,n,asNumeric(p))
zstar<-rep(mpfr(0,precBits = precBits),K-1)
zstar=c(zstar)
for(i in 1:(K-1)){
  zstar[i]=zu.mpfr(sample,i+1,precBits)
}
  
end_time <- Sys.time()
end_time-start.time

The timing part needs "Time difference of 6.207007 secs" on my end. I am wondering if I can make it better. I need it faster because I need to run the simulation with thousands of iterations. I am not sure if the slowness is caused by mpfr but I do need it to enforce accuracy.