在人间漂流

厚德 求真 励学 笃行


  • 首页

  • 关于

  • 标签

  • 分类

  • 归档

  • 公益404

spark tips

发表于 2019-09-01 | 分类于 Distributed System , Bigdata Technology | 0 comments
阅读次数:
  |   字数统计: 710(字)   |   阅读时长: 3(分)

spark tips

jar包上传到hdfs

  • 分布式集群上跑spark时每次都要将${SPARK_HOME}/jars/下所有的jar包上传到分布式系统中供调用,为了避免每次的上传时间可以将这些jar包存放到hdfs中,==注意:我在2.\的版本中将所有jar包打包成.archive文件上传到hdfs,但是这样做导致spark在yarn上启动失败,后来尝试将这些jar包分散上传到hdfs的文件夹下就能正常在yarn上启动spark了。==具体做法如下:
1
2
3
4
将${SPARK_HOME}/jars下的所有jar包上传到hdfs中:
hadoop fs -put jars/* /spark
在${SPARK_HOME}/conf/spark-defaults.conf中添加配置项:
spark.yarn.jars hdfs://localhost:9000/spark/jars/*
  • 如果要使用hbase则需要将hbase的lib目录下hbase开头的jar包和guava-{version}.jar、protobuf-java-{version}.jar上传到hdfs上的同一目录jars,hbase.xml复制一份到hadoop的配置文件目录etc/hadoop下

spark on yarn

1. 调参

1.1 client mode

  • client mode下,spark的driver与applicationmaster是分离的,要分别设置内存大小。任务提交的同时driver就开始运行了,所以driver内存不可以在程序里使用SparkConf类设置,可通过spark-defaults.conf或者CLI下—driver-memory配置,而applicationmaster的内存大小没有限制可以在SparkConf里面通过spark.yarn.am.memory来配置

1.2 cluster mode

  • cluster mode下,spark的driver和applicationmaster在同一个container里面,可以通过在SparkConf,spark-defaults.conf或者CLI下用spark.yarn.am.memory来对这个container内存大小进行配置

1.3 CLI参数设置

spark在yarn上运行能使用的container数目是由核数vcores和executor内存大小共同决定的,由于hadoop不像spark,hadoop可以设置一个结点的总内存和核数,所以当cpu或者内存资源超过设置的额度时会限制container到尽可能大的数目

  • —num-executors设置executor数目
  • —executor-memory设置executor内存
  • —executor-cores设置executor核数
  • —conf spark.yarn.am.memory=2g设置application master内存大小
  • —conf spark.yarn.am.cores=2设置application master核数

1.4 运行方式

  • cluster
1
spark-submit --class <yourclass> --master yarn --deploy-mode cluster --conf spark.yarn.am.memory=2g --executor-memory 2g --executor-cores 2 <yourjar> <args>
  • client
1
spark-submit --class <yourclass> --master yarn --deploy-mode client --driver-meemoty 2g --executor-memory 2g --executor-cores 2 <yourjar> <args>

spark standalone

1.调参

  • spark无需像mapreduce那样给每个节点设置总的内存限制,spark总内存使用会自动扩展。

  • spark集群总核数无需配置,自动识别集群中所有机器的虚拟内核数

  • spark的executor数目由集群中的总核数和exexutor的核数决定:

  • executor.memory不等于RDD可用的memory size,因为JVM等也需要占据一定内存资源,差不多只能用一半,driver.memory也一样

  • 运行方式:

    1
    spark-submit --class <yourclass> --master spark://slave103:7077 --conf spark.driver.memery=6g --conf spark.executor.memory=2g --conf spark.executor.cores=2 <yourjar> <args>
---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址:
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢!

TensorBoard tutorial

发表于 2019-09-01 | 分类于 Machine Learning , Tensorflow | 0 comments
阅读次数:
  |   字数统计: 45(字)   |   阅读时长: 1(分)
参考链接:
  1. google关于tensorboard的官方视频教程
  2. 视频教程对应代码
  3. 视频配套博客:TENSORBOARD可视化
  4. tensorflow官方网站关于tensorboard projector简介
  5. my tensorboard notebook
---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址:
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢!

后缀自动机(SAM)的实现与应用

发表于 2019-08-16 | 分类于 Algorithm , ACM Solution | 0 comments
阅读次数:
  |   字数统计: 5k(字)   |   阅读时长: 25(分)

原理

后缀自动机(Suffix Automaton,简称SAM)是一个能解决许多字符串相关问题的有力的数据结构。举几个例子

  • 一个字符串是否出现在另一个字符串中
  • 一个字符串包含的各不相同的子串数目
  • 查询一个字符串在另一个字符串中的出现次数

  • 一个字符串在另一个字符串中第一次出现位置或者所有出现位置

  • 求两个或多个字符串的最长公共子串

上面这几个例子在构建好SAM之后都可以在线性时间内解决

具体的原理可以参考这篇文章1,里面不仅给出了算法流程还给出了算法的正确性证明.

1. https://oi-wiki.org/string/sam/ ↩

实现

首先我们假设字符集大小为常数。如果字符集大小不是常数,SAM 的时间复杂度就不是线性的。从一个结点出发的转移存储在支持快速查询和插入的平衡树(map)中。因此如果我们记$\sum$为字符集,$|\sum|$为字符集大小,则算法的渐进时间复杂度为$O(nlog|\sum|)$,空间复杂度为$O(n)$。然而如果字符集足够小,可以不写平衡树,以空间换时间将每个结点的转移存储为长度为$|\sum|$的数组(用于快速查询)和链表(用于快速遍历所有可用关键字)。这样算法的时间复杂度为$O(n)$,空间复杂度为$O(n|\sum|)$。

这两种转移方式我都实现了一遍,并且在SAM基础上解决了前面列举的问题,下面是具体实现的代码

平衡树(map)版本

展开代码
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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
#include <bits/stdc++.h>
using namespace std;

// 后缀自动机,可以在O(n)复杂度内创建一个字符串的SAM
// 最多2*n - 1个状态节点(n >= 2),最多3*n - 4个转换(n >= 3)
namespace SAM {
​ struct state {
​ // SAM必须变量
​ int len, link; // len表示状态包含的后缀字符串中的最长串长度,link为后缀链接len(link[v]) + 1 = minlen(v)
​ std::map<char, int> next;

// 用于求子串出现次数的变量
int cnt; // 子串出现次数

// firstpos用于求子串第一次出现位置,配合is_clone和inv_link可用于求子串所有出现位置
int firstpos; // 子串第一次出现位置
bool is_clone; // 标志节点是否是复制的
vector<int> inv_link; // 节点的后缀引用列表,该节点对应串是列表内节点的后缀

// ml和nl用于求多字符串的最长公共子串
int nl; // 当前字符串与状态节点对应最长后缀的公共子串长
int ml; // 所有字符串与状态节点对应最长后缀的公共子串长中的最小长度
}; // 状态节点,一个状态对应一个endpoint集合

const int MAXLEN = 100000; // 字符串长度限制
state st[MAXLEN * 2];
int sz, last;

// 初始化
void sam_init() {
for (int i = 0; i < sz; i++) {
st[i].next.clear();
st[i].inv_link.clear();
}
st[0].len = 0;
st[0].link = -1;
sz = 1;
last = 0;
}

// 向sam中插入字符
void sam_extend(char c) {
int cur = sz++;
st[cur].len = st[last].len + 1;
st[cur].ml = st[cur].len;
st[cur].nl = 0;
st[cur].firstpos = st[cur].len - 1;
st[cur].is_clone = false;
st[cur].cnt = 1;
int p = last;
while (p != -1 && !st[p].next.count(c)) {
st[p].next[c] = cur;
p = st[p].link;
}
if (p == -1) {
st[cur].link = 0;
} else {
int q = st[p].next[c];
if (st[p].len + 1 == st[q].len) {
st[cur].link = q;
} else {
int clone = sz++;
st[clone].len = st[p].len + 1;
st[clone].ml = st[clone].len;
st[clone].nl = 0;
st[clone].is_clone = true;
st[clone].firstpos = st[q].firstpos;
st[clone].cnt = 0;
st[clone].next = st[q].next;
st[clone].link = st[q].link;
while (p != -1 && st[p].next.count(c) && st[p].next[c] == q) {
st[p].next[c] = clone;
p = st[p].link;
}
st[q].link = st[cur].link = clone;
}
}
last = cur;
}
}

// 获取子串数目
bool vis[SAM::MAXLEN * 2];
int get_substr_num(int idx) {
​ int res = (idx == 0 ? 0 : SAM::st[idx].len - SAM::st[SAM::st[idx].link].len);
​ map<char, int>::iterator iter = SAM::st[idx].next.begin();
​ for (; iter != SAM::st[idx].next.end(); iter++) {
​ if (!vis[iter->second]) res += get_substr_num(iter->second);
​ }
​ vis[idx] = true;
​ return res;
}

// 子串查询,判断子串s是否包含在SAM对应的字符串中
bool lookup(char *s) {
​ int idx = 0;
​ for (int i = 0; s[i]; i++) {
​ if (SAM::st[idx].next.count(s[i]))
​ idx = SAM::st[idx].next[s[i]];
​ else return false;
​ }
​ return true;
}

// 路径计数(可用于计算子串数目)
int d[SAM::MAXLEN * 2]; // d[i]表示从节点i出发的路径数目
int path_count(int idx) {
​ // 从根节点出发的路径数目即为子串数目
​ d[idx] = 0;
​ map<char, int>::iterator iter = SAM::st[idx].next.begin();
​ for (; iter != SAM::st[idx].next.end(); iter++) {
​ d[idx] += 1 + path_count(iter->second);
​ }
​ return d[idx];
}

// 查找字典序第k大的子串
bool find_kth_substr(int k, string &res, int idx) {
​ map<char, int>::iterator iter = SAM::st[idx].next.begin();
​ for (; iter != SAM::st[idx].next.end(); iter++) {
​ int t = 1 + d[iter->second];
​ if (k > t) k -= t;
​ else {
​ res += iter->first;
​ if (k > 1) find_kth_substr(k - 1, res, iter->second);
​ return true;
​ }
​ }
​ return false;
}

// 查询子串出现次数
int csort[SAM::MAXLEN + 1]; // 用于对len计数排序
int idsort[2 * SAM::MAXLEN]; // 存储按len排好序的索引id
int query_count(char *s, int idx) {
std::map<char, int>::iterator iter = SAM::st[idx].next.begin();
for (; iter != SAM::st[idx].next.end(); iter++) {
​ if (iter->first == s[0]) {
​ return !s[1] ? SAM::st[iter->second].cnt : query_count(s + 1, iter->second);
​ }
}
return 0;
}

// 查询子串第一次出现位置,返回对应节点索引号,不存在则返回-1
int query_firstpos(char *s, int idx) {
​ std::map<char, int>::iterator iter = SAM::st[idx].next.begin();
​ for (; iter != SAM::st[idx].next.end(); iter++) {
​ if (iter->first == s[0]) {
​ return !s[1] ? iter->second : query_firstpos(s + 1, iter->second);
​ }
​ }
​ return -1;
}

// 查询子串所有出现位置
void query_allpos(int v, int plen, vector<int> &res) {
​ if (!SAM::st[v].is_clone) res.push_back(SAM::st[v].firstpos - plen + 1);
​ for (size_t i = 0; i < SAM::st[v].inv_link.size(); i++) query_allpos(SAM::st[v].inv_link[i], plen, res);
}

// 最长公共子串
string LCS(string t) {
​ int v = 0, l = 0, best = 0, bestpos = 0;
​ for (size_t i = 0; i < t.size(); i++) {
​ while (v && !SAM::st[v].next.count(t[i])) {
​ v = SAM::st[v].link;
​ l = SAM::st[v].len;
​ }
​ if (SAM::st[v].next.count(t[i])) {
​ v = SAM::st[v].next[t[i]];
​ l++;
​ }
​ if (l > best) {
​ best = l;
​ bestpos = i;
​ }
​ }
​ return t.substr(bestpos - best + 1, best);
}

// 多个字符串的最长公共前缀
// 依据SAM的状态节点id在SAM中提取指定长度的子串,存在则返回真
// 并将结果写到lcs中否则返回false
bool get_substr(int idx, int sid, size_t len, string &lcs) {
​ if (idx != 0 && idx == sid) return true;
​ map<char, int>::iterator iter = SAM::st[idx].next.begin();
​ for (; iter != SAM::st[idx].next.end(); iter++) {
​ if (get_substr(iter->second, sid, len, lcs)) {
​ lcs += iter->first;
​ if (idx != 0) return true;
​ else {
​ if (lcs.size() == len) {
​ reverse(lcs.begin(), lcs.end());
​ return true;
​ } else {
​ lcs.clear();
​ }
​ }
​ }
​ }
​ return false;
}

// 时间复杂度为各字符串长总和
void LCS_Mul(int n, string &lcs) {
​ int ans = 0, nid = 0;
​ string s;
​ while (n--) {
​ cin >> s;
​ int tnl = 0, p = 0;
​ for (size_t j = 0; j < s.size(); j++) {
​ if (SAM::st[p].next.count(s[j])) {
​ tnl++;
​ p = SAM::st[p].next[s[j]];
​ } else {
​ while (p && !SAM::st[p].next.count(s[j])) p = SAM::st[p].link;
​ if (SAM::st[p].next.count(s[j])) {
​ tnl = SAM::st[p].len + 1;
​ p = SAM::st[p].next[s[j]];
​ } else {
​ tnl = 0;
​ }
​ }
​ SAM::st[p].nl = max(SAM::st[p].nl, tnl);
​ }
​ for (int j = SAM::sz - 1; j > 0; j--) {
​ p = idsort[j];
​ if (SAM::st[p].nl < SAM::st[p].ml) SAM::st[p].ml = SAM::st[p].nl;
​ if (SAM::st[p].link && SAM::st[SAM::st[p].link].nl < SAM::st[p].nl)
​ SAM::st[SAM::st[p].link].nl = SAM::st[p].nl;
​ SAM::st[p].nl = 0;
​ }
​ }
​ for (int i = 1; i < SAM::sz; i++) {
​ if (SAM::st[i].ml > ans) {
​ ans = SAM::st[i].ml;
​ nid = i;
​ }
​ }
​ cout << ans << endl;
​ get_substr(0, nid, ans, lcs);
}


char instr[SAM::MAXLEN];
int main() {
​ cin >> instr;
​ SAM::sam_init();
​ for (int i = 0; instr[i]; i++) SAM::sam_extend(instr[i]);
​ cout << "SAM节点数: " << SAM::sz << endl; // SAM节点数

// 获取子串总数
memset(vis, false, sizeof(vis));
cout << get_substr_num(0) << " " << path_count(0) << endl;

// 查询串s是否在母串中
char s[100];
cout << "查询子串: ";
cin >> s;
if (lookup(s)) cout << "exists!" << endl;
else cout << "not exists!" << endl;

// 查找字典序第k大的子串
string substr_k = "";
cout << "查询第k大子串:";
int k;
cin >> k;
if (find_kth_substr(k, substr_k, 0))
cout << "第" << k << "大子串: " << substr_k << endl;
else
cout << "不存在," << k << "已超过" << instr << "子串总数" << endl;

// 查询子串出现次数
cout << "查询子串出现次数: ";
cin >> s;
// 对len计数排序
int m = 0;
memset(csort, 0, sizeof(csort));
for (int i = 1; i < SAM::sz; i++) {
csort[SAM::st[i].len]++;
m = max(m, SAM::st[i].len);
}
for (int i = 2; i <= m; i++) csort[i] += csort[i - 1];
for (int i = 1; i < SAM::sz; i++) idsort[csort[SAM::st[i].len]--] = i;
for (int i = SAM::sz - 1; i > 0; i--) {
SAM::st[SAM::st[idsort[i]].link].cnt += SAM::st[idsort[i]].cnt;
}
cout << "子串数目: " << query_count(s, 0) << endl;

// 查询子串第一次出现位置,不存在则返回-1,否则返回节点索引对应到第一次出现结束位置
cout << "查询子串第一次出现位置: ";
cin >> s;
int fid = query_firstpos(s, 0);
if (fid == -1)
cout << "子串不存在!" << endl;
else
cout << SAM::st[fid].firstpos - int(strlen(s)) + 1 << endl;

// 查询子串的所有出现位置(依赖于子串第一次出现节点)
cout << "查询子串所有出现位置: ";
cin >> s;
fid = query_firstpos(s, 0);
if (fid == -1) {
cout << "子串不存在!" << endl;
} else {
for (int i = 1; i < SAM::sz; i++) SAM::st[SAM::st[i].link].inv_link.push_back(i);
vector<int> occur_pos; // 所有出现的起始位置存放在occur_pos中
query_allpos(fid, strlen(s), occur_pos);
for (size_t i = 0; i < occur_pos.size(); i++) cout << occur_pos[i] << " ";
cout << endl;
}

// 两个字符串的最长公共子串
cout << "两个字符串的最长公共子串: ";
string s1;
cin >> s1;
cout << LCS(s1) << endl;

// 多个字符串的最长公共子串
cout << "多个字符串的最长公共子串: ";
int n;
cin >> n;
// 对len计数排序
int lens = 0; // len类别数
memset(csort, 0, sizeof(csort));
for (int i = 1; i < SAM::sz; i++) {
csort[SAM::st[i].len]++;
lens = max(lens, SAM::st[i].len);
}
for (int i = 2; i <= lens; i++) csort[i] += csort[i - 1];
for (int i = 1; i < SAM::sz; i++) idsort[csort[SAM::st[i].len]--] = i;
string lcs;
LCS_Mul(n, lcs);
cout << "多个子串的最长公共子串为" << lcs << ", 长度" << lcs.size() << endl;
return 0;
}

数组(array)版本

展开代码
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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
#include <bits/stdc++.h>
using namespace std;

// 后缀自动机,可以在O(n)复杂度内创建一个字符串的SAM
// 最多2*n - 1个状态节点(n >= 2),最多3*n - 4个转换(n >= 3)
namespace SAM {
struct state {
// SAM必须变量
int len, link; // len表示状态包含的后缀字符串中的最长串长度,link为后缀链接len(link[v]) + 1 = minlen(v)
int trans[26];

// 用于求子串出现次数的变量
int cnt; // 子串出现次数

// firstpos用于求子串第一次出现位置,配合is_clone和inv_link可用于求子串所有出现位置
int firstpos; // 子串第一次出现位置
bool is_clone; // 标志节点是否是复制的
vector<int> inv_link; // 节点的后缀引用列表,该节点对应串是列表内节点的后缀

// ml和nl用于求多字符串的最长公共子串
int nl; // 当前字符串与状态节点对应最长后缀的公共子串长
int ml; // 所有字符串与状态节点对应最长后缀的公共子串长中的最小长度
}; // 状态节点,一个状态对应一个endpoint集合

const int MAXLEN = 100000; // 字符串长度限制
state st[MAXLEN * 2];
char alpha[26]; // 小写字母表,限制输入字符串仅含小写字母
int sz, last;

// 初始化
void sam_init() {
for (int i = 0; i < sz; i++) {
memset(st[i].trans, 0 , sizeof(st[i].trans));
st[i].inv_link.clear();
}
for (int i = 0; i < 26; i++) {
alpha[i] = i + 'a';
}
st[0].len = 0;
st[0].link = -1;
sz = 1;
last = 0;
}

// 向sam中插入字符
void sam_extend(char ch) {
int c = ch - 'a';
int cur = sz++;
st[cur].len = st[last].len + 1;
st[cur].ml = st[cur].len;
st[cur].nl = 0;
st[cur].firstpos = st[cur].len - 1;
st[cur].is_clone = false;
st[cur].cnt = 1;
int p = last;
while (p != -1 && !st[p].trans[c]) {
st[p].trans[c] = cur;
p = st[p].link;
}
if (p == -1) {
st[cur].link = 0;
} else {
int q = st[p].trans[c];
if (st[p].len + 1 == st[q].len) {
st[cur].link = q;
} else {
int clone = sz++;
st[clone].len = st[p].len + 1;
st[clone].ml = st[clone].len;
st[clone].nl = 0;
st[clone].is_clone = true;
st[clone].firstpos = st[q].firstpos;
st[clone].cnt = 0;
memcpy(st[clone].trans, st[q].trans, sizeof(st[q].trans));
st[clone].link = st[q].link;
while (p != -1 && st[p].trans[c] == q) {
st[p].trans[c] = clone;
p = st[p].link;
}
st[q].link = st[cur].link = clone;
}
}
last = cur;
}
}

// 获取子串数目
bool vis[SAM::MAXLEN * 2];
int get_substr_num(int idx) {
int res = (idx == 0 ? 0 : SAM::st[idx].len - SAM::st[SAM::st[idx].link].len);
for (int i = 0; i < 26; i++) {
if (!SAM::st[idx].trans[i]) continue;
if (!vis[SAM::st[idx].trans[i]]) res += get_substr_num(SAM::st[idx].trans[i]);
}
vis[idx] = true;
return res;
}

// 子串查询,判断子串s是否包含在SAM对应的字符串中
bool lookup(char *s) {
int idx = 0;
for (int i = 0; s[i]; i++) {
if (SAM::st[idx].trans[s[i] - 'a'])
idx = SAM::st[idx].trans[s[i] - 'a'];
else return false;
}
return true;
}

// 路径计数(可用于计算子串数目)
int d[SAM::MAXLEN * 2]; // d[i]表示从节点i出发的路径数目
int path_count(int idx) {
// 从根节点出发的路径数目即为子串数目
d[idx] = 0;
for (int i = 0; i < 26; i++) {
if (SAM::st[idx].trans[i])
d[idx] += 1 + path_count(SAM::st[idx].trans[i]);
}
return d[idx];
}

// 查找字典序第k大的子串
bool find_kth_substr(int k, string &res, int idx) {
for (int i = 0; i < 26; i++) {
if (!SAM::st[idx].trans[i]) continue;
int t = 1 + d[SAM::st[idx].trans[i]];
if (k > t) k -= t;
else {
res += SAM::alpha[i];
if (k > 1) find_kth_substr(k - 1, res, SAM::st[idx].trans[i]);
return true;
}
}
return false;
}

// 查询子串出现次数
int csort[SAM::MAXLEN + 1]; // 用于对len计数排序
int idsort[2 * SAM::MAXLEN]; // 存储按len排好序的索引id
int query_count(char *s, int idx) {
for (int i = 0; i < 26; i++) {
if (!SAM::st[idx].trans[i]) continue;
if (SAM::alpha[i] == s[0]) {
return !s[1] ? SAM::st[SAM::st[idx].trans[i]].cnt : query_count(s + 1, SAM::st[idx].trans[i]);
}
}
return 0;
}

// 查询子串第一次出现位置,返回对应节点索引号,不存在则返回-1
int query_firstpos(char *s, int idx) {
for (int i = 0; i < 26; i++) {
if (!SAM::st[idx].trans[i]) continue;
if (SAM::alpha[i] == s[0]) {
return !s[1] ? SAM::st[idx].trans[i] : query_firstpos(s + 1, SAM::st[idx].trans[i]);
}
}
return -1;
}

// 查询子串所有出现位置
void query_allpos(int v, int plen, vector<int> &res) {
if (!SAM::st[v].is_clone) res.push_back(SAM::st[v].firstpos - plen + 1);
for (size_t i = 0; i < SAM::st[v].inv_link.size(); i++) query_allpos(SAM::st[v].inv_link[i], plen, res);
}

// 最长公共子串
string LCS(string t) {
int v = 0, l = 0, best = 0, bestpos = 0;
for (size_t i = 0; i < t.size(); i++) {
while (v && !SAM::st[v].trans[t[i] - 'a']) {
v = SAM::st[v].link;
l = SAM::st[v].len;
}
if (SAM::st[v].trans[t[i] - 'a']) {
v = SAM::st[v].trans[t[i] - 'a'];
l++;
}
if (l > best) {
best = l;
bestpos = i;
}
}
return t.substr(bestpos - best + 1, best);
}

// 多个字符串的最长公共前缀
// 依据SAM的状态节点id在SAM中提取指定长度的子串,存在则返回真
// 并将结果写到lcs中否则返回false
bool get_substr(int idx, int sid, size_t len, string &lcs) {
if (idx != 0 && idx == sid) return true;
for (int i = 0; i < 26; i++) {
if (!SAM::st[idx].trans[i]) continue;
if (get_substr(SAM::st[idx].trans[i], sid, len, lcs)) {
lcs += SAM::alpha[i];
if (idx != 0) return true;
else {
if (lcs.size() == len) {
reverse(lcs.begin(), lcs.end());
return true;
} else {
lcs.clear();
}
}
}
}
return false;
}

// 时间复杂度为各字符串长总和
void LCS_Mul(int n, string &lcs) {
int ans = 0, nid = 0;
string s;
while (n--) {
cin >> s;
int tnl = 0, p = 0;
for (size_t j = 0; j < s.size(); j++) {
if (SAM::st[p].trans[s[j] - 'a']) {
tnl++;
p = SAM::st[p].trans[s[j] - 'a'];
} else {
while (p && !SAM::st[p].trans[s[j] - 'a']) p = SAM::st[p].link;
if (SAM::st[p].trans[s[j] - 'a']) {
tnl = SAM::st[p].len + 1;
p = SAM::st[p].trans[s[j] - 'a'];
} else {
tnl = 0;
}
}
SAM::st[p].nl = max(SAM::st[p].nl, tnl);
}
for (int j = SAM::sz - 1; j > 0; j--) {
p = idsort[j];
if (SAM::st[p].nl < SAM::st[p].ml) SAM::st[p].ml = SAM::st[p].nl;
if (SAM::st[p].link && SAM::st[SAM::st[p].link].nl < SAM::st[p].nl)
SAM::st[SAM::st[p].link].nl = SAM::st[p].nl;
SAM::st[p].nl = 0;
}
}
for (int i = 1; i < SAM::sz; i++) {
if (SAM::st[i].ml > ans) {
ans = SAM::st[i].ml;
nid = i;
}
}
cout << ans << endl;
get_substr(0, nid, ans, lcs);
}


char instr[SAM::MAXLEN];
int main() {
cin >> instr; // 仅限小写字母
SAM::sam_init();
for (int i = 0; instr[i]; i++) SAM::sam_extend(instr[i]);
cout << "SAM节点数: " << SAM::sz << endl; // SAM节点数

// 获取子串总数
memset(vis, false, sizeof(vis));
cout << get_substr_num(0) << " " << path_count(0) << endl;

// 查询串s是否在母串中
char s[100];
cout << "查询子串: ";
cin >> s;
if (lookup(s)) cout << "exists!" << endl;
else cout << "not exists!" << endl;

// 查找字典序第k大的子串
string substr_k = "";
cout << "查询第k大子串:";
int k;
cin >> k;
if (find_kth_substr(k, substr_k, 0))
cout << "第" << k << "大子串: " << substr_k << endl;
else
cout << "不存在," << k << "已超过" << instr << "子串总数" << endl;

// 查询子串出现次数
cout << "查询子串出现次数: ";
cin >> s;
// 对len计数排序
int m = 0;
memset(csort, 0, sizeof(csort));
for (int i = 1; i < SAM::sz; i++) {
csort[SAM::st[i].len]++;
m = max(m, SAM::st[i].len);
}
for (int i = 2; i <= m; i++) csort[i] += csort[i - 1];
for (int i = 1; i < SAM::sz; i++) idsort[csort[SAM::st[i].len]--] = i;
for (int i = SAM::sz - 1; i > 0; i--) {
SAM::st[SAM::st[idsort[i]].link].cnt += SAM::st[idsort[i]].cnt;
}
cout << "子串数目: " << query_count(s, 0) << endl;

// 查询子串第一次出现位置,不存在则返回-1,否则返回节点索引对应到第一次出现结束位置
cout << "查询子串第一次出现位置: ";
cin >> s;
int fid = query_firstpos(s, 0);
if (fid == -1)
cout << "子串不存在!" << endl;
else
cout << SAM::st[fid].firstpos - int(strlen(s)) + 1 << endl;

// 查询子串的所有出现位置(依赖于子串第一次出现节点)
cout << "查询子串所有出现位置: ";
cin >> s;
fid = query_firstpos(s, 0);
if (fid == -1) {
cout << "子串不存在!" << endl;
} else {
for (int i = 1; i < SAM::sz; i++) SAM::st[SAM::st[i].link].inv_link.push_back(i);
vector<int> occur_pos; // 所有出现的起始位置存放在occur_pos中
query_allpos(fid, strlen(s), occur_pos);
for (size_t i = 0; i < occur_pos.size(); i++) cout << occur_pos[i] << " ";
cout << endl;
}

// 两个字符串的最长公共子串
cout << "两个字符串的最长公共子串: ";
string s1;
cin >> s1;
cout << LCS(s1) << endl;

// 多个字符串的最长公共子串
cout << "多个字符串的最长公共子串: ";
int n;
cin >> n;
// 对len计数排序
int lens = 0; // len类别数
memset(csort, 0, sizeof(csort));
for (int i = 1; i < SAM::sz; i++) {
csort[SAM::st[i].len]++;
lens = max(lens, SAM::st[i].len);
}
for (int i = 2; i <= lens; i++) csort[i] += csort[i - 1];
for (int i = 1; i < SAM::sz; i++) idsort[csort[SAM::st[i].len]--] = i;
string lcs;
LCS_Mul(n, lcs);
cout << "多个子串的最长公共子串为" << lcs << ", 长度" << lcs.size() << endl;
return 0;
}
---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址:
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢!

Word2Vec原理详解

发表于 2019-08-14 | 分类于 Research , Paper | 0 comments
阅读次数:
  |   字数统计: 1.3k(字)   |   阅读时长: 6(分)
paper
  • Efficient Estimation of Word Representations in
    Vector Space

  • Distributed Representations of Words and Phrases
    and their Compositionality

简介

Word2Vec分为CBOW模型(Continuous Bag-of-Words Model)和Skip-gram(Continuous Skip-gram Model)模型,如下图所示

当训练文本集规模较小时适合使用CBOW模型,当训练文本集规模较大时适合使用Skip-gram模型

训练方法

Hierarchical Softmax

方法一:基于Huffman树的Hierarchical Softmax

Negative Sampling

方法二:作者Mikolov提出的Negative Sampling(NEG),该方法是Noise Contrastive Estimation(NCE)的一个简化版本

CBOW模型

在CBOW模型中,已知词$w$的上下文$Context(w)$,需要预测$w$,因此对于给定的$Context(w)$,词$w$就是一个正样本,其它词就是负样本了。假定现在已经选好了一个关于$w$的负样本子集$NEG(w) \ne \emptyset$,且对$\forall\ \widetilde{w} \in \mathcal{D}$,定义

令$\sigma(x)=\frac{1}{1+\mathrm{e}^{x}}$,语料库$C$,概率生成函数如下

其中

可将$p(u|Context(w))$写成整体表达式如下

从而$g(w)$可简写为

最终的目标函数对G取对数求最大似然。从形式上看,最大化$g(w)$相当于最大化,同时最小化所有的,使得增大正样本概率的同时降低负样本的概率,这正是我们所希望看到的结果

为方便梯度推导,将(1.6)式花括号中的内容简记为$\mathcal{L}(w,u)$,即

利用梯度上升法对$\mathcal{L}$进行优化求最大值,$\mathcal{L}(w,u)$分别对$x_{w}和\theta^{u}$求偏导

于是$\theta^{u}$的更新公式可写成

接下来考虑$\mathcal{L}(w,u)$对的偏导,由于式(1.7)中和的对称性可得

于是,利用式(1.10)的结果可得$v(\widetilde{w}),\widetilde{w} \in Context(w)$的更新公式为

伪代码如下

Skip-gram模型

结合SKip-gram和CBOW模型的区别,将优化目标函数由原来的$G=\prod_{w \in C}g(w)$改写为

这里,$\prod_{u \in Context(w)}g(u)$表示对于一个给定的样本$(w,Context(w))$,我们希望最大化的量,$g(u)$类似于上一节$g(w)$,定义为

其中$NEG(u)$表示处理词u时生成的负样本子集,条件概率

(2.3)式写成整体表达式为

同样,我们取G的对数,最终的目标函数为

与CBOW模型一样利用梯度上升法对$\mathcal{L}$进行优化求最大值,将式(2.5)中花括号中的内容简记为$\mathcal{L}(w,u,z)$,然后将$\mathcal{L}(w,u,z)$分别对$v(w)$和$\theta^{z}$求偏导,具体求解过程省略,可参考文末链接参考1中的5.2节

参考

  1. * word2vec 中的数学原理详解
  2. Word2Vec教程 - Skip-Gram模型
  3. Word2Vec教程(2)- Negative Sampling
  4. Hierarchical Softmax(视频)
  5. Noise Contrastive Estimation(paper)
---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址:
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢!

三色旗问题

发表于 2019-08-14 | 分类于 Algorithm | 0 comments
阅读次数:
  |   字数统计: 411(字)   |   阅读时长: 1(分)

问题描述

假设有一条绳子,上面有红、白、蓝三种颜色的旗子,起初绳子上的旗子的颜色并没有顺序,你希望将之分类并排列为蓝、白、红的顺序,要如何移动次数最少,注意你只能在绳子上进行这个动作,并且一次只能调换两个旗子。

解决思路

快排的思想,分别设三个标记wflag,rflag,bflag,初始bflag和wflag置首,rflag置末,然后分以下几个步骤解决:

  1. rflag向左移动到第一个非红色旗处停下,bflag向前移动到第一个非蓝色旗停下,wflag向前移动超越bflag并在非白色旗处停下
  2. wflag处如果为红旗则与rflag处的旗交换位置,如果为蓝旗则与bflag处的旗交换位置,交换位置后移动更新三个标志位的位置
  3. 重复1~2步直到wflag > rflag

具体实现

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
#pragma GCC optimize("O2")
#include <bits/stdc++.h>
using namespace std;

//char color[] = {'r','w','b','w','w','b','r','b','w','r','\0'};
char color[] = {'r','r','r','r','r','b', '\0'};
int n = strlen(color);
int wflag = 0, bflag = 0, rflag = n - 1;

void move() {
while (bflag < n && color[bflag] == 'b') bflag++;
while (wflag < n && (wflag < bflag || color[wflag] == 'w')) wflag++;
while (rflag >= 0 && color[rflag] == 'r') rflag--;
}

int main()
{
int res = 0;
cout << color << " " << n << endl;
move();
while (wflag < n && wflag <= rflag) {
if (color[wflag] == 'r') {
swap(color[wflag], color[rflag]);
res++;
} else if (color[wflag] == 'b') {
swap(color[wflag], color[bflag]);
res++;
}
move();
}
cout << color << endl;
cout << "至少交换 " << res << " 次" << endl;
return 0;
}
---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址:
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢!

皮尔逊相关系数绝对值小于1的证明

发表于 2019-08-09 | 分类于 Math | 0 comments
阅读次数:
  |   字数统计: 517(字)   |   阅读时长: 2(分)

皮尔逊相关系数的介绍可参考维基百科:皮尔逊积矩相关系数

在统计学中,皮尔逊积矩相关系数(Pearson product-moment correlation coefficient,又称作 PPMCC或PCCs, 文章中常用r或Pearson’s r表示)用于度量两个变量X和Y之间的相关程度(线性相关)

两个变量的总体相关系数定义为两个变量之间的协方差与标准差的商

估算样本的协方差和标准差,可得到样本相关系数

r 亦可由样本点()的标准分数均值估算,得到与上式等价的表达式

其中分别是样本的标准分数、样本均值和样本标准差

总体和样本皮尔逊相关系数都小于或者等于1,可以由柯西-施瓦茨不等式证明,关于柯西-施瓦茨不等式的定义与证明具体可参考维基百科:柯西-施瓦茨不等式,下面将直接利用柯西-施瓦茨不等式来证明皮尔逊相关系数的绝对值小于1

证明:

要证明式(1.2)中的r绝对值小于等于1即证明小于等于1,设式(1.2)中的,代入有

又由柯西-施瓦茨不等式知

由此可知式(1.4)中的,即,证毕!

---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址:
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢!

各类概率分布函数tutorial

发表于 2019-08-08 | 分类于 Math | 0 comments
阅读次数:
  |   字数统计: 1.1k(字)   |   阅读时长: 4(分)

泊松分布、指数分布与伽马分布简介及其关系

泊松分布:

指数分布:

伽马分布: 先看伽马函数,伽马函数又称第二类欧拉积分,函数式为,对伽马函数做一个变形即可得到

式(1.1)去除积分即为伽马分布的概率密度函数

令,代入积分公式(1.1)中则可得到新的概率密度函数

其中为shape parameter,主要决定了分布曲线的形状,为scale parameter,主要决定曲线有多陡,越大则曲线越陡。

伽马分布是为解决阶乘数列的插值而产生的,伽马函数使得任意正实数都有阶乘,伽马函数有几条特性如下:

其中,

物理解释:

泊松分布: 单位时间t内独立事件发生次数n的概率分布

指数分布: 独立事件的时间间隔的概率分布

伽马分布: 个独立事件的时间间隔的概率分布

参考链接:
  • 伽马分布,指数分布,泊松分布的关系

    这篇文章先讲述泊松分布,然后由泊松分布引出指数分布,说明了泊松分布与指数分布的关系,最后说明了伽马分布与指数分布的关系.

  • LDA-math-神奇的Gamma函数(3)

    这篇文章由二项分布和Beta分布的关系式,借用二项分布的极限是泊松分布的性质转换得到泊松分布与伽马分布的关系式,以及泊松分布的概率累积函数与伽马分布的概率累积函数互补的关系。

Beta分布简介

Beta分布: 首先来看Beta函数,Beta函数又称贝塔函数或者第一类欧拉积分,它的函数式如下

对Beta函数作一个变形即可得到

去除式(2.1)的积分即可得到Beta分布的概率密度函数

参考链接:
  • LDA-math-认识Beta/Dirichlet分布(1)

    这篇文章由掷骰子游戏从物理上解释了Beta分布的由来

  • LDA-math-认识Beta/Dirichlet分布(2)

    这篇文章在上一篇文章的基础上通过更新掷骰子的游戏规则引出了Beta分布与二项分布的共轭分布关系,共轭的意思就是,数据符合二项分布的时候,参数的先验分布和后验分布都能保持Beta 分布的形式。我们知道贝叶斯的参数估计过程是

    对应到beta分布和二项分布的共轭

Dirichlet分布简介

Multinomial分布: 中文名多项式分布

Dirichlet分布: 中文名狄利克雷分布

参考链接:
  • LDA-math-认识Beta/Dirichlet分布(3)

    这篇文章由掷骰子的游戏引出Dirichlet分布,并且通过改变游戏规则引出Dirichlet分布与Multinomial分布的关系,但是在推导共轭关系时定义有些许错误要注意

  • LDA基础知识—-Dirichlet 分布

    这是优酷上讲解LDA的Dirichlet分布的视频

---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址:
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢!

多项式相关算法集成

发表于 2019-07-24 | 分类于 Algorithm | 0 comments
阅读次数:
  |   字数统计: 2.1k(字)   |   阅读时长: 11(分)

原理

快速数论变换(FNT)

原理同FFT一样,用与弧度具有同样性质的一组值替代弧度进行计算,具体性质如下:

FFT 中能在 O(nlog⁡n) 时间内变换用到了单位根 的什么性质:

  1. $\omega_n^n = 1$

  2. , , ⋯, 是互不相同的,这样带入计算出来的点值才可以用来还原出系数

  3. , ,这使得在按照指数奇偶分类之后能够把带入的值也减半使得问题规模能够减半

        这点保证了能够使用相同的方法进行逆变换得到系数表示

参考:

  1. 从多项式乘法到快速傅里叶变换

  2. 快速数论变换 (NTT)

  3. FNT用到的各种素数

多项式逆元

参考:

  1. 多项式求逆元

多项式除法及求模

参考:

  1. 多项式除法及求模

多项式多点求值与拉格朗日插值

参考:

  1. 多项式的多点求值与快速插值
  2. 多项式多点求值和插值
  3. 维基百科:拉格朗日插值法

实现

展开代码
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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
/*
* 多项式相关算法模板集合
*/
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
// 2281701377是一个原根为3的指数,平方刚好不会爆long long
// 另外一个常用的998244353 = 119ll * (1 << 23) + 1
// 1004535809=479⋅221+1 加起来刚好不会爆 int 也不错
const ll mod_v = 17ll * (1 << 27) + 1;
const int MAXL = 18, N = 1 << MAXL, M = 1 << MAXL; // MAXL最大取mov_v的2次幂
struct polynomial {
​ ll coef[N]; // size设为多项式系数个数的两倍
​ ll a[N], b[N], d[N], r[N];
​ ll xcoor[N], ycoor[N], v[N], mcoef[N]; // size设为点个数的两倍
​ vector<vector<ll> > poly_divisor;
​
​ static ll mod_pow(ll x,ll n,ll m) {
​ ll ans = 1;
​ while(n > 0){
​ if(n & 1)
​ ans = ans*x%m;
​ x = x*x%m;
​ n >>= 1;
​ }
​ return ans;
​ }
​
​ struct FastNumberTheoreticTransform {
​ ll omega[N], omegaInverse[N];
​ int range;
​
​ // 初始化频率
​ void init (const int& n) {
​ range = n;
​ ll base = mod_pow(3, (mod_v - 1) / n, mod_v);
​ ll inv_base = mod_pow(base, mod_v - 2, mod_v);
​ omega[0] = omegaInverse[0] = 1;
​ for (int i = 1; i < n; ++i) {
​ omega[i] = omega[i - 1] * base % mod_v;
​ omegaInverse[i] = omegaInverse[i - 1] * inv_base % mod_v;
​ }
​ }
​
​ // Cooley-Tukey算法:O(n*logn)
​ void transform (ll *a, const ll *omega, const int &n) {
​ for (int i = 0, j = 0; i < n; ++i) {
​ if (i > j) std::swap (a[i], a[j]);
​ for(int l = n >> 1; ( j ^= l ) < l; l >>= 1);
​ }
​
​ for (int l = 2; l <= n; l <<= 1) {
​ int m = l / 2;
​ for (ll *p = a; p != a + n; p += l) {
​ for (int i = 0; i < m; ++i) {
​ ll t = omega[range / l * i] * p[m + i] % mod_v;
​ p[m + i] = (p[i] - t + mod_v) % mod_v;
​ p[i] = (p[i] + t) % mod_v;
​ }
​ }
​ }
​ }
​

// 时域转频域
void dft (ll *a, const int& n) {
transform(a, omega, n);
}

// 频域转时域
void idft (ll *a, const int& n) {
transform(a, omegaInverse, n);
for (int i = 0; i < n; ++i) a[i] = a[i] * mod_pow(n, mod_v - 2, mod_v) % mod_v;
}
} fnt ;

// 与模比较转换值
ll mod_trans(ll v) {
return abs(v) <= mod_v / 2 ? v : (v < 0 ? v + mod_v : v - mod_v);
}

// 二分求\prod{x-xi}的所有二分子多项式的系数
void binary_subpoly(int l, int r, int idx) {
if (l == r - 1) {
poly_divisor[idx].push_back(-xcoor[l]);
poly_divisor[idx].push_back(1);
} else {
int lidx = (idx << 1) + 1, ridx = lidx + 1;
binary_subpoly(l, (l + r) / 2, lidx);
binary_subpoly((l + r) / 2, r, ridx);
int t = poly_divisor[lidx].size() + poly_divisor[ridx].size() - 1;
int p = 1;
while(p < t) p <<= 1;
copy(poly_divisor[lidx].begin(), poly_divisor[lidx].end(), a);
fill(a + poly_divisor[lidx].size(), a + p, 0);
copy(poly_divisor[ridx].begin(), poly_divisor[ridx].end(), b);
fill(b + poly_divisor[ridx].size(), b + p, 0);
fnt.dft(a, p);
fnt.dft(b, p);
for (int i = 0; i < p; i++) a[i] *= b[i];
fnt.idft(a, p);
for (int i = 0; i < t; i++) poly_divisor[idx].push_back(mod_trans(a[i]));
}
}

// 模x^deg,a为要求逆元的多项式系数,结果存放在b[0~deg]中
// T(deg) = T(deg/2) + deg*log(deg),复杂度O(deg*log(deg))
void polynomial_inverse(int deg, ll* a, ll* b, ll* tmp) {
if(deg == 1) {
b[0] = mod_pow(a[0], mod_v - 2, mod_v);
} else {
polynomial_inverse((deg + 1) >> 1, a, b, tmp);
int p = 1;
while(p < (deg << 1) - 1) p <<= 1;
copy(a, a + deg, tmp);
fill(tmp + deg, tmp + p, 0);
fill(b + ((deg + 1) >> 1), b + p, 0);
//fnt.init(p);
fnt.dft(tmp, p);
fnt.dft(b, p);
for(int i = 0; i != p; ++i) {
b[i] = (2 - tmp[i] * b[i] % mod_v) * b[i] % mod_v;
if(b[i] < 0) b[i] += mod_v;
}
fnt.idft(b, p);
fill(b + deg, b + p, 0);
}
}

// A = D*B + R,A为n项n-1次幂,B为m项m-1次幂,D为n-m+1项n-m次幂,R为m-1项m-2次幂
// 要求a,b中系数以低次到多次顺序排列
// n >= m,复杂度O(n*logn);n < m,复杂度O(n)
int polynomial_division(int n, int m, ll *A, ll *B, ll *D, ll *R) {
if (n < m) {
copy(A, A + n, R);
return n;
} else {
static ll A0[N], B0[N], tmp[N]; //数组太大会爆栈,添加到全局区

int p = 1, t = n - m + 1;
while(p < (t << 1) - 1) p <<= 1;

fill(A0, A0 + p, 0);
reverse_copy(B, B + m, A0);
polynomial_inverse(t, A0, B0, tmp);
fill(B0 + t, B0 + p, 0);
fnt.dft(B0, p);

reverse_copy(A, A + n, A0);
fill(A0 + t, A0 + p, 0);
fnt.dft(A0, p);

for(int i = 0; i != p; ++i)
A0[i] = A0[i] * B0[i] % mod_v;
fnt.idft(A0, p);
reverse(A0, A0 + t);
copy(A0, A0 + t, D);

for(p = 1; p < n; p <<= 1);
fill(A0 + t, A0 + p, 0);
fnt.dft(A0, p);
copy(B, B + m, B0);
fill(B0 + m, B0 + p, 0);
fnt.dft(B0, p);
for(int i = 0; i != p; ++i)
A0[i] = A0[i] * B0[i] % mod_v;
fnt.idft(A0, p);
for(int i = 0; i != m - 1; ++i)
R[i] = ((A[i] - A0[i]) % mod_v + mod_v) % mod_v;
//fill(R + m - 1, R + p, 0);
return m - 1;
}
}

// 多项式的点值计算
// l和r为存储要求的点数组的左右边界(左闭右开), idx为除数多项式索引(初始0)
// polycoef为用于计算点的多项式系数,num为其系数个数
// 设多项式项数x=num,点数y=r-l,n=max(x,y),复杂度O(n(logn)^2)
void polynomial_calculator(int l, int r, int idx, int num, ll *polycoef) {
int mid = (l + r) / 2;
int lidx = (idx << 1) + 1, ridx = lidx + 1;
int lsize = poly_divisor[lidx].size(), rsize = poly_divisor[ridx].size();
ll *lmod_poly = new ll[lsize - 1], *rmod_poly = new ll[rsize - 1];
copy(poly_divisor[lidx].begin(), poly_divisor[lidx].end(), a);
copy(poly_divisor[ridx].begin(), poly_divisor[ridx].end(), b);
int lplen = polynomial_division(num, lsize, polycoef, a, d, lmod_poly);
int rplen = polynomial_division(num, rsize, polycoef, b, d, rmod_poly);
if (l == mid - 1) {
v[l] = mod_trans(lmod_poly[0]);
} else {
polynomial_calculator(l, (l + r) / 2, lidx, lplen, lmod_poly);
}
if (r == mid + 1) {
v[(l + r) / 2] = mod_trans(rmod_poly[0]);
} else {
polynomial_calculator((l + r) / 2, r, ridx, rplen, rmod_poly);
}
delete []lmod_poly;
delete []rmod_poly;
}

// 拉格朗日插值:二分+快速数论变换
// l和r为存储要求的点数组的左右边界(左闭右开), idx为由点二分构造出多项式的索引(初始0)
// polycoef为插值得到的多项式结果
// 设点个数为n,复杂度O(n*(logn)^2),结果polycoef为n-1次多项式
void polynomial_interpolate(int l, int r, int idx, ll *polycoef) {
if (l == r - 1) {
polycoef[0] = ycoor[l] * mod_pow(v[l], mod_v - 2, mod_v) % mod_v;
} else {
int mid = (l + r) >> 1;
int lidx = (idx << 1) + 1, ridx = lidx + 1;
int sz = poly_divisor[idx].size() - 1;
int lsize = poly_divisor[lidx].size(), rsize = poly_divisor[ridx].size();
int p = 1;
while (p < sz) p <<= 1;
ll *leftpoly = new ll[p], *rightpoly = new ll[p];
polynomial_interpolate(l, mid, lidx, leftpoly);
polynomial_interpolate(mid, r, ridx, rightpoly);
copy(poly_divisor[lidx].begin(), poly_divisor[lidx].end(), a);
copy(poly_divisor[ridx].begin(), poly_divisor[ridx].end(), b);
fill(leftpoly + lsize - 1, leftpoly + p, 0);
fill(rightpoly + rsize - 1, rightpoly + p, 0);
fill(a + lsize, a + p, 0);
fill(b + rsize, b + p, 0);
fnt.dft(leftpoly, p);
fnt.dft(b, p);
for (int i = 0; i < p; i++) leftpoly[i] = leftpoly[i] * b[i] % mod_v;
fnt.idft(leftpoly, p);
fnt.dft(rightpoly, p);
fnt.dft(a, p);
for (int i = 0; i < p; i++) rightpoly[i] = rightpoly[i] * a[i] % mod_v;
fnt.idft(rightpoly, p);
for (int i = 0; i < sz; i++) polycoef[i] = mod_trans((leftpoly[i] + rightpoly[i]) % mod_v);
delete []leftpoly;
delete []rightpoly;
}
}

// 初始化点个数到二分子多项式个数
// 调用polynomial_calculator和polynomial_interpolate前调用
void init(int vnum) {
int vnum2 = 1;
while (vnum2 < vnum) vnum2 <<= 1;
for (int i = 0; i < 2 * vnum2 - 1; i++) poly_divisor.push_back(vector<ll>());
}
} poly;

应用

多项式快速乘

展开代码 ​```c++ int main() { ios_base::sync_with_stdio(false); poly.fnt.init(1 << MAXL); int n, m; cin >> n >> m; for (int i = 0; i < n; i++) cin >> poly.a[i]; for (int i = 0; i < m; i++) cin >> poly.b[i]; cout << "a*b之后的真实系数:" << endl; for (int i = 0; i < n; i++) for (int j = 0; j < m; j++) poly.r[i + j] += poly.a[i] * poly.b[j]; for (int i = 0; i < n + m - 1; i++) cout << poly.r[i] << " "; cout << endl; cout << "利用fft计算得出的系数:" << endl; int p = 1; while (p < n + m - 1) p <<= 1; // 只要p大于多项式结果中的最大次幂即可 poly.fnt.dft(poly.a, p); poly.fnt.dft(poly.b, p); for (int i = 0; i < p; i++) { poly.d[i] = poly.a[i] * poly.b[i] % mod_v; } poly.fnt.idft(poly.d, p); for (int i = 0; i < n + m - 1; i++) cout << poly.d[i] << " "; return 0; }
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
</details>

<div align="center">
<img src="/images/poly/poly5.png">
</div>

### 多项式逆元

<details>
<summary>展开代码</summary>
```c++
int main() {
ios_base::sync_with_stdio(false);
poly.fnt.init(1 << MAXL);
int n;
cin >> n;
for(int i = 0; i != n; ++i)
cin >> poly.a[i];
int m;
cin >> m; // 输入模x^m
poly.polynomial_inverse(m, poly.a, poly.b, poly.d);
cout << "inverse: " << endl;
for(int i = 0; i != m; ++i)
cout << (poly.b[i] + mod_v) % mod_v << " ";
cout << endl;
cout << "a*b相乘取模验证取模后的前m项系数:" << endl;
memset(poly.d, 0, sizeof(poly.d));
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) {
poly.d[i + j] = (poly.d[i + j] + poly.a[i] * poly.b[j] % mod_v)% mod_v;
}
}
cout << m << endl;
for (int i = 0; i < m; i++) {
cout << poly.d[i] << " ";
}
return 0;
}

输入数据得到输出后的结果如下图

多项式除法以及取模

展开代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
int main() {
ios_base::sync_with_stdio(false);
poly.fnt.init(1 << MAXL);
int n, m;
cin >> n >> m;
for(int i = 0; i != n; ++i)
cin >> poly.a[i]; // 0次幂系数开始输入,缺失的幂系数输入0
for (int i = 0; i < m; i++)
cin >> poly.b[i]; // 0次幂系数开始输入
int rlen = poly.polynomial_division(n, m, poly.a, poly.b, poly.d, poly.r);
cout << "输出商多项式:" << endl;
for (int i = 0; i < n - m + 1; i++)
cout << poly.d[i] << " ";
cout << endl;
cout << "输出余数多项式:" << endl;
for (int i = 0; i < rlen; i++)
cout << poly.r[i] << " ";
return 0;
}

多项式多点求值

展开代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
int main() {
ios_base::sync_with_stdio(false);
poly.fnt.init(1 << MAXL);
int n, vnum;
// 输入要计算的点
cin >> n;
for (int i = 0; i < n; i++) {
cin >> poly.coef[i];
}
cin >> vnum;
for (int i = 0; i < vnum; i++) {
cin >> poly.xcoor[i];
}
// 二分求子多项式
poly.init(vnum);
poly.binary_subpoly(0, vnum / 2, 1);
poly.binary_subpoly(vnum / 2, vnum, 2);
// 将输入点代入插值得到的多项式中进行验证,输出计算得到的结果
cout << "计算结果:" << endl;
poly.polynomial_calculator(0, vnum, 0, n, poly.coef);
for (int i = 0; i < vnum; i++) cout << poly.v[i] << " ";
cout << endl;
}

输入数据后得到结果如下

拉格朗日快速插值计算

展开代码
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
int main() {
// 多项式系数默认低次到高次排列
ios_base::sync_with_stdio(false);
poly.fnt.init(1 << MAXL);
// 多项式插值
int vnum;
// 输入要计算的点
cin >> vnum;
for (int i = 0; i < vnum; i++) {
cin >> poly.xcoor[i] >> poly.ycoor[i];
}
// 二分求子多项式
poly.init(vnum);
poly.binary_subpoly(0, vnum, 0);
for (unsigned int i = 1; i < poly.poly_divisor[0].size(); i++) {
poly.mcoef[i - 1] = poly.poly_divisor[0][i] * i % mod_v;
}
// 遍历i计算所有\sum_{j!=i}{xi-xj}
poly.polynomial_calculator(0, vnum, 0, poly.poly_divisor[0].size() - 1, poly.mcoef);
// 拉格朗日插值计算多项式
poly.polynomial_interpolate(0, vnum, 0, poly.coef);
// 输出插值得到的多项式系数
for (int i = 0; i < vnum; i++) {
cout << poly.coef[i] << " ";
}
cout << endl;
// 将输入点代入插值得到的多项式中进行验证,输出计算得到的结果
poly.polynomial_calculator(0, vnum, 0, vnum, poly.coef);
for (int i = 0; i < vnum; i++) cout << poly.v[i] << " ";
cout << endl;
return 0;
}

输入数据得到的结果如下

---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址:
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢!

快速傅里叶变换原理及其实现

发表于 2019-07-18 | 分类于 Algorithm | 0 comments
阅读次数:
  |   字数统计: 1.1k(字)   |   阅读时长: 6(分)

原理

参考:

如何给文科生解释傅里叶变换

一小时学会快速傅里叶变换(Fast Fourier Transform)

实现

  • 使用C++内置的complex

    展开代码
    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
    #include <bits/stdc++.h>
    using namespace std;

    const double PI = acos(-1);
    const int MAXL = 18, N = 1 << MAXL; // 序列周期大小

    struct FastFourierTransform {
    ​ complex<double> omega[N], omegaInverse[N];
    ​ int range;

    // 初始化频率
    // 可以只进行一次初始化,只要大于可能出现的最高次幂即可,后面的转换就可以复用初始化
    void init (const int& n) {
    range = n;
    for (int i = 0; i < n; ++i) {
    omega[i] = complex<double>(cos(2 * PI / n * i), sin( 2 * PI / n * i )) ;
    omegaInverse[i] = std::conj(omega[i]);
    }
    }

    // Cooley-Tukey算法:O(n*logn)
    void transform (complex<double> *a, const complex<double> *omega, const int &n) {
    for (int i = 0, j = 0; i < n; ++i) {
    if (i > j) std::swap (a[i], a[j]);
    for(int l = n >> 1; ( j ^= l ) < l; l >>= 1);
    }

    for (int l = 2; l <= n; l <<= 1) {
    int m = l / 2;
    for (complex<double> *p = a; p != a + n; p += l) {
    for (int i = 0; i < m; ++i) {
    complex<double> t = omega[range / l * i] * p[m + i];
    p[m + i] = p[i] - t;
    p[i] += t;
    }
    }
    }
    }

    // 时域转频域
    void dft (complex<double> *a, const int& n) {
    transform(a, omega, n);
    }

    // 频域转时域
    void idft (complex<double> *a, const int& n) {
    transform(a, omegaInverse, n);
    for (int i = 0; i < n; ++i) a[i] /= n;
    }
    } fft ;
  • 使用自定义的complex,两者在精度上有细微差别但影响不大

    展开代码
    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
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    #include <bits/stdc++.h>
    using namespace std;

    const double PI = acos(-1);
    const int MAXL = 18, N = 1 << MAXL; // 序列周期大小

    struct Complex {
    double r, i ;
    Complex ( ) { }
    Complex ( double r, double i ) : r ( r ), i ( i ) { }
    inline void real ( const double &x ) { r = x ; }
    inline double real ( ) { return r ; }
    inline void imag(const double &x) { i = x; }
    inline double imag() { return i; }
    inline Complex operator + ( const Complex& rhs ) const {
    return Complex ( r + rhs.r, i + rhs.i ) ;
    }
    inline Complex operator - ( const Complex& rhs ) const {
    return Complex ( r - rhs.r, i - rhs.i ) ;
    }
    inline Complex operator * ( const Complex& rhs ) const {
    return Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ;
    }
    inline void operator /= ( const double& x ) {
    r /= x, i /= x ;
    }
    inline void operator *= ( const Complex& rhs ) {
    *this = Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ;
    }
    inline void operator += ( const Complex& rhs ) {
    r += rhs.r, i += rhs.i ;
    }
    inline Complex conj ( ) {
    return Complex ( r, -i ) ;
    }
    friend ostream& operator << (ostream& out, const Complex& cx) {
    out << "(" << cx.r << ", " << cx.i << ")" << endl;
    return out;
    }
    friend istream& operator >> (istream& in, Complex& cx) {
    in >> cx.r >> cx.i;
    return in;
    }
    } ;


    struct FastFourierTransform {
    Complex omega [N], omegaInverse [N] ;
    int range;

    // 初始化
    void init ( const int& n ) {
    range = n;
    for ( int i = 0 ; i < n ; ++ i ) {
    omega [i] = Complex ( cos ( 2 * PI / n * i), sin ( 2 * PI / n * i ) ) ;
    omegaInverse [i] = omega [i].conj ( ) ;
    }
    }

    // 迭代式二分转换
    void transform ( Complex *a, const int& n, const Complex* omega ) {
    for ( int i = 0, j = 0 ; i < n ; ++ i ) {
    if ( i > j ) std :: swap ( a [i], a [j] ) ;
    for( int l = n >> 1 ; ( j ^= l ) < l ; l >>= 1 ) ;
    }

    for ( int l = 2 ; l <= n ; l <<= 1 ) {
    int m = l / 2;
    for ( Complex *p = a ; p != a + n ; p += l ) {
    for ( int i = 0 ; i < m ; ++ i ) {
    Complex t = omega [range / l * i] * p [m + i] ;
    p [m + i] = p [i] - t ;
    p [i] += t ;
    }
    }
    }
    }

    // 时域转频域
    void dft ( Complex *a, const int& n ) {
    transform ( a, n, omega ) ;
    }

    // 频域转时域
    void idft ( Complex *a, const int& n ) {
    transform ( a, n, omegaInverse ) ;
    for ( int i = 0 ; i < n ; ++ i ) a [i] /= n ;
    }
    } fft ;

应用

  • 加速多项式计算O(nlogn)

    展开代码
    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
    int main() {
    // 应用到加速多项式乘法
    complex<double> a[500], b[500];
    int n, m;
    cin >> n >> m;
    for (int i = 0; i < n; i++) cin >> a[i];
    for (int i = 0; i < m; i++) cin >> b[i];

    // a*b之后的真实系数
    int r[20] = {0};
    for (int i = 0; i < n; i++)
    for (int j = 0; j < m; j++)
    r[i + j] += a[i].real() * b[j].real();
    cout << "真实系数:";
    for (int i = 0; i < n + m - 1; i++)
    cout << r[i] << " ";
    cout << endl;

    // 利用fft计算得出的系数
    int p = 1;
    while (p < n + m - 1) p <<= 1; // 只要p大于多项式结果中的最大次幂即可
    fft.init(p);
    fft.dft(a, p);
    fft.dft(b, p);
    complex<double> c[500];
    for (int i = 0; i < p; i++) {
    c[i] = a[i] * b[i];
    }
    fft.idft(c, p);
    // 考虑到精度问题取实部+0.5向下取整
    cout << "fft计算得出的系数:";
    for (int i = 0; i < n + m - 1; i++) cout << int(c[i].real() + 0.5) << " ";
    return 0;
    }

    输入一组数据得到输出结果如下

---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址:
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢!

动态规划题集

发表于 2019-07-18 | 分类于 Algorithm , ACM Solution | 0 comments
阅读次数:
  |   字数统计: 553(字)   |   阅读时长: 2(分)

01背包问题

HDU-3466

题目: HDU-3346. Proud Merchants

分析:  首先对于任意两个物品i和j,如果我们只选他们中的一个或全不选,我们所需要的初始金钱是相同的.现在如果我们对于它们两个都选,初始金钱就跟选择的顺序有关了.比如:

3 5 6

5 10 5

两个物品,如果先1后2,初始需要13,如果先2后1,初始需要10.

对于任意两个物品i和j,如果我们先i后j初始:

金钱= Q[i]+Q[j]-(Q[i]-p[i]) = p[i]+Q[j]

同理可得先j后i,初始:

金钱= Q[j]+Q[i]-(Q[j]-p[j]) = p[j]+Q[i]

所以对于i和j,假设我们手上有M金钱,我们肯定先以(Q-P)值从大到小开始选择,这样我才能保证对于同一种选择方案(某些选,某些不选),我们所需要的初始金钱最小.或者换一种说法,给定一定的金额,如果按照(Q-P)值从大到小的顺序进行选择,可以获得尽可能多的物品.

solution:  01背包+贪心,关键在于找到要进行贪心的值

implementation:

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
/*******************************************************************************
Author: liujian
Email: brooksj@foxmail.com
Date: 2019-07-18 10:51
File: hdu_3466.cpp
Last modified: 2019-07-18 10:51
*******************************************************************************/

#pragma GCC optimize("O2")
#include <bits/stdc++.h>
using namespace std;
const int maxn = 5 * 1e2 + 1;
const int maxm = 5 * 1e3 + 1;

struct item {
int p, q, v;
item(int _p = 0, int _q = 0, int _v = 0):p(_p),q(_q),v(_v){}
friend bool operator < (const item &a, const item &b) {
return (a.q - a.p) < (b.q - b.p);
}
};

int dp[maxm], sum[maxn];
item items[maxn];

int main() {
ios_base::sync_with_stdio(false);
int n, m;
while (cin >> n >> m) {
sum[0] = 0;
for (int i = 1; i <= n; i++) {
cin >> items[i].p >> items[i].q >> items[i].v;
}
sort(items + 1, items + n + 1);
for (int i = 1; i <= n; i++) {
sum[i] = sum[i - 1] + items[i].p;
}
for (int i = 0; i <= m; i++) {
dp[i] = 0;
}
for (int i = 1; i <= n; i++) {
int lowerb = max(max(m - sum[n] + sum[i], items[i].p), items[i].q);
for (int j = m; j >= lowerb; j--) {
dp[j] = max(dp[j], dp[j - items[i].p] + items[i].v);
}
}
cout << dp[m] << endl;
}
return 0;
}
---------------- The End ----------------

作者: brooksjay
联系邮箱: jaypark@smail.nju.edu.cn
本文地址:
本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
转载请注明出处, 谢谢!

<i class="fa fa-angle-left"></i>1234<i class="fa fa-angle-right"></i>
brooksj

brooksj

我听说虫洞可以穿梭时空

31 日志
15 分类
40 标签
GitHub E-Mail
Links
  • 博客旧址
  • 苏剑林's blog
  • Karpathy's Blog
  • Ruder's Blog
  • Edward Ma's Blog
  • Miskcoo's Blog
© 2023 brooksj
由 Hexo 强力驱动
|
主题 — NexT.Mist v5.1.4