Boost :

hinnant_at_[hidden]

---
Interface
---
I've followed the advice of John Dibling in the CUJ article except that 
the functor doesn't receive the entire range, but only the first K 
elements.  For example consider:
char array[] = {'a', 'b', 'c', 'd', 'e'};
unsigned size_array = sizeof(array)/sizeof(array[0]);
const int k = 3;
print p = std::for_each_combination(array, array + size_array, k, 
print());
This calls the print() functor (say p) like:
     p(array, array+k);
for array permuted such that [array, array+k) will hold each 
combination of 3 items taken from the total of 5 items.  And after the 
for_each_combination call, array will be in its original order.  While 
always permuting the elements into positions [0, k) may at first seem 
limiting in how the for_each_* algorithms can be implemented, it offers 
interesting possibilities for the functors.
If a functor knows the total length of the sequence (say N), and 
receives [array, array+k) as a subsequence representing the current 
combination/permutation, that functor can then call for_each_* on the 
remaining sequence [array+k, array+N).  This allows for a very complex 
analysis where one considers all combinations (or permutations) of the 
first k elements, and then simultaneously considers permutations of the 
remaining elements for each of the outer combinations.
One could use this capability to search for solutions to variations of 
the traveling salesman problem.  Instead of minimizing the route of 1 
salesman over N cities.  You might want to minimize two salesmen's 
routes over the same set of cities, each city being visited by one and 
only one salesman.  For example you might consider:
For each circular-reversible permutation of 10 cities considered 5 at a 
time:
     Compute distance around the 5 city route and add to that:
     For each circular-reversible permutation of the remaining 5 cities 
considered 5 at a time:
        Compute distance around the 5 city route.
        Record minimum combined distance and routes.
Naturally this technique can easily be generalized to M salesmen 
covering N cites, N/M cites covered by each salesman.  Unfortunately I 
currently lack the circular-reversible permutation algorithm and have 
to instead use the twice-as-expensive circular permutation algorithm.
The three algorithms I do have are here:
#include <iterator>
#include <algorithm>
namespace test
{
namespace detail
{
template<class BidirectionalIterator, class Function, class Size>
Function
permute(BidirectionalIterator first, BidirectionalIterator mid,
           BidirectionalIterator active, BidirectionalIterator last,
           Size k, Function f, Size n)
{
     if (k == 0)
         f(first, mid);
     else
     {
         using std::swap;
         BidirectionalIterator activep1 = active;
         ++activep1;
         BidirectionalIterator p = active;
         ++p;
         for (Size i = 0; i < n; ++i, ++p)
         {
             if (k == 1)
                 f(first, mid);
             else
                 f = permute<BidirectionalIterator,Function,Size>(first, 
mid, activep1, last, k-1, f , n-1);
             if (i != n-1)
                 swap(*active, *p);
         }
         std::rotate(active, activep1, last);
     }
     return f;
}
template<class BidirectionalIterator, class Function, class Size>
Function
combine(BidirectionalIterator first, BidirectionalIterator mid,
            BidirectionalIterator active, BidirectionalIterator last,
            Size k, Function f, Size n)
{
     if (k == 0)
         f(first, mid);
     else
     {
         Size e = n - k + 1;
         Size nn = n-1;
         BidirectionalIterator l2 = last;
         BidirectionalIterator ap1 = active;
         ++ap1;
         BidirectionalIterator apk = active;
         std::advance(apk, k);
         for (Size i = 0; i < e; ++i, --nn, --l2)
         {
             if (k == 1)
                 f(first, mid);
             else
                 f = combine<BidirectionalIterator,Function,Size>(first, 
mid, ap1, l2, k-1, f, nn);
             std::rotate(active, i < e-1 ? ap1 : apk, last);
         }
     }
     return f;
}
}  // detail
template<class BidirectionalIterator, class Function, class Size>
Function
for_each_permutation(BidirectionalIterator first, BidirectionalIterator 
last, Size k, Function f)
{
     typedef typename 
std::iterator_traits<BidirectionalIterator>::difference_type 
difference_type;
     Size n = static_cast<Size>(std::distance(first,last));
     if (k < 1 || n == 0)
         return f;
     if (k > n)
         k = n;
     BidirectionalIterator mid = first;
     std::advance(mid, k);
     return detail::permute<BidirectionalIterator,Function,Size>(first, 
mid, first, last, k, f, n);
}
template<class BidirectionalIterator, class Function, class Size>
Function
for_each_circular_permutation(BidirectionalIterator first, 
BidirectionalIterator last, Size k, Function f)
{
     typedef typename 
std::iterator_traits<BidirectionalIterator>::difference_type 
difference_type;
     Size n = static_cast<Size>(std::distance(first,last));
     if (k < 1 || n == 0)
         return f;
     if (k > n)
         k = n;
     BidirectionalIterator mid = first;
     std::advance(mid, k);
     BidirectionalIterator f2 = first;
     ++f2;
     for (Size i = 0, e = n - k + 1; i < e; ++i, --n, ++first, ++f2, 
++mid)
         f = detail::permute<BidirectionalIterator,Function,Size>(first, 
mid, f2, last, k-1, f, n-1);
     return f;
}
template<class BidirectionalIterator, class Function, class Size>
Function
for_each_combination(BidirectionalIterator first, BidirectionalIterator 
last, Size k, Function f)
{
     typedef typename 
std::iterator_traits<BidirectionalIterator>::difference_type 
difference_type;
     Size n = static_cast<Size>(std::distance(first,last));
     if (k < 1)
         return f;
     if (k > n)
         k = n;
     BidirectionalIterator mid = first;
     std::advance(mid, k);
     return detail::combine<BidirectionalIterator,Function,Size>(first, 
mid, first, last, k, f, n);
}
}  // test
And here's a demo that exercises all 3 algorithms:
// demo
#include <iostream>
#include <iterator>
#include <algorithm>
class print
{
public:
     print() : i_(0) {}
     unsigned get() const {return i_;}
     template <class It>
     void operator()(It first, It last)
     {
         typedef typename std::iterator_traits<It>::value_type 
value_type;
         std::ostream_iterator<value_type> out(std::cout, " ");
         std::cout << ++i_ << "\t: ";
         std::copy(first, last, out);
         std::cout << '\n';
     }
private:
     unsigned i_;
};
int main()
{
     char array[] = {'a', 'b', 'c', 'd', 'e'};
     unsigned size_array = sizeof(array)/sizeof(array[0]);
     const int k = 3;
     std::cout << "Consider " << size_array << " items " << k << " at a 
time from this sequence\n";
     print()(array, array + size_array);
     {
     std::cout << "Permutations\n";
     print p = test::for_each_permutation(array, array + size_array, k, 
print());
     std::cout << "--- Total of " << p.get() << " calls ---\n";
     std::cout << "Sequence left in orginal order\n";
     print()(array, array + size_array);
     std::cout << '\n';
     }
     {
     std::cout << "Circular Permutations\n";
     print p = test::for_each_circular_permutation(array, array + 
size_array, k, print());
     std::cout << "--- Total of " << p.get() << " calls ---\n";
     std::cout << "Sequence left in orginal order\n";
     print()(array, array + size_array);
     std::cout << '\n';
     }
     {
     std::cout << "Combinations\n";
     print p = test::for_each_combination(array, array + size_array, k, 
print());
     std::cout << "--- Total of " << p.get() << " calls ---\n";
     std::cout << "Sequence left in orginal order\n";
     print()(array, array + size_array);
     std::cout << '\n';
     }
}
The output should be:
Consider 5 items 3 at a time from this sequence
1   : a b c d e
Permutations
1   : a b c
2   : a b d
3   : a b e
4   : a c b
5   : a c d
6   : a c e
7   : a d b
8   : a d c
9   : a d e
10  : a e b
11  : a e c
12  : a e d
13  : b a c
14  : b a d
15  : b a e
16  : b c a
17  : b c d
18  : b c e
19  : b d a
20  : b d c
21  : b d e
22  : b e a
23  : b e c
24  : b e d
25  : c a b
26  : c a d
27  : c a e
28  : c b a
29  : c b d
30  : c b e
31  : c d a
32  : c d b
33  : c d e
34  : c e a
35  : c e b
36  : c e d
37  : d a b
38  : d a c
39  : d a e
40  : d b a
41  : d b c
42  : d b e
43  : d c a
44  : d c b
45  : d c e
46  : d e a
47  : d e b
48  : d e c
49  : e a b
50  : e a c
51  : e a d
52  : e b a
53  : e b c
54  : e b d
55  : e c a
56  : e c b
57  : e c d
58  : e d a
59  : e d b
60  : e d c
--- Total of 60 calls ---
Sequence left in orginal order
1   : a b c d e
Circular Permutations
1   : a b c
2   : a b d
3   : a b e
4   : a c b
5   : a c d
6   : a c e
7   : a d b
8   : a d c
9   : a d e
10  : a e b
11  : a e c
12  : a e d
13  : b c d
14  : b c e
15  : b d c
16  : b d e
17  : b e c
18  : b e d
19  : c d e
20  : c e d
--- Total of 20 calls ---
Sequence left in orginal order
1   : a b c d e
Combinations
1   : a b c
2   : a b d
3   : a b e
4   : a c d
5   : a c e
6   : a d e
7   : b c d
8   : b c e
9   : b d e
10  : c d e
--- Total of 10 calls ---
Sequence left in orginal order
1   : a b c d e
Note that the functor (print() in this case) can be stateful, and that 
state is accumulated from call to call, and returned back to the 
client.  In this case the state is a simple line number counter.
Anyway, I offer these up on the chance they might be useful, and on the 
chance that someone might help complete the library to compute 
"reversible" permutations, which may indeed be quite useful (as in the 
traveling salesman problem, and variations of it).
-Howard