题目来源:https://leetcode.com/problems/random-pick-with-blacklist/description/

标记难度:Hard

提交次数:3/6

代码效率:

  • 二分:10.06%
  • 重映射:45.67%

题意

给定区间[0, N)和一个该区间内的“黑名单”,要求以均匀的概率返回该区间内非黑名单的数字,且调用rand()的次数尽量少。

分析

二分

我感觉这是一道很妙的题。看到题之后,我的第一反应是:假定黑名单的长度为M,那么我们只需要随机[0, N - M)范围内的数字,然后把它们映射到[0, N)区间内就可以了。问题是怎么映射呢?

于是我想出了这样一种做法:随机得到一个数字之后,计算黑名单中有多少个数小于等于这个数,然后把这个数量加到这个数上,返回结果。然后我就把代码交上去了,于是,很快我就发现这个做法错得有点离谱,因为它完全没有保证返回的数不在黑名单里。反例如N = 4, blacklist = [0, 1],随机得到0时,上述做法会返回1。

于是我尝试从另一个角度来思考这个问题。我们显然可以显式地构建一个映射:顺序枚举[0, N - M)范围内的数,维护一个指针,指向下一个能够被映射的数,遇到黑名单中的数则跳过。这个做法是正确的,但考虑到1 <= N <= 1000000000的数据量,显然不可能把整个映射表都建立出来,然后去查表。

于是这个思路被否决了,我又回头去考虑怎么在线计算出准确的映射这个问题。我想了想:考虑“黑名单中有多少个数小于等于当前的数”这个思路其实是正确的,问题在于,被考虑的不应该是[0, N - M)中的数,而应该是[0, N)中的数。考察这个例子,N = 9, balcklist = [2, 3, 5],令blacklist_before(i)表示黑名单中<=i的数的个数:

[0, N)中的数i blacklist_before(i) i - blacklist_before(i) 映射到[0, N - M)
0 0 0 0
1 0 1 1
2 1 1 -
3 2 1 -
4 2 2 2
5 3 2 -
6 3 3 3
7 3 4 4
8 3 5 5

可以观察得到一些非常有趣的性质:

  • i - blacklist_before(i)是非单调递增的,因为它代表的是[0, i]区间内不属于黑名单内的数的个数
  • i - blacklist_before(i)没有递增时,表示i是一个黑名单内的数字
  • i[0, N - M)区间中应该映射到的数与i - blacklist_before(i)两列很接近(从意义上来说也是)

所以可以从这个逆向映射的角度考虑,用二分的方法解决问题:假定在[0, N - M)中随机得到了y,我们需要找到满足x - blacklist_before(x) == yxlower_bound。然后就可以写了。


在讨论区里我看到了类似的做法,但是思考的角度不太一样。我们可以把从区间[0, N - M)中生成随机数r的情形这样分类:[1]

  • r[0, B[0])区间内,可以直接返回r
  • r[B[0], B[1] - 1)区间内,应返回r + 1
  • ...
  • r[B[i]-i, B[i+1]-(i+1))区间内,应返回r + i + 1。注意到r + i + 1位于[B[i] + 1, B[i+1])区间内,因此这样做是安全的。

因此可以在B[i] - (i+1)数组中进行二分查找。

这种做法是在经过处理的blacklist数组上进行二分查找,而我是在[0, N)区间上直接查找,所以比我的做法更压缩一些。

HashMap重映射法

我之前考虑过建立整个映射表,但因为数据量而放弃了。然而,事实上,完全不需要把整个映射表都建立出来;或者不如说,完全不需要按顺序建立映射,只要保证所有数字都可以被随机取到即可。所以可以采取这样的一种做法:将黑名单中的元素分成两类,一类在[0, N - M)区间内,一类在[N - M, N)区间内。然后在[N - M, N)区间内为[0, N - M)区间内的黑名单元素顺序寻找映射(我看到了从前向后[2]和从后向前[3]两种方法,不过本质上是相同的),同时注意跳过那些同样在黑名单内的元素。

上述从后向前找映射的一个图示例子:N=10, blacklist=[3, 5, 8, 9],将3和5映射为7和6。

从后向前映射

代码

二分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
class Solution {
private:
int N;
vector<int> blacklist;
unordered_set<int> blackSet;
random_device rd;
mt19937 gen;
uniform_int_distribution<> dist;

int findMap(int y) {
int l = 0, r = N - 1;
// need to find lower_bound
// 二分真是一件Tricky的事情
while (l < r) {
int mid = (l + r) / 2;
int smaller = upper_bound(blacklist.begin(), blacklist.end(), mid) - blacklist.begin();
if (mid - smaller < y)
l = mid + 1;
else
r = mid;
}

return l;
}

public:
Solution(int N, vector<int> blacklist): gen(rd()), dist(0, N - blacklist.size() - 1) {
this->N = N;
this->blacklist = blacklist;
sort(this->blacklist.begin(), this->blacklist.end());
for (int b: blacklist)
blackSet.insert(b);
}

int pick() {
int r = dist(gen);
int m = findMap(r);
// cout << r << ' ' << m << endl;
// 事实证明,保证lower_bound之后,找到的数必然不是黑名单中的数
/* while (blackSet.find(m) != blackSet.end())
m++; */
return m;
}
};

/**
* Your Solution object will be instantiated and called as such:
* Solution obj = new Solution(N, blacklist);
* int param_1 = obj.pick();
*/
// 这个问题很有趣。
// 既然从N个数的区间中屏蔽了B.length个点,那么就是从N - B.length的区间中随机选点
// 然后换算一下。
// 但是换算方法看起来是个很大的问题……

重映射

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
class Solution {
private:
int N;
vector<int> blacklist;
unordered_set<int> blackSet;
unordered_map<int, int> blackMap;
random_device rd;
mt19937 gen;
uniform_int_distribution<> dist;

public:
Solution(int N, vector<int> blacklist): gen(rd()), dist(0, N - blacklist.size() - 1) {
this->N = N;
this->blacklist = blacklist;
sort(this->blacklist.begin(), this->blacklist.end());
for (int b: blacklist)
blackSet.insert(b);

int M = blacklist.size();
int mapTo = N - M;
for (int b: this->blacklist) {
if (b >= N - M)
break;
while (blackSet.find(mapTo) != blackSet.end())
mapTo++;
blackMap[b] = mapTo++;
}
}

int pick() {
int r = dist(gen);
if (blackSet.find(r) != blackSet.end())
return blackMap[r];
return r;
}
};

/**
* Your Solution object will be instantiated and called as such:
* Solution obj = new Solution(N, blacklist);
* int param_1 = obj.pick();
*/

  1. Simple Java solution with Binary Search

  2. Java O(B) constructor and O(1) pick, HashMap

  3. Java O(B) / O(1), HashMap

题目来源:https://leetcode.com/problems/generate-random-point-in-a-circle/description/

标记难度:Medium

提交次数:1/1

代码效率:

  • 极坐标法:64.06%
  • 拒绝取样:97.20%

题意

给定平面上一个圆的圆心位置和半径,从圆中以均匀的概率随机选取点。

分析

拒绝取样

其实我的第一反应是用拒绝取样(Rejection Sampling)的思路来做:首先从这个圆的与坐标轴平行的外切正方形中均匀随机选取点,然后判断点是否位于圆中;如果不在,重新生成一个新的点,再次进行判断;否则直接返回。

直觉上来说,拒绝取样显然是正确的;不过我们可以用一种稍微更加形式化的方法来描述。(以下内容参考了拒绝采样(reject sampling)的简单认识,非常直观形象。)

下图是一个随机变量的密度函数曲线,试问如何获得这个随机变量的样本呢?

如果你像我一样,已经把概率论与数理统计统统还给数学老师了,那么提示一下,概率密度函数(PDF)是累积分布函数(CDF)的导数,反映的是概率的“密集程度”。你可以设想一根极细的无穷长的金属杆,概率密度相当于杆上各点的质量密度。由于上述原因,概率密度函数曲线下的面积表示的就是取值概率。(参考了应该如何理解概率分布函数和概率密度函数?

概率密度函数曲线

我们首先用一个矩形将这个密度曲线套起来,把密度曲线框在一个矩形里,如下图所示:

被框起来的概率密度函数曲线

然后,向这个矩形里随机投点。随机投点意味着在矩形这块区域内,这些点是满足均匀分布的。投了大概10000个点,如下面这个样子:

投点之后的分布情况

显然,有的点落在了密度曲线下侧,有的点落在了密度曲线的上侧。上侧的点用绿色来表示,下侧的点用蓝色来表示,如下图:

落在上侧和下侧的点

只保留密度曲线下侧的点,即蓝色的点:

保留下侧的点

到这里,提一个问题:在密度曲线以下的这块区域里,这些点满足什么分布?均匀分布!这是拒绝采样最关键的部分,搞个矩形、向矩形里投点等等,所做的一切都是为了获得一个密度曲线所围成区域的均匀分布。只要能获得这样一个在密度曲线下满足均匀分布的样本,我们就可以获得与该密度曲线相匹配的随机变量的采样样本。方法是,只需把每个蓝点的横坐标提取出来,这些横坐标所构成的样本就是我们的目标样本。下图左侧,是按照以上方法获得的一个样本的直方图以及核密度估计,下图右侧,是开始的密度曲线。

统计结果

Reject sampling solution中给出了一种简单明了的用拒绝采样实现此题的方法。

极坐标法

但是我觉得这个方法的逼格不够高。于是我心想,我们怎么表示一个圆内的点呢?那自然就是用极坐标法,半径+角度。那么我们随机一个半径,再随机一个角度,最后换算成直角坐标系里的坐标,不就可以了嘛!

……只是听上去可以而已……

交上去之后我发现,一共有8个测试样例,最后一个我总是过不去。我的随机性是不是真的写出了bug?我实在是想不出来错在哪里了,于是去看了看题解区——原来有这么多人和我有一样的困惑。是的,上述思路确实有问题。

当我写出上述代码的时候,我心里想的是:先用角度随机选择一条半径,然后在这条半径上随机选择一个点。但问题是我们不能把圆这样看成是很多半径的集合——或者说,这条半径上不同位置的点的密集程度是不一样的。显然距离圆心更远的一端被选择的概率应该更大。

在圆上随机选择

直接对角度和半径都进行随机取样会产生这样的后果,靠近圆心的点被选择的概率偏大(Uniform random points in a circle using polar coordinates):

左侧是错误的随机结果,右侧是正确的结果

我们重新考虑一下半径的问题。以半径$r$作为随机变量,则随机点落在$r$范围内的概率分布为

$$\text{CDF}(r) = \frac{r^2}{R^2}$$

$$\text{PDF}(r) = \frac{d}{dr} \text{CDF}(r) = \frac{2r}{R^2}$$

[0, 1]均匀分布对$\text{CDF}(r)$进行随机,则有$r = R \sqrt{\text{rand}()}$。

C++中的数学运算

使用std::uniform_real_distribution生成随机实数比直接用rand()好得多,至少stackoverflow上是这么说的。具体使用方法如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include <iostream>
#include <iomanip>
#include <string>
#include <map>
#include <random>

int main()
{
std::random_device rd;

//
// Engines
//
std::mt19937 e2(rd());
//std::knuth_b e2(rd());
//std::default_random_engine e2(rd()) ;

//
// Distribtuions
//
std::uniform_real_distribution<> dist(0, 10);
//std::normal_distribution<> dist(2, 2);
//std::student_t_distribution<> dist(5);
//std::poisson_distribution<> dist(2);
//std::extreme_value_distribution<> dist(0,2);

double rand = dist(e2));

return 0;
}

目前C++中并没有提供什么特别高级的获得PI值的方法,仍然是用C中提供的M_PI宏:

1
2
3
#define _USE_MATH_DEFINES  // 根据实际情况;参见https://stackoverflow.com/questions/1727881/how-to-use-the-pi-constant-in-c
#include <math.h>
M_PI

std::cosstd::sin函数使用的是弧度。好久不写,都忘光了。

代码

拒绝取样

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
class Solution {
private:
double radius, x_center, y_center, x_leftdown, y_leftdown;
std::random_device rd;
std::mt19937 gen;
std::uniform_real_distribution<> dis;

public:
Solution(double radius, double x_center, double y_center): gen(rd()), dis(0, 1) {
this->radius = radius;
this->x_center = x_center;
this->y_center = y_center;
x_leftdown = x_center - radius;
y_leftdown = y_center - radius;
}

vector<double> randPoint() {
double x, y;
do {
x = x_leftdown + 2 * radius * dis(gen);
y = y_leftdown + 2 * radius * dis(gen);
} while ( (x - x_center) * (x - x_center) + (y - y_center) * (y - y_center) > radius * radius);

return {x, y};
}
};

极坐标

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
class Solution {
private:
double radius, x_center, y_center;
std::random_device rd;
std::mt19937 gen;
std::uniform_real_distribution<> randDeg, randRadius;

public:
Solution(double radius, double x_center, double y_center): gen(rd()), randDeg(0, 2 * M_PI), randRadius(0, 1) {
this->radius = radius;
this->x_center = x_center;
this->y_center = y_center;
}

vector<double> randPoint() {
double r = sqrt(randRadius(gen)) * radius; // !
double deg = randDeg(gen);
double x = x_center + r * cos(deg);
double y = y_center + r * sin(deg);

return {x, y};
}
};

题目来源:https://leetcode.com/problems/linked-list-random-node/description/

标记难度:Medium

提交次数:3/3

代码效率:

  • 平凡的解法:71.27%
  • 平凡但更好的写法:99.11%
  • 水塘抽样法:70.97%

题意

从一个链表中随机取元素。

分析

平凡的思路

一种非常直接的想法是把链表当成数组,先遍历一遍链表,得到链表的总长度,然后每次按下标取随机数,在链表中顺序查找对应的数。此时预处理的时间复杂度为O(N),每次获取随机数的平摊时间复杂度为O(N)。当然,如果直接新开一个O(N)的数组(vector),把链表缓存下来,然后直接按下标寻址,则每次获取随机数的时间复杂度会下降到O(1)

虽然我们平时常用的生成随机数的搭配是srand(time(0))rand(),但它们其实不是很符合伪随机数的要求,一部分原因是rand()生成的随机数的大小最大为RAND_MAX。(参见RAND_MAX)既然我们写的是C++,那么更标准的方法是采用std::uniform_int_distribution使用方法如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#include <random>
#include <iostream>

int main()
{
std::random_device rd; // 将用于为随机数引擎获得种子
std::mt19937 gen(rd()); // 以播种标准 mersenne_twister_engine
std::uniform_int_distribution<> dis(1, 6);

for (int n=0; n<10; ++n)
// 用 dis 变换 gen 所生成的随机 unsigned int 到 [1, 6] 中的 int
std::cout << dis(gen) << ' ';
std::cout << '\n';
}

我写了相应的代码,但实在是太麻烦了,平时还是用rand()方便啊!

水塘抽样法

至于水塘抽样法——它的优点是不需要提前扫描一遍以获得整个链表的大小,缺点是每次都需要遍历整个链表并取多次随机数来进行抽样,虽然时间复杂度仍然是O(N),但隐含的常数变大了。我认为它适合的场景是一次取多个随机数,而不是这种情况。

我总是记不住水塘抽样法的原理,不如在这里把维基上的最简单的例子复述一遍:假定你想要从网易云音乐的全体曲库(显然原文不是网易云音乐,我就随便举个我熟悉例子)里随机取出10首乐曲,保存在“我喜欢的音乐”歌单中,但是你并不知道网易云音乐的曲库有多大。此时,你可以采用这种方法:

  • 顺序浏览所有乐曲。
  • 将前10首乐曲保存到歌单中。
  • 浏览到第$i$($i > 10$)首乐曲的时候,以$\frac{10}{i}$的概率保存这首新的乐曲(同时删除一首已保存的旧的曲子;其中每首旧乐曲被删除的概率都是$\frac{1}{10}$),也即以$1 - \frac{10}{i}$的概率删除这首乐曲。
  • 不断重复上一过程,直到遍历完所有乐曲。

如果网易云音乐里一共只有10首乐曲,那么显然每首乐曲被保存下来的概率都是1。如果有11首乐曲,则在浏览到第11首乐曲的时候,我们保存这首乐曲的概率为$\frac{10}{11}$;对于之前的旧乐曲,其中某一首被保存下来的概率为

$$P(\text{新乐曲没有被选择保存}) + P(\text{新乐曲被选择保存}) \cdot P(\text{这一首没有被选择替换掉}) \\
= \frac{1}{11} + \frac{10}{11} \cdot \frac{9}{10} = \frac{10}{11}$$

也就是说,到目前为止,对于这11首乐曲,每一首仍然处于歌单中的概率都是$\frac{10}{11}$,这是符合我们的要求的。

在浏览到第12首乐曲的时候,我们保存这首乐曲的概率为$\frac{10}{12}$;而对于之前的那11首乐曲,它们仍然处于歌单中的概率为

$$P(\text{在只看到过11首曲子的的时候,它们还在歌单里}) \cdot (P(\text{第12首乐曲被直接丢弃了}) \\ + P(\text{第12首乐曲被选择保留}) \cdot P(\text{这一首没有被选择替换掉})) \\
= \frac{10}{11} (\frac{2}{12} + \frac{10}{12} \cdot \frac{9}{10})$$

也就是说,现在每一首乐曲仍然处于歌单中的概率都是$\frac{10}{12}$,这是符合常理的。

严谨的证明需要用到数学归纳法,不过形式和上述思考过程非常类似,所以我就懒得再抄一遍了。

我觉得这是一种非常精妙的incremental的算法。

代码

平凡的解法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
/**
* Definition for singly-linked list.
* struct ListNode {
* int val;
* ListNode *next;
* ListNode(int x) : val(x), next(NULL) {}
* };
*/
class Solution {
int length;
ListNode* head;
public:
/** @param head The linked list's head.
Note that the head is guaranteed to be not null, so it contains at least one node. */
Solution(ListNode* head) {
length = 0;
ListNode* p = head;
this->head = head;
while (p != NULL) {
length++;
p = p->next;
}
srand(time(0));
}

/** Returns a random node's value. */
int getRandom() {
int x = rand() % length;
int i = 0;
ListNode* p = head;
while (i < x) {
p = p->next;
i++;
}
return p->val;
}
};

/**
* Your Solution object will be instantiated and called as such:
* Solution obj = new Solution(head);
* int param_1 = obj.getRandom();
*/

平凡但更好的写法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
/**
* Definition for singly-linked list.
* struct ListNode {
* int val;
* ListNode *next;
* ListNode(int x) : val(x), next(NULL) {}
* };
*/
class Solution {
int length;
ListNode* head;
random_device rd; // 将用于为随机数引擎获得种子
mt19937* gen; // 以播种标准 mersenne_twister_engine
uniform_int_distribution<>* dis;
public:
/** @param head The linked list's head.
Note that the head is guaranteed to be not null, so it contains at least one node. */
Solution(ListNode* head) {
length = 0;
ListNode* p = head;
this->head = head;
while (p != NULL) {
length++;
p = p->next;
}
gen = new mt19937(rd());
dis = new uniform_int_distribution<>(0, length-1);
}

/** Returns a random node's value. */
int getRandom() {
int x = (*dis)(*gen);
int i = 0;
ListNode* p = head;
while (i < x) {
p = p->next;
i++;
}
return p->val;
}
};

/**
* Your Solution object will be instantiated and called as such:
* Solution obj = new Solution(head);
* int param_1 = obj.getRandom();
*/

水塘抽样法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
/**
* Definition for singly-linked list.
* struct ListNode {
* int val;
* ListNode *next;
* ListNode(int x) : val(x), next(NULL) {}
* };
*/
class Solution {
ListNode* head;
public:
/** @param head The linked list's head.
Note that the head is guaranteed to be not null, so it contains at least one node. */
Solution(ListNode* head) {
this->head = head;
}

/** Returns a random node's value. */
int getRandom() {
// 进行水塘抽样
int ans = head->val;
int i = 2;
ListNode* p = head->next;
while (p != NULL) {
int x = rand() % i;
if (x == 0)
ans = p->val;
i++;
p = p->next;
}
return ans;
}
};

/**
* Your Solution object will be instantiated and called as such:
* Solution obj = new Solution(head);
* int param_1 = obj.getRandom();
*/