# Shuffle product of two lists

I want to do the following: I have two lists {a_1,...a_n}, {b_1,...,b_n} and I would like to build now all shuffles out of this. This means all unions of these lists while still keeping the individual order of the parent lists.

Example

{a,b}, {c,d}


{{c,d,a,b},{c,a,d,b},{c,a,b,d},{a,c,d,b},{a,c,b,d},{a,b,c,d}}

Is there any easy way to do this?

• A friendly helper, I'm glad you like my answer, but I always recommend waiting 24 hours before Accepting an answer, to let everyone around the world have a chance to reply. You may like other answers better if you give them a chance to happen. :-) Feb 3, 2014 at 10:37
• Ok, I'll wait :) Feb 3, 2014 at 10:43

## Recursion

This is top-level code and therefore unlikely to be as efficient as a compiled solution, but the recursive algorithm should have reasonable computational complexity:

f[u : {a_, x___}, v : {b_, y___}, c___] := f[{x}, v, c, a] ~Join~ f[u, {y}, c, b]

f[{x___}, {y___}, c___] := {{c, x, y}}  (* rule for empty-set termination *)


Now:

f[{a, b}, {c, d}]

{{a, b, c, d}, {a, c, b, d}, {a, c, d, b}, {c, a, b, d}, {c, a, d, b}, {c, d, a, b}}


### Proof of result

A question was raised regarding the validity of the result of this function. You can graphically demonstrate that this function yields the correct result (as I understand it) for length 2,4 lists with this:

ArrayPlot @ f[{Pink, Red}, {1, 2, 3, 4}]


• Pink always precedes Red
• The gray values are always in order
• No shuffles are missing
• No shuffles are duplicated

I must conclude that the referenced paper is in error, or the original question diverges from the definition of "shuffle product" therein.

## Extension to multiple lists

While it is possible to extend a two-list shuffle product to multiple lists using Fold it may be of interest to do it directly. This code is slow however and should not be used where performance matters.

f2[in_, out___] :=
Join @@ ReplaceList[in, {x___, {a_, b___}, y___} :> f2[{x, {b}, y}, out, a]]

f2[{{} ..}, out__] := {{out}}

f2[{{1, 2}, {Cyan}, {LightRed, Pink, Red}}] // Transpose // ArrayPlot


## Alternative method

Edit: I had a function g which used Permutations on a binary list to generate all the base orderings, but I did not fill these orderings efficiently. rasher used the same start but came up with a clever and fast way to fill those orderings. I am replacing this section of my answer with a refactored version of his code; credit to him for making the Permutations approach competitive. (If you are interested in my original, slow g see the edit history.)

Here is a refactoring of rasher's shufflem function; I shall call mine shuffleW.

Edit 2017: Now significantly faster after moving or eliminating several Transpose operations.

shuffleW[s1_, s2_] :=
Module[{p, tp, ord},
p = Permutations @ Join[1 & /@ s1, 0 & /@ s2]\[Transpose];
tp = BitXor[p, 1];
ord = Accumulate[p] p + (Accumulate[tp] + Length[s1]) tp;
Outer[Part, {Join[s1, s2]}, ord, 1][[1]]\[Transpose]
]

shuffleW[{a, b}, {c, d}]

{{a, b, c, d}, {a, c, b, d}, {a, c, d, b}, {c, a, b, d}, {c, a, d, b}, {c, d, a, b}}


## Timings

Timings of these functions as well as Simon's shuffles and rasher's shufflem, performed in v10.1.

I am also including rasher's unnamed but very fast penultimate Ordering method. I shall call it shuffleO and format it for readability.

shuffleO[s1_, s2_] :=
With[{all = Join[s1, s2]},
Join[1 & /@ s1, 0 & /@ s2]
// Permutations
// Map[Ordering @ Ordering @ # &]
// Flatten
// all[[#]] &
// Partition[#, Length[all]] &
]

s1 = CharacterRange["a", "k"];
s2 = Range @ 11;

f[s1, s2]        // Length // RepeatedTiming
shuffles[s1, s2] // Length // RepeatedTiming
shufflem[s1, s2] // Length // RepeatedTiming
shuffleO[s1, s2] // Length // RepeatedTiming
shuffleW[s1, s2] // Length // RepeatedTiming

{2.6523, 705432}

{2.0618, 705432}

{1.455, 705432}

{0.78, 705432}

{0.759, 705432}


rasher's methods are the fastest in this test. I find my f the most readable as the mechanism of its action is directly visible. Take your pick. :-)

### Disparate list lengths

Yi Wang posted a nice clean method that has some interesting properties. Specifically it is the best performing solution so far in the case of input lists of significantly disparate length, but they must be fed in the correct order or the performance is magnitudes worse. First in the symmetric test above:

Fold[insertElem, {{s2, 0}}, s1][[All, 1]] // Length // RepeatedTiming

{2.67, 705432}


Now with disparate lists:

s1 = Range[80];
s2 = {"a", "b", "c"};

Fold[insertElem, {{s2, 0}}, s1][[All, 1]] // Length // RepeatedTiming

{6.25, 91881}


Swap the lists and try again (so that insertElem is Folded over the shorter list):

{s2, s1} = {s1, s2};

Fold[insertElem, {{s2, 0}}, s1][[All, 1]] // Length // RepeatedTiming

{0.121, 91881}


This is significantly faster than all the others (order makes little difference to these):

f[s1, s2]        // Length // RepeatedTiming
shuffles[s1, s2] // Length // RepeatedTiming
shufflem[s1, s2] // Length // RepeatedTiming
shuffleO[s1, s2] // Length // RepeatedTiming
shuffleW[s1, s2] // Length // RepeatedTiming

{0.529, 91881}

{0.5637, 91881}

{0.562, 91881}

{0.4135, 91881}

{0.4356, 91881}

• Very nice, that's what I had in mind. Thanks! Feb 3, 2014 at 10:33
• @Mr.Wizard: Curious what my real entry does in your environment...
– ciao
Feb 4, 2014 at 9:00
• @rasher That looks very interesting. I'm working on another post but I'll run timings a bit later. (If I forget, poke me.) I see yours is also based on Permutations of a binary list; I was not happy with the way I filled in those values. I hope you worked some magic! :-) Feb 4, 2014 at 9:04
• @Mr.Wizard: Yes, I scratched my beard for a bit on that! I'm sure there's still something left in it, was pleased with initial results - and of course, it extends to more than two lists fairly easily, I think. Look forward to your results!
– ciao
Feb 4, 2014 at 9:06
• I picked this one, since - as Mr. Wizard stated and I agree - the mechanism of its action is directly visible. I don't need performance at this moment so much :) Feb 4, 2014 at 9:44

Probably not too bad performance-wise, haven't tested though:

x = {a, b};
y = {c, d};

n = Length[x] + Length[y];
px = Subsets[Range[n], {Length[x]}];
py = Reverse[Subsets[Range[n], {Length[y]}]];

Normal /@ MapThread[SparseArray[{#1 -> x, #2 -> y}] &, {px, py}]

(* {{a, b, c, d}, {a, c, b, d}, {a, c, d, b}, {c, a, b, d}, {c, a, d, b}, {c, d, a, b}} *)


Update

This has better performance:

shuffles[x_, y_] := Module[{n, px, py, z, xy},
n = Length /@ {x, y};
{px, py} = Subsets[Range @ Tr @ n, {#}] & /@ n;
py = Reverse @ py;
xy = z = Join[x, y];
(z[[#]] = xy; z) & /@ Join[px, py, 2]
]

shuffles[{a, b}, {c, d}]
(* {{a, b, c, d}, {a, c, b, d}, {a, c, d, b}, {c, a, b, d}, {c, a, d, b}, {c, d, a, b}} *)


Tweaking other people's code

This is a refactoring of Mr Wizard's refactoring of rasher's code. The primary changes are to Flatten the ordering matrix so that the whole output is extracted with a single Part call (and subsequently partitioned), and to reduce the number of Transpose operations. On my PC this gives about a 30% speed gain over Mr Wizard's version:

shuffleSWR[s1_, s2_] := Module[{p, tp, accf, ord, ss},
p = Transpose @ Permutations @ Join[1 & /@ s1, 0 & /@ s2];
tp = BitXor[p, 1];
ord = Accumulate[p] p + (Accumulate[tp] + Length[s1]) tp;
ss = Join[s1, s2];
ss[[Flatten @ Transpose @ ord]] ~Partition~ Length[ss]]

• I refactored shuffles. If you don't like the changes just revert it. Feb 3, 2014 at 11:56

Something different, impractical yet amusing, with a hint of mathematical insouciance...

s1 = {a, b, c}
s2 = {f, g, h}

sall = Join[s1, s2];
ln = Length[sall];

s = Solve[
Max[sall] == ln && Min[sall] == 1 && Less[Sequence @@ s1] &&
Less[Sequence @@ s2] && Unequal[Sequence @@ sall], sall, Integers];

Range@ln /. MapAt[Reverse, s, {All, All}]


And, for my actual entry into the horse race, I offer the following. Seems to beat the current record holder pretty handily in my limited tests.

shufflem[s1_, s2_] := Module[{ss, sr, s1l, s2l, p, pp, tp, pres},

{ss, sr} = {Join[s1, s2], Range@((s1l = Length@s1) + (s2l = Length@s2))};

p = Permutations[Join[ConstantArray[1, s1l], ConstantArray[0, s2l]]];

pp = (Accumulate@Transpose[p] // Transpose) p;

tp = Clip[p, {1, 0}, {1, 0}];

pres = pp +
Clip[(Accumulate@Transpose[tp] // Transpose) tp + s1l, {s1l + 1, Max[sr]}, {0, 0}];

ss[[#]] & /@ pres
]


As requested by Mr. Wizard, an explanation for readers.

Pondering this over dinner, I imagined it as one list "bubbling" through the other. The pattern that came to mind was obvious: a binary counter kind of thing, where the 1 and 0 represent positions in the shuffled list.

Having distinct values makes the permutation operation reasonable, i.e., just naively permuting the whole "deck" of the lists and the filtering results in mostly wasted work, since most will have elements violating the order rule, while having only representations of some element in some list means we'll never see such transpositions in the permutation.

With that in mind, I build a list of possible positions in the shuffled list (the range of total length in sr), along with making a joined version of the lists (ss).

I then build (in p) a list with 1 and 0, in quantities matching the lengths of each original list, and generate the permutations. So now I have a big list of all possible permutations of the final shuffled list, with 1's where an element from the first list goes, and a 0 for those from the second.

We want the elements to fill these 1/0 slots to be in the order of the originals, that is, element 1,2,3... for each.

So I accumulate over each permutation (I double transpose for performance), leaving me with say something like {1,0,1,1,0,0} going to {1,1,2,3,3,3}. I'm only interested in the items from the first list on the first pass, so I multiply the accumulants by their respective permutation, e.g., in the current example we'd get {1,0,2,3,0,0} for that one permutation, meaning from the first list, slots 1,3, and 4 get filled from the list elements 1,2 and 3 of the joined list.

The second pass fills the 0 spots from the second list, by doing the same thing with some hand-waving to get the values into the second half of the joined list: I clip the permutations to turn 1->0 and 0->1, do the same machination as the last, but I add an offset (the length of the first list prefix in the joined list). I clip out noise values, ending up with the list of permutations with elements looking like {1,4,2,3,5,6}, etc.

Finally, I map over that result of permutation lists indexing them into the joined list, so for example, if the joined list were {a,b,c,d,e,f}, the example permutation element would result in {a,d,b,c,e,f} being pulled.

Vectorizing the work makes it snappy - the vast majority of the time overall is in mapping over the work to get the final result.

Lastly, a surprisingly quick (I thought it would be much slower than it turns out to be) method that avoids any mathematical manipulation:

Module[{jl = Flatten[{#1, #2}]},
Partition[
jl[[Flatten[
Ordering /@
Ordering /@
Permutations[
Join[ConstantArray[1, Length[#1]],
ConstantArray[0, Length[#2]]]]]]], Length[jl]]] &[s1, s2]


And this use of ordering is even faster due to the relatively small size of input lists, turns out mapping over them is quicker, resulting in overall time that seems competitive with the refactored*refactored version of my original. Interesting:

Partition[ Join[#1, #2][[Flatten[ Ordering[Ordering[#]] & /@
Permutations[Flatten[{1 & /@ #1, 0 & /@ #2}]]]]],
Length[#1] + Length[#2]] &[s1, s2]


Finally, the above generalized to allow more lists (timings on a netbook that was near bursting into flames):

lists = {{a, b, c, d}, {e, f, g, h, i}, {j, k, l, m, n, o}}

res2 = Partition[ Flatten[##][[Flatten[ Ordering[Ordering[#]] & /@
Permutations[Flatten[MapIndexed[(idx = #2[[1]]; idx & /@ #1) &, ##]]]]]],
Length[Flatten[##]]] &[lists]; // Timing

Length[res2]

res2 // Short

(* {4.539629, Null} *)

(* 630630 *)

(* {{a,b,c,d,e,f,g,h,i,j,k,l,m,n,o},{a,b,c,d,e,f,g,h,j,i,k,l,m,n,o},<<630627>>,{j,k,l,m,n,o,e,f,g,h,i,a,b,c,d}} *)


Any of the methods become pretty useless with large lists and/or large number of lists:

e.g, two lists of 40 elements, were one able to generate permutations that adhere to the ordering rule at a rate of 1,000,0000 per second would far exceed the age of the universe before completing.

The following generates a random permutation that adheres to the ordering rules for such cases:

randomOrderedShuffle[lists__] :=
Join[lists][[Ordering@Ordering@RandomSample@(Join @@ ConstantArray @@@
Transpose[{Range@Length@{lists}, Length /@ {lists}}])]]

(* The following example has 3,644,153,415,887,633,116,359,073,848,179,365,185,734,400  "valid" permutations... *)

randomOrderedShuffle[Range[10], Range[11, 20], Range[21, 30], Range[31, 40], Range[41, 50], Range[51, 60]] // Timing

(*  {0., {1, 21, 41, 11, 22, 12, 31, 23, 51, 2, 3, 52, 42, 13, 14, 53, 43, 15, 32, 4, 54, 16, 33, 17, 34, 55, 5, 44, 45, 56, 6, 35, 57, 36, 7, 46, 8, 37, 58, 24, 59, 47, 25, 9, 10, 60, 26, 48, 18, 27, 19, 38, 49, 28, 39, 29, 40, 30, 20, 50}}  *)

• :) ... it does only work with unassigned variables, though. Feb 4, 2014 at 7:37
• Yes, it was meant as a goof. Adding my real entry now...
– ciao
Feb 4, 2014 at 8:58
• This is the approach I wanted to take, but I couldn't see a fast way to compute pres. The Accumulate trick is very clever. Excellent work! Feb 4, 2014 at 10:07
• @SimonWoods: Appreciate the comment! Yes, as I commented to Mr.W, I was scratching around for a bit on that. I think there's more tweaking ;-)
– ciao
Feb 4, 2014 at 10:36
• Regarding the last point, that the Map is the slowest part, you can save some time by doing ss[[Flatten@pres]] ~Partition~ Length[ss] instead Feb 4, 2014 at 11:38

What about inserting the first list into the second list at every allowed position (where allowed means preserving the order of the first list)?

lst1 = {a, b};
lst2 = {c, d};

insertElem[lst_, x_] :=
Join @@ (Table[{Insert[#, x, i], i}, {i, #2 + 1, Length@# + 1}] & @@@ lst)

Fold[insertElem, {{lst2, 0}}, lst1][[All, 1]]

{{a, b, c, d}, {a, c, b, d}, {a, c, d, b}, {c, a, b, d}, {c, a, d, b}, {c, d, a, b}}


I was told not to use the slow Insert. But here generating permutations is much harder. Thus Insert may not be too bad.

EDIT: Another worse attempt

This is not as efficient as the above one. But just for fun, one can generate the positions of lst2 in the result, then generate positions of lst1~Join~lst2 in the result. Finally, permute lst1~Join~lst2 into the final result:

comb[lenIn_, lenOut_] := Module[{i},
i[0] = 0;
Flatten[#, lenOut - 1] & @ With[{
elem = i@# & /@ Range@lenOut,
ranges = Sequence @@ ({i[#], i[# - 1], lenIn - 1} & /@ Range@lenOut)},
Table[elem, ranges]]]

(* e.g. comb[3,2] gives {{0, 0}, {0, 1}, {0, 2}, {1, 1}, {1, 2}, {2, 2}} *)
(* One can also write comb as follows, with even worse performance *)
(* comb[lenIn_, lenOut_] := Select[Tuples[Range[0, lenIn - 1], lenOut], OrderedQ] *)

shufflep[lst1_, lst2_] := Module[{pos2, pos},
pos2 = Range@Length@lst2 + # & /@ comb[Length@lst1 + 1, Length@lst2];
pos = Complement[Range[Length@lst1 + Length@lst2], #]~Join~# & /@ pos2;
Permute[lst1~Join~lst2, InversePermutation@#] & /@ pos]

• That is quite a clean method, and I see the possibility to make this code a bit cleaner still. May I edit your post with these changes? If you don't like them you can simply "revert" the edit to undo them. Feb 4, 2014 at 10:08
• Sure. You are welcome to edit! Feb 4, 2014 at 10:09
• Thanks a lot @Mr.Wizard! This is indeed cleaner! Feb 4, 2014 at 10:13
• Notes: (1) If the lists that are assembled are short Insert is not expensive; this is a good method. (2) Because you Fold over lst1, lst1 should be the shorter of the two lists if there is one; this will be much faster than the other way around. Feb 4, 2014 at 10:33
• Please see the timings I added to my answer. In some cases your method is far and away the fastest method yet posted, when applied correctly. Well done. Feb 4, 2014 at 10:51

Quite horribly inefficient:

l = {{a, b}, {c, d}};

Intersection @@ (Cases[Permutations[Flatten[l]],Riffle[#, ___, {1, -1, 2}]] & /@ l)


{{a, b, c, d}, {a, c, b, d}, {a, c, d, b}, {c, a, b, d}, {c, a, d, b}, {c, d, a, b}}

... and it only gets worse for the general case with $n$ lists:

l = {{a, b}, {c, d}, {e, f}};

Intersection @@ (Cases[Permutations[Flatten[l]],Riffle[#, ___, {1, -1, 2}]] & /@ l)


{{a, b, c, d, e, f}, {a, b, c, e, d, f}, {a, b, c, e, f, d}, {a, b, e, c, d, f}, {a, b, e, c, f, d}, {a, b, e, f, c, d}, {a, c, b, d, e, f}, {a, c, b, e, d, f}, {a, c, b, e, f, d}, {a, c, d, b, e, f}, {a, c, d, e, b, f}, {a, c, d, e, f, b}, {a, c, e, b, d, f}, {a, c, e, b, f, d}, {a, c, e, d, b, f}, {a, c, e, d, f, b}, {a, c, e, f, b, d}, {a, c, e, f, d, b}, {a, e, b, c, d, f}, {a, e, b, c, f, d}, {a, e, b, f, c, d}, {a, e, c, b, d, f}, {a, e, c, b, f, d}, {a, e, c, d, b, f}, {a, e, c, d, f, b}, {a, e, c, f, b, d}, {a, e, c, f, d, b}, {a, e, f, b, c, d}, {a, e, f, c, b, d}, {a, e, f, c, d, b}, {c, a, b, d, e, f}, {c, a, b, e, d, f}, {c, a, b, e, f, d}, {c, a, d, b, e, f}, {c, a, d, e, b, f}, {c, a, d, e, f, b}, {c, a, e, b, d, f}, {c, a, e, b, f, d}, {c, a, e, d, b, f}, {c, a, e, d, f, b}, {c, a, e, f, b, d}, {c, a, e, f, d, b}, {c, d, a, b, e, f}, {c, d, a, e, b, f}, {c, d, a, e, f, b}, {c, d, e, a, b, f}, {c, d, e, a, f, b}, {c, d, e, f, a, b}, {c, e, a, b, d, f}, {c, e, a, b, f, d}, {c, e, a, d, b, f}, {c, e, a, d, f, b}, {c, e, a, f, b, d}, {c, e, a, f, d, b}, {c, e, d, a, b, f}, {c, e, d, a, f, b}, {c, e, d, f, a, b}, {c, e, f, a, b, d}, {c, e, f, a, d, b}, {c, e, f, d, a, b}, {e, a, b, c, d, f}, {e, a, b, c, f, d}, {e, a, b, f, c, d}, {e, a, c, b, d, f}, {e, a, c, b, f, d}, {e, a, c, d, b, f}, {e, a, c, d, f, b}, {e, a, c, f, b, d}, {e, a, c, f, d, b}, {e, a, f, b, c, d}, {e, a, f, c, b, d}, {e, a, f, c, d, b}, {e, c, a, b, d, f}, {e, c, a, b, f, d}, {e, c, a, d, b, f}, {e, c, a, d, f, b}, {e, c, a, f, b, d}, {e, c, a, f, d, b}, {e, c, d, a, b, f}, {e, c, d, a, f, b}, {e, c, d, f, a, b}, {e, c, f, a, b, d}, {e, c, f, a, d, b}, {e, c, f, d, a, b}, {e, f, a, b, c, d}, {e, f, a, c, b, d}, {e, f, a, c, d, b}, {e, f, c, a, b, d}, {e, f, c, a, d, b}, {e, f, c, d, a, b}}

• @David Relative efficiency may be overrated but computational complexity is not. This fails on the latter. Feb 7, 2014 at 11:35
• @Mr.Wizard O(no!) Feb 9, 2014 at 19:53
• @YvesKlett LOL -- I hope that was written with a smile on your face as I imagine. :-) Feb 9, 2014 at 20:05
• @Mr.Wizard indeed :D Feb 9, 2014 at 20:51

One more way!

list1 = {a, b};
list2 = {c, d};
list = Range[2 Length@list1];
first = Select[Permutations[list, {(Length@list1)}], Sort[#] === # &];
second = Complement[list, #] & /@ first;
Normal[SparseArray[#]] & /@
(Join @@@ (Transpose@{Thread[Rule[#, list1]] & /@ first,Thread[Rule[#, list2]] & /@ second}))


{{a, b, c, d}, {a, c, b, d}, {a, c, d, b}, {c, a, b, d}, {c, a, d, b}, {c, d, a, b}}

I'm late to arrive. How about a one-liner using pattern matching:

Cases[Cases[Permutations[{a,b,c,d}], {___, a, ___, b, ___}], {___, c, ___, d, ___}]


Automated form:

MrShuffles[lis1_, lis2_] := Cases[ Cases[Permutations[Join[lis1, lis2]],
Riffle[lis1, ___, {1, -1, 2}]], Riffle[lis2, ___, {1, -1, 2}]]


Example:

MrShuffles[{a, b}, {c, d}]


{{a, b, c, d}, {a, c, b, d}, {a, c, d, b}, {c, a, b, d}, {c, a, d, b}, {c, d, a, b}}

MrShuffles[{a, b, c}, {d, e}]


{{a, b, c, d, e}, {a, b, d, c, e}, {a, b, d, e, c}, {a, d, b, c, e}, {a, d, b, e, c}, {a, d, e, b, c}, {d, a, b, c, e}, {d, a, b, e, c}, {d, a, e, b, c}, {d, e, a, b, c}}

shuffle[list1_, list2_] :=
Module[{i = 0, len1 = Length[list1], len2 = Length[list2]},
Fold[Insert[#1, list2[[Mod[++i, len2, 1]]], #2] &, list1, #] & /@
Subsets[Range[len1 + len2], {len2}]]


• @Mr.Wizard Well,I'm little blush to see you do this one by one for me,I will edit those image link in future. :)
– yode
Mar 13, 2017 at 1:41
• I was testing to see that this method works properly. See my update to meta.mathematica.stackexchange.com/a/2113/121 -- you can do the rest. Mar 13, 2017 at 1:46