随机算法:Random_shuffle 的一种实现

Fifnmar edited 4 年,10 月前

之前学习过 Linear Conguential Generator, LCG 的实现,觉得很有意思。今天又看到了,并学习了一番著名的「洗牌算法」。

Linear Congruential Generator

LCG 使用一个递归定义:
$$
X_n = (aX_{n-1} + c) \mod m
$$

C++ 标准库在 <random> 中提供了相应 Linear Congruential Engine:

template<
    typename UIntType,  // Seed's type
    UIntType a,
    UIntType c,
    UIntType m
> class linear_congruential_engine;

通过选取合适的 a, c, 和 m,可以使得这个数列表现出很强的随机性。比如,可以使用 Mersenne primes $2^{31} -1$ 或 $2^{61} - 1$ 作 a。

std::linear_congruential_engine<
    unsigned, 2689234, 12345, 15892213> engine(some_seed = 20200202);

已经有一些广泛使用的好用的参数,可以直接拿来用。1969 年,Lewis, Goodman 和 Miller 发现了十分高效的参数:
$$
\begin{cases}
a &=& 16807 \
c &=& 0 \
m &=& 2^{31} -1
\end{cases}
$$
使用这组参数的 LCG,即 $X_n = 16807X_{n-1}\mod {2^{31}-1}$ 又被称作 Minimal Standard。1993 年,Park, Miller 和 Stockmeyer 声称发现了更好的参数 a,这就是「Newer Minimum Standard」:
$$
\begin{cases}a &=& 48271 \c &=& 0 \m &=& 2^{31} -1\end{cases}
$$
Minimum Standard 已经可以很好地满足大部分使用情况。标准已经把上面两个设立为 typedef,分别是 minstd_rand0minstd_rand,可以直接使用。标准库的 default_random_engine 一般就是它们中的一个。

Random Shuffle

有时会用到随机打乱数组的功能,一个朴素的想法是不断生成随机索引对,把这两个随机索引对对应的元素交换。这样做,一方面可能有重复交换且有的元素没有被交换过;另一方面效率太低。

这里就用到洗牌算法。其基本思想是,把数组分成「已打乱的部分」和「未打乱的部分」即剩余部分,每次都从剩余的元素中随机抽一个放到已打乱的部分中。这样,每个元素被选中的概率就是相等的。

代码:

#define current_time                                                           \
    (std::chrono::system_clock::now().time_since_epoch().count())

template <typename Iter> void random_shuffle(Iter const bg, Iter const ed) {
    std::default_random_engine e(current_time);
    for (auto i = ed - bg; i > 0; --i)
        std::swap(bg[i - 1], bg[e() % i]);
}

这里吐糟一下,C++ 的时间库怎么和别的库那么另类啊。

下面是几次测试把序列 $1,2,3…50$ 打乱的结果:

Case 0:
Seed: 1580715449358169800
27 43 46 15 40 33 16 12 36 39 38 21 49 30 35 41 1 13 32 25 2 9 11 44 10
38 21 49 30 35 41 1 13 32 25 2 9 11 44 10 29 7 22 14 47 4 5 37 17 24

Case 1:
Seed: 1580715449360163800
45 48 28 7 6 38 8 9 47 49 16 15 18 17 20 4 12 34 41 29 10 13 26 39 43
16 15 18 17 20 4 12 34 41 29 10 13 26 39 43 35 40 50 21 19 37 23 27 11 36

Case 2:
Seed: 1580715449362159300
44 28 34 26 35 31 12 50 14 10 25 43 11 48 3 40 18 19 27 33 42 29 45 2 32
25 43 11 48 3 40 18 19 27 33 42 29 45 2 32 41 22 8 6 1 37 9 15 5 21

Case 3:
Seed: 1580715449364153800
3 42 2 20 21 41 31 1 17 19 49 34 47 45 12 37 14 23 29 32 50 46 33 35 48
49 34 47 45 12 37 14 23 29 32 50 46 33 35 48 40 30 10 36 39 4 28 16 5 22

Case 4:
Seed: 1580715449365184200
22 8 37 32 23 41 44 7 29 45 17 30 36 43 1 6 24 49 19 5 46 34 20 14 2
17 30 36 43 1 6 24 49 19 5 46 34 20 14 2 9 48 28 50 10 26 47 33 11 16

Case 5:
Seed: 1580715449367146100
15 27 32 17 22 42 31 14 38 24 13 46 30 6 41 5 12 25 49 34 33 37 45 29 23
13 46 30 6 41 5 12 25 49 34 33 37 45 29 23 2 50 1 48 28 3 26 19 21 43

Case 6:
Seed: 1580715449369141200
26 44 21 24 23 41 2 43 46 17 48 32 27 13 18 34 38 30 9 22 14 28 8 15 1
48 32 27 13 18 34 38 30 9 22 14 28 8 15 1 7 37 33 10 16 42 19 12 31 35

Case 7:
Seed: 1580715449370175100
20 25 19 39 30 41 43 44 15 17 50 7 45 8 3 13 2 6 35 5 33 32 36 49 47
50 7 45 8 3 13 2 6 35 5 33 32 36 49 47 48 21 27 31 23 34 16 1 37 29

Case 8:
Seed: 1580715449372132800
38 4 44 2 37 21 41 32 18 11 28 26 9 39 48 8 45 3 1 20 50 25 22 31 6
28 26 9 39 48 8 45 3 1 20 50 25 22 31 6 14 43 23 34 19 10 5 7 24 16

Case 9:
Seed: 1580715449374136900
13 35 3 27 30 40 8 28 16 26 41 33 36 21 12 43 47 29 1 48 15 20 5 24 10
41 33 36 21 12 43 47 29 1 48 15 20 5 24 10 39 19 45 46 11 4 25 18 2 42

Case 10:
Seed: 1580715449376154700
30 28 39 34 41 6 10 15 9 19 13 27 22 46 36 47 25 8 49 24 16 3 18 7 40
13 27 22 46 36 47 25 8 49 24 16 3 18 7 40 31 42 43 32 20 23 17 4 37 38

Case 11:
Seed: 1580715449378148100
17 16 41 20 18 3 47 35 11 22 7 37 49 36 27 40 1 25 13 9 21 5 19 33 8
7 37 49 36 27 40 1 25 13 9 21 5 19 33 8 31 6 34 42 28 14 43 2 48 26

Case 12:
Seed: 1580715449380146400
36 7 50 23 24 13 48 41 5 15 46 10 32 20 38 18 47 43 9 37 33 31 14 26 29
46 10 32 20 38 18 47 43 9 37 33 31 14 26 29 3 4 40 19 17 42 1 45 49 11

Case 13:
Seed: 1580715449382107200
23 46 38 37 2 6 7 41 27 1 11 33 5 48 50 36 22 39 30 12 14 43 31 8 44
11 33 5 48 50 36 22 39 30 12 14 43 31 8 44 19 34 45 13 15 10 26 4 35 24

Case 14:
Seed: 1580715449384100100
22 44 14 8 33 48 4 32 34 38 11 47 18 16 43 24 41 20 30 6 49 17 3 25 1
11 47 18 16 43 24 41 20 30 6 49 17 3 25 1 21 29 37 50 9 45 7 26 2 19

Case 15:
Seed: 1580715449386095000
39 49 25 44 14 27 29 28 23 12 34 7 33 9 1 3 32 24 6 13 35 50 26 47 11
34 7 33 9 1 3 32 24 6 13 35 50 26 47 11 48 4 17 22 31 8 40 37 20 43

Case 16:
Seed: 1580715449388119200
16 40 35 27 17 14 18 41 44 5 50 28 22 15 2 12 29 46 32 9 6 4 31 21 26
50 28 22 15 2 12 29 46 32 9 6 4 31 21 26 42 10 19 43 3 33 1 25 48 23

Case 17:
Seed: 1580715449390084700
20 26 18 28 12 2 23 42 9 46 29 40 16 50 24 32 27 37 10 38 43 22 39 14 49
29 40 16 50 24 32 27 37 10 38 43 22 39 14 49 6 44 1 19 7 21 5 31 36 33

Case 18:
Seed: 1580715449392079100
12 3 45 26 49 24 15 25 8 2 46 36 38 31 20 30 5 19 37 34 41 42 9 35 39
46 36 38 31 20 30 5 19 37 34 41 42 9 35 39 50 11 48 40 29 28 33 17 21 13

Case 19:
Seed: 1580715449394073800
30 8 23 33 28 21 24 27 20 6 5 26 38 35 11 25 48 14 7 46 2 10 39 3 13
5 26 38 35 11 25 48 14 7 46 2 10 39 3 13 37 1 19 45 47 41 16 42 4 12

这个算法的复杂度很容易看出是 $\Theta(n)$ ,而且因为每个元素都得安排,这已经是最低的复杂度了,而且我的 Overhead 也可以说非常小,应该算是很够用。

Churn-out 出来一份测试代码:

void crude_test() {
    constexpr int n = 50;
    std::ofstream out("temp.txt");

    for (int i = 0; i < 20; ++i) {
        out << "Case " << i << ":\n" << "Seed: "
            << std::chrono::system_clock::now().time_since_epoch().count()
            << '\n';

        int a[n];
        std::iota(a, a + n, 1);
        random_shuffle(a, a + n);

        constexpr int line_size = 25;
        for (int i = 0; i < n / line_size; ++i) {
            out << a[i * 10];
            for (int j = 1; j < line_size; ++j)
                out << ' ' << a[i * 10 + j];
            out << '\n';
        }
        out << '\n';

        Sleep(1); // 1ms;
    }
}

进行压力测试:

Scale: 32767 numbers,
Time cost: 2 Time Units.

Scale: 65533 numbers,
Time cost: 5 Time Units.

Scale: 131065 numbers,
Time cost: 11 Time Units.

Scale: 262129 numbers,
Time cost: 24 Time Units.

Scale: 524257 numbers,
Time cost: 53 Time Units.

Scale: 1048513 numbers,
Time cost: 123 Time Units.

Scale: 2097025 numbers,
Time cost: 261 Time Units.

Scale: 4194049 numbers,
Time cost: 811 Time Units.

Scale: 8388097 numbers,
Time cost: 1820 Time Units.

Scale: 16776193 numbers,
Time cost: 4429 Time Units.

Scale: 33552385 numbers,
Time cost: 10442 Time Units.

Scale: 67104769 numbers,
Time cost: 22807 Time Units.

Scale: 134209537 numbers,
Time cost: 48565 Time Units.

Scale: 268419073 numbers,
Time cost: 99549 Time Units.

Scale: 536838145 numbers,
Time cost: 203518 Time Units.

大体上是符合一次曲线的。不过随着规模上升,时间消耗的系数也有上涨,不知道为什么。

Comments

Fifnmar

Orz,测试代码贴错了,是这个

void stress_test() {
    std::ofstream out("Random_shuffle_stress_testing.txt");

    for (int scale = (1 << 15) - 1; scale<std::numeric_limits<int>::max()>> 2;
         scale <<= 1, --scale) {
        out << "Scale: " << scale << " numbers,\n";
        std::vector<int> vec(scale);
        std::iota(vec.begin(), vec.end(), 1);

        auto counter = std::chrono::steady_clock::now();
        random_shuffle(vec.begin(), vec.end());
        out << "Time cost: "
            << (std::chrono::steady_clock::now() - counter).count() / 100000
            << " Time Units.\n\n";
    }
}