Random number generator that returns unique 64-bit numbers in sorted order

Question

(I have asked the question at stackoverflow.com here but maybe it's better to ask here for the statistical part. Feel free to correct my question, fix the tags, redirect me...)

I need a generator for many (up to one trillion, $10^{12}$) unique, pseudo random 64-bit numbers, that is range $[0..2^{64})$. It should be as "fair" as possible, have the same statistical distribution as if generated using a regular pseudo-RNG and then sorted. The problem is that sorting $10^{12}$ numbers is slow. The use case is replicating a test that was run for BBHash (in the paper, 4.5 Indexing a trillion keys).

The simplest (wrong, in my view) solution, which is not random, is to generate an evenly spaced sequence like 0, 1000, 2000,... . Looking at the source code, this seems to be about what the BBHash test does. Another solution is to iterate over all numbers, and with probability about ${10^{12}}/{2^{64}}$ return the number, otherwise skip it. But this is too slow, and will not return exactly that many results (and adjusting the probability while you go sounds like the distribution would not be fully correct).

I think there are three approaches:

Random Gaps

Maybe the solution in "An Efficient Algorithm for Sequential Random Sampling" would work, but it is slow because for each iteration, many floating point numbers need to be calculated. Also, the sum of such gaps, after $10^{12}$ , may not land where it should due to rounding (I wouldn't want to try that).

Randomly sized ranges with fixed number of entries

With a random number generator, split the range $[0..2^{64})$ into one million randomly sized sub-ranges that each should contain one million entries. But how, and which distribution to use. Then generate one million numbers in that range using a pseudo-RNG.

Random number of entries for each fixed (sub-)range

Create a tree: For each bit level (starting with the most significant bit), recursively generate a random number of how many entries should have the bit at that level set to 0, using the normal distribution. The remaining entries have the bit at this level set to 1. At each recursion level, this will narrow the range by about half. Stop for example when there are less than 1 million entries, and then switch to using an in-memory pseudo-RNG and sort those numbers (or use a bit field).

I can calculate the probability that, for $x$ coin flips, more than $y$ are heads, using the normal distribution (with a CDF approximation). And then pick a random number and do a binary search. Is that the correct, best (easiest, fastest) way? Here is the code that I use so far (Java code, but the algorithm and formulas should be relatively easy to understand even if you don't know Java):

public static void main(String... args) {
    Random r = new Random();
    Iterator<Long> it = randomSequence(r, 10, 32);
    while(it.hasNext()) {
        System.out.println(it.next());
    }
}

/**
 * Random sequence generator.
 *
 * @param r the random generator
 * @param size the number of entries to generate
 * @param shift the number of bits of the result
 * @return the iterator
 */
static Iterator<Long> randomSequence(final Random r, final long size, final int shift) {
    if (size < 5) {
        // small lists are generated using a regular hash set
        TreeSet<Long> set = new TreeSet<Long>();
        while (set.size() < size) {
            set.add(r.nextLong() & ((2L << shift) - 1));
        }
        return set.iterator();
    }
    // large lists are created recursively
    return new Iterator<Long>() {
        long remaining = size, zeros = randomHalf(r, size);
        Iterator<Long> lowBits0 = randomSequence(r, zeros, shift - 1);
        Iterator<Long> lowBits1;
        @Override
        public boolean hasNext() {
            return remaining > 0;
        }
        @Override
        public Long next() {
            remaining--;
            if (lowBits0.hasNext()) {
                return lowBits0.next();
            }
            if (lowBits1 == null) {
                lowBits1 = randomSequence(r, size - zeros, shift - 1);
            }
            return (1L << shift) + lowBits1.next();
        }
    };
}

/**
 * Get the number of entries that are supposed to be below the half,
 * according to the probability theory. For example, for a number of coin
 * flips, how many are heads.
 *
 * @param r the random generator
 * @param samples the total number of entries
 * @return the number of entries that should be used for one half
 */
static long randomHalf(Random r, long samples) {
    long low = 0, high = samples;
    double x = r.nextDouble();
    while (low + 1 < high) {
        long mid = (low + high) / 2;
        double p = probabilityBucketAtMost(samples, mid);
        if (x > p) {
            low = mid;
        } else {
            high = mid;
        }
    }
    return (low + high) / 2;
}

static double probabilityBucketAtMost(long flips, long heads) {
    // https://www.fourmilab.ch/rpkp/experiments/statistics.html
    long x = heads;
    long n = flips;
    double variance = Math.sqrt(n/4);
    // mean
    long mu = n / 2;
    // https://en.wikipedia.org/wiki/Normal_distribution
    // Numerical approximations for the normal CDF
    // the probability that the value of a standard normal random variable X is <= x
    return phi((x - mu) / variance);
}

static double phi(double x) {
    return 0.5 * (1 + Math.signum(x) * Math.sqrt(1 - Math.exp(-2 * x * x / Math.PI)));
}

Answer 1

I doubt you need such precision, but it's an interesting challenge to achieve it anyway, especially with minimal additional software. Can these trillion 64-bit integers be efficiently created without an extended arithmetic library, for instance?

Divide-and-conquer is a good strategy, provided it ultimately results in massively parallelizable code. (Waiting for a single CPU to output a trillion numbers can be frustrating...) The trick will be to divide the problem up in an efficient way.

One method suggested by these considerations is to bin the population $1,2,\ldots, 2^{64}=N$ into $n$ equal-size bins. We might as well take $n$ to be a power of two. A good choice (for double-precision floating point arithmetic) is $n \approx 2^{24}$. The numbers of elements in the sample $x$ falling within the bins constitute one draw from a multinomial distribution: these are the counts you would get if you put $10^{12}$ balls randomly into $n$ bins.

I will discuss drawing from this (huge) multinomial distribution in a bit. For now, let $k=(k_1,k_2, \ldots, k_n)$ be the vector of these multinomial counts. For each $i$, draw $k_i$ values from the set $1, 2, \ldots, N/n$ without replacement. Their quantity $N/n$ is sufficiently small that we don't have to be clever about this: just draw them and sort them. (You will be sorting around $10^{12}/2^{24}\approx 60000$ numbers each time, not the entire list of a trillion numbers.) All this is routine, requiring no extended precision work. This is the parallelizable part: you could assign blocks of this vector $k$ to different processors.

Upon adding $(i-1)N/n$ to the values drawn for bin $i$, we will have a sorted array of $k_1$ integers between $1 + (i-1)N/n$ and $iN/N$, inclusive. Stringing out all $k_1+k_2+\cdots+k_n=10^{12}$ such integers gives the sorted list we desire. This addition can be done with bit arithmetic: just concatenate the $24$ bits for the binary representation of $i$ to the $40$ bits for the binary representation of the sample values.

Finally, there are various ways to obtain the multinomial vector $k$. The value of $S=10^{12}$ is so large it tends to overflow standard Multinomial generators. One obvious way to get around this is to generate $n$ independent Poisson$(S/n - \epsilon)$ variates where $\epsilon$ is a smallish number selected to assure, with high probability, that the sum of those variates will not exceed $S$. (A value of about $3\sqrt{S}$ will work around $99.8\%$ of the time. A Google search turns up an abstract that sounds very much like this approach: see Brown and Bromberg, An Efficient Two-Stage Procedure for Generating Random Variates from the Multinomial Distribution. TAS Aug 1984 vol. 38, no.3.) Then recursively draw from a Multinomial distribution designed to sum to the much smaller value of $S - \sum_{i=1}^n k_i$. (Once the target sum is sufficiently small, you might just as well sample with replacement from $1, 2, \ldots, n$ and add those counts to $k$. The R code at the end of this post implements this approach. It works.)

To summarize, the algorithm is

1. Draw the $n$ counts $k$ from a multinomial distribution. They sum to $S$.
2. For each $k_i$, sample the integers $1, 2, \ldots, N/n$ without replacement and sort them, giving $x_{i1} \lt x_{i2} \lt \cdots \lt x_{ik_i}$.
3. Form a 64-bit integer by concatenating each $(i,x_{ij})$ pair. Output these by looping over the $k_i$.

Asymptotically--fwiw--this algorithm requires $O(n + S\log(S/n))$ computational effort. Since $n\ll S$ and $\log(S/n) \approx 1$, this is practically $O(S)$ with a decently small multiplier--and a great deal faster than any $O(N)$ algorithm might hope to be.

Using R, I can typically accomplish $(1)$ in a couple of seconds; $(2)$ will take a couple of days; and $(3)$ depends heavily on the system's I/O capabilities, but is about as fast as one can imagine. With a few thousand processors deployed in the cloud, potentially the entire sorted list can be produced in a minute.

The following R code embodies this algorithm for the purposes of testing it. For testing, it is configured to process only the first $B=2^8$ entries of $k$ rather than all $2^{24}$ of them. (It is capable of doing the entire operation as-is if you have eight terabytes of RAM and a few days to wait :-). If you write the output for $k_i$ to disk at the time you generate it, instead of keeping it in RAM, you can get away with just a few hundred MB of RAM for the whole operation--but you still need a large disk.) The code needed to parallelize it depends on what system you use, but ought to be straightforward.

mult <- 2^40  # Bin size
n <- 2^24     # Eventually, 2^64/mult
N <- n * mult # Eventually, 2^64
size <- 10^12 # Amount of output. Eventually, 10^12
#
# Generate multinomial counts summing to `size`.  Store them in `k`.
#
system.time({
  s <- size
  k <- rep(0, n)
  while (s != 0) {
    if (s > 1e5) {
      s0 <- ceiling(s - 3*sqrt(s))
      cat("Generating", s0, "Poisson values...\n")
      dk <- rpois(n, s0/n)
      sdk <- sum(as.numeric(dk))
      if (sdk <= s) {
        k <- k + dk
        s <- s - sdk
      }
    } else {
      cat("Sequentially generating", s, "Poisson values.\n")
      u <- table(sample.int(n, s))
      q <- as.numeric(names(u))
      k[q] <- k[q] + u
      s <- 0
    }
  }
}) # 2-3 seconds typically.
#
# Generate random values within each bin according to its count.
# This can be deployed in parallel across P processors, assigning a block
# of n/P to each processor.  Total computing time is a few days, but with
# P large, it can be brought down to a minute or so.
#
# Output `i` is a 2 X `size` array of integers; each column represents a
# long integer.  These will need to be stored somewhere.
#
B <- 2^8 # Block length
system.time({
  j <- which(k > 0)[1:B]
  i <- sapply(j, function(l) rbind(l-1, sort(sample.int(mult, k[l]))))
  i <- matrix(unlist(i), nrow=2)
})
#------------------------------------------------------------------------------#
#
# Convert `i` to integers and plot.
#
x <- i[1,] * mult + i[2, ]
par(mfrow=c(1,3))
hist(k)
hist(x, breaks=seq(0, mult*B, length.out=64))
hist(diff(x), freq=FALSE, breaks=64)
curve(dexp(x, rate=size/N), add=TRUE, col="Red", lwd=2)
par(mfrow=c(1,1))

At the left is the histogram of $n=2^{24}$ multinomial counts. The middle is a histogram of the sequence of integers generated from the first $B=2^8$ entries in $k$. At right is a histogram of their differences, showing (a) they are positive (that is, the output truly is sorted and unique) and (b) they follow the expected Exponential distribution (superimposed in red).

Answer 2

You want $10^{12}$ numbers out of a range of $2^{64}$ possible numbers. Every number in the range has a $\frac{10^{12}}{2^{64}}$ chance of being selected. I think you can look at it like a coin toss but not a 50-50 chance. You can loop over the full $2^{64}$ range, at each iteration get a random binomial with the correct probability of success, $\frac{10^{12}}{2^{64}}$ in your case. If 1, then return that number of the iteration.

This can be implemented via a generator so you do not need to store all of the random numbers. You loop over all possibilities - the random numbers are in sequence. You call a random binomial with the same probability of success for all possible numbers - it is fair.

Given probability, you may not get exactly $10^{12}$ numbers. But you should get close when run over $2^{64}$ trials.

If you need exactly $10^{12}$ numbers, I can think of storing the first set of random numbers then count how many there are. If too many, call a random number generator based on the index of each number to delete that number. If too few, run the original loop again with a different seed, and insert into the original list until you get the right count - though this places a bias into the lower numbers in the possible range.