30
$\begingroup$

I wonder what is the fastest method to generate derangements?

Both the Combinatorica function and Martin Ender's answer to Permutations of lists of fixed even numbers are based on filtering the output of Permutations. Let's compare them.

s = Range @ 9;

Needs["Combinatorica`"] // Quiet
Derangements[s] // Length // RepeatedTiming

Select[Permutations[s], FreeQ[s - #, 0] &] // Length // RepeatedTiming

(* {5.056, 133496} *)

(* {1.172, 133496} *)

Martin's code improves handily on the package code. Can we do better?

Can we generate these directly and avoid filtering entirely?

$\endgroup$
  • 4
    $\begingroup$ This is rather unexpected. $\endgroup$ – Martin Ender Mar 23 '17 at 9:19
  • $\begingroup$ Check this 2002 MathGroup thread in particular I think the suggestion by Hartmut Wolf. $\endgroup$ – Daniel Lichtblau Mar 23 '17 at 14:11
  • 2
    $\begingroup$ @MartinEnder Do you mean that the Combinatorica function is slow? That package was written when Mathematica had fewer choices, and the package seems to have been written by mathematicians rather than expert programmers, by which I mean the algorithms are often advanced but the implementations can be less than optimal. I learned long ago (my Project Euler days) to reimplement functions when speed was critical. $\endgroup$ – Mr.Wizard Mar 23 '17 at 14:39
  • 2
    $\begingroup$ I should add that I misunderstood the intent here. The reference link I gave was for generating random derangements quickly, not for generating all derangements quickly. $\endgroup$ – Daniel Lichtblau Mar 23 '17 at 15:26
  • 2
    $\begingroup$ @yode I seem to recall that question being asked before on this site. As far as I can recall that is not directly possible, however you can load the package and export the function you need. Would you like me to post an example of that, if I don't find an existing Q&A that describes this? $\endgroup$ – Mr.Wizard Mar 25 '17 at 2:17
27
$\begingroup$

Chunks of derangements

Since I've already written library link code generating permutations, generating derangements requires just few tweaks:

/* derangements.c */
#include "WolframLibrary.h"

DLLEXPORT mint WolframLibrary_getVersion() {
    return WolframLibraryVersion;
}
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) {
    return LIBRARY_NO_ERROR;
}
DLLEXPORT void WolframLibrary_uninitialize(WolframLibraryData libData) {}

DLLEXPORT int nextDerangementsChunk(
        WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res
) {
    /* Values tensor: list of integers in original order. */
    MTensor valuesT = MArgument_getMTensor(Args[0]);
    /* Actual data of values tensor. */
    mint* values = libData->MTensor_getIntegerData(valuesT);
    /* Number of elements in list. */
    mint n = libData->MTensor_getDimensions(valuesT)[0];

    /* Ordered values tensor: list of integers in non-increasing order. */
    MTensor orderedValuesT = MArgument_getMTensor(Args[1]);
    /* Actual data of ordered values tensor. */
    mint* orderedValues = libData->MTensor_getIntegerData(orderedValuesT);

    /* `stateT` tensor: `{next1, next2, ..., head, ref}`. */
    MTensor stateT = MArgument_getMTensor(Args[2]);
    /*
     * First `n` elements of `next` array contain indices of next nodes
     * in emulated linked list. Other elements of `stateT` tensor are used
     * only through direct pointers.
     */
    mint* next = libData->MTensor_getIntegerData(stateT);
    /* Pointer to index of head node. */
    mint* head = next + n;
    /* Pointer to index of reference node. */
    mint* ref = head + 1;

    /* Number of permutations in returned chunk. */
    mint chunkSize = MArgument_getInteger(Args[3]);
    /* Dimensions of returned `chunk` tensor. */
    mint chunkDims[2] = {chunkSize, n};
    /* 2 dimentional tensor with chunk of permutations to be returned. */
    MTensor chunkT;
    libData->MTensor_new(MType_Integer, 2, chunkDims, &chunkT);
    /* Actual data of the chunk tensor. */
    mint* chunk = libData->MTensor_getIntegerData(chunkT);

    mint i;
    for (i = 0; i < chunkSize; i++) {
        /*
         * Based on:
         * Aaron Williams. 2009. Loopless generation of multiset permutations
         * using a constant number of variables by prefix shifts.
         * http://webhome.csc.uvic.ca/~haron/CoolMulti.pdf
         */
        mint afterRef = next[*ref];
        mint localRef;
        if (next[afterRef] < n && orderedValues[*ref] >= orderedValues[next[afterRef]]) {
            localRef = afterRef;
        } else {
            localRef = *ref;
        }
        mint newHead = next[localRef];

        next[localRef] = next[newHead];
        next[newHead] = *head;

        if (orderedValues[newHead] < orderedValues[*head]) {
            *ref = newHead;
        }
        *head = newHead;

        /* Populate i-th permutation in chunk. */
        mint j, index;
        for (j = 0, index = *head; j < n; j++) {
            if (orderedValues[index] == values[j]) {
                /*
                 * This is not a derangement. Decrement index so that i-th place
                 * will be populated with next permutation.
                 */
                i--;
                break;
            }
            chunk[i*n + j] = orderedValues[index];
            index = next[index];
        }
    }

    /* Return control over state tensor back to Wolfram Language. */
    libData->MTensor_disown(stateT);

    /* Set chunk tensor as returned value of LibraryFunction. */
    MArgument_setMTensor(Res, chunkT);

    return LIBRARY_NO_ERROR;
}

Save above code in derangements.c file in same directory as current notebook, or paste it as a string, instead of {"derangements.c"}, as first argument of CreateLibrary in code below. Pass, in "CompileOptions", appropriate optimization flags for your compiler, the ones below are for GCC.

Needs@"CCompilerDriver`"
SetDirectory@NotebookDirectory[];
CreateLibrary[{"derangements.c"}, "derangements"(*,
    "CompileOptions" -> "-Wall -march=native -O3"*)
]
nextDerangementsChunk = LibraryFunctionLoad[%, "nextDerangementsChunk",
    {{Integer, 1, "Constant"}, {Integer, 1, "Constant"}, {Integer, 1, "Shared"}, Integer},
    {Integer, 2}
]

nextDerangementsChunk accepts four arguments: list of integers for which we want to generate derangements, list of same integers but in non-increasing order, list representing "state" of generator, and number of derangements in returned chunk. "Generator state" is described more precisely in my permutations post.

As a usage example let's generate derangements of {2, 1, 4, 1, 3} in two 5-element, and one 2-element chunks:

values = {2, 1, 4, 1, 3};
ordered = Reverse@Sort@values;
state = Join[Range@Length@values, {0, Length@values - 2}];
nextDerangementsChunk[values, ordered, state, 5]
nextDerangementsChunk[values, ordered, state, 5]
nextDerangementsChunk[values, ordered, state, 2]
(* {{1, 4, 3, 2, 1}, {3, 4, 1, 2, 1}, {4, 3, 1, 2, 1}, {1, 4, 1, 3, 2}, {1, 3, 1, 4, 2}} *)
(* {{1, 4, 2, 3, 1}, {4, 2, 1, 3, 1}, {1, 3, 2, 4, 1}, {1, 2, 3, 4, 1}, {3, 2, 1, 4, 1}} *)
(* {{1, 3, 1, 2, 4}, {1, 2, 1, 3, 4}} *)

Currently nextDerangementsChunk does no checks of given arguments, passing inconsistent arguments can lead to infinite loop, or kernel crash.


Number of derangements

Above algorithm requires explicit number of expected derangements, so we need to calculate in advance how many derangements, of our list, are there.

In general number of derangements is given by certain integral of product of Laguerre polynomials. For list of unique elements there's a built-in function that gives number of derangements: Subfactorial.

We'll use Subfactorial function for mentioned special case and Laguerre polynomials in general:

multiSubfactorial = With[{tallied = Tally@#},
    If[tallied === {{1, Length@#}},
        Subfactorial@Length@#
    (* else *),
        With[
            {coeffs = Block[{x}, 
                CoefficientList[Times @@ (LaguerreL[#1, x]^#2 & @@@ tallied), x]
            ]},
            Abs@Total[Factorial@Range[0, Length@coeffs - 1] coeffs]
        ]
    ]
]&;

All derangements

derangements // ClearAll
derangements[empty:_[]] := {empty}
derangements[_[_]] = {};
derangements[list_List /; VectorQ[Unevaluated@list, IntegerQ]] :=
    With[{n = Length@list},
        nextDerangementsChunk[
            list,
            Reverse@Sort@list,
            Join[Range@n, {0, n - 2}],
            multiSubfactorial@Tally[list][[All, 2]]
        ]
    ]
derangements[expr_ /; Not@AtomQ@Unevaluated@expr] :=
    With[{n = Length@expr, list = List @@ expr},
    With[{tallied = Sort@Tally@list},
    With[{unique = Head@expr @@ tallied[[All, 1]]},
        unique[[#]] & /@ nextDerangementsChunk[
            Lookup[PositionIndex@tallied[[All, 1]], list][[All, 1]],
            Flatten@Reverse@
                MapIndexed[ConstantArray[First@#2, Last@#1]&, tallied],
            Join[Range@n, {0, n - 2}],
            multiSubfactorial@tallied[[All, 2]]
        ]
    ]]]

Check that it generates same derangements as other methods for integer lists:

And @@ (Function[s, Sort@derangements@s === Sort@Select[Permutations@s, FreeQ[s - #, 0] &]] /@ Join @@ (Tuples[Range@#, #] & /@ Range@6))
(* True *)

and symbolic lists:

ClearAll[f]
And @@ (Function[s, Sort@derangements@s === Sort@Select[Permutations@s, FreeQ[s - #, 0] &]] /@ Join @@ (Tuples[f /@ Range@#, #] & /@ Range@6))
(* True *)

Benchmarks

For list of unique integers, from OP, derangements is ten times faster than Pick:

s = Range@9;
(res1 = Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1]&@Permutations[s]) // MaxMemoryUsed // RepeatedTiming
(res2 = derangements@s) // MaxMemoryUsed // RepeatedTiming
Sort@res1 === Sort@res2
(* {0.052, 78385160} *)
(* {0.0043, 9613720} *)
(* True *)

Speed and memory usage difference is bigger for multisets with multiple duplicates where ratio of derangements to permutations can be much lower than 1/E.

s = Join[ConstantArray[1, 6], Range[2, 7]];
(res1 = Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1] &@Permutations[s]) // MaxMemoryUsed // RepeatedTiming
(res2 = derangements@s) // MaxMemoryUsed // RepeatedTiming
Sort@res1 === Sort@res2
(* {0.13, 191603344} *)
(* {0.0054,   70728} *)
(* True *)

s = Join[ConstantArray[1, 7], ConstantArray[2, 5], Range[3, 5]];
(res1 = Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1] &@Permutations[s]) // MaxMemoryUsed // RepeatedTiming
(res2 = derangements@s) // MaxMemoryUsed // RepeatedTiming
Sort@res1 === Sort@res2
(* {0.518, 778380768} *)
(* {0.016,    182984} *)
(* True *)
$\endgroup$
  • 1
    $\begingroup$ Pretty awesome! $\endgroup$ – Jacob Akkerboom Mar 23 '17 at 20:40
  • $\begingroup$ I still don't have a C compiler installed. Any chance I could get a Windows x64 binary from you? $\endgroup$ – Mr.Wizard Mar 23 '17 at 23:26
  • $\begingroup$ @Mr.Wizard Sorry, but I don't have access to Windows. $\endgroup$ – jkuczm Mar 24 '17 at 22:55
  • 1
    $\begingroup$ @JacobAkkerboom It's actually a "brute force" approach. It does generate all permutations, just doesn't store invalid ones in memory. It's fast because it's written in C, so algorithmic incompetence is well hidden behind low-levelness of implementation. $\endgroup$ – jkuczm Mar 24 '17 at 23:01
23
$\begingroup$

This is the fastest method I have come up with:

s = Range @ 9;

Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1] & @ Permutations[s] // 
  Length // RepeatedTiming
{0.0408, 133496}
$\endgroup$
  • 4
    $\begingroup$ such a neat way...brought me knowledge (and joy after a bad day) +1:) $\endgroup$ – ubpdqn Mar 23 '17 at 11:00
  • $\begingroup$ @ubpdqn Thank you. :-) $\endgroup$ – Mr.Wizard Mar 23 '17 at 13:28
  • $\begingroup$ Have you tried a non-discarding approach and found it to be slower, or did you just look into faster discarding approaches so far? $\endgroup$ – Martin Ender Mar 23 '17 at 13:39
  • $\begingroup$ @MartinEnder I only considered faster discarding (filtering) methods. I don't even know how one would directly generate these without discarding, hence Can we generate these directly and avoid filtering entirely? which you'll note I did not self-answer. $\endgroup$ – Mr.Wizard Mar 23 '17 at 14:44
  • $\begingroup$ @Mr.Wizard Right, I just wasn't sure whether "this is the fastest method I have come up with" included attempting non-filtering methods or not. :) I think it should be possible to generate them directly recursively, based on the recurrence of the subfactorial numbers, but I'm not sure how good I am with writing fast Mathematica code. I'll give it a try and see what I come up with. $\endgroup$ – Martin Ender Mar 23 '17 at 14:46
17
$\begingroup$

Here is one way to generate them directly: it is based on a way to generate all permutations but discards invalid ones early:

derangements[{}, ___] = {{}};
derangements[list_List, orig_List] := 
 Union @@ 
   (Prepend[#] /@ derangements[DeleteCases[list, #, 1, 1], Rest@orig] &) /@ 
     DeleteCases[list, First@orig]
derangements[list_List] := derangements[list, list]

Basically, we generate them recursively, keeping track of the original list (or what's left of it) to make sure that don't use an element in the same position that it originally appeared in.

For each element i in the input list that isn't equal to the currently forbidden one, (DeleteCases[list, First @ orig]), we recursively call derangements on list with that element removed and the remainder of orig, and prepend i to each of the resulting derangements.

This in itself is quite slow, and comes somewhere between Combinatorica and my original filtering approach:

s = Range@9;
Needs["Combinatorica`"] // Quiet
Derangements[s] // Length // AbsoluteTiming
Select[Permutations[s], FreeQ[s - #, 0] &] // Length // AbsoluteTiming
Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1] &@Permutations[s] // Length // AbsoluteTiming
derangements[s] // Length // AbsoluteTiming
{4.91098, 133496}
{1.09925, 133496}
{0.0509919, 133496}
{2.44123, 133496}

However, since this is recursive, it's possible to memoise this approach, after which it's only 3-4 times slower than Mr. Wizard's filtering approach:

derangementsMemo[{}, ___] = {{}};
derangementsMemo[list_List, orig_List] := derangementsMemo[list, orig] = 
 Union @@ 
   (Prepend[#] /@ derangementsMemo[DeleteCases[list, #, 1, 1], Rest@orig] &) /@ 
     DeleteCases[list, First@orig]
derangementsMemo[list_List] := derangementsMemo[list, list]

derangementsMemo[s] // Length // AbsoluteTiming
{0.179871, 133496}

That said, I'm not sure how feasible memoisation of this kind of combinatorial function is in the long run. You might want to clean out the cache after each use of derangementsMemo.

I'm sure some more experienced Mathematica users would be able to optimise this approach even further.

I believe there must be also at least one other strategy for generating the derangements recursively, based on the recurrence of the subfactorial numbers:

$$ !n = (n-1)(!(n-1) + !(n-2)) $$

However, I haven't yet tried to wrap my head around it, but will report back if I do.


Packing

Mr. Wizard's solution is fast because it utilizes functions optimized for packed arrays. Above recursive solutions can, to some extend, also use packed arrays.

To take advantage of packing we must make sure that calls ending recursion return packed arrays, and replace unpacking Prepend[#] /@ ... with PadLeft. Adding some cosmetic changes and localization of memoisation we get:

ClearAll[pad, derangementsInternal, derangementsPacked]
pad[_, _]@{} = {};
pad[i_, x_]@l_ := PadLeft[l, {Length@l, i}, x];

derangementsInternal[1] =
  If[Last@#1 === Last@#2, {}, Developer`ToPackedArray@{#1}]&;
derangementsInternal[i_][list_, orig_] := derangementsInternal[i][list, orig] =
  Join @@
    (pad[i, #]@derangementsInternal[i - 1][DeleteCases[list, #, 1, 1], orig]&) /@
      Complement[list, orig[[-i ;; -i]]]

derangementsPacked[{}] = {{}};
derangementsPacked[list_List] := Internal`InheritedBlock[{derangementsInternal},
  derangementsInternal[Length@list][list, list]
]

which is faster and more memory efficient than Pick-based solution:

s = Range@9;
(res1 = Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1] &@Permutations[s]) // MaxMemoryUsed // RepeatedTiming
(res2 = derangementsPacked@s) // MaxMemoryUsed // RepeatedTiming
res1 === res2
{0.051, 78384288}
{0.032, 34599000}
True

Advantage of this approach, over filtering ones, grows with number of duplicates:

s = Join[ConstantArray[1, 7], ConstantArray[2, 5], Range[3, 5]];
(res1 = Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1] &@Permutations[s]) // MaxMemoryUsed // RepeatedTiming
(res2 = derangementsPacked@s) // MaxMemoryUsed // RepeatedTiming
res1 === res2
{0.63, 778380768}
{0.0019,  824784}
True

For large enough number of duplicates it can beat library link filtering, which for above example is ten times slower: {0.016, 182984}.

$\endgroup$
  • $\begingroup$ Very nice. The memoized code is more competitive that I thought it might be. $\endgroup$ – Mr.Wizard Mar 23 '17 at 23:21
  • $\begingroup$ Very nice (sorry, the system only let's me upvote once). I had thought about the "other" strategy of recursion you have mentioned and to me it builds upon some kind of graph-representation: first element -> another element -> another element. Once the cycle is completed (last element chooses the first), the problem reduces to adding all derangements for what elements are left. I will try to post this eventually but it is too slow yet... $\endgroup$ – gwr Mar 27 '17 at 13:41
9
$\begingroup$

Here is a straightforward compiled implementation of Knuth's "Algorithm X" for lexicographically generating restricted permutations, specialized to the derangement case:

derange = Compile[{{list, _Integer, 1}},
                  Module[{n = Length[list], a, db, k, l, p, q, t, u},
                         a = Range[n]; l = Append[a, 0];
                         u = Table[0, {n}];
                         db = Internal`Bag[Most[{1}]]; k = 1;

                         Label[2]; p = 0; q = l[[1]];

                         Label[3]; a[[k]] = q;
                         If[a[[k]] == k, Goto[5]];

                         If[k == n,
                            t = list[[a]]; Internal`StuffBag[db, t, 1],
                            u[[k]] = p; l[[p + 1]] = l[[q + 1]]];
                         If[k == n, Goto[6]];
                         k++; Goto[2];

                         Label[5]; p = q; q = l[[p + 1]];
                         If[q != 0, Goto[3]];

                         Label[6]; k--;
                         If[k == 0, Goto[7]];
                         p = u[[k]]; q = a[[k]]; l[[p + 1]] = q;
                         Goto[5];

                         Label[7]; Partition[Internal`BagPart[db, All], n]],
                  RuntimeOptions -> "Speed"]

It is slightly slower than the Wizard's method, and (gasp!) uses Label[]/Goto[], but does not generate unneeded permutations. Rewriting this not to use Label[]/Goto[] ought to be a stiff challenge.

wizardDerange[s_List] :=
      Pick[#, Unitize[Times @@ (Transpose[#] - s)], 1] & @ Permutations[s]

RepeatedTiming[MaxMemoryUsed[wizardDerange[Range[9]]]]
   {0.063, 78384288}

RepeatedTiming[MaxMemoryUsed[derange[Range[9]]]]
   {0.11, 39655464}
$\endgroup$
  • 1
    $\begingroup$ You've been warned. $\endgroup$ – Martin Ender Mar 24 '17 at 8:42
  • 1
    $\begingroup$ Indeed: "A number of rules have been discovered, violation of which will either seriously impair or totally destroy the intellectual manageability of the program. […] Examples are the exclusion of goto-statements and of procedures with more than one output parameter." - Dijkstra . Just kidding, I just wanted to show off I paid attention to your link, my own code is beyond sinful. $\endgroup$ – Jacob Akkerboom Mar 24 '17 at 9:15
  • 1
    $\begingroup$ (Yes, I know what Dijkstra said, you guys. :P) I tried figuring out how to de-Goto[] "Algorithm X", but I gave up. Maybe somebody better than me can try. $\endgroup$ – J. M. is away Mar 24 '17 at 9:20
  • 1
    $\begingroup$ @J.M.: The canonical way to de-Goto[] this algorithm would be to recognize that you're creating a simple state machine, declare symbols for the states, Switch on the current state within a loop. It's debatable how much this increases the readability/maintainability. :-) $\endgroup$ – wchargin Mar 25 '17 at 16:15
8
$\begingroup$

Recursive Derangements

Here is my approach at doing this recursively:

ClearAll[pickSlot];
pickSlot[ elem_Integer, {slot_List}, taken_List ] := Sow @ Flatten[ { taken, elem }]
pickSlot[ elem_Integer, slots__List , taken_List ] := With[
   {
       remaining = DeleteCases[ Rest@slots, elem, 2 ] (* remove element from other lists *)
   },
   Scan[
       pickSlot[
           #,
           remaining, 
           {taken, elem }  (* append the picked element to results *)
       ] &,
       First @ remaining
   ]
]

ClearAll[derangements];
derangements[ s_?VectorQ ] := Module[
    {
        $slots 
    },

    $slots =  MapIndexed[
        DeleteCases[ #1, First@ #2] &, 
        Map[ s &, s ]
    ];

    Flatten[
        Rest @ Reap @ Scan[
            pickSlot[#, $slots, {}] &,
            First @ $slots
        ],
        2
    ]
]

With[
    { s = Range @ 9 },
    RepeatedTiming @ Length @ derangements[s]
]

{1.68, 133 496}

How it works

The basic strategy is to start out with a list $slots of possible positions for each element, where we remove the original position from each list.

We then use the recursice function pickSlot[ elem, slots, taken ], which will pick elem from the first list in slots. It will append this elemnt to the list of choices (taken) and adjust the list of slots, e.g. the first list is dropped and the chosen element is removed from the other lists.

$\endgroup$
7
$\begingroup$

Below we define a function findDerangements using Compile. I it's not as fast as jkuczm's answer, but I guess I'll post it anyway.

store[x_] :=
  (
   res[[nRes]] = x;
   nRes++
   );
fbRuleDef[list_, pos_] :=
 Block[{res},
  res = 0;
  Do[
   If[
    list[[ii]] == 1
    ,
    res = ii;
    Break[];
    ]
   ,
   {ii, pos + 1, Length@list}
   ];
  res
  ]
heldBlock3 =
  Hold@
    Block[
     {min}
     ,
     min = fbRuleDef[leftOvers, 0];
     If[
      min == depth 
      ,
      If[
       depth == nn,
       {0, 0}
       ,
       {1, fbRuleDef[leftOvers, min]}
       ]
      ,
      {1, min}
      ]
     ] /. DownValues@fbRuleDef;
heldBlock4 =
  Hold@
    Block[
     {min}
     ,
     min = fbRuleDef[leftOvers, cur];
     If[
      min > 0,
      If[
       min == depth 
       ,
       min = fbRuleDef[leftOvers, min];
       If[
        min > 0
        ,
        {1, min}
        ,
        {0, 0}
        ]
       ,
       {1, min}
       ]
      ,
      {0, 0}
      ]
     ] /. DownValues[fbRuleDef];
feasAndMinInl4 =
  Function[
    Null,
    Compile[
     {{leftOvers, _Integer, 1}, {depth, _Integer}, {cur, _Integer}, { 
       nn, _Integer}}
     ,
     #
     ,
     CompilationTarget -> "C"
     ],
    HoldAll
    ] @@ heldBlock4;
feasAndMinInl3 =
  Function[Null,
    Compile[
     {{leftOvers, _Integer, 1}, {depth, _Integer}, {nn, _Integer}}
     ,
     #
     ,
     CompilationTarget -> "C"
     ]
    ,
    HoldAll
    ] @@ heldBlock3;
SetDelayed @@ {
  Unevaluated@feasAndMinInlRuleDef3[leftOvers_, depth_, nn_],
  Unevaluated @@ heldBlock3
  }
SetDelayed @@ {
  Unevaluated@feasAndMinInlRuleDef4[leftOvers_, depth_, cur_, nn_],
  Unevaluated @@ heldBlock4
  }
heldBlockInlRef =
  Hold@
    Block[{depth, iters, test, leftOvers, dir, resLen, res, nRes, 
      feasibleLeftOversExist, min, tempTens},
     depth = 1;
     iters = Table[0, nn];
     test = True;
     leftOvers = Table[1, nn];
     dir = 1;
     resLen = Subfactorial[nn];
     res = Table[0, {resLen}, {nn}];
     nRes = 1;
     feasibleLeftOversExist = 0;

     While[
      test
      ,
      If[
       dir == 0 ,
       tempTens = 
        feasAndMinInlRuleDef4[leftOvers, depth, iters[[depth]], nn];
       feasibleLeftOversExist = tempTens[[1]]; min = tempTens[[2]] ;
       If[
        feasibleLeftOversExist == 1
        ,
        iters[[depth]] = min;
        leftOvers[[iters[[depth]]]] = 0;
        dir = 1; depth++;
        ,
        If[
         depth > 1,

         dir = 0; depth--;
         leftOvers[[iters[[depth]]]] = 1;
         ,
         test = False
         ]
        ]
       ,
       (*dir \[Equal]1*)

       tempTens = feasAndMinInlRuleDef3[leftOvers, depth, nn];
       feasibleLeftOversExist = tempTens[[1]]; min = tempTens[[2]] ;
       If[
        feasibleLeftOversExist == 1
        ,
        If[
         depth == nn
         ,
         iters[[depth]] = min;
         store[iters];
         dir = 0; depth--;
         leftOvers[[iters[[depth]]]] = 1;
         ,
         iters[[depth]] = min;
         leftOvers[[iters[[depth]]  ]] = 0;
         dir = 1; depth++
         ]
        ,
        dir = 0; depth--;
        leftOvers[[iters[[depth]]]] = 1;
        ]
       ]
      ];
     res
     ] /. 
   Join[ DownValues[store], DownValues[feasAndMinInlRuleDef3], 
    DownValues[feasAndMinInlRuleDef4]];

findDerangements =
  Function[
    Null,
    Compile[
     {{nn, _Integer}}
     ,
     #
     ,
     CompilationTarget -> "C"
     ,
     CompilationOptions -> {"ExpressionOptimization" -> True}
     ,
     RuntimeOptions -> "Speed"

     ]
    , HoldAll
    ] @@ heldBlockInlRef;

I get the following timings, slightly disappointing for me, but oh well :)

s = Range@9;
(res1 = Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1] &@
     Permutations[s]) // MaxMemoryUsed // RepeatedTiming
(res2 = derangements@s) // MaxMemoryUsed // RepeatedTiming    
(res3=findDerangements[9])//MaxMemoryUsed//RepeatedTiming
Sort@res1===Sort@res2=== res3
{0.043,78384288}
{0.0051,9612408}
{0.044,28836192}
True
$\endgroup$
5
$\begingroup$

Cycles!

I'm arriving late to this party, but thought I would post my alternative anyway :) I doubt it will be able to beat the speed of the best of the other answers, but my version has some other benefits, such as short code, no Compileing needed to get decent speed, and transparent relation to the recursive equation describing subfactorials.

In Mathematica, the documentation typically refers to two ways of representing permutations: simple permutation lists, as all other answers to this question have been dealing with, and permutation Cycles. The two appear together in the documentation in a similar fashion to the way lists and SparseArrays do.

The Cycles of a derangement is simply a Cycles where all positions 1,...,n are present and no cycle has Length < 2 (no element is sent to itself). The different derangements will differ in

  1. how many cycles there are / the length of the cycles
  2. the indices that appear in the same cycle
  3. the order within a cycle up to full rotations

Cycles does some autosorting such that a) indices inside a cycle are rotated such that the smallest element is the first one, and b) cycles are sorted on their first element.

Recursion of subfactorials

The formula $$!n = (n-1)(!(n-1) + !(n-2))$$ is quite easily understood when considering Cycles. Consider the derangements of Range[2] and Range[3]. Since cycles must have Length >= 2, the only possibilities are

Cycles[{{1,2}}]
Cycles[{{1,2,3}}], Cycles[{{1,3,2}}]

We can create the derangements of Range[4] from these in the following way. First, we can add the element 4 to the second line. Since 4 is the largest element, it cannot be the first element in a cycle, but apart from this, it can appear anywhere. If n = 4, that is n-1 places where we can put the 4, immediately giving the $$(n-1)(!(n-1))$$ contribution! What are we missing? Well, placing only one element cannot change the number of cycles, since cycles like {4} is not allowed in a derangement. But we can add the cycle {3,4} to the first line, and get

Cycles[{{1,2}, {3,4}}]

In addition, we can swap 3 and any smaller index, and this would give new cycles not captured above, here

Cycles[{{3,2}, {1,4}}], Cycles[{{1,3}, {2,4}}]

which is autosorted into

Cycles[{{1,4}, {2,3}}], Cycles[{{1,3}, {2,4}}]

On the other hand, swapping 4 and something lower will just either reproduce results from adding 4 to the derangements of Range[3] or reproduce the swapping of 3, like here:

Cycles[{{4,2}, {3,1}}], Cycles[{{1,4}, {3,2}}] ==
Cycles[{{1,3}, {2,4}}], Cycles[{{1,4}, {2,3}}]

The number of ways we can swap a 3 with something smaller is n-2, and if we don't swap, we have all the n-1 terms we need from $$(n-1)(!(n-2))$$

Implementation

Here is the code (finally!) :

deranCycles[2] = {Cycles[{{1,2}}]};
deranCycles[3] = {Cycles[{{1,2,3}}], Cycles[{{1,3,2}}]};

deranCycles[n_] := deranCycles[n] =
Block[{prev, prevprev, ans1, ans2},
 prevprev = deranCycles[n - 2];
 prev = deranCycles[n - 1];

 ans1 = Flatten[
   Table[Insert[cyc, n, {1, i, j}], {cyc, prev}, {i, 
     Length@cyc[[1]]}, {j, 2, 1 + Length@cyc[[1, i]]}]];

 ans2 = Table[Insert[cyc, {n - 1, n}, {1, -1}], {cyc, prevprev}];
 ans2 = With[{rules = Table[{i -> n - 1, n - 1 -> i}, {i, n - 2}]},
   Flatten[{ans2,
     Table[cyc /. rules, {cyc, ans2}]}]
   ];
 Flatten[{ans1, ans2}]
 ]

It is quite fast, but not amazing:

wiz[n_] := 
 With[{s = Range@n}, 
  Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1] &@Permutations[s]];
AbsoluteTiming[MaxMemoryUsed[wiz[9]]]
AbsoluteTiming[MaxMemoryUsed[deranCycles[9]]]

{0.140247, 104 511 152}

{0.516991, 42 172 288}

and they are the same:

Sort[PermutationList /@ deranCycles[9]] == Sort[wiz[9]]

True

Note that we have not Compiled, none of the lists are packed, and we are using Insert and even ReplaceRepeated! So I'm guessing there is room for improvement here, but I don't have time to try it right now. If you think of something, I'd love to know!

$\endgroup$
  • $\begingroup$ Wonderful answer! And it picks up very knowledgable at what my comment at Martin Ender's answer was meant to say.(+1) $\endgroup$ – gwr Mar 29 '17 at 9:39
  • $\begingroup$ @gwr Thank you! $\endgroup$ – Marius Ladegård Meyer Mar 29 '17 at 10:17
3
$\begingroup$

Chunks of (partial) derangements

A library link function generating (partial) derangements directly. Conceptually this is iterative version of recursive algorithm from Martin Ender's answer.

It is based on special purpose submultiset data structure, that allows efficient iteration over subsets of given multiset, and over elements of those subsets. submultiset supports O(1) removing and restoring of multiset elements and keeps numeric index uniquely identifying current subset among all possible subsets of multiset.

Algorithm fills preallocated result array while iterating over subsets of multiset. It caches starting and ending row for each subset, using mentioned unique index, so that those results can be copied, without recalculation, if subset is encountered again, which is a lightweight form of memoization.

Mutable state of submultiset is kept in single array of integers which can be easily exchanged between C and Mathematica to generate subsequent chunks of required derangements.

/* derangements.c */
#include <limits.h>
#include <stdbool.h>
#include "WolframLibrary.h"


/* Define macro holding maximal value of `mint` type. */
#ifndef MINT_MAX
#ifdef MINT_32
#define MINT_MAX INT_MAX
#else
#define MINT_MAX LONG_MAX
#endif
#endif


typedef struct {
    /*
     * Array with mutable state of submultiset:
     * {
     *   removed_count,
     *   subset_index,
     *   head, next_0, next_1, ..., next_{all_unique_count},
     *   mult_0, mult_1, ..., mult_{all_unique_count},
     *   before_removed_0, before_removed_1, ..., before_removed_{all_count-1},
     *   removed_0, removed_1, ..., removed_{all_count-1}
     * }
     */
    mint * state;
    /* Number of elements in full multiset. */
    mint const all_count;
    /* Number of unique values in full multiset. */
    mint const all_unique_count;
    /* Array of unique values of full multiset. */
    mint const * const values;
    /* Array of subset index multipliers for consecutive unique elements of full multiset. */
    mint const * const index_multipliers;
} submultiset;


/* "Private methods" for accessing of various parts of submultiset state. */
/* Return pointer to number of elements removed from full multiset. */
static inline mint * submultiset_removed_count_ptr(submultiset const * const s) { return s->state; }
/* Return pointer to index uniquely identifying current submultiset of full multiset. */
static inline mint * submultiset_index_ptr(submultiset const * const s) { return submultiset_removed_count_ptr(s) + 1; }
/* Return array of indices of "head" and "next" nodes of current submultiset. */
static inline mint * submultiset_next(submultiset const * const s) { return submultiset_index_ptr(s) + 1; }
/* Return array of multiplicities of elements in current submultiset. */
static inline mint * submultiset_multiplicities(submultiset const * const s) { return submultiset_next(s) + s->all_unique_count + 1; }
/* Return array of indices of elements before subsequent elements removed from full multiset. */
static inline mint * submultiset_before_removed(submultiset const * const s) { return submultiset_multiplicities(s) + s->all_unique_count; }
/* Return array of indices of subsequent elements removed from full multiset. */
static inline mint * submultiset_removed(submultiset const * const s) { return submultiset_before_removed(s) + s->all_count - 1; }
/* Return pointer to index of element before current element of current submultiset. */
static inline mint * submultiset_before_current_ptr(submultiset const * const s) { return submultiset_before_removed(s) + *submultiset_removed_count_ptr(s); }
/* Return pointer to index of current element of current submultiset. */
static inline mint * submultiset_current_ptr(submultiset const * const s) { return submultiset_removed(s) + *submultiset_removed_count_ptr(s); }


/* "Accessors" */
/* Return size of state of submultiset. */
inline size_t * submultiset_state_size(submultiset const * const s) { return (sizeof *(s->state)) * (submultiset_removed(s) + s->all_count - 1 - s->state); }
/* Return number of elements removed from full multiset. */
inline mint submultiset_removed_count(submultiset const * const s) { return *submultiset_removed_count_ptr(s); }
/* Return index uniquely identifying current submultiset of full multiset. */
inline mint submultiset_index(submultiset const * const s) { return *submultiset_index_ptr(s); }
/* Return number of all subsets of full multiset. */
inline mint submultiset_subsets_number(submultiset const * const s) { return s->index_multipliers[s->all_unique_count]; }
/* Return value of current element in current submultiset. */
inline mint submultiset_current_value(submultiset const * const s) { return s->values[*submultiset_current_ptr(s)]; }
/* Return value of element removed as `removed_number` in sequence of removals leading to current submultiset. */
inline mint submultiset_removed_value(submultiset const * const s, mint const removed_number) { return s->values[submultiset_removed(s)[removed_number]]; }
/* Return value of first element in current submultiset. */
inline mint submultiset_first_value(submultiset const * const s) { return s->values[submultiset_next(s)[0]]; }
/* Return `true` if iteration over submultisets finished, return `false` otherwise. */
inline bool submultiset_is_depleted(submultiset const * const s) { return submultiset_removed_count(s) == -1; }

/* "Mutators" */
/* Change submultiset state to indicate that iteration over submultisets finished. */
inline void submultiset_deplete(submultiset * const s) { *submultiset_removed_count_ptr(s) = -1; }

/* Set next element as current, if next element exists. */
inline void submultiset_advance(submultiset * const s) {
    if (*submultiset_current_ptr(s) != -1) {
        *submultiset_before_current_ptr(s) = *submultiset_current_ptr(s);
        *submultiset_current_ptr(s) = submultiset_next(s)[*submultiset_current_ptr(s) + 1];
    }
}

/* Advance current element until it's different than given `value` or last element is reached. */
inline void submultiset_skip(submultiset * const s, mint const value) {
    while (s->values[*submultiset_current_ptr(s)] == value && *submultiset_current_ptr(s) != -1) {
        *submultiset_before_current_ptr(s) = *submultiset_current_ptr(s);
        *submultiset_current_ptr(s) = submultiset_next(s)[*submultiset_current_ptr(s) + 1];
    }
}

/* Set first element of current submultiset as current element. */
inline void submultiset_reset_current(submultiset * const s) {
    *submultiset_before_current_ptr(s) = -1;
    *submultiset_current_ptr(s) = submultiset_next(s)[0];
}

/* If current element exists remove it from given submultiset and return `true`, otherwise return `false`. */
inline bool submultiset_remove_current(submultiset * const s) {
    if (*submultiset_current_ptr(s) != -1) {
        mint const curr = *submultiset_current_ptr(s);

        /*
         * If there's more than one element with same value as current element decrease its multiplicity,
         * otherwise remove node from "linked list".
         */
        if (submultiset_multiplicities(s)[curr] > 1) {
            --(submultiset_multiplicities(s)[curr]);
        } else {
            submultiset_next(s)[*submultiset_before_current_ptr(s) + 1] = submultiset_next(s)[curr + 1];
        }
        *submultiset_index_ptr(s) += s->index_multipliers[curr];
        ++(*submultiset_removed_count_ptr(s));

        return true;
    }

    return false;
}

/* If removed elements exist restore last removed element and return `true`, otherwise return `false`. */
inline bool submultiset_restore(submultiset * const s) {
    --(*submultiset_removed_count_ptr(s));
    if (submultiset_removed_count(s) >= 0) {
        mint const curr = *submultiset_current_ptr(s);
        mint const prev = *submultiset_before_current_ptr(s);

        if (submultiset_next(s)[prev + 1] == submultiset_next(s)[curr + 1]) {
            submultiset_next(s)[prev + 1] = curr;
        } else {
            ++(submultiset_multiplicities(s)[curr]);
        }
        *submultiset_index_ptr(s) -= s->index_multipliers[curr];

        return true;
    }

    return false;
}


/* Helper functions. */
inline void minus_one_fill(mint * const arr, mint const len) {
    for (mint i = 0; i < len; ++i) {
        arr[i] = -1;
    }
}


DLLEXPORT mint WolframLibrary_getVersion() { return WolframLibraryVersion; }
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) { return LIBRARY_NO_ERROR; }
DLLEXPORT void WolframLibrary_uninitialize(WolframLibraryData libData) {}

DLLEXPORT int nextDerangementsChunk(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
    /* Values tensor: list of integers in original order. */
    MTensor const valuesT = MArgument_getMTensor(Args[0]);
    /* Unique values tensor: duplicate free list of integers in original order. */
    MTensor const uniqueT = MArgument_getMTensor(Args[1]);
    /* Tensor of subset index multipliers for consecutive elements of multiset. */
    MTensor const indexMultipliersT = MArgument_getMTensor(Args[2]);
    /* Tensor containing state of used submultiset. */
    MTensor const stateT = MArgument_getMTensor(Args[3]);
    /* Maximal size of returned chunk. */
    mint chunkSize = MArgument_getInteger(Args[4]);

    /* Number of elements in list. */
    mint const n = libData->MTensor_getDimensions(valuesT)[0];
    /* Actual data of values tensor. */
    mint * const values = libData->MTensor_getIntegerData(valuesT);

    submultiset * const s = &(submultiset) {
        .state = libData->MTensor_getIntegerData(stateT),
        .all_count = n,
        .all_unique_count = libData->MTensor_getDimensions(uniqueT)[0],
        .values = libData->MTensor_getIntegerData(uniqueT),
        .index_multipliers = libData->MTensor_getIntegerData(indexMultipliersT)
    };

    /* Tensor that will hold chunk of derangements to be returned. */
    MTensor chunkT;
    /* Actual data of chunk tensor. */
    mint * chunk;

    if (chunkSize == 0 || submultiset_is_depleted(s)) {
        /* Return empty chunk. Since it was requested, or generator was already depleted. */
        libData->MTensor_disown(stateT);
        mint const chunkDims[2] = {0, 0};
        int err = libData->MTensor_new(MType_Integer, 2, chunkDims, &chunkT);
        if (err) { return err; }
        MArgument_setMTensor(Res, chunkT);
        return LIBRARY_NO_ERROR;
    }

    switch (n) {
        case 0: {
            mint const chunkDims[2] = {1, 0};
            int err = libData->MTensor_new(MType_Integer, 2, chunkDims, &chunkT);
            if (err) {
                libData->MTensor_disown(stateT);
                return err;
            }
            submultiset_deplete(s);
            libData->MTensor_disown(stateT);
            MArgument_setMTensor(Res, chunkT);
            return LIBRARY_NO_ERROR;
        }
        case 1: {
            if (submultiset_first_value(s) == values[0]) {
                mint const chunkDims[2] = {0, 0};
                int err = libData->MTensor_new(MType_Integer, 2, chunkDims, &chunkT);
                if (err) {
                    libData->MTensor_disown(stateT);
                    return err;
                }
            } else {
                mint const chunkDims[2] = {1, 1};
                int err = libData->MTensor_new(MType_Integer, 2, chunkDims, &chunkT);
                if (err) {
                    libData->MTensor_disown(stateT);
                    return err;
                }
                chunk = libData->MTensor_getIntegerData(chunkT);
                chunk[0] = submultiset_first_value(s);
            }
            submultiset_deplete(s);
            libData->MTensor_disown(stateT);
            MArgument_setMTensor(Res, chunkT);
            return LIBRARY_NO_ERROR;
        }
    }

    /*
     * Allocate single chunk of memory for all needed temporary arrays.
     * `n - 1` for `first_rows`, `subsets_number` for `first_rows_cache` and `last_rows_cache`.
     * For automatic `chunkSize` calculation, instead of memory for `last_rows_cache` we need memory for copy of submultiset state.
     */
    size_t pool_size = (sizeof (mint)) * (n + 2 * submultiset_subsets_number(s));
    if (chunkSize < 0 && submultiset_state_size(s) > (sizeof (mint)) * submultiset_subsets_number(s)) {
        pool_size += submultiset_state_size(s) - (sizeof (mint)) * submultiset_subsets_number(s);
    }
    mint * const pool = malloc(pool_size);
    if (!pool) {
        libData->MTensor_disown(stateT);
        return LIBRARY_MEMORY_ERROR;
    }

    if (chunkSize < 0) {
        /* Negative `chunkSize` means that size of chunk should be automatically calculated. */

        /* Array of indices of first rows filled in current iteration for subsequent columns. */
        mint * const first_rows = pool;
        minus_one_fill(first_rows, n);

        /* Array of cached numbers of rows filled for subsequent submultisets. */
        mint * const rows_num_cache = first_rows + n;
        minus_one_fill(rows_num_cache, submultiset_subsets_number(s));

        /* Copy initial generator state, so that it can be restored after calculating chunk size. */
        mint * const state_copy = rows_num_cache + submultiset_subsets_number(s);
        memcpy(state_copy, s->state, submultiset_state_size(s));
        mint * const old_state = s->state;
        s->state = state_copy;

        mint row = 0;

        for (;;) {
            /* Elements same as value of current clumn can't be used in column of derangement. */
            submultiset_skip(s, values[submultiset_removed_count(s)]);

            if (submultiset_remove_current(s)) {
                first_rows[submultiset_removed_count(s)] = row;

                if (rows_num_cache[submultiset_index(s)] != -1) {
                    /* Current multiset was already encountered use cached number of derangements. */
                    if (row > MINT_MAX - rows_num_cache[submultiset_index(s)]) {
                        free(pool);
                        libData->MTensor_disown(stateT);
                        return LIBRARY_MEMORY_ERROR;
                    }
                    row += rows_num_cache[submultiset_index(s)];
                } else {
                    /* Not using cache. */
                    if (submultiset_removed_count(s) == n - 1) {
                        /* One element remains in submultiset. */
                        if (submultiset_first_value(s) != values[n - 1]) {
                            if (row >= MINT_MAX) {
                                free(pool);
                                libData->MTensor_disown(stateT);
                                return LIBRARY_MEMORY_ERROR;
                            }
                            ++row;
                        }
                    } else {
                        /* Proceed to next column. */
                        submultiset_reset_current(s);
                        continue;
                    }
                }
            }

            if (first_rows[submultiset_removed_count(s)] == -1) {
                first_rows[submultiset_removed_count(s)] = 0;
            } else {
                rows_num_cache[submultiset_index(s)] = row - first_rows[submultiset_removed_count(s)];
            }

            if (submultiset_restore(s)) {
                submultiset_advance(s);
            } else {
                break;
            }
        }

        chunkSize = row;

        /* Restore submultiset state. */
        s->state = old_state;
    }

    /* Array of indices of first rows filled in current iteration for subsequent columns. */
    mint * const first_rows = pool;
    minus_one_fill(first_rows, n);

    /* Array of cached indices of first rows filled for subsequent submultisets. */
    mint * const first_rows_cache = first_rows + n;
    minus_one_fill(first_rows_cache, submultiset_subsets_number(s));

    /* Array of cached indices of last rows filled for subsequent submultisets, no need for initialization. */
    mint * const last_rows_cache = first_rows_cache + submultiset_subsets_number(s);

    /* Index of currently filled row in chunk. */
    mint row = 0;

    mint const chunkDims[2] = {chunkSize, n};
    int err = libData->MTensor_new(MType_Integer, 2, chunkDims, &chunkT);
    if (err) {
        free(pool);
        libData->MTensor_disown(stateT);
        return err;
    }
    /* Actual data of chunk tensor. */
    chunk = libData->MTensor_getIntegerData(chunkT);

    for (;;) {
        /* Elements same as value of current clumn can't be used in column of derangement. */
        submultiset_skip(s, values[submultiset_removed_count(s)]);

        if (submultiset_remove_current(s)) {
            first_rows[submultiset_removed_count(s)] = row;

            if (submultiset_removed_count(s) == n - 1) {
                /* One element remains in submultiset. Take care of last element and go back to previous collumn. */
                if (submultiset_first_value(s) != values[n - 1]) {
                    chunk[(row + 1)*n - 1] = submultiset_first_value(s);
                    if (row == chunkSize - 1) {
                        break;
                    }
                    ++row;
                }
            } else {
                /*
                 * If derangements for current multiset were cached and all of them will fit in current chunk copy cached derangements.
                 * If there are more derangements for current multiset than there's space left in chunk ignore cache,
                 * since proper state, from "middle" of cached derangements, for next chunk, must be prepared.
                 */
                mint requiredSize = last_rows_cache[submultiset_index(s)] - first_rows_cache[submultiset_index(s)] + row + 1;
                if (first_rows_cache[submultiset_index(s)] != -1 && requiredSize <= chunkSize) {
                    mint offset = (row - first_rows_cache[submultiset_index(s)]) * n;
                    for (mint i = first_rows_cache[submultiset_index(s)] * n + submultiset_removed_count(s); i <= last_rows_cache[submultiset_index(s)] * n + submultiset_removed_count(s); i += n) {
                        memcpy(chunk + i + offset, chunk + i, (sizeof *chunk) * (n - submultiset_removed_count(s)));
                    }
                    if (requiredSize == chunkSize) {
                        row = requiredSize - 1;
                        break;
                    }
                    row = requiredSize;
                } else {
                    /* Not using cache, proceed to next column. */
                    submultiset_reset_current(s);
                    continue;
                }
            }
        }

        /*
         * If first row for current column is -1, it means that it was not changed after initialization,
         * and current submultiset was not started in this chunk.
         * In that case don't cache it since it'll be incomplete.
         * Otherwise chache index of first and last row of its derangements.
         */
        if (first_rows[submultiset_removed_count(s)] == -1) {
            first_rows[submultiset_removed_count(s)] = 0;
        } else {
            first_rows_cache[submultiset_index(s)] = first_rows[submultiset_removed_count(s)];
            last_rows_cache[submultiset_index(s)] = row - 1;
        }
        if (submultiset_restore(s)) {
            /*
             * Fill current column of chunk, from first to last row of this iteration,
             * with current element of submultiset.
             */
            for (mint i = first_rows[submultiset_removed_count(s) + 1] * n + submultiset_removed_count(s); i < row * n + submultiset_removed_count(s); i += n) {
                chunk[i] = submultiset_current_value(s);
            }

            submultiset_advance(s);
        } else {
            --row;
            break;
        }
    }

    if (submultiset_removed_count(s) > 0) {
        submultiset_restore(s);
        /*
         * Main loop finished, without depleting all submultisets,
         * which means that some collumns might be not fully filled.
         * Fill remaining empty columns in all rows up to last row.
         */
        for (mint col = submultiset_removed_count(s); col >= 0; --col) {
            if (first_rows[col+1] == -1) {
                first_rows[col+1] = 0;
            }
            for (mint i = first_rows[col+1] * n + col; i <= row * n + col; i += n) {
                chunk[i] = submultiset_removed_value(s, col);
            }
        }

        submultiset_advance(s);
    }
    free(pool);
    libData->MTensor_disown(stateT);

    if (chunkSize != row + 1) {
        /* Not all chunk rows were filled, retun shrinked chunk. */
        chunkSize = row + 1;
        MTensor shrinkedChunkT;
        mint const chunkDims[2] = {chunkSize, n};
        int err = libData->MTensor_new(MType_Integer, 2, chunkDims, &shrinkedChunkT);
        if (err) {
            libData->MTensor_free(chunkT);
            return err;
        }
        memcpy(libData->MTensor_getIntegerData(shrinkedChunkT), chunk, (sizeof *chunk) * chunkSize * n);
        libData->MTensor_free(chunkT);
        MArgument_setMTensor(Res, shrinkedChunkT);
        return LIBRARY_NO_ERROR;
    }

    MArgument_setMTensor(Res, chunkT);
    return LIBRARY_NO_ERROR;
}

Save above code in derangements.c file in same directory as current notebook, or paste it as a string, instead of {"derangements.c"}, as first argument of CreateLibrary in code below. Pass, in "CompileOptions", appropriate flags for your compiler, the ones below are for GCC. Above code uses some C99 features, so make sure to compile with C99 support.

Needs@"CCompilerDriver`"
SetDirectory@NotebookDirectory[];
lib = CreateLibrary[{"derangements.c"}, "derangements"(*,
  "CompileOptions" -> "-std=c99 -Wall -march=native -O3"*)
]
nextDerangementsChunk = LibraryFunctionLoad[lib, "nextDerangementsChunk",
  {{Integer, 1, "Constant"}, {Integer, 1, "Constant"}, {Integer, 1, "Constant"}, {Integer, 1, "Shared"}, Integer},
  {Integer, 2}
]

initialState = Function[{n, mult}, With[{m = Max[(n - 1), 0]}, Join[
  ConstantArray[0, 3],
  Range[Length@mult - 1],
  {-1},
  mult,
  ConstantArray[-1, m],
  ConstantArray[0, m]
]]];

nextDerangementsChunk accepts five arguments.

  1. List of "reference values" containing elements that can't be used in subsequent places of derangement. For ordinary derangements this should be original list of values for which derangements are calculated, for partial derangements - elements for which fixed points are allowed should be replaced with element not present in values.

  2. List of unique values of multiset which derangements we're generating.

  3. List of "index multipliers", used to calculate index of submultiset, obtained by multiplying subsequent multiplicities of values in multiset.

  4. List representing "state" of generator. Initial state can be generated using initialState[n, mult] where n is number of elements in multiset which derangements we're calculating, and mult is list of multiplicities of subsequent values in this multiset. -1 as first element of state means that generator is depleted.

  5. Maximal number of derangements in returned chunk. -1 means all remaining derangements.

Basic examples

As a usage example let's print partial derangements of {2, 1, 4, 1, 3}, with allowed fixed points in 1-st and 3-rd position, in chunks of at most 5 elements:

values = {2, 1, 4, 1, 3};
{unique, mult} = Transpose@Tally@values;
state = initialState[Length@values, mult];
indexMultipliers = FoldList[Times, 1, mult + 1];
refValues = ReplacePart[values, {{1}, {3}} -> Min@values - 1];
While[First@state > -1,
  Print@nextDerangementsChunk[refValues, unique, indexMultipliers, state, 5]
]
(* {{2,4,1,3,1},{2,3,1,4,1},{1,2,1,3,4},{1,2,4,3,1},{1,2,3,4,1}} *)
(* {{1,4,2,3,1},{1,4,1,3,2},{1,4,3,2,1},{1,3,2,4,1},{1,3,1,2,4}} *)
(* {{1,3,1,4,2},{1,3,4,2,1},{4,2,1,3,1},{4,3,1,2,1},{3,2,1,4,1}} *)
(* {{3,4,1,2,1}} *)

All (partial) derangements

derangements // ClearAll
derangements[empty:_[], Repeated[_, {0, 1}]] := {empty}
derangements[expr_ /; Not@AtomQ@expr, Optional[pos_List /; VectorQ[pos, IntegerQ], {}]] :=
  Module[{packed, unique, mult, list, uniqueVal, posInd, postprocess},
    packed = Developer`ToPackedArray@expr;
    If[Developer`PackedArrayQ@packed && VectorQ[packed, IntegerQ],
      {unique, mult} = Transpose@Tally@packed;
      postprocess = Identity;
    (* else *),
      list = List @@ expr;
      {uniqueVal, mult} = Transpose@Tally@list;
      posInd = PositionIndex@uniqueVal;
      packed = Lookup[posInd, list][[All, 1]] // Developer`ToPackedArray;
      unique = Values[posInd][[All, 1]];
      uniqueVal = Head@expr @@ uniqueVal;
      postprocess = Map[uniqueVal[[#]]&]
    ];
    nextDerangementsChunk[
      ReplacePart[packed, Transpose@{pos} -> (Min@packed - 1)],
      unique,
      FoldList[Times, 1, mult + 1],
      initialState[Length@packed, mult],
      -1
    ] // postprocess
  ]

Basic examples

Partial derangements of f[a, c, b, a] expression, with allowed fixed points at first and penultimate positions:

derangements[f[a, c, b, a], {1, -2}]
(* {f[a, a, c, b], f[a, a, b, c], f[a, b, a, c], f[c, a, a, b], f[b, a, a, c]} *)

Benchmarks

We'll compare above function with Pick-based solution and derangementsPacked.

List of unique integers, from OP:

s = Range@9;
(res1 = Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1]&@Permutations[s]) // MaxMemoryUsed // RepeatedTiming
(res2 = derangementsPacked@s) // MaxMemoryUsed // RepeatedTiming
(res3 = derangements@s) // MaxMemoryUsed // RepeatedTiming
res1 === res2 === res3
(* {0.050, 78384440} *)
(* {0.031, 34618920} *)
(* {0.0031, 9618368} *)
(* True *)

while library link filtering gives {0.0043, 9613720}.

As it's a non-filtering approach its advantage grows with number of duplicates:

s = Join[ConstantArray[1, 7], ConstantArray[2, 5], Range[3, 5]];
(res1 = Pick[#, Unitize[Times @@ (#\[Transpose] - s)], 1]&@Permutations[s]) // MaxMemoryUsed // RepeatedTiming
(res2 = derangementsPacked@s) // MaxMemoryUsed // RepeatedTiming
(res3 = derangements@s) // MaxMemoryUsed // RepeatedTiming
res1 === res2 === res3
(* {0.5481, 778380768} *)
(* {0.0018,    824784} *)
(* {0.0000702, 187920} *)
(* True *)

library link filtering gives {0.016, 182984}.

$\endgroup$
  • $\begingroup$ Thank you yet again, and I am glad that you (apparently) found my question so inspirational! $\endgroup$ – Mr.Wizard Jun 5 '17 at 0:06
  • $\begingroup$ Is there a case where this method is inferior to your earlier answer? $\endgroup$ – Mr.Wizard Jun 5 '17 at 0:10
  • 1
    $\begingroup$ @Mr.Wizard Only inferiority that comes to my mind is that this version requires additional memory for caching. It seems that memory allocated and freed inside LibraryFunction is not measured by MaxMemoryUsed, so it's not easy to detect from inside Mathematica. This memory is independent of requested chunk size and equal to 2 times size of machine integer times last element of indexMultipliers. If you're generating all derangements, then cache size is negligible compared to size of output, but if you're generating small chunks for some larger set, then cache size might be dominating. $\endgroup$ – jkuczm Jun 5 '17 at 19:44
  • 1
    $\begingroup$ @Mr.Wizard e.g. for Range@27 even if you're generating chunk containing only one derangement, you need something of order of 2GB for caching. This could be changed by adding additional argument allowing manual switch-off of caching, or by switching it off automatically if there's not enough memory, but this is not implemented currently. $\endgroup$ – jkuczm Jun 5 '17 at 19:46
1
$\begingroup$

Very slow, but I got curious,

derangements[n_Integer] := Module[{list = Range[n], f, g},
  f[{k_, s_}] := 
   Map[{k + 1, Append[s, {#}]} &, Delete[list, Append[s, {k + 1}]]];
  g[ql_] := Flatten[Map[f, ql], 1];
  Flatten /@ Nest[g, {{0, {}}}, Length[list]][[All, 2]]
  ]

derangements[9] // Length // RepeatedTiming
(* {2.60, 133496} *)
$\endgroup$

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.