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.