nthPermutation (#5068)

add the function nthPermutation that permutates a given range in place n
permutations in O(1).

This is O(n - 1) steps faster than calling nextPermutation n times.

std/algorithm/sorting.d

@@ -35,6 +35,9 @@ $(T2 nextEvenPermutation,
 $(T2 nextPermutation,
         Computes the next lexicographically greater permutation of a range
         in-place.)
+$(T2 nthPermutation,
+        Computes the nth permutation of a range
+        in-place.)
 $(T2 partialSort,
         If `a = [5, 4, 3, 2, 1]`, then `partialSort(a, 3)` leaves
         `a[0 .. 3] = [1, 2, 3]`.
@@ -4466,3 +4469,175 @@ shapes. Here's a non-trivial example:
     }
     assert(n == 60);
 }
+
+/** Permutes `range` into the `perm` permutation.
+The algorithm has a constant runtime complexity with respect to the number of
+permutations created.
+Due to the number of unique values of `ulong` only the first 21 elements of
+`range` can be permuted. The rest of the range will therefore not be
+permuted.
+This algorithm uses the $(HTTP en.wikipedia.org/wiki/Lehmer_code, Lehmer
+Code).
+
+The algorithm works as follows:
+$(D_CODE
+    auto pem = [4,0,4,1,0,0,0]; // permutation 2982 in factorial
+    auto src = [0,1,2,3,4,5,6]; // the range to permutate
+
+    auto i = 0;                    // range index
+    // range index iterates pem and src in sync
+    // pem[i] + i is used as index into src
+    // first src[pem[i] + i] is stored in t
+    auto t = 4;                    // tmp value
+    src = [0,1,2,3,n,5,6];
+
+    // then the values between i and pem[i] + i are moved one
+    // to the right
+    src = [n,0,1,2,3,5,6];
+    // at last t is inserted into position i
+    src = [4,0,1,2,3,5,6];
+    // finally i is incremented
+    ++i;
+
+    // this process is repeated while i < pem.length
+
+    t = 0;
+    src = [4,n,1,2,3,5,6];
+    src = [4,0,1,2,3,5,6];
+    ++i;
+    t = 6;
+    src = [4,0,1,2,3,5,n];
+    src = [4,0,n,1,2,3,5];
+    src = [4,0,6,1,2,3,5];
+)
+
+Returns:
+    The permuted range.
+
+Params:
+    range = The Range to permute. The original ordering will be lost.
+    perm = The permutation to permutate `range` to.
+*/
+auto ref Range nthPermutation(Range)
+                             (auto ref Range range, const ulong perm)
+if (isRandomAccessRange!Range && hasLength!Range)
+{
+    if (!nthPermutationImpl(range, perm))
+    {
+        throw new Exception(
+            "The range to permutate must not have less"
+            ~ " elements than the factorial number has digits");
+    }
+
+    return range;
+}
+
+///
+pure @safe unittest
+{
+    auto src = [0, 1, 2, 3, 4, 5, 6];
+    auto rslt = [4, 0, 6, 2, 1, 3, 5];
+
+    src = nthPermutation(src, 2982);
+    assert(src == rslt);
+}
+
+/**
+Returns: `true` in case the permutation worked, `false` in case `perm` had
+    more digits in the factorial number system than range had elements.
+    This case must not occur as this would lead to out of range accesses.
+*/
+bool nthPermutationImpl(Range)
+                       (auto ref Range range, ulong perm)
+if (isRandomAccessRange!Range && hasLength!Range)
+{
+    import std.range.primitives : ElementType;
+    import std.numeric : decimalToFactorial;
+
+    // ulong.max has 21 digits in the factorial number system
+    ubyte[21] fac;
+    size_t idx = decimalToFactorial(perm, fac);
+
+    if (idx > range.length)
+    {
+        return false;
+    }
+
+    ElementType!Range tmp;
+    size_t i = 0;
+
+    for (; i < idx; ++i)
+    {
+        size_t re = fac[i];
+        tmp = range[re + i];
+        for (size_t j = re + i; j > i; --j)
+        {
+            range[j] = range[j - 1];
+        }
+        range[i] = tmp;
+    }
+
+    return true;
+}
+
+///
+pure @safe unittest
+{
+    auto src = [0, 1, 2, 3, 4, 5, 6];
+    auto rslt = [4, 0, 6, 2, 1, 3, 5];
+
+    bool worked = nthPermutationImpl(src, 2982);
+    assert(worked);
+    assert(src == rslt);
+}
+
+pure @safe unittest
+{
+    auto rslt = [4, 0, 6, 2, 1, 3, 5];
+
+    auto src = nthPermutation([0, 1, 2, 3, 4, 5, 6], 2982);
+    assert(src == rslt);
+}
+
+pure @safe unittest
+{
+    auto src = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
+    auto rslt = [4, 0, 6, 2, 1, 3, 5, 7, 8, 9, 10];
+
+    src = nthPermutation(src, 2982);
+    assert(src == rslt);
+}
+
+pure @safe unittest
+{
+    import std.exception : assertThrown;
+
+    auto src = [0, 1, 2, 3];
+
+    assertThrown(nthPermutation(src, 2982));
+}
+
+pure @safe unittest
+{
+    import std.internal.test.dummyrange;
+    import std.meta : AliasSeq;
+
+    auto src = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
+    auto rsl = [4, 0, 6, 2, 1, 3, 5, 7, 8, 9, 10];
+
+    foreach (T; AliasSeq!(
+            DummyRange!(ReturnBy.Reference, Length.Yes, RangeType.Random, int[]),
+            DummyRange!(ReturnBy.Value, Length.Yes, RangeType.Random, int[])))
+    {
+        static assert(isRandomAccessRange!(T));
+        static assert(hasLength!(T));
+        auto dr = T(src.dup);
+        dr = nthPermutation(dr, 2982);
+
+        int idx;
+        foreach (it; dr)
+        {
+            assert(it == rsl[idx++]);
+        }
+    }
+}

std/numeric.d

@@ -3335,6 +3335,89 @@ C swapRealImag(C)(C input)
     return C(input.im, input.re);
 }
 
+/** This function transforms `decimal` value into a value in the factorial number
+system stored in `fac`.
+
+A factorial number is constructed as:
+$(D fac[0] * 0! + fac[1] * 1! + ... fac[20] * 20!)
+
+Params:
+    decimal = The decimal value to convert into the factorial number system.
+    fac = The array to store the factorial number. The array is of size 21 as
+        `ulong.max` requires 21 digits in the factorial number system.
+Returns:
+    A variable storing the number of digits of the factorial number stored in
+    `fac`.
+*/
+size_t decimalToFactorial(ulong decimal, ref ubyte[21] fac)
+        @safe pure nothrow @nogc
+{
+    import std.algorithm.mutation : reverse;
+    size_t idx;
+
+    for (ulong i = 1; decimal != 0; ++i)
+    {
+        auto temp = decimal % i;
+        decimal /= i;
+        fac[idx++] = cast(ubyte)(temp);
+    }
+
+    if (idx == 0)
+    {
+        fac[idx++] = cast(ubyte) 0;
+    }
+
+    reverse(fac[0 .. idx]);
+
+    // first digit of the number in factorial will always be zero
+    assert(fac[idx - 1] == 0);
+
+    return idx;
+}
+
+///
+@safe pure @nogc unittest
+{
+    ubyte[21] fac;
+    size_t idx = decimalToFactorial(2982, fac);
+
+    assert(fac[0] == 4);
+    assert(fac[1] == 0);
+    assert(fac[2] == 4);
+    assert(fac[3] == 1);
+    assert(fac[4] == 0);
+    assert(fac[5] == 0);
+    assert(fac[6] == 0);
+}
+
+@safe pure unittest
+{
+    ubyte[21] fac;
+    size_t idx = decimalToFactorial(0UL, fac);
+    assert(idx == 1);
+    assert(fac[0] == 0);
+
+    fac[] = 0;
+    idx = 0;
+    idx = decimalToFactorial(ulong.max, fac);
+    assert(idx == 21);
+    auto t = [7, 11, 12, 4, 3, 15, 3, 5, 3, 5, 0, 8, 3, 5, 0, 0, 0, 2, 1, 1, 0];
+    foreach (i, it; fac[0 .. 21])
+    {
+        assert(it == t[i]);
+    }
+
+    fac[] = 0;
+    idx = decimalToFactorial(2982, fac);
+
+    assert(idx == 7);
+    t = [4, 0, 4, 1, 0, 0, 0];
+    foreach (i, it; fac[0 .. idx])
+    {
+        assert(it == t[i]);
+    }
+}
+
 private:
 // The reasons I couldn't use std.algorithm were b/c its stride length isn't
 // modifiable on the fly and because range has grown some performance hacks