Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.


{a,b}, {c,d}

Is there any easy way to do this?

share|improve this question
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. :-) – Mr.Wizard Feb 3 '14 at 10:37
Ok, I'll wait :) – A friendly helper Feb 3 '14 at 10:43

6 Answers 6

up vote 15 down vote accepted


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 *)


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}]

enter image description here

  • 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 original question diverges from the definition of "shuffle product" therein.

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.

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

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 of these functions as well as Simon's shuffles and rasher's shufflem, performed in v7:

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

f[s1, s2]        // Length // Timing
shuffles[s1, s2] // Length // Timing
shufflem[s1, s2] // Length // Timing
shuffleW[s1, s2] // Length // Timing
{1.997, 705432}

{1.576, 705432}

{0.843, 705432}

{0.624, 705432}

rasher's method (and my refactor of it) is 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 // Timing
{2.091, 705432}

Now with disparate lists:

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

Fold[insertElem, {{s2, 0}}, s1][[All, 1]] // Length // Timing
{5.164, 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 // Timing
{0.094, 91881}

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

f[s1, s2]        // Length // Timing
shuffles[s1, s2] // Length // Timing
shufflem[s1, s2] // Length // Timing
shuffleW[s1, s2] // Length // Timing
{0.39, 91881}

{0.436, 91881}

{0.406, 91881}

{0.359, 91881}
share|improve this answer
Very nice, that's what I had in mind. Thanks! – A friendly helper Feb 3 '14 at 10:33
@Mr.Wizard: Curious what my real entry does in your environment... – ciao Feb 4 '14 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! :-) – Mr.Wizard Feb 4 '14 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 '14 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 :) – A friendly helper Feb 4 '14 at 9:44

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}}

share|improve this answer
+1 Efficiency is overrated. – David Carraher Feb 6 '14 at 9:28
@David Relative efficiency may be overrated but computational complexity is not. This fails on the latter. – Mr.Wizard Feb 7 '14 at 11:35
@Mr.Wizard O(no!) – Yves Klett Feb 9 '14 at 19:53
@YvesKlett LOL -- I hope that was written with a smile on your face as I imagine. :-) – Mr.Wizard Feb 9 '14 at 20:05
@Mr.Wizard indeed :D – Yves Klett Feb 9 '14 at 20:51

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}} *)


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]]
share|improve this answer
I refactored shuffles. If you don't like the changes just revert it. – Mr.Wizard Feb 3 '14 at 11:56

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}}

share|improve this answer

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}]},
       Ordering /@ 
        Ordering /@ 
          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


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:

enter image description here

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}}  *)
share|improve this answer
:) ... it does only work with unassigned variables, though. – Yves Klett Feb 4 '14 at 7:37
Yes, it was meant as a goof. Adding my real entry now... – ciao Feb 4 '14 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! – Simon Woods Feb 4 '14 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 '14 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 – Simon Woods Feb 4 '14 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]
share|improve this answer
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. – Mr.Wizard Feb 4 '14 at 10:08
Sure. You are welcome to edit! – Yi Wang Feb 4 '14 at 10:09
Thanks a lot @Mr.Wizard! This is indeed cleaner! – Yi Wang Feb 4 '14 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. – Mr.Wizard Feb 4 '14 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. – Mr.Wizard Feb 4 '14 at 10:51

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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