/* 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]]
]
]]]