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.
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.
List of unique values of multiset which derangements we're generating.
List of "index multipliers", used to calculate index of submultiset, obtained by multiplying subsequent multiplicities of values in multiset.
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.
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}.
