next_permutation and cartesian product for ranges?

[bearophile]

Andrei Alexandrescu:

> Yes, we need next_permutation. It will be up to you to convince 
> the Grand Inquisition Committee of Reviewers that 
> next_even_permutation is necessary and/or sufficient.

I don't like the design of C++ STL here. Instead of a 
next_permutation(), what about ranges that yield all 
permutations, all adjacent permutations, all combinations?


Permutations:

import std.algorithm, std.conv, std.traits;

struct Permutations(bool doCopy=true, T) if (isMutable!T) {
     private immutable size_t num;
     private T[] items;
     private uint[31] indexes;
     private ulong tot;

     this (/*in*/ T[] items) /*pure*/ nothrow
     in {
         static enum string L = text(indexes.length); // impure
         assert(items.length >= 0 && items.length <= 
indexes.length,
                "Permutations: items.length must be >= 0 && < " ~ 
L);
     } body {
         static ulong factorial(in uint n) pure nothrow {
             ulong result = 1;
             foreach (i; 2 .. n + 1)
                 result *= i;
             return result;
         }

         this.num = items.length;
         this.items = items.dup;
         foreach (i; 0 .. cast(typeof(indexes[0]))this.num)
             this.indexes[i] = i;
         this.tot = factorial(this.num);
     }

     @property T[] front() pure nothrow {
         static if (doCopy)
             //return items.dup; // not nothrow
             return items ~ []; // slower
         else
             return items;
     }

     @property bool empty() const pure nothrow {
         return tot == 0;
     }

     void popFront() pure nothrow {
         tot--;
         if (tot > 0) {
             size_t j = num - 2;

             while (indexes[j] > indexes[j + 1])
                 j--;
             size_t k = num - 1;
             while (indexes[j] > indexes[k])
                 k--;
             swap(indexes[k], indexes[j]);
             swap(items[k], items[j]);

             size_t r = num - 1;
             size_t s = j + 1;
             while (r > s) {
                 swap(indexes[s], indexes[r]);
                 swap(items[s], items[r]);
                 r--;
                 s++;
             }
         }
     }
}

Permutations!(doCopy,T) permutations(bool doCopy=true, T)
                                     (/*in*/ T[] items)
/*pure*/ nothrow if (isMutable!T) {
     return Permutations!(doCopy, T)(items);
} unittest {
     import std.bigint;
     foreach (p; permutations([BigInt(1), BigInt(2), BigInt(3)]))
         assert((p[0] + 1) > 0);
}

version (permutations2_main) {
     void main() {
         import std.stdio, std.bigint;
         foreach (p; permutations!false([1, 2, 3]))
             writeln(p);
         alias BigInt B;
         foreach (p; permutations!false([B(1), B(2), B(3)])) {}
     }
}

----------------------------------

Adjacent permutations by swapping (to be converted in a range):

import std.algorithm, std.array, std.typecons, std.range;

struct Spermutations {
     const int n;
     alias Tuple!(int[], int) TResult;

     int opApply(int delegate(ref TResult) dg) {
         int result;
         TResult aux;

         int sign = 1;
         alias Tuple!(int, int) Pair;
         auto p = iota(n).map!(i => Pair(i, i ? -1 : 0))().array();

         aux = tuple(p.map!(pp => pp[0])().array(), sign);
         result = dg(aux); if (result) goto END;

         while (p.canFind!(pp => pp[1])()) {
             // Failed using std.algorithm here, too much complex
             auto largest = Pair(-100, -100);
             int i1 = -1;
             foreach (i, pp; p)
                 if (pp[1]) {
                     if (pp[0] > largest[0]) {
                         i1 = i;
                         largest = pp;
                     }
                 }
             int n1 = largest[0], d1 = largest[1];

             sign *= -1;
             int i2;
             if (d1 == -1) {
                 i2 = i1 - 1;
                 swap(p[i1], p[i2]);
                 if (i2 == 0 || p[i2 - 1][0] > n1)
                     p[i2][1] = 0;
             } else if (d1 == 1) {
                 i2 = i1 + 1;
                 swap(p[i1], p[i2]);
                 if (i2 == n - 1 || p[i2 + 1][0] > n1)
                     p[i2][1] = 0;
             }
             aux = tuple(p.map!(pp => pp[0])().array(), sign);
             result = dg(aux); if (result) goto END;

             foreach (i3, ref pp; p) {
                 auto n3 = pp[0], d3 = pp[1];
                 if (n3 > n1)
                     pp[1] = (i3 < i2) ? 1 : -1;
             }
         }

         END: return result;
     }
}

void main() {
     import std.stdio;
     foreach (n; [3, 4]) {
         writefln("\nPermutations and sign of %d items", n);
         foreach (tp; Spermutations(n))
             writefln("Perm: %s Sign: %2d", tp.tupleof);
     }
}

----------------------------------

Combinations (to be converted in a range):

ulong binomial(long n, long k) pure nothrow
in {
     assert(n > 0, "binomial: n must be > 0.");
} body {
     if (k < 0 || k > n)
         return 0;
     if (k > (n / 2))
         k = n - k;
     ulong result = 1;
     foreach (ulong d; 1 .. k + 1) {
         result *= n;
         n--;
         result /= d;
     }
     return result;
}


struct Combinations(T, bool copy=true) {
     // Algorithm by Knuth, Pre-fascicle 3A, draft of
     // section 7.2.1.3: "Generating all partitions".
     T[] items;
     int k;
     size_t len = -1; // computed lazily

     this(in T[] items, in int k)
     in {
         assert(items.length, "combinations: items can't be 
empty.");
     } body {
         this.items = items.dup;
         this.k = k;
     }

     @property size_t length() /*logic_const*/ {
         if (len == -1) // set cache
             len = cast(size_t)binomial(items.length, k);
         return len;
     }

     int opApply(int delegate(ref T[]) dg) {
         if (k == items.length)
             return dg(items);         // yield items

         auto outarr = new T[k];
         if (k == 0)
             return dg(outarr);        // yield []

         if (k < 0 || k > items.length)
             return 0;                 // yield nothing

         int result, x;
         immutable n = items.length;
         auto c = new uint[k + 3]; // c[0] isn'k used

         foreach (j; 1 .. k + 1)
             c[j] = j - 1;
         c[k + 1] = n;
         c[k + 2] = 0;
         int j = k;

         while (true) {
             // The following lines equal to:
             //int pos;
             //foreach (i; 1 .. k +1)
             //    outarr[pos++] = items[c[i]];
             auto outarr_ptr = outarr.ptr;
             auto c_ptr = &(c[1]);
             auto c_ptrkp1 = &(c[k + 1]);
             while (c_ptr != c_ptrkp1)
                 *outarr_ptr++ = items[*c_ptr++];


             static if (copy) {
                 auto outarr2 = outarr.dup;
                 result = dg(outarr2); // yield outarr2
             } else {
                 result = dg(outarr); // yield outarr
             }

             if (j > 0) {
                 x = j;
                 c[j] = x;
                 j--;
                 continue;
             }

             if ((c[1] + 1) < c[2]) {
                 c[1]++;
                 continue;
             } else
                 j = 2;

             while (true) {
                 c[j - 1] = j - 2;
                 x = c[j] + 1;
                 if (x == c[j + 1])
                     j++;
                 else
                     break;
             }

             if (j > k)
                 return result; // End

             c[j] = x;
             j--;
         }
     }
}

Combinations!(T,copy) combinations(bool copy=true, T)
                                   (in T[] items, in int k)
in {
     assert(items.length, "combinations: items can't be empty.");
} body {
     return Combinations!(T, copy)(items, k);
}

----------------------------------

Bye,
bearophile