凸优化用于一类dp优化的问题
可以把一个$O(n)$的复杂度优化成$O(\log n)$
如果我们最后的dp数组在某一维上是凸壳就可以用凸壳来维护
在这种问题中我们如果不管这一维的影响可以很快地求出极值
我们就可以用凸优化的思想来操作
假设我们的dp数组是$F(x)$并且满足$F(x)$的导数是单调函数
我们就可以对$F$函数进行一定的修改使得$F’(x)=F(x)+kx$函数的极值落在我们需要求的那一个点上
于是我们可以对斜率k进行二分，来求出最合适的那一个值
具体来说，这道题里面我们需要求的是树上k+1条点不相交路径的最大权值和
我们发现如果设$F(x)$表示选出x条不相交路径的最大权值，最后的F函数是一个上凸壳
然后我们就可以考虑二分出斜率的偏移量$w$
考虑它的实际意义，就是说每次多选出一条路径就需要加上额外的$w$代价，这样dp我们只需要记录最值和对应的路径条数就可以了
至于是要让路径条数最大化还是最小化因题而异
其实挺简单的算法，多做几道题就会了
在一些神奇的问题中，光是维护点分树（动态点分治）是完全没有办法解决问题的，例如这道题，我们发现查询的区间是和树的形态无关的，这样的问题就不能单纯维护点分树，需要对点分树进行可持久化
关于这个可以去看我上古时期的blog，如果这个都看不懂就直接退群吧。。。。
听起来在点分治上加上了动态两个字，其实就是把点分治中分治中心用树型关系表示出来
在统计到一个点的信息的时候，我们就直接从这个点出发，往点分树父亲节点上跳，比如我们需要查询到u节点距离小于等k的点的个数，我们就先看当前分治中心u管辖区域内有多少个点和他的距离小于等于k，然后我们跳到prt[u]，查询在prt[u]的子树里面有多少个距离小于等于k - dis(prt[u], u)的节点，但是我们发现u的管辖区域内的点存在重复计算，所以我们需要把这一部分用相同的方法容斥掉，就可以计算答案了
修改和查询十分类似，不赘讲（江老著名梗）了
可持久化点分树，就拿上面哪个cf题做例子吧，用问题作为切入点可能会好一些
首先我们可以发现这道题的查询是和树无关的，如果只有单点的查询是很快可以在点分树上处理出来的
并且我们知道查询是可以转化成前缀和查询的，于是就有了维护可持久化点分树的思路
因为平常的点分我们是从下到上更新，这次需要变成从上到下的插入，并且在点分树上添加出来一条新的链，所以在考虑继承上一个点分树的状态的时候，我们需要在插入的时候继承不改变的子树的信息，因此我们需要先将原来的树进行三度化处理，这样每个分治重心的儿子最多只有三个，就可以处理了
三度化的技巧主要在于理解左兄弟右儿子的思路，把所有的兄弟用长度为0的边串起来，这样在新的树上的距离就是在原树上的距离了，并且因为串两个兄弟只需要1个点，所以新的树的总点数也是$O(n)$级别的，不影响复杂度
具体代码实现长成这样：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
28void rebuild(int u, int fa) {
int last = 0, tmp = 0;
for (auto cur : edge[u]) {
int v = cur.first, w = cur.second;
if (v == fa) continue;
++tmp;
if (tmp == 1) {
addedge(u, v, w);
addedge(v, u, w);
last = u;
} else if (tmp == ((signed) edge[u].size()) - (u != 1)) {
addedge(last, v, w);
addedge(v, last, w);
} else {
m++;
addedge(last, m, 0);
addedge(m, last, 0);
last = m;
addedge(m, v, w);
addedge(v, m, w);
}
}
for (auto cur : edge[u]) {
int v = cur.first;
if (v == fa) continue;
rebuild(v, u);
}
}
其中edge是原树的边，addedge是新树的边
然后因为我们需要在点分树上更新一直到节点u（当前的更新节点）
所以我们需要记录一下在点分树上u的每个祖先v是prt[v]的哪一个儿子
这个可以在第一次点分树预处理出来
然后为了保证复杂度（不要多一个log）
最好用$O(n\log n)$预处理，$O(1)$查询的LCA
就是在欧拉序上进行处理
1 | namespace LCA { |
然后我们考虑普通的点分树在维护什么东西
首先是用来记录信息的值和用来进行容斥的值，因为这道题只需要记录距离，所以我们用两个值来表示就可以了，如果需要更复杂的维护就需要用到高级的数据结构了
我们查询到一个点的距离和的时候首先从根往下走，假设当前节点是v，查询节点是u，那么v除了u所在子树的所有点到u的所有距离就是除了u子树所有节点到v的距离和加上除了u所有子树的大小乘上u到v的距离
于是我们就可以进行容斥了
点分树部分的代码是这样的：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
namespace Tree_Devide {
struct Node {
int ch[3], id, siz;
ll sum[2];
} s[N * LOG];
int siz[N], F[N], dep[N], siz_all, rt, cnt = 0;
int root[N];
bool vis[N];
vector<int> pre[N];
void getroot(int u, int fa) {
siz[u] = 1;
F[u] = 0;
for (int i = head[u]; i; i = E[i].nxt) {
int v = E[i].v;
if (v == fa || vis[v]) continue;
getroot(v, u);
siz[u] += siz[v];
F[u] = max(F[u], siz[v]);
}
F[u] = max(F[u], siz_all - siz[u]);
if (F[u] < F[rt]) rt = u;
}
void devide(int u) {
vis[u] = 1;
s[u].id = u;
for (int i = head[u]; i; i = E[i].nxt) {
int v = E[i].v;
if (vis[v]) continue;
F[rt = 0] = siz_all = siz[v];
getroot(v, 0);
dep[rt] = dep[u] + 1;
for (int j : pre[u]) {
pre[rt].push_back(j);
}
for (int j = 0; j < 3; j++) {
if (s[u].ch[j]) continue;
s[u].ch[j] = rt;
pre[rt].push_back(j);
devide(rt);
break;
}
}
}
void insert(int &t, int las, int u, int fa) {
t = ++cnt;
s[t] = s[las];
s[t].siz++;
s[t].sum[0] += getdis(s[t].id, u);
if (fa) {
s[t].sum[1] += getdis(s[fa].id, u);
}
if (s[t].id == u) return;
insert(s[t].ch[pre[u][dep[s[t].id]]], s[las].ch[pre[u][dep[s[t].id]]], u, t);
}
ll query(int t, int u) {
ll res = 0;
for (int i = 0; i < dep[u]; i++) {
int cur = s[t].ch[pre[u][i]];
res += (s[t].siz - s[cur].siz) * getdis(s[t].id, u) + s[t].sum[0] - s[cur].sum[1];
t = cur;
}
res += s[t].sum[0];
return res;
}
void build() {
F[rt = 0] = siz_all = cnt = m;
getroot(1, 0);
root[0] = rt;
devide(rt);
for (int i = 1; i <= n; i++) {
insert(root[i], root[i - 1], a[i], 0);
}
}
}
然后放一下cf757g的完整代码：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#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int, int> pi;
const int N = 4e5 + 10;
const int LOG = 20;
struct Edge {
int v, w, nxt;
} E[N << 1];
int head[N], tot = 0;
int n, m, q, a[N];
vector<pi> edge[N];
void addedge(int u, int v, int w) {
E[++tot] = (Edge) {v, w, head[u]};
head[u] = tot;
}
void rebuild(int u, int fa) {
int last = 0, tmp = 0;
for (auto cur : edge[u]) {
int v = cur.first, w = cur.second;
if (v == fa) continue;
++tmp;
if (tmp == 1) {
addedge(u, v, w);
addedge(v, u, w);
last = u;
} else if (tmp == ((signed) edge[u].size()) - (u != 1)) {
addedge(last, v, w);
addedge(v, last, w);
} else {
m++;
addedge(last, m, 0);
addedge(m, last, 0);
last = m;
addedge(m, v, w);
addedge(v, m, w);
}
}
for (auto cur : edge[u]) {
int v = cur.first;
if (v == fa) continue;
rebuild(v, u);
}
}
namespace LCA {
int dep[N], st[N << 1][LOG], dfn[N], Log[N << 1], ind = 0;
ll dis[N];
void dfs(int u, int fa) {
st[dfn[u] = ++ind][0] = u;
for (int i = head[u]; i; i = E[i].nxt) {
int v = E[i].v;
if (v == fa) continue;
dep[v] = dep[u] + 1;
dis[v] = dis[u] + E[i].w;
dfs(v, u);
st[++ind][0] = u;
}
}
int checkmin(int x, int y) {
return dep[x] < dep[y] ? x : y;
}
void ST() {
dfs(1, 0);
for (int i = 2; i <= ind; i++) {
Log[i] = Log[i >> 1] + 1;
}
for (int j = 1; j < LOG; j++) {
for (int i = 1; i + (1 << j) - 1 <= ind; i++) {
st[i][j] = checkmin(st[i][j - 1], st[i + (1 << (j - 1))][j - 1]);
}
}
}
int lca(int x, int y) {
x = dfn[x], y = dfn[y];
if (x > y) swap(x, y);
int k = Log[y - x + 1];
return checkmin(st[x][k], st[y - (1 << k) + 1][k]);
}
ll getdis(int x, int y) {
return dis[x] + dis[y] - dis[lca(x, y)] * 2;
}
}
using LCA::getdis;
namespace Tree_Devide {
struct Node {
int ch[3], id, siz;
ll sum[2];
} s[N * LOG];
int siz[N], F[N], dep[N], siz_all, rt, cnt = 0;
int root[N];
bool vis[N];
vector<int> pre[N];
void getroot(int u, int fa) {
siz[u] = 1;
F[u] = 0;
for (int i = head[u]; i; i = E[i].nxt) {
int v = E[i].v;
if (v == fa || vis[v]) continue;
getroot(v, u);
siz[u] += siz[v];
F[u] = max(F[u], siz[v]);
}
F[u] = max(F[u], siz_all - siz[u]);
if (F[u] < F[rt]) rt = u;
}
void devide(int u) {
vis[u] = 1;
s[u].id = u;
for (int i = head[u]; i; i = E[i].nxt) {
int v = E[i].v;
if (vis[v]) continue;
F[rt = 0] = siz_all = siz[v];
getroot(v, 0);
dep[rt] = dep[u] + 1;
for (int j : pre[u]) {
pre[rt].push_back(j);
}
for (int j = 0; j < 3; j++) {
if (s[u].ch[j]) continue;
s[u].ch[j] = rt;
pre[rt].push_back(j);
devide(rt);
break;
}
}
}
void insert(int &t, int las, int u, int fa) {
t = ++cnt;
s[t] = s[las];
s[t].siz++;
s[t].sum[0] += getdis(s[t].id, u);
if (fa) {
s[t].sum[1] += getdis(s[fa].id, u);
}
if (s[t].id == u) return;
insert(s[t].ch[pre[u][dep[s[t].id]]], s[las].ch[pre[u][dep[s[t].id]]], u, t);
}
ll query(int t, int u) {
ll res = 0;
for (int i = 0; i < dep[u]; i++) {
int cur = s[t].ch[pre[u][i]];
res += (s[t].siz - s[cur].siz) * getdis(s[t].id, u) + s[t].sum[0] - s[cur].sum[1];
t = cur;
}
res += s[t].sum[0];
return res;
}
void build() {
F[rt = 0] = siz_all = cnt = m;
getroot(1, 0);
root[0] = rt;
devide(rt);
for (int i = 1; i <= n; i++) {
insert(root[i], root[i - 1], a[i], 0);
}
}
}
using Tree_Devide::root;
using Tree_Devide::build;
using Tree_Devide::insert;
using Tree_Devide::query;
int main() {
scanf("%d %d", &n, &q);
m = n;
for (int i = 1; i <= n; i++) {
scanf("%d", &a[i]);
}
for (int i = 1; i < n; i++) {
int u, v, w;
scanf("%d %d %d", &u, &v, &w);
edge[u].push_back(make_pair(v, w));
edge[v].push_back(make_pair(u, w));
}
rebuild(1, 0);
LCA::ST();
build();
ll lastans = 0;
while (q--) {
int op, l, r, u;
scanf("%d", &op);
if (op == 1) {
scanf("%d %d %d", &l, &r, &u);
l ^= lastans;
r ^= lastans;
u ^= lastans;
lastans = query(root[r], u) - query(root[l - 1], u);
printf("%lld\n", lastans);
lastans %= (1 << 30);
} else {
scanf("%d", &u);
u ^= lastans;
swap(a[u], a[u + 1]);
insert(root[u], root[u - 1], a[u], 0);
}
}
return 0;
}
min25 筛法可以解决一类积性函数前缀和的问题
这样的函数需要满足两个个条件才可以用Min25筛法：
我们先考虑算出所有质数的$f(x)$函数的和，也就是素数和减去素数个数
首先我们假设需要求出
即小于等于n的素数和
考虑函数$g(n,y)$表示小于等于n的所有数里面是质数或者最小质因子大于$pri_j$的数的和
那么最后需要求出的答案就是$g(n,Psiz)$
考虑g函数怎么转移
当$pri_j^2>n$的时候，显然不会有新的贡献，于是
否则的话我们把最小质因子是$pri_j$的贡献去掉
因为最小质因子小于$pri_j$的情况我们已经统计过了，所以需要把这些情况减掉
我们可以发现这里的f需要满足是完全积性函数
我们现在算出来了所有质数的$f(x)$的和
那么就可以开始算$f(x)$的前缀和了
设$S(n, j)$表示小于等于n的数里面最小质因子大于等于$pri_j$的所有数的函数和
首先考虑所有质数的贡献
然后枚举和数的最小质因子和其次幂来进行统计
这样的算法复杂度是$\frac{n^{\frac{3}{4}}}{\log(n)}$
一般来说可能比杜教快一点？
适用范围比杜教大一些吧
挺好用的一个算法
$\begin{bmatrix} n \\ m \end{bmatrix}$表示把$n$个数划分成$m$个环排列的方案数
枚举最后一个点的位置
枚举最后一个环排列的大小
$\begin{Bmatrix} n \\ m \end{Bmatrix}$表示把$n$个数划分成$m$个无序集合的方案数
枚举最后一个点的位置
枚举最后一个集合的大小
证明如下：
比较简单，直接考虑怎么算出$x^k$项的系数就可以了，省略
也可以转换成组合数和阶乘的形式
证明如下：
将公式2代入公式1
得到：
交换和号得：
显然当$n = m$的时候$\sum_{k = m} ^ n (-1)^{n - k}\begin{bmatrix} n \\ k \end{bmatrix} \begin{Bmatrix} k \\ m \end{Bmatrix}$是1
所以有：
因为对于任意的x上面的式子都成立，所以可以得到当$n\not= m$的时候$\sum_{k = m} ^ n (-1)^{n - k}\begin{bmatrix} n \\ k \end{bmatrix} \begin{Bmatrix} k \\ m \end{Bmatrix}$是0
证明如下：
将公式1代入公式2
得到：
交换和号得：
显然当$n = m$的时候$\sum_{k = m} ^ n \begin{Bmatrix} n \\ k \end{Bmatrix}\begin{bmatrix} k \\ m \end{bmatrix} (-1)^{k - m}$是1
所以有：
因为对于任意的x上面的式子都成立，所以可以得到当$n\not= m$的时候$\sum_{k = m} ^ n \begin{Bmatrix} n \\ k \end{Bmatrix}\begin{bmatrix} k \\ m \end{bmatrix} (-1)^{k - m}$是0
ps:
当然$(-1)^{k - m}$和$(-1)^{n - k}$没有任何区别
因为$n - m$是偶数的时候，没有影响
$n - m$是奇数的时候，相当于全体取负，不影响结论
证明如下：
因为一个环排列对应着一个置换，然后一个置换又可以差分成很多的环排列，所以假设前$n$个数分成了$k$个环排列，并且这k个环排列中有$m$个不变，剩下的$k - m$个和新的环排列组合起来变成一个完整的环排列，也就是第$m+1$个环排列
证明如下：
假设不和当前数在同一个集合里面的数是k个，所以这k个数会分成m个集合，然后再用组合数算一算和当前数在一个集合的数的方案数
其实就是第一类斯特林数的递推式2的化简版本，组合意义一模一样
证明如下：
考虑第k+1个数是第一个被放进第m+1个集合中的数
所以前k个数会被放进m个集合中，并且还有$(m+1)^{n - k}$个数可以随便放
对于：
有：
反之亦然
证明可以直接带入然后用公式3和公式4进行理解
假设需要求第一类斯特林数第n行
首先有：
就是下降幂换成上升幂，过程不变是是没有符号了，和公式1同理
那么也就是说$x^{\overline{n}}$是第n行斯特林数的生成函数
然后可以用$O(n\log^2 n)$的时间分治fft来算
也可以用$O(n\log n)$的倍增算法
然后$f_n(x+n)$实际上就是把$f_n(x)$的各项$x^k$换成$(x+n)^k$
可以用二项式系数展开并且用$O(n\log n)$时间卷积求出
实际上就是第二类斯特林数的特殊表达形式
发现没有其实这是一个卷积：
然后可以直接卷积算出来了
本文所讲所有的多项式运算全部基于vector的使用
1 | typedef vector<int> Poly; |
自认为多项式最重要的思想就是倍增了
其次就是分治
建议先看看fft残疾人手册和ntt残疾人手册
$C(n)=A(n)+B(n)$
就是系数相加，非常简单的实现
1 | Poly add(Poly a, Poly b) { |
$C(n)=A(n)-B(n)$
系数相减，依旧简单的实现
1 | Poly sub(Poly a, Poly b) { |
$C(n)=\sum_{i=0}^nA(i)B(n-i)$
就是ntt
需要预处理原根的次幂
交换位置的时候动态处理就好，因为在倍增过程中ntt的长度不一样
1 | int w[2][N]; |
对于多项式$A$，我们要求出$B$
满足$A\ast B=1(\mod x^{n+1})$
这里的$\mod x^n$是保留前$n+1$个系数的意思，不然$B(n)$会有无穷项
我们假设已经求出了$B’$满足$A\ast B’=1(\mod x^{\lceil\frac{n}{2}\rceil})$
那么因为存在$A\ast B=1(\mod x^{\lceil\frac{n}{2}\rceil})$
所以$A\ast (B-B’)=0(\mod x^{\lceil\frac{n}{2}\rceil})$
又因为A不是0，所以$(B-B’)=0(\mod x^{\lceil\frac{n}{2}\rceil})$
两边同时平方（模数也要平方）
$B^2-2\ast B\ast B’+B’^2=0(\mod x^n)$
因为有个$B^2$不好计算，所以把等式两边同时乘上$A$
得到$B-2\ast B’+A\ast B’^2=0(\mod x^n)$
即$B=2\ast B’-A\ast B’^2(\mod x^n)$
这样就可以递归求解了
1 | Poly inv(Poly a, int n) { |
对于n次多项式$A$和m次多项式$B$，构造出小于n-m次多项式$C$和小于m次的多项式$D$
满足$A=B*C+D$
我们考虑一下多项式$A$的系数反转形式记做$RevA$
所以$RevA(x)=x^nA(\frac{1}{x})$
我们把原来的式子左右同时乘上$x^n$
得到$A(\frac{1}{x})\ast x^n=B(\frac{1}{x})\ast C(\frac{1}{x})\ast x^n+D(\frac{1}{x})\ast x^n$
所以$RevA(x)=RevB(x)\ast RevC(x)+x^{n-m+1}RevD(x)$
在$\mod x^{n-m+1}$的意义下发现$RevD$就没了，然而并不影响$RevA,RevB,RevC$
所以我们只需要求出$RevB$在$\mod x^{n-m+1}$意义下的逆元，然后就可以算出$RevC$和$RevD$
1 | Poly rev(Poly a) { |
因为多项式的形式很简单，所以求导直接可以模拟出来
1 | Poly deri(Poly a) { |
积分也直接模拟吧。。
1 | Poly inte(Poly a) { |
对于n次多项式$A$我们需要求出n次多项式$B$满足$B=\ln(A)$
两边同时求导
$B’=\frac{A’}{A}$
所以把$A$的导和逆乘起来然后积分回去就可以了。。
1 | Poly ln(Poly a) { |
对于n次多项式$A$我们需要求出来n次多项式$B$满足$B=e^A$
所以$\ln(B)=A $
然后我们构造一个函数$F(B)=\ln(B)-A$
我们要求函数$F(A)=0$的时候的解
根据泰勒展开和牛顿迭代：
$f(x)=f(x_0)+f’(x_0)(x-x_0)$
我们可以得到迭代公式：$x=x_0-\frac{f(x_0)}{f’(x_0)}$
对应到求exp上面来：$B=B_0-\frac{F(B_0)}{F’(B_0)}$
展开得到：$B=B_0-\frac{\ln(B_0)-A}{\frac{1}{B_0}}=B_0(1-\ln(B_0)+A)$
然后就可以根据上面的这个式子进行倍增，因为每一次计算长度会变成原来的两倍
1 | Poly exp(Poly a, int n) { |
给定n次多项式$A$，求$B$满足$B^2=A(\mod x^{n+1})$
和求exp的思路比较类似，同样用到了牛顿迭代法
构造函数$F(B)=B^2-A$，求$F(B)=0$的解
带入迭代公式：$x=x_0-\frac{f(x_0)}{f’(x_0)}$
得到$B=B_0-\frac{B_0^2-A}{2B_0}=\frac{B_0^2+A}{2B_0}=\frac{1}{2}\ast (B_0+\frac{A}{B_0})$
依旧是倍增进行求解
1 | Poly sqrt(Poly a, int n) { |
数论函数是莫比乌斯反演和杜教筛的前置知识
非常简单的定义：$e(n)=[n=1]$
也很简单的定义：$1(n)=1$
依旧简单的定义：$id(n)=n$
莫比乌斯函数$\mu(n)$，当n有平凡因子的时候值是$0$，否则就是$(-1)^{质因子个数}$，即：
莫比乌斯函数的性质：
欧拉函数$\phi(n)$的定义是$[1,n]$中和n互质的数的个数
假设集合$P$包含了所有n的质因子，那么有：
欧拉函数的性质：
对于两个数论函数$f(n),g(n)$，定义Dirichlet卷积为：
并且Dirichlet卷积满足交换律、结合律和分配率，即：
由Dirichlet卷积我们可以推导出：
如果函数$f,g$都是积性函数，那么$f*g$也是积性函数
如果我们无法快速（在$n\log(n)$以内）求出$f\ast g$
那么我们就可以用$O(n\log((n))$的时间求出我们想要的$f\ast g $
可以这样处理：
1 | typedef vector<int> Poly; |
如果两个函数$f, g $，满足
那么有：
也可以表示成：
但是后面这种形式应该不太常用
证明很简单，我们其实就是要证明：
$f=g\ast1$和$g=f \ast\mu$是等价的
证明其实很简单，第一个式子两边同时卷上$\mu$或者第二个式子两边同时卷上$1$
注意在Mobius反演的时候有几种常用技巧：
杜教筛，很方便，但是有其局限性
一般情况是我们需要计算一个积性函数$f$的前缀和：
于是我们可以构造一个有良好特殊性质的积性函数$g$，使得$f\ast g$可以快速算出并具有一定性质
然后我们可以快速算出：
把上面的式子交换一下和号得到：
就变成了：
也就是说：
所以：
我们发现上面的这个式子右边的部分，一个可以快速计算，一个可以根号分块递归成子问题
注意我们需要用线性筛预处理一部分前缀和大概$O(n^{\frac{2}{3}})$个，这样时间效率是$O(n^{\frac{2}{3}})$的
常数比较大，主要是需要map来对ll值和前缀和进行映射
定义神树集是一个二维平面上的点集，使得存在至少一条直线恰好经过两个点
求$x\in[0,n],y\in[0,n]$的神树集个数$\mod 12345678$
首先有一个结论：只有点数小于2或者所有点共线的点集不是神树集
那我们就把这一部分答案容斥掉就可以了
首先我们可以枚举两个端点的横纵坐标差然后进行统计（这里没有统计平行坐标轴的）
其中$(n+1-x)(n+1-y)$是可以选择的左边端点的数量
然后我们枚举gcd
然后把后面拆开按照g的次数分类
变成了
其中
因为我们有：
所以可以对$G(n)$进行化简，把x和y的贡献分别来看，得到：
同理，我们也可以对$H(n)$进行化简
至于有关$i=1$的细节问题，请自行自习思考
然后发现$\sum_{i=1}^n\phi(i),\sum_{i=1}^ni\phi(i),\sum_{i=1}^ni^2\phi(i) $都是可以用杜教筛算出来的
然后$\sum_{i=1}^n 2^i, \sum_{i=1}^n 2^ii,\sum_{i=1}^n 2^i i^2$都是可以用做差法在$O(\log n)/O(1)$的时间内算出来的
到这里我们已经完成了一大部分
剩下的就是对还没有讨论的情况进行分析
首先因为我们要算的是不合法情况，所以应该去掉被包含合法的两个点的个数
因此我们需要算出不和坐标轴平行的两个点点集的个数
所有点对个数是$\frac{(n+1)^2((n+1)^2-1)}{2}$，去掉与坐标轴平行的个数是$2(n+1)\frac{(n+1)n}{2}=(n+1)^2n$
其次因为在与坐标轴平行的情况中，点数大于2的情况我们没有减去
这一部分的情况实际上我们只用考虑$2(n+1)$个单独的行就可以了
对于每一行我们用所有的方案数减去小于3的方案数
就是$2^{n+1}-\frac{n(n+1)}{2}-(n+1)-1$
最后是所有方案数中点数小于2的也不合法，一共有$(n+1)^2+1$种
然后就直接用总的方案数$2^{(n+1)^2}$种减去不合法的就可以了
注意一下$2^{(n+1)^2}$快速幂的时候指数不能直接取mod
并且还有一堆恶心的地方需要分类讨论，因为模数真的很奇葩
1 | #include<bits/stdc++.h> |
建议先阅读FFT的残疾人手册
首先fft的计算方式是在实数域下进行运算的
所以难免会有精度误差，但是如果我们需要在模的意义下快速求多项式的卷积FFT就失灵了
因此NTT就出现了
我们需要一个在模意义下成立的具有单位根性质的东西
那么单位根的性质我们用到了哪一些呢？
模数Mod的原根g的定义是需要满足：
而模数Mod定要满足是质数且可以表示成$Mod = k\cdot2^r+1$的形式
在这种情况下令$g_n=g^k\pmod {Mod}$
就可以使得：
这样的话，我们就可以用原根来代替单位根进行模意义下的快速变换了
主要过程和fft没有本质区别
FFT/NTT可以把本来复杂度是$\mathcal O(n^2)$的多项式乘法优化成$\mathcal O(n\log n)$的复杂度
在比较高等级的比赛中会用吧
在电脑上存一个板子也挺方便的
一个$n$次多项式
的点值表示法为
其实就是把$x_0,x_1,…,x_{n-1}$代入$A(x)$中得到的$y_0,y_1,…,y_{n-1}$
一个大小为$n$点值表示法可以唯一确定一个$n$次多项式
因为点值表示法实际上就是n个方程
这样是可以唯一确定系数集合$a_0,a_1,…,a_{n-1}$的
我们把多项式$A(x),B(x)$在$x_0,x_1,…,x_{2n-2}$处的点值分别求出来，得到
把这些点值乘起来变成
就可以还原成多项式$A(x)*B(x)$
给没有复数基础知识的OIer普及一下：$i=\sqrt{-1}$，是一个常数
然后形如$cos(\theta)+i\cdot sin(\theta)$的复数有一个很好的性质
$(cos(\theta)+i\cdot sin(\theta))^k=cos(k\theta)+i\cdot sin(k\theta)$
这个东西可以用数学归纳法暴力展开证明
n次单位根是形如$cos(\frac{2\pi k}{n})+i\cdot sin(\frac{2\pi k}{n})$的复数
我们记$\omega_{n}=cos(\frac{2\pi}{n})+i\cdot sin(\frac{2\pi}{n})$
然后有$\omega_n^k=cos(\frac{2\pi k}{n})+i\cdot sin(\frac{2\pi k}{n})$
所以说n次单位根其实就是把单位元分成了n等份
根据这个性质很容易得到
也就是说存在：
根据三角函数的周期性可以得到一个式子：
结合$(7)(8)$可以得到：
还有一个很好的性质：
第$(9)$和$(10)$个式子都是我们要用到重要性质
DFT就是所谓的离散傅里叶变换（Discrete Fourier Transform, DFT）
因为朴素的多项式转化成点值的复杂度是$n^2$的，很慢
所以考虑怎么把这个东西优化下来
首先DFT的优化用到了分治的思想
根据第$(9)$个式子，我们来猜想一些性质
多项式：
把A的所有项按照x指数的奇偶性来分类变成
令：
所以有：
那么$A(x)$就变成了$F(x)$在$x^2$处的点值和$G(x)$在$x^2$处的点值和x的积之和
在取x集合的时候令：
就可以得到：
这样的话我们只需要计算出$F(x)$和$G(x)$的值就可以算出多项式$A(x)$的点值了
而因为每次我们递归的时候都会缩小一半的计算范围，只不过需要预先把n补到2的次幂
所以这样的时间复杂度是：
我们做完了DFT，现在需要用点值还原出多项式
也就是说要从这个有n个等式的线性方程组：
解出$a_0,a_1,…a_{n-1}$的值。
我们把上面的线性方程组写成矩阵形式，发现是：
也就是说我们要求出矩阵
的逆矩阵$P’$。
在这里直接给出逆矩阵$P’$并证明一定成立
那么对于矩阵$Q=PP’$，有$Q_{i,j}=\sum_{k=0}^{n-1}P_{i,k}\cdot P’_{k,j}$
即$Q_{i,j}=\sum_{k=0}^{n-1}\frac{1}{n}\omega_n^{-ik}\cdot \omega_n^{kj}=\frac{1}{n}\sum_{k=0}^{n-1}(\omega_n^{j-i})^k $
$i\not=j$：$Q_{i,j}=\frac{1}{n}\sum_{k=0}^{n-1}(\omega_n^{i-j})^k=\frac{1}{n}\cdot \frac{(\omega_n^{i-j})^n-1}{\omega_n^{i-j}-1}$
因为公式$(8)$，可以得到$(\omega_n^{i-j})^n=1$，所以$Q_{i,j}=0$
$i=j$：$Q_{i,j}=1$
所以$P$和$P’$互为逆矩阵
也就是说现在要做的就是求出$P’$和点值矩阵的乘积，即：
再来观察一下原来的DFT过程中我们做的事，就是求出了
那这个时候我们只需要把$A(\omega_n^0) ,A(\omega_n^1) ,\dots,A(\omega_n^{n-1})$当作点值，并且把DFT过程中的$\omega_n^{k}$替换成$\omega_n^{-k}$就可以了
dd只要理解了最初的迭代过程应该就问题不大了
就是简单的模拟出分治过程
omiga数组最好是可以预处理出来，但是这样的实现方式依旧有很大的常数
1 | void fft(int n, Complex *p, int pos, int step, Complex *omg) { |
我们考虑每一次的分组递归
第一次显然是按照最后一个二进制位的奇偶性来分组的
进入下一层迭代后就发现变成了一个$\frac{n}{2}$次多项式，可以看做所有的指数都除以了2，参考$(14)$
所以最后一位就没有意义了，下一次排序就是按照倒数第二位的奇偶性来进行分组的
发现本质就是把二进制反转进行交换，很奇妙的
可以这样预处理：
1 | for (int i = 0; i < n; i++) { |
所以就可以用迭代来进行运算啦
1 | void transform(Complex f[N], Complex w[N]) { |
总的板子
1 | #include<bits/stdc++.h> |
当DP的转移可以用转移矩阵表示出来的时候，这个DP一般都是可以用矩阵快速幂来进行优化的
而且矩阵快速幂优化DP的特征非常的明显，一般是有一维状态数特别大，一维状态数特别小
然后建立转移矩阵的小技巧就是$trans_{i,j}$表示从状态i转移到状态j的贡献
然后矩阵快速幂还可以把乘法变成加法，加法变成取min/max，可以证明这样还是满足结合律的，实际的dp含义也不会被改变
是一种优化转移的很好的方式
当然，矩阵维护DP也可以扩展到动态DP，在这里跟题目无关就不讲了
数据结构优化DP
一般都是直接在数据结构里面维护每个DP值（和那个位置上对应的贡献），这个时候需要把移动指针产生的贡献全部拆分开来，然后如果贡献的叠加范围是存在某个特殊的性质（比如区间，或者某个数的倍数）就可以选择相应的数据结构进行维护
数据结构的使用非常灵活，需要多加练习才可以熟练掌握
一种比较简单的实现形式，通常是有比较显然的决策单调性
通常有两种方式：
数据结构中只需要存点，然后可以直接每次取维护的队首/栈顶来进行更新，维护简单
我们发现每个决策点的最优贡献是一个区间
那么我们就维护一个很多个区间组成的队列
这个队列每次加入的时候就可以进行比较并维护
具体实现是这样的：
3.1.2的单调队列还可以用来维护四边形不等式优化的DP
斜率优化的一个典型形式就是$dp_{i}=dp_{j}+a_i*b_j+c_j$
这时候因为出现了$a_i*b_j$，所以如果可以优化一般都是斜率优化（大概是因为其他的优化方式都不涉及乘除）
然后就可以非常无脑的把这个式子用一次函数式表示出来
比如上面的这个式子经过移项变成$-a_i*b_j+dp_{i}=dp_j+c_j$
观察一下这个式子发现如果令
那么原来的式子实际上就变成了$y=kx+b$
我们显然要求的就是截距b的某个极值
这样我们就可以把一个决策点j转化成二维平面上的一个点$(b_j,dp_j+c_j)$
我们每次要用斜率$-a_i$去截出一个最大的$dp_i$
很显然这个最优决策点一定存在于凸壳上面，所以直接用单调队列/单调栈维护出凸壳就可以了
然后根据题目看需不需要在凸壳上二分出最优的决策点就可以了
至于是上凸壳还是下凸壳，就可以根据k的单调性确定了
但是在x坐标$b_j$不单调的时候就只有使用3.2.2的方法了
有的时候我们不方便找到决策点但是可以知道，如果$i$的决策点是$p_i$，那么对于$j<i<k$有$p_j\le p_i\le p_k$
那我们直接对于mid包里找到$p_{min}$然后递归成区间求解就可以了
这里对应斜率优化的时候$x$坐标不单调的情况
因为x不单调所以我们没有办法（其实也有办法）动态维护凸壳所以我们再用分治解决这个问题
首先选择一个mid，然后递归求解左边区间的问题，这个时候我们就已经知道了左边区间的所有的DP值
然后于是就可以把左边的所有点按照x归并排序之后建成凸壳并用右边的点进行询问，询问完了就可以继续递归右子区间进行解题了
分治复杂度是$nlogn$的
看起来是最简单的结构，但是用到了最复杂的思想
一般是需要证明一个这样的式子：
$w_{i,j+1}+w_{i+1,j}\le w_{i, j}+ w_{i+1,j+1}$
所以就可以得到
$w_{i+1,j}-w_{i+1,j+1}\le w_{i,j}-w_{i,j+1}$
$w_{i,j+1}-w_{i+1,j+1}\le w_{i,j}-w_{i+1,j}$
推导出的通式就是：$对于i<i’\le j<j’,有w_{i,j}+w_{i’,j’}≤w_{i’,j}+w_{i,j’} $
然后对于一个常见的转移方程$dp_i=dp_j+w_{i,j}$
设i的决策点是$p_i$，那么i-1的决策点$p_{i-1}$一定满足：
任意的$k<p_{i-1}$
有$dp_{p_{i-1}}+w_{p_{i-1},i-1}\le dp_{k}+w_{k,i-1}$
即$w_{p_{i-1},i-1}-w_{k,i-1}\le dp_{k}-dp_{p_{i-1}}$
如果w满足平行四边形不等式，那么也一定有：
$w_{p_{i-1},i}-w_{k,i}\le dp_{k}-dp_{p_{i-1}}$
所以有$dp_{p_{i-1}}+w_{p_{i-1},i} \le dp_{k}+w_{k,i}$
就证明了单调性，因为这个东西也是区间覆盖
所以还是可以用单调队列进行维护的
完结撒花