All Permutations - non-recursive

The technique shown by @spacerat is the factorial-number (factoradic) mapping: an index is decomposed into mixed-radix digits that tell you which remaining element to pick at each position. That is why you can jump straight to permutation number a without generating earlier permutations. The tradeoff is cost per query — naive removal with erase is O(n) so the full generation is O(n^2) per permutation.

Practical improvements: use a clear name (avoid default), accept the input as std::vector<T> (as suggested by @CppBuilder2006), precompute factorials once and range-check the requested index (or throw). Beware factorial overflow: 64-bit indexes cover up to 20! — beyond that you need big integers (Boost.Multiprecision, etc.). Also heed @mrnutty: prefer well-tested components for production; this technique is useful when you must index permutations directly, not to wholesale-replace std::next_permutation.

A faster selection step avoids repeated O(n) removals. Use a Fenwick (binary indexed) tree to represent which slots are still available and find the k-th free slot in O(log n), yielding O(n log n) total per permutation. Example (concise):

#include <vector>
#include <cstdint>
#include <limits>
#include <stdexcept>

template<typename Int = int>
struct Fenwick {
    int n; std::vector<Int> bit;
    Fenwick(int n_=0): n(n_), bit(n_+1,0) {}
    void add(int idx, Int val){ for(++idx; idx<=n; idx += idx & -idx) bit[idx] += val; }
    int find_kth(Int k){ // 0-based k-th one
        int idx = 0, mask = 1;
        while ((mask << 1) <= n) mask <<= 1;
        for (; mask; mask >>= 1){
            int nxt = idx + mask;
            if (nxt <= n && bit[nxt] <= k){ k -= bit[nxt]; idx = nxt; }
        }
        return idx; // 0-based position
    }
};

template<typename T>
std::vector<T> permutation_by_index(const std::vector<T>& items, uint64_t index){
    size_t n = items.size();
    std::vector<uint64_t> fact(n+1,1);
    for(size_t i=1;i<=n;++i){
        if (fact[i-1] > std::numeric_limits<uint64_t>::max() / i)
            throw std::overflow_error("factorial overflow");
        fact[i] = fact[i-1]*i;
    }
    if (index >= fact[n]) throw std::out_of_range("index");
    Fenwick<> fw((int)n);
    for(int i=0;i<(int)n;++i) fw.add(i,1);
    std::vector<T> out; out.reserve(n);
    uint64_t rem = index;
    for(size_t pos=n; pos>0; --pos){
        uint64_t f = fact[pos-1];
        int digit = (int)(rem / f);
        rem %= f;
        int ix = fw.find_kth(digit);
        out.push_back(items[ix]);
        fw.add(ix, -1);
    }
    return out;
}

For CUDA (answering @chtulu): do the index->permutation work per thread, precompute factorials in constant memory, avoid std::vector on device, and use plain arrays. For small n a simple linear scan for the k-th unused element is easiest on GPU; for larger n implement a per-thread BIT or do host-side generation. If lexicographic order is not required, the O(1) swap-remove method (replace removed slot with last element) is the fastest on massively parallel hardware.