My goal is to do some large iteration Monte Carlo simulation of different card games that require on the order of 1 billion random decks to get a reliably accurate answer. The bulk of the computer time is actually shuffling/preparing the decks, not searching for the "hit" criteria I am looking for.

I am using an string function rich interpreted language. I have bit functions available that work up to 32 bit at a time. I am only interested in the ranks of the 52 card deck, not the suits. So I am looking for a fast algorithm that will give me 4 of each of the 13 ranks for a full 52 card deck. It doesn't matter which order the suits of a particular rank come in since we are ignoring the suits (other than we are enforcing that there only be 4 of each suit for each rank).

My first attempt was to just generate a normal deck with the ranks and suits "encoded" using cards numbered 1 thru 52 with an array of flags telling which cards were already seen. Each random number generated, (from 1 to 52), was for determining the next card (if there was no collision). This algorithm worked but was kinda slow, only producing about 6000 full decks per second on a 3.06 Ghz dual core computer (using only about 50% of the CPU).

My 2nd attempt was just to tweak the first attempt by only picking the first 48 cards (since the "seen" ratio is very high for the last 4 cards and likely many collisions will occur, slowing it down), then choosing a random number from 1 to 24 to handle the permutations of the last 4 "missing" cards which are known from simply scanning the flag array and finding the 4 that are marked "false" meaning they haven't been seen yet. This was roughly 50% faster as I was able to generate about 9000 full decks per second vs. about 6000 previously.

For my 3rd attempt, I am using a 32 bit random number to place 4 or 8 cards at a time from the unshuffled deck into a built up shuffled deck.

Any help/ideas would be greatly appreciated.

  • $\begingroup$ To be super clearer, is it correct that what you want is a random sequence of 52 numbers consisting of 4 1's, 4 2's, 4 3's, ..., and 4 13's? $\endgroup$ – John L. Oct 7 '18 at 18:52
  • $\begingroup$ Are you familiar with the Fisher–Yates shuffle? $\endgroup$ – Yuval Filmus Oct 7 '18 at 18:56
  • $\begingroup$ I will try Fisher-Yates and time it. Yes I am famiar with it I just didn't try it yet. Yes I only want four 1s, four 2s... four 13s (Kings). I heard some people say the Fisher-Yates shuffle algorithm is slow but I will try it and compare. $\endgroup$ – David James Oct 7 '18 at 19:03
  • 3
    $\begingroup$ It seems like the easiest suggestion is to switch from whatever programming language you are using to a faster one like C. $\endgroup$ – Yuval Filmus Oct 8 '18 at 0:47
  • 1
    $\begingroup$ FWIW, I was able to get around 1 million shuffles per second on a straightfoward C implementation of Fisher-Yates without any optimizations (other than "-O3") on an Intel Core i5-3740. $\endgroup$ – mhum Oct 17 '18 at 21:31

The questioner, David James wrote on the occasion of bounty.

Is there a way to speed this shuffle up being that there are repeated ranks? For example, suppose we only had 2 ranks but still 52 cards, could that be sped up over 52 distinct cards? if so, then how can we speed up 13 different ranks with 4 of each rank as in this problem? It seems like the rank only deck of 52 has many fewer possible distinct decks so there should be a way to shuffle them quicker. For example, suppose the deck starts out in rank order such as AAAA22223333... It won't matter if our random # generator (RNG) picks any of the first 4 cards, the result will be the same (a rank A card). I am hoping to use this to our advantage rather than just using Fisher-Yates shuffle as if we are dealing with 52 distinct cards.

Driven by David James's insight and insistence that there should be some way to speed this shuffle up being that there are repeated ranks, this new answer will be focusing whether it is possible to and how to take advantage of the repeated ranks. Other methods to improve the performance, which might be much more effective, powerful and randomness-correct are not considered.

Yes, there is a way to take advantage of the repeated ranks. Please go to my Java playground provided by repl.it and hit the button "run". You will see output like the following.

[6, 10, 3, 8, 7, 10, 3, 3, 5, 13, 2, 2, 4, 5, 4, 4, 1, 9, 11, 11, 6, 9, 10, 10, 5, 12, 13, 3, 1, 7, 2, 4, 8, 2, 3, 13, 8, 12, 6, 9, 4, 13, 9, 12, 5, 7, 1, 7, 11, 1, 3, 6]
    warm-up one: 532197/second
first 4 skipped: 552181/second
    warm-up two: 528262/second
    David-Apass: 656598/second
          Naive: 528541/second
    Improvement: 24.2

Your actual output might vary greatly because of many factors. What is important is the last three lines of output. "David-Apass: 656598/second" means 656598 decks per a second by David-Apass shuffle, a variation of Fisher-Yates shuffle that takes advantages of the equal ranks. "Naive: 528541/second" means 528541 decks per a second by a common implementation of the Fisher–Yates shuffle. The line of comparison "Improvement: 24.2%" shows the speed-up of David-Apass shuffle over the common Fisher-Yates shuffle. The nature or the mechanism of David-Apass shuffle apart might have been discovered many years ago. However, I have not found any reference to it yet. All I can claim is that I have discovered it independently, thanks to David James' insight and insistence.

** The following is the accepted answer to the version of the question at https://cs.stackexchange.com/revisions/98273/5, which wass before the bounty.**

It is not a trivial task to create 1 billion random decks, although it is not a daunting task any more because of today's fast computer.

Firstly, let us talk about the question proper. What would be a fast algorithm?

It looks like that the standard algorithm to shuffle cards, Fisher–Yates shuffle is a pretty good choice in term of speed. It is certainly one of the easiest to implement as well. Here is the complete pseudocode.

for t from 1 to 1000000000 do
    let arr=[1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10,11,11,11,11,12,12,12,12,13,13,13,13]
    for i from 51 down to 1 do
         let j  be a random integer such that 0 ≤ j ≤ i
         exchange arr[j] and arr[i]
    output arr

How can we make this algorithm run faster?

The very first rule about performance improvement in real-life task is measure, measure and measure. A very rough measurement in my Python code reveals that the most of the time is spent in the generation of random integers. This is not surprising once you take a close look at the algorithm above, assuming outputting arr is fast enough. In fact, that observation is probably true in any reasonably good algorithm to create a random deck of cards. So the next attack should be how to produce the random number faster or how to reduce the numbers of random numbers need to be generated. If a user has recorded billions of random numbers of similar kind, then fetching them from the records might be a faster way. Or use a faster random number generator. On the other hand, I cannot see how to reduce the number of needed random numbers without affecting the randomness of the generated deck.

Note that manipulation of string is much slower than manipulation of numbers in almost all programming languages, it is preferred to operate the cards in the form of numbers. Since there are a billion of decks, the storage becomes an issue. One of the most efficient ways to encode a deck is to encode them as a 52-digit number in base 13. However, it is not be very convenient to encode to decode from that big number. Another way is, as the questioner mentioned to me in chatroom, to encode each card in a nibble, i.e., four bits, similar to BCD, which is compact enough while much easier to manipulate. Also note that compression cannot probably do much to reduce the storage of this bunch of random numbers.

As Yuval mentioned, an obvious suggestion is to take the advantage of a fast programming language such as C. Let us do not consider assembly languages yet.

Let me conclude this answer with the strategy taken by the questioner, showing the art of applying compute science to real life situations may lead to very different directions.

I decided to just store the 1 billion shuffled decks in a disk file (about 26GB for 1 billion decks), then I can just reread it when I want to do other analysis. Generating it once will isolate the speed reduction and then the analysis part should "rip" (relatively speaking).

  • $\begingroup$ One slight tweak to the Fisher-Yates shuffle algorithm is to not shuffle in place. Have a 2nd array to hold the shuffled cards so then you don't need to do a 3 step swap, you can just do it in 2 assignments. $\endgroup$ – David James Oct 8 '18 at 5:12
  • $\begingroup$ Your version of Fisher–Yates shuffle is interesting. Although I doubt it will be significant in language like C, where swapping is fast because of locality while allocation might be slow, it seems helping significantly in your situation. $\endgroup$ – John L. Oct 8 '18 at 6:29
  • $\begingroup$ Also what is very interesting about this particular task is that the # of possible distinct decks is quite large at $52!/(4!)^{13}$ which is about $9.2$ x $10^{49}$. 1 billion decks in only magnitude 9 which means even if we sort the decks, the difference between sorted decks will still require most of the original 26 bytes to encode, although there would be some storage savings. It seems in this case, the simplest way is to just store 2 cards per byte, each in base 16 (such as 13 and an 8 ranks would be 16*13+8) so I would store Chr(216) for those. 26 billion bytes is actually about 24.2GB. $\endgroup$ – David James Oct 8 '18 at 7:01
  • $\begingroup$ 24.2 GB is a reasonable size data file by todays standards since most hard disks can easily handle that. If I ever need more than 1 billion decks, then I would consider packing them tighter, sorting the decks, and using some compression. However I would still expect something like 10 billion decks to take up about 200GB of space so they would still be about 20 bytes per deck. $\endgroup$ – David James Oct 8 '18 at 11:50
  • $\begingroup$ Unfortunately, what algorithm is fast in one language might not be the fastest in some other language depending on many factors. I think it is important to try different algorithms in the same programming environment to get a feel for which is the fastest and to gauge the relative speed of them, but also to try other languages to see if there is a performance enhancement (maybe better code optimization in some language for example). When I save the 1 billion random decks to a disk file and then reread them and process them, I will report back the speed enhancement here. $\endgroup$ – David James Oct 9 '18 at 16:31

I'm pretty sure there's no way to take any significant advantage of the fact that you don't care about the suits. The most effective solution is likely to just generate random permutations of the integers from 0 to 51, and map these to card ranks by dividing by 4 and rounding down (which can be done with a simple bit shift, since 4 is a power of 2).

As Apass.Jack notes in their answer, for generating the permutations you almost certainly want to use some form of the Fisher–Yates–Durstenfeld–Knuth shuffle, since it avoids the inefficiency of the sample-and-reject-duplicates method when only a few valid choices remain. Since you're generating permutations of consecutive integers, you may want to consider using the "inside-out" variant of the shuffle (probably first described by Knuth in 1994):

for i from 0 to n − 1 do
    j ← random integer such that 0 ≤ j ≤ i
    if j ≠ i
        a[i] ← a[j]
    a[j] ← i

The main advantage of this algorithm over the "standard" FYDK shuffle is that you don't need to initialize the results array a before shuffling it, since it gets dynamically initialized during the shuffle.

Nonetheless, you may want to implement both a standard FYDK shuffle and the inside-out variant and benchmark them, since their different memory access patterns can affect real-world performance in ways that are hard to predict in advance.

In any case, at this point the biggest remaining opportunities for optimization are likely to be:

  1. using random numbers more efficiently,
  2. using a faster random number generator, and
  3. switching to a lower-level language.

Of these three, I'd actually recommend starting with #3; this is exactly the kind of problem where using a low-level compiled language like C or C++ (or even a mid-level language with a good JIT compiler, like Java) can provide significant performance benefits with relatively little cost in development effort.

Once you've ported your code to a low-level language, you may want to consider optimizing your choice of RNG. There are some pretty fast PRNGs out there (my personal favorite being Bob Jenkins' "flea" RNG, as reviewed e.g. here), but to properly realize their benefits, you need to also minimize the overhead cost of calling the generator (ideally, letting your compiler inline the PRNG code directly into your shuffling loop).

Finally, it's worth noting that most PRNGs generate something like 32 or more bits of (pseudo)randomness per call, whereas a Fisher-Yates shuffle of a 52 card deck only uses less than 6 bits per iteration. Thus, there could certainly be room for optimization by splitting RNG outputs into smaller pieces and reusing them for multiple iterations of the shuffling loop.

That said, I'd personally rather avoid this option, at least until I had exhausted all other possible avenues for optimization. One problem is that, even if your RNG provides perfectly random outputs, it's easy to introduce subtle bugs in the splitting process that can bias your results in ways that may not be so easy to detect. Another problem is that trying to squeeze every last drop of entropy out of your (P)RNG output is a good way to magnify any non-random biases that output may have. Personally, I would feel safer using a fast but possibly not quite perfect PRNG to generate lots of extra random bits, and converting them to random numbers in the required range in a simple and relatively foolproof way, than trying to use up every last bit of my RNG output and relying on that RNG and my conversion code having no subtle weaknesses.

Still, there are some simple optimizations you can safely make here. For example, you can trivially save one RNG call (and some other unnecessary work) per shuffle by modifying the code above to look like this:

a[0] ← 0
for i from 1 to n − 1 do
    j ← random integer such that 0 ≤ j ≤ i
    if j ≠ i
        a[i] ← a[j]
    a[j] ← i

The only difference between this code and the version I copied from Wikipedia above is that I skipped the first iteration, where we always had i = j = 0, and instead just initialized a[0] before the loop. While such little optimizations won't have much effect alone, they're pretty easy and safe and can add up.

In fact, since you said don't care about the card suits, I guess we can save a few more RNG calls:

for i from 0 to 3 do
    a[i] ← 0
for i from 4 to n − 1 do
    j ← random integer such that 0 ≤ j ≤ i
    if j ≠ i
        a[i] ← a[j]
    a[j] ← ⌊i / 4⌋

This generates the card ranks directly, by doing the division by 4 inside the initialize-and-shuffle loop, and splits the first four iterations of that loop (where all elements of the array are still identical) into a separate loop. This way, you'll only need 48 RNG calls to shuffle 52 cards. That's almost an 8% speedup (assuming that random number generation dominates the time cost), which might be significant just by itself.

Ps. FWIW, while precalculating the shuffled decks and storing them on disk is certainly an option worth considering, I very much doubt it is optimal. With a fast shuffle and a fast RNG, both implemented in a low-level language, shuffling the decks dynamically should be orders of magnitude faster than reading them from even a very fast solid-state disk.

Addendum: There's a potential advantage to reshuffling your previously shuffled decks, rather than generating a new deck for each shuffle: it can smooth out possible biases in your shuffling process and give you a more uniform distribution of shuffles.

Basically, if you start with a new in-order deck for each shuffle (or generate one during the shuffle, like the code above does), then you're repeatedly sampling decks from the same distribution. If your RNG (or the way you convert its output into array indices) is not perfect, then this distribution may be biased.

In particular, if you're using a PRNG with less than $\log_2(52!) \approx 225.6$ bits of internal state, then you're guaranteed to have some decks that you can never generate by shuffling an in-order deck, simply because there will be more possible decks than there are possible initial states for your PRNG at the beginning of the shuffle. If you don't care about the suits, the corresponding threshold drops to $\log_2(52! \mathbin/ 4!^{13}) \approx 166$ bits. But of course, in either case, you'd need at least twice as many bits just to be reasonably sure that all possible decks can be generated, and even more to guarantee a reasonable level of uniformity.

On the other hand, if you keep reshuffling a previously shuffled deck, then you're effectively performing a (pseudo)random walk on the set of possible permutations, with the possible transitions defined by the relation "$A$ can be shuffled into $B$ by the PRNG". As long as this random walk is ergodic (and it should usually be, even for pretty crappy RNGs), it should converge towards a uniform distribution on the set of all possible permutations (since the transition kernel is also invariant with respect to the semigroup action).

Thus, if you could keep shuffling the same deck over and over, even with a poor-quality PRNG, you should in principle have a pretty high chance of eventually encountering each possible shuffle equally often. Of course, the order in which you'd encounter these shuffles might not be fully random, and (if you never reseeded the PRNG from a truly random source) there would still be a small probability of getting caught in a short cycle that visits only some of the possible shuffles. And of course, with a 52 card deck (with or without suits), you'd never actually be able to test this within your lifetime on any computer you could possibly build.

But still, in general, a good rule of thumb is that reshuffling the same deck over and over will tend to smooth out any possible biases in the shuffling process, compared to starting from an unshuffled deck each time. In fact, that's pretty much the reason why the often rather imperfect shuffling done by humans using a physical deck of cards is still sufficient to avoid generating noticeably biased decks, and why card players are advised to be extra careful in shuffling a new deck before using it for the first time.

Fortunately, both the standard Fisher-Yates shuffle and the inside-out version described above can be easily made to reshuffle an existing deck. In fact, the standard version does that by default, while for the inside-out variant the changes are minor. Specifically, assuming that a is already a (possibly shuffled) permutation of n cards, either of the following loops will serve to reshuffle it:

// standard FYDK shuffle (descending order)
for i from n - 1 down to 1 do
    j ← random integer such that 0 ≤ j ≤ i
    if j ≠ i: swap a[i] and a[j]

// "inside-out" version (Knuth 1994)
for i from 1 up to n − 1 do
    j ← random integer such that 0 ≤ j ≤ i
    if j ≠ i: swap a[i] and a[j]

Either of these will require $n-1$ RNG calls to (re)shuffle an $n$ card deck; unfortunately, in this case, I really don't see any practical way to save any more RNG calls even if some of the cards in the deck are known to be equivalent.

Addendum 2: If your language / CPU can do fast 256-bit math (specifically, division and remainder), and you have access to a high-quality RNG than can generate 256-bit numbers quickly (such as using two calls to AES_CTR_DRBG, which should be pretty fast on x86 CPUs that support the AES-NI instruction set, if your crypto library makes use of it), then you could do something like this:

r ← random 256-bit integer (i.e. 0 ≤ r < 2^256)

// optional bias elimination via rejection sampling
while r < 2^256 mod 52!
    r ← random 256-bit integer (i.e. 0 ≤ r < 2^256)

// inside-out FYDK shuffle-and-initialize
a[0] ← 0
for i from 1 to 51 do
    j ← r mod (i + 1)
    r ← ⌊r / (i + 1)⌋
    if j ≠ i: a[i] ← a[j]
    a[j] ← i

where hopefully either your compiler is smart enough to fuse the modular reduction and division of r by i + 1 into a single operation, or your language directly provides a combined "divmod" operator or function for it.

As noted in the comments above, the bias elimination by rejection sampling is not strictly required; omitting it just means that about 30% of all possible shuffles are very slightly more likely (about $1 + \frac{52!}{2^{256}} \approx 1.0000000007$ times as likely, to be exact) to occur than the rest. It's up to you to decide whether that bias is something you're willing to accept. If you do decide to use the bias elimination, do precalculate the constant 2^256 mod 52!, as calculating it from scratch is likely to be almost as expensive as the shuffle itself.

If you don't have access to efficient 256-bit math, you can always split the shuffling loop into smaller pieces. For example, $20!$, $32!\mathbin/20!$, $43!\mathbin/32!$ and $52!\mathbin/43!$ are all less than $2^{64}$, so using 64-bit math you could do something like this:

s ← 1
a[0] ← 0
for t in (20, 32, 43, 52) do
    m ← 2^64 mod (t! / s!)
    r ← random integer such that m ≤ r < 2^64
    for i from s to t - 1 do
        j ← r mod (i + 1)
        r ← ⌊r / (i + 1)⌋
        if j ≠ i: a[i] ← a[j]
        a[j] ← i
    s ← t

where I'm assuming that the code used to generate a random integer such that m ≤ r < 2^64 internally consumes 64-bit random numbers and does any rejection sampling needed to ensure that the fall into the desired range. Of course, in practice you'd probably want to unroll the outer loop and to precalculate the m values for each iteration, or at least to store them in a look-up table.

If you're limited to 32-bit math (or, in particular, if you only have a 64 → 32 + 32 bit divmod operation), you can instead iterate over t in (12, 19, 25, 30, 35, 40, 44, 48, 52) and set m to 2^32 mod (t! / s!). Note that this requires nine iterations of the outer loop instead of eight, since there will be some unused entropy left in r at the end of each iteration. In principle we could recycle this entropy, but that would complicate the code even more, so I'll leave that as an exercise. ;-)

In any case, whichever method you choose, you will definitely want to test your code very throughly before using it. I have not done so myself for the pseudocode examples above, so they might contain bugs, and porting them from pseudocode to your actual programming language could easily introduce more.

Unfortunately, shuffling is exactly the kind of task where it's very easy to make subtle mistakes that are hard to catch in testing due to the inherent stochasticity of the task. Some useful tests I'd recommend include:

  • shuffling a small deck of, say, 3 to 5 cards many times, and checking that each possible shuffle occurs approximately equally often (e.g. using Pearson's $\chi^2$ test);
  • shuffling a full 52 card deck many times (starting with an ordered deck each time, to maximize any observable bias), recording the position of some specific card in the deck, and checking that those positions are approximately uniformly distributed (again using e.g. a $\chi^2$ test; repeat this for each of the 52 cards);
  • same as above, but instead recording the distribution of cards occurring in a given position in the shuffled deck;
  • same as the two tests above, but instead record the difference in the positions of two specific cards (modulo 52) or the difference in the card values at two specific positions in the shuffled deck (note that in both cases there will be 51 possible differences, since 0 is not a possible result);
  • same as above, but instead of just recording the difference, record each observed pair of positions or card values (note that you'll need more shuffles per test for this, since there will now be 52 × 51 possible outcomes);
  • count the number of disjoint cycles in the permutation generated by the shuffle, and check that the distribution of this count over many shuffles matches the expected distribution given by the Stirling numbers of the first kind (in particular, this test should catch any off-by-one errors leading to an accidental implementation of Sattolo's algorithm instead of the standard FYDK shuffle);
  • other statistical randomness tests, especially any tests that would directly measure some specific property that you both care about in your actual simulation and can predict the theoretically expected distribution of.
  • $\begingroup$ Thanks for your detailed answer. I am somewhat surprised that no optimization can be made if we ignore the suits (other than we know there must be 4 of each rank). Yes I agree in a low level language like compiled c, generating the random decks is likely faster than storing them on disk and rereading them, however, in some slow interpreted languages, it is actually faster to read them in from a disk file. Also, one way around the "sample and reject duplicates method", as you refer to it, is to stop using that method when only a few choices remain and use a 2nd random number for what's left. $\endgroup$ – David James Oct 20 '18 at 8:43
  • $\begingroup$ Continuation of above comment... For example, in my code, I would select 48 random cards using the random number generator from 1 to 52 (inclusive) and an array of flags keeping track which cards were already chosen, but then for the last 4 cards, I would pick a random number between 1 and 24 (inclusive) and then have code to place the last 4 cards in all possible 24 permutations of them. I suppose I could take it 1 step father but that would then require a large if or case statement cuz the last 5 cards would have 120 possible permutations which is a lot of code just for the last 5 cards. $\endgroup$ – David James Oct 20 '18 at 8:46
  • $\begingroup$ @DavidJames: That's certainly possible; in fact, on modern CPUs, I suspect it would be even faster to just use plain old arithmetic (rather than a lookup table or a big case statement) to decode your mod-24 (or mod-120) random number into the last four or five indices for the shuffle. But my general intuition is still that it's best to keep your RNG fast and your shuffling loop simple. $\endgroup$ – Ilmari Karonen Oct 20 '18 at 11:53
  • $\begingroup$ ... That said, your comment did inspire me to come up with a simple optimization that lets you shuffle a previously unshuffled 52 card deck with only 48 RNG calls, assuming that you don't care about the suits, which I've edited into my answer above. But I also added (a rather lengthy, as it turns out) addendum about why always starting from a new, unshuffled deck may not be the best idea in practice, since it's more likely to expose any RNG biases than reshuffling the same deck over and over. $\endgroup$ – Ilmari Karonen Oct 20 '18 at 11:53

Just about the random number generation: You will need 52 billion random numbers, so the period of a pseudo-random generator must far exceed 52 billion.

I assume you don't need cryptograhically secure random numbers, so you can use a linear congruential random number generator, or some variation of Marsaglia's xorshift generator.

IBM implemented a linear congruential RNG taking one processor cycle per random number on a POWER processor; it should be portable to any processor with a FPU with fused multiply-add.

A 64-bit xorshift generator can probably be vectorised with 256 bit vectors to create four random numbers at a time, and the code to extract a number from 0 to 51 / 0 to 50 / 0 to 49 etc can be vectorised as well. Replicate this four times, and you will get 16 shuffled sequences at a time very, very fast indeed. Ten million sequences per second should be not very hard.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.