如何用 shuffle 排序

· · 休闲·娱乐

看到这个标题，相信大家都会想到猴子排序。假如我们观察猴子排序的最后一次打乱过程，发现这就相当于随机数生成器给出了一个完美的随机数序列，恰好使其打乱后变得有序，这启发我们思考：能否设计一个特殊的随机数生成器，使得用初始序列构造后再传入 std::shuffle 后可以恰好是序列升序排序。

答案是肯定的，那这魔法一般的操作该如何实现呢，请听我一步步讲解。

首先，我们需要能够自定义一个随机数生成器，要求传入 std::shuffle 后能够正常运行。实际上，这件事的要求并不高，只需要在类中定义 result_type 类型，result_type operator()() 用于生成随机数，static constexpr result_type min() 和 static constexpr result_type max() 分别表示可以生成的最大值和最小值即可。

举个例子：

class R{
public:
    using result_type=uint32_t;
    R():state(1){}
    R(result_type seed):state(seed){}
    result_type operator()(){
        state^=state<<13;
        state^=state>>17;
        state^=state<<5;
        return state;
    }
    static constexpr result_type min(){return 0;}
    static constexpr result_type max(){return numeric_limits<result_type>::max();}
private:
    result_type state;
};

这样，我们就实现了一个使用 shift-xor 算法的随机数生成器。

接下来，我们分析一下 std::shuffle 的具体实现（经过部分美化修改）。

template<typename _RandomAccessIterator, typename _UniformRandomNumberGenerator>
voidshuffle(_RandomAccessIterator __first, _RandomAccessIterator __last, _UniformRandomNumberGenerator&& __g){
    if (__first == __last) return;

    typedef typename iterator_traits<_RandomAccessIterator>::difference_type _DistanceType;
    typedef typename std::make_unsigned<_DistanceType>::type __ud_type;
    typedef typename std::uniform_int_distribution<__ud_type> __distr_type;
    typedef typename __distr_type::param_type __p_type;
    typedef typename remove_reference<_UniformRandomNumberGenerator>::type _Gen;
    typedef typename common_type<typename _Gen::result_type, __ud_type>::type __uc_type;

    const __uc_type __urngrange = __g.max() - __g.min();
    const __uc_type __urange = __uc_type(__last - __first);

    if (__urngrange / __urange >= __urange){
        _RandomAccessIterator __i = __first + 1;
        if ((__urange % 2) == 0){
            __distr_type __d{0, 1};
            std::iter_swap(__i++, __first + __d(__g));
        }
        while (__i != __last){
            const __uc_type __swap_range = __uc_type(__i - __first) + 1;
            const pair<__uc_type, __uc_type> __pospos = __gen_two_uniform_ints(__swap_range, __swap_range + 1, __g);
            std::iter_swap(__i++, __first + __pospos.first);
            std::iter_swap(__i++, __first + __pospos.second);
        }
        return;
    }
    __distr_type __d;
    for (_RandomAccessIterator __i = __first + 1; __i != __last; ++__i){
        std::iter_swap(__i, __first + __d(__g, __p_type(0, __i - __first)));
    }
}

是不是很难理解，别急，让我们一步步分析。首先，我们先不管那些 typedef，并假设 __urngrange / __urange >= __urange 不成立，于是代码进入到最后的循环中。其中 __d 是一个 std::uniform_int_distribution 的随机数传输器，循环的含义是对于 i=1,2,3,\dots,n-1，选择一个 j=0,1,2,\dots,i，并交换 a_i 和 a_j，这就是著名的 Knuth 洗牌算法。

接下来再分析 __urngrange / __urange >= __urange 成立时时怎么样的。先看 __gen_two_uniform_ints 的实现。

template<typename _IntType, typename _UniformRandomBitGenerator>
pair<_IntType, _IntType> __gen_two_uniform_ints(_IntType __b0, _IntType __b1, _UniformRandomBitGenerator&& __g){
    _IntType __x = uniform_int_distribution<_IntType>{0, (__b0 * __b1) - 1}(__g);
    return std::make_pair(__x / __b1, __x % __b1);
}

可以发现，这就是把一个随机数拆成两个用，这两个随机数分别在 [0,__b0) 和 [0,__b1) 范围内。

我们回到上面的代码中，发现 if 内外的算法其实是相同的，if 判断里的代码也只是起到了在区间长度较小时提升算法效率的作用。

std::shuffle 的原理分析完了，接下来我们考虑如何让 std::uniform_int_distribution 生成我们想要的数。

std::uniform_int_distribution 的源代码极度繁琐，故这里只给出最重要的部分简化后的代码。

template<typename _Wp, typename _Urbg, typename _Up>
static _Up _S_nd(_Urbg& __g, _Up __range){
    using _Up_traits = __gnu_cxx::__int_traits<_Up>;
    using _Wp_traits = __gnu_cxx::__int_traits<_Wp>;

    // reference: Fast Random Integer Generation in an Interval
    // ACM Transactions on Modeling and Computer Simulation 29 (1), 2019
    // https://arxiv.org/abs/1805.10941
    _Wp __product = _Wp(__g()) * _Wp(__range);
    _Up __low = _Up(__product);
    if (__low < __range){
        _Up __threshold = -__range % __range;
        while (__low < __threshold){
            __product = _Wp(__g()) * _Wp(__range);
            __low = _Up(__product);
        }
    }
    return __product >> _Up_traits::__digits;
}

template<typename _IntType>
template<typename _UniformRandomBitGenerator>
typename uniform_int_distribution<_IntType>::result_type
uniform_int_distribution<_IntType>::
operator()(_UniformRandomBitGenerator& __urng, const param_type& __param){
    typedef typename _UniformRandomBitGenerator::result_type _Gresult_type;
    typedef typename make_unsigned<result_type>::type __utype;
    typedef typename common_type<_Gresult_type, __utype>::type __uctype;

    constexpr __uctype __urngmin = _UniformRandomBitGenerator::min();
    constexpr __uctype __urngmax = _UniformRandomBitGenerator::max();
    constexpr __uctype __urngrange = __urngmax - __urngmin;

    const __uctype __urange = __uctype(__param.b()) - __uctype(__param.a());

    __uctype __ret;

    const __uctype __uerange = __urange + 1;

    __UINT32_TYPE__ __u32erange = __uerange;
    __ret = _S_nd<__UINT64_TYPE__>(__urng, __u32erange);
    return __ret + __param.a();
}

_S_nd 函数使用拒绝采样法生成 [0,__range) 内无偏的随机数，我们要想办法让我们的生成器生成的数从不会被拒绝并返回我们想要的值。

若不被拒绝，则它返回的数是 \left\lfloor\dfrac{r\times range}{2^{32}}\right\rfloor，其中 r 为随机数生成器生成的数字。我们可以选择让随机数生成器生成 \left\lceil\dfrac{x\times2^{32}}{range}\right\rceil 使得它返回 x。为了避免被拒绝，于是可以在返回值上加 1 来保证 low\ge range。

现在，我们终于完成了所有准备工作，可以开始设计随机数生成器了。

首先，我们考虑如何算出每个数需要和谁交换，这非常简单，只需要先偷偷排序，记录每个数需要到达我位置，再从后往前构造，从排序后的状态开始判断每个数的交换对象。

算出操作序列后再判断一下 __urngrange / __urange >= __urange 是否成立，若成立则把相邻两个数合并即可，注意要严格按照 stl 源码实现。

为了让这个生成器更像一个真的随机数生成器，可以在里面加上一个 std::mt19937，若前面用于排序的“随机数”生成完毕后就调用真正的生成器来生成。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;

class bogo_sort{
public:
    using result_type=uint32_t;
    bogo_sort(const vector<int> &a):p(0),r(reduce(a.begin(),a.end(),0,[](int x,int y){return x^y;})){
        int n=a.size();
        vector<pair<int,int>> b(n);
        for (int i=0;i<n;i++) b[i]={a[i],i};
        sort(b.begin(),b.end());
        vector<int> id(n),pos(n);
        for (int i=0;i<n;i++) id[i]=b[i].second,pos[id[i]]=i;
        gen.resize(n-1);
        for (int i=n-1;i>=1;i--){
            int p=pos[i];
            gen[i-1].first=p;
            int t=id[i];
            swap(id[i],id[p]);
            pos[i]=i,pos[t]=p;
        }
        for (int i=0;i<n-1;i++) gen[i].second=i+1;
        if (max()/n>=result_type(n)){
            int i=((n&1)==0),j=i;
            for (;i<n-2;i+=2,j++){
                result_type x=gen[i].second,y=gen[i+1].second;
                gen[j]={gen[i].first*(y+1)+gen[i+1].first,(x+1)*(y+1)-1};
            }
            gen.resize(j);
        }
    }
    result_type operator()(){
        if (p==int(gen.size())) return r();
        p++;
        return ((ull(gen[p-1].first)<<32)+gen[p-1].second)/(gen[p-1].second+1)+1;
    }
    static constexpr result_type min(){return 0;}
    static constexpr result_type max(){return numeric_limits<result_type>::max();}
private:
    vector<pair<result_type,result_type>> gen;
    int p;
    mt19937 r;
};

int main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout.tie(nullptr);
    int n;
    cin>>n;
    vector<int> a(n);
    for (auto &i:a) cin>>i;
    bogo_sort rnd(a);
    shuffle(a.begin(),a.end(),rnd);
    for (auto i:a) cout<<i<<" ";
    return 0;
}

提交记录。