Brewer method unequal probabilities speed up

up_brewer.c

#include <R.h>
#include <Rmath.h>

void C_up_brewer(double* pik, int* N, double* sb) {
    double n = 0.0;
    double a;
    double u;
    int j;

    // Initialize sb and calculate sum of pik (n)
    for (int i = 0; i < *N; i++) {
	sb[i] = 0;
        n += pik[i];
    }

    // Random seed
    GetRNGstate();

    // Main loop
    for (int i = 1; i <= (int)n; i++) {
        double p_sum = 0.0;
        double* p = (double*) malloc(*N * sizeof(double));

        // Calculate a
        a = 0.0;
        for (int k = 0; k < *N; k++) {
            a += pik[k] * sb[k];
        }

        // Calculate p
        for (int k = 0; k < *N; k++) {
	    p[k] = (1 - sb[k]) * pik[k] * ((n - a) - pik[k]) / ((n - a) - pik[k] * (n - i + 1));
            p_sum += p[k];
        }

        // Generate a random number from a uniform distribution
        u = runif(0, 1);

        // Find the first index j such that p[j] - u > 0
        double cumsum = 0.0;
        for (j = 0; j < *N; j++) {
            cumsum += p[j]/p_sum;
            if (cumsum - u > 0) {
                break;
            }
        }

        // Update sb
        sb[j] = 1;

        // Free the dynamically allocated memory for p
        free(p);
    }
    PutRNGstate();
}

up_brewer.R

up_brewer <- function(pik, eps = 1e-06) {
  if (any(is.na(pik)))
    stop("there are missing values in the pik vector")
  list <- pik > eps & pik < 1 - eps
  pikb <- pik[list]
  s <- pik
  N <- as.integer(length(pikb))
  sb <- vector(mode = "double", length = N)
  r <- .C("C_up_brewer", pikb, N, result = sb)
  s[list] <- r$result
  s
}

In linux, you can compile up_brewer.c from R using

system("rm up_brewer.o up_brewer.so")
system("R CMD SHLIB -lRmath up_brewer.c")

And don't forget to load .so file before running up_brewer function in R.
dyn.load("up_brewer.so")