I want a fast rank only shuffling algorithm for a full 52 card deck
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:


using random numbers more efficiently,
using a faster random number generator, and
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 &lt; 2^256)

// optional bias elimination via rejection sampling
while r &lt; 2^256 mod 52!
    r ← random 256-bit integer (i.e. 0 ≤ r &lt; 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 &lt; 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 &lt; 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 &rarr; 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 &times; 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.