Skip to main content
replaced http://mathematica.stackexchange.com/ with https://mathematica.stackexchange.com/
Source Link

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

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 postmy permutations post.

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

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.

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

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.

Fix `derangements` for empty expressions.
Source Link
jkuczm
  • 15.1k
  • 2
  • 54
  • 85
derangements // ClearAll
derangements[_[]]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]]
        ]
    ]]]
derangements // ClearAll
derangements[_[]] = 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]]
        ]
    ]]]
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]]
        ]
    ]]]
Add example of generation in chunks, and few words about number of derangements.
Source Link
jkuczm
  • 15.1k
  • 2
  • 54
  • 85
/* 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];
    
    /* ValuesOrdered values tensor: multisetlist of integers in non-increasing order. */
    MTensor orderedValuesT = MArgument_getMTensor(Args[1]);
    /* Actual data of orderedValuesordered 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;
}
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]
        ]
    ]
]&;
multiSubfact = If[MemberQ[#, Except@1],
    With[
        {coeffs = Block[{x},
            CoefficientList[Times @@ (LaguerreL[#, x] & /@ #), x]
        ]},
        Abs@Total[Factorial@Range[0, Length@coeffs - 1] coeffs]
    ]
(* else *),
    Subfactorial@Length@#
]&;

derangements // ClearAll
derangements[_[]] = derangements[_[_]] = {};
derangements[list_List /; VectorQ[Unevaluated@list, IntegerQ]] :=
    With[{n = Length@list},
        nextDerangementsChunk[
            list,
            Reverse@Sort@list,
            Join[Range@n, {0, n - 2}],
            multiSubfact@Tally[list][[AllmultiSubfactorial@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}],
            multiSubfact@tallied[[AllmultiSubfactorial@tallied[[All, 2]]
        ]
    ]]]
/* 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];
    
    /* Values tensor: multiset of integers in non-increasing order. */
    MTensor orderedValuesT = MArgument_getMTensor(Args[1]);
    /* Actual data of orderedValues 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]) {
                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;
}
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}
]
multiSubfact = If[MemberQ[#, Except@1],
    With[
        {coeffs = Block[{x},
            CoefficientList[Times @@ (LaguerreL[#, x] & /@ #), x]
        ]},
        Abs@Total[Factorial@Range[0, Length@coeffs - 1] coeffs]
    ]
(* else *),
    Subfactorial@Length@#
]&;

derangements // ClearAll
derangements[_[]] = derangements[_[_]] = {};
derangements[list_List /; VectorQ[Unevaluated@list, IntegerQ]] :=
    With[{n = Length@list},
        nextDerangementsChunk[
            list,
            Reverse@Sort@list,
            Join[Range@n, {0, n - 2}],
            multiSubfact@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}],
            multiSubfact@tallied[[All, 2]]
        ]
    ]]]
/* 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;
}
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]
        ]
    ]
]&;
derangements // ClearAll
derangements[_[]] = 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]]
        ]
    ]]]
Source Link
jkuczm
  • 15.1k
  • 2
  • 54
  • 85
Loading