【做题笔记】多项式/FFT/NTT

HDU1402 - A * B Problem Plus

题目链接

大数乘法是多项式的基础应用，其原理是将多项式$f(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+\cdots +a_n{x}^n$中的$x=10$，然后让大数的每一位当做对应的$10^y$的系数，然后进行多项式乘法运算，以降低大数乘法的时间复杂度（高精度乘法时间复杂度是$O(n^2)$的，多项式乘法是$O(n\log n)$，$n$是位数）。

在进行多项式乘法之后，我们得到了新的多项式$h(x)=c_0+c_1\times 10+c_2\times 10^2+\cdots +c_m\times 10^m$。
这里的$c_0,c_1,\cdots ,c_m$由于在乘法时会发生进位，所以每一位可能会大于$9$，此时要将后面的位的值挪到下一位上去，让这一位保证在$9$以内。

然后从高位到低位找到第一个不是$0$的位，然后逐位输出即可。

如果最终结果是$0$，要特别判断。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod = (1ll << 47) * 7 * 4451 + 1;
const LL g = 3;LL mul(LL x, LL y) {return (x * y - (LL)(x / (long double)mod * y + 1e-3) * mod + mod) % mod;
}template<class T> T power(T a, LL b) {T res = 1;for (; b; b >>= 1) {if (b & 1) res = mul(res, a);a = mul(a, a);}return res;
}LL rev[4000005];void change(LL *y, int len) {for (int i = 0; i < len; ++i) {if (i < rev[i]) swap(y[i], y[rev[i]]);}
}void DFT(LL *y, int len, int on) {change(y, len);for (LL h = 2; h <= len; h <<= 1) {LL wn = power(g, (mod - 1) / h);if (on == -1) wn = power(wn, mod - 2); for (int j = 0; j < len; j += h) {LL w = 1ll;for (int k = j; k < j + h / 2; ++k) {LL u = y[k] % mod;LL t = mul(w, y[k + h / 2]);y[k] = (u + t) % mod;y[k + h / 2] = (u - t + mod) % mod;w = mul(w, wn);}}}
}void NTT(LL *x1, LL *x2, LL len) {for (int i = 0; i < len; ++i) {rev[i] = rev[i >> 1] >> 1;if (i & 1) rev[i] |= len >> 1;}DFT(x1, len, 1);DFT(x2, len, 1);for (int i = 0; i < len; ++i) {x1[i] = mul(x1[i], x2[i]);}DFT(x1, len, -1);LL inv = power(len, mod - 2);for (int i = 0; i < len; ++i) {x1[i] = mul(x1[i], inv);}
}string a, b;
LL x1[2000005], x2[2000005];void main2() {cin >> b;int n = a.length(), m = b.length();int len = 1;while (len <= (n + m) * 2 + 1) len <<= 1;int pt = n - 1;for (char c: a) {x1[pt--] = c - '0';}for (int i = n; i < len + 300; ++i) {x1[i] = 0;}pt = m - 1;for (char c: b) {x2[pt--] = c - '0';}for (int i = m; i < len; ++i) {x2[i] = 0;}NTT(x1, x2, len);for (int i = 0; i < len; ++i) {if (x1[i] > 9) {x1[i + 1] += (x1[i] / 10);x1[i] %= 10;}}pt = len;while (x1[pt] > 0) {if (x1[pt] > 9) {x1[pt + 1] += (x1[pt] / 10);x1[pt++] %= 10;}}while (pt >= 0 and x1[pt] == 0) --pt;if (pt == -1) {cout << 0 << '\n';return;}for (int i = pt; i >= 0; --i) {cout << x1[i];}cout << '\n';
}int main() {ios::sync_with_stdio(false);cin.tie(0); cout.tie(0);LL _ = 1;
//	cin >> _;while (cin >> a) main2();return 0;
}

HDU4609 - 3-idiots

题目链接

要求能够成三角形的概率，我们可以通过可以构成三角形的情况总数/总的情况数来得到。其中总的情况数就是$\left(\genfrac{}{}{0ex}{}{n}{3}\right)$，重点是求解可以构成三角形的情况。

首先我们可以维护一个多项式，$x_i$的系数表示长度为$i$的线段有多少条。那么以样例$1$为例，得到的多项式的系数就是：$[0,1,0,2,1]$。

接下来，我们将这个多项式自己与自己相乘，也就是$[0,1,0,2,1]\times [0,1,0,2,1]$，得到的多项式的系数是：$[0,0,1,0,4,2,4,4,1]$。

这个多项式表示的含义是从第一个$[0,1,0,2,1]$取出一个边，从第二个$[0,1,0,2,1]$取出一个边：
取出两个边的长度和为$0$的有$0$种；
取出两个边的长度和为$1$的有$0$种；
取出两个边的长度和为$2$的有$1$种；
取出两个边的长度和为$3$的有$0$种；
取出两个边的长度和为$4$的有$4$种；
取出两个边的长度和为$5$的有$2$种；
取出两个边的长度和为$6$的有$4$种；
取出两个边的长度和为$7$的有$4$种；
取出两个边的长度和为$8$的有$1$种。

我们现在要利用得到的这个信息来求解答案。

我们先来看一看我们的三条边形成三角形，有怎样的条件：设三角形三条边是$a,b,c$，且$a<b<c$。于是有$a+b>c$。

首先，在我们得到的$[0,0,1,0,4,2,4,4,1]$中，存在同一条边被采用两次的情况，首先我们要遍历每一条边，将长度之和为其二倍的情况$-1$。原始的$4$条边的边长分别是$[1,3,3,4]$，所以我们要减去两条边长度为$[2,6,6,8]$的一种方案。于是变成：$[0,0,0,0,4,2,2,4,0]$。

然后，由于我们两边的边是一样的，左边取$a$右面取$b$，左边取$b$右边取$a$，这两种其实应当是一种方案。所以我们应当将每一个可能的两边长度的可能情况减半。于是变成$[0,0,0,0,2,1,1,2,0]$。

如果我们设最终答案为$total$，将卷积后得到的结果求一个前缀和叫做$pre$，$n$条边最长的那条边为$mx$。先将所有边从小到大进行排序。

接下来依次枚举$n$条边，对于第$i$条边，我们设其长度为$c_i$，并认定其为最长边的情况，求解这种情况时的可行方案数。首先，根据两边之和大于第三边，我们剩下两条边的长度之和应当大于$c_i$，所以从刚刚卷积和处理后得到的信息中，将长度之和为$[c_i+1,\infty ]$的可能情况加到答案中。

那么这一步就是$total\leftarrow total+pre[2mx]-pre[c_i]$。

然后根据我们刚才所设，当前边我们令其为最大边，所以一切包含一个小于$c_i$的边、包含一个大于$c_i$的边都是不合法的，但是因为这两条边的长度和一定大于$c_i$，所以被算进了$total$里面，所以要从$total$里面减掉。我们考虑这样的可能发生情况是$(n-i)\times (i-1)$，其中$n-i$是比$c_i$长的边的数量，$i-1$是比$c_i$短的数量和。

这一步就是$total\leftarrow total-(n-i)\times (i-1)$。

还有一种情况，就是两条边都取了比$c_i$长的，那就是$\left(\genfrac{}{}{0ex}{}{n-i}{2}\right)$种情况，因为是从比$c_i$长的$n-i$个边里任意挑选两个边。

这一步就是$total\leftarrow total-\left(\genfrac{}{}{0ex}{}{n-i}{2}\right)$。

还有一种情况，就是一条边取了自身，另一条边取了其他，存在$n-1$种情况。这里因为之前已经除以$2$，所以不用考虑从哪里搬来的情况。

这一步就是$total\leftarrow total-(n-1)$。

将$n$条边对$total$的贡献全部计算一遍之后，最后的结果就是所有合法的情况数了。

最后的答案就是$\frac{total}{\left(\genfrac{}{}{0ex}{}{n}{3}\right)}$。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 100005;
const LL mod = 998244353;
const LL g = 3, gi = 332748118;
const double PI = acos(-1.0);struct Complex {double x, y;Complex(double _x = 0.0, double _y = 0.0) {x = _x; y = _y;}Complex operator - (const Complex &b) const {return Complex(x - b.x, y - b.y);}Complex operator + (const Complex &b) const {return Complex(x + b.x, y + b.y);}Complex operator * (const Complex &b) const {return Complex(x * b.x - y * b.y, x * b.y + y * b.x);}
};LL rev[4000005];void change(Complex *y, int len) {for (int i = 0; i < len; ++i) {if (i < rev[i]) swap(y[i], y[rev[i]]);}
}void FFT(Complex *y, int len, int on) {change(y, len);for (int h = 2; h <= len; h <<= 1) {Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));for (int j = 0; j < len; j += h) {Complex w(1, 0);for (int k = j; k < j + h / 2; ++k) {Complex u = y[k];Complex t = w * y[k + h / 2];y[k] = u + t;y[k + h / 2] = u - t;w = w * wn;}}}if (on == -1) {for (int i = 0; i < len; ++i) {y[i].x = y[i].x / len + 0.5;}}
}LL n, len = 289999;
LL a[N], x[N * 3], pre[N * 2];
Complex F[N * 2];void main2() {cin >> n;LL mx = 0;memset(x, 0, sizeof(x));for (int i = 1; i <= n; ++i) {cin >> a[i];mx = max(mx, a[i]);++x[a[i]];}len = 1;while (len <= mx * 2 + 1) len <<= 1;for (int i = 0; i <= mx; ++i) {F[i].x = (double)x[i];F[i].y = 0;}for (int i = mx + 1; i < len; ++i) {F[i].x = F[i].y = 0;}for (int i = 0; i < len; ++i) {rev[i] = rev[i >> 1] >> 1;if (i & 1) rev[i] |= len >> 1;}FFT(F, len, 1);for (int i = 0; i < len; ++i) {F[i] = (F[i] * F[i]);}FFT(F, len, -1);for (int i = 0; i < len; ++i) {x[i] = (LL)F[i].x;}for (int i = 1; i <= n; ++i) {--x[a[i] * 2];}for (int i = 1; i <= mx * 2; ++i) {x[i] /= 2;}pre[0] = 0;for (int i = 1; i <= mx * 2; ++i) {pre[i] = pre[i - 1] + x[i];}LL total = 0;for (LL i = 1; i <= n; ++i) {total += (pre[mx * 2] - pre[a[i]]);total -= ((n - i) * (i - 1));total -= (n - 1);total -= ((n - i) * (n - i - 1) / 2);}LL all = n * (n - 1) * (n - 2) / 6;double ans = (double)total / (double)all;cout << fixed << setprecision(7) << ans << '\n';
}int main() {ios::sync_with_stdio(false);cin.tie(0); cout.tie(0);LL _ = 1;cin >> _;while (_--) main2();return 0;
}

CF954I - Yet Another String Matching Problem

题目链接

首先考虑，如果$S$串和$T$串长度相同，该怎么做。我们发现，我们的目的是消除差异性，即对于上下每一对字母，他们最后都要变成同一个字母。

以abcd和ddcb为例，那么上下每一对依次是a-d，b-d，c-c，d-b。我们将每一个字母设为一个点，将每一对当做一个边连起来，连完之后长这个样子：

会发现几个字母会连成一个连通图：

我们将四个字母分为了两个部分：$c$自己是一个部分，$abd$三个字母是一个部分。我们发现，连的边指在原来的两个串中对应的同一个位置的两个字母，他们最终要变成同一个字母。而每一次操作要同时改变两个串中的这个字母，所以我们不难推出，被我们连上边形成连通块的这些字母，他们最后必然要变成同一个字母。我们就让他变成他们之中的一个字母，那么我们发现，其实我们的操作数，就是我们连的边的数量，因为我们每连一次边，就相当于要进行一次题目中的操作，变换一个字母。

这样我们用并查集来实现：对于每一个位置，判断这个位置的两个字母是否相同，如果相同那就不需要我们动了；如果不同，就需要判断这两个字母现在是否在一个连通块，如果在，则不管（已经变成同样的字母了），不在同一个连通块就用并查集合并一下，然后答案$+1$。

因为我们的字符集只有前$6$个字母$abcdef$，所以合并操作这类的都是常数级别。那么，在两个串长度相等的情况下，得到答案的时间复杂度是$O(m)$的。

但是我们原题中$S$串会比较长，我们需要枚举每一个长度为$|T|$的子串，每个子串都要求一次答案，于是时间复杂度变成$O(nm)$的了，直接爆炸。

我们来尝试进行一波优化，看看哪里是最消耗时间的地方。

我们注意到，我们的字符串长度的规模是$1.25\times 10^5$，而我们的字符集大小只有$6$，意味着我们同一个位置可能的字符对只有$6\times 36$种。又由于我们不考虑上下相同的字符串，所以我们需要考虑的字符对情况只有$36-6=30$种。对于一次子串和$T$串的求解过程中，如果串过长，那么相同的字符对会出现很多次，而其中只有一次是有用的，在第一次的时候我们对这两个字母进行连边，从第二次开始再遇到这个字符对，根本就是毫无用处了。所以，如果我们能够判断出，每一个子串跟$T$串匹配的过程中，所有位里是否包含这样的字符对就好了。

怎么来做？先说结论：用到卷积。

先来看一个关于卷积的性质。假设我们有两个多项式：
$f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$
$g(x)=b_0+b_{1}x+b_2{x}^2+\cdots +b_n{x}^n$

试求$f(x)\cdot g(x)$。
我们发现，$f(x)\cdot g(x)=a_0{b}_0+(a_0{b}_1+a_1{b}_0)x+(a_0{b}_2+a_1{b}_1+a_2{b}_0)x^2+(a_0{b}_3+a_1{b}_2+a_2{b}_1+a_3{b}_0)x^3+\cdots$

发现了什么？最后得到的多项式系数，$x^k$的系数是$\sum _{i=0}^k{a}_i{b}_{k-i}$，你会发现他们由所有$i+j=k$的$a_i{b}_j$相加得到，换句话说，组成$x^k$的系数的所有$a_i{b}_j$，都满足他们的下标和$i+j$是定值，且定值为$k$。

然后看看我们需要的是什么？事实上，我们可以总结出我们现阶段的问题是：从$S$串的第$i$位作为第一个字符的长度为$|T|$的子串和$T$串匹配的过程中，字符$c_1$和字符$c_2$是否出现。

在其中的一个固定了$c_1$和$c_2$的子问题中，我们设一个数组$a_i$，表示$S$串的第$i$位是否为$c_1$。如果$a_i=1$，说明是$c_1$；如果$a_i=0$，说明第$i$不是$c_1$。同样的方式设一个数组$b_i$，表示$T$串第$i$位是否为$c_2$。

那么，$T$串的第一个字母和$S$串第$i$个字母对上的时候，字符对$(c1,c2)$出现的条件是什么？如果$\sum _{k=0}^{|T|-1}{a}_{i+k}{b}_k>0$，说明存在这样的字符对$(c_1,c_2)$。
因为这种情况下，$S$串的第$i+k$位和$T$串的第$k$是对上的，只有$a_{i+k}=b_k=1$成立，我们才能说这样的字符对存在。

但是怎么求呢？我们现在有了$a_i,b_i$，但是因为$\sum _{k=0}^{|T|-1}{a}_{i+k}{b}_k$中$a,b$的下标加起来不是一个常值。如果是一个常值的话，我们就可以把$a_i$数组和$b_i$数组当做两个多项式的系数，然后我们进行多项式的乘法，这时$a_i{b}_j$自然就是我们两个多项式乘在一起之后$x^{i+j}$的系数中的一部分。
我们发现，我们要求的这个求和公式的形式等同于$\sum _{k=0}^{|T|-1}{a}_i{b}_j$，如果$i+j$是定值$C$，那么我们所求的求和公式的结果就是多项式乘在一起之后$x^C$的系数。

观察我们想求的公式$\sum _{k=0}^{|T|-1}{a}_{i+k}{b}_k$，我们发现，如果$b_k$中的$k$变成$-k$，哪怕是加一点常数，就完美了。因为这样两个变量$k$就抵消掉了，剩下的和就是一个定值，就可以通过两个多项式相乘之后的系数求得。

把$b_k$中的$k$变成$-k$很简单，我们可以直接将$T$串翻转，这样我们原本的$b_k$就变成了$b_{|T|-k-1}$，我们原本要求的求和公式$\sum _{k=0}^{|T|-1}{a}_{i+k}{b}_k$也变成了$\sum _{k=0}^{|T|-1}{a}_{i+k}{b}_{|T|-k-1}$。

这个时候我们$a,b$的下标之和就是$|T|+i-1$。这里面$|T|$是$T$的长度，是定值；$i$是当前子问题中$T$串第一个字符跟$S$串的哪个字符对齐，在子问题中也是一个定值。于是$|T|$就是一个定值。这样的话，我们就可以利用我们刚刚推出的结论来求解这个求和式了。将两个数组$a,b$当做多项式系数后，进行多项式乘法，得到的结果多项式中，$x^{|T|+i-1}$的系数，就是我们所求的结果。如果这个结果大于$0$，则说明存在这样的字符对$(c_1,c_2)$。

设$f[i][c_1][c_2]$为在$T$串第一个字符跟$S$串的哪个字符对齐时是否存在字符对$(c_1,c_2)$。如果存在，其值为$1$，否则为$0$。

所以整体的流程就是：

1. 对于字符$c_1$，我们先看$S$串中哪些位置是$c_1$。设数组$a$来表示这个信息。如果$a_i=1$，说明字符串$S$的第$i$位是$c_1$；如果$a_i=0$，说明字符串的第$i$位不是$c_1$。
2. 翻转字符串$T$。
3. 对于字符$c_2$，我们先看$T$串中哪些位置是$c_2$。设数组$b$来表示这个信息。如果$b_i=1$，说明字符串$T$的第$i$位是$c_2$；如果$b_i=0$，说明字符串的第$i$位不是$c_2$。
4. 把得到的$a,b$数组当做两个多项式的系数，进行卷积。
5. 得到的结果多项式的系数中，如果$0\le |T|+i-1<n$，那么只要$x^{|T|+i-1}$的系数大于$0$，则可以认为$f[i][c_1][c_2]=1$。如果这一项的系数为$0$，则$f[i][c_1][c_2]=0$。

这样的话，对于每一对$c_1,c_2$，我们都这样操作一遍，就可以得到全部的$f$数组的信息了，即：我们求得了从$S$串的第$i$位作为第一个字符的长度为$|T|$的子串和$T$串匹配的过程中，字符$c_1$和字符$c_2$是否出现。接下来用这个信息，对于每一个$i$的位置，看看所有点对的情况，然后用上面说到的方法用并查集合并一下，就得到了每一个位置的答案。

试试分析现在的时间复杂度。多项式乘法采用FFT或者NTT的话，一次DFT/IDFT的时间复杂度是$O(n\log n)$。我们一共进行了$30$个两个字母不同的字符对的查询，每一次查询进行了$3$次DFT/IDFT，也就是说光是在多项式乘法上，时间复杂度就是$O(30\times 3n\log n)$。经过计算发现这个时间复杂度略微有点大。主要原因是$30\times 3$这个常数太大了。

每一个点对都要进行$3$次DFT，这三次$DFT$分别做了：

1. 将$S$串的$a$数组进行DFT
2. 将翻转后的$T$串的$b$数组进行DFT
3. 将得到的结果进行IDFT。

但是，由于我们的字符集是有限大小的，我们每一个字母作为字符对的第一个字母时，都要求一次$a$数组，进行一次DFT，这是重复操作。$b$数组也是同理。所以我们可以预处理出$S$串所有字符的$a$数组的DFT后的结果，预处理出翻转后$T$串所有字符的$b$数组的DFT后的结果。后续求的时候，直接将预处理出的两个数组相乘，然后进行一次IDFT即可。这样，我们进行DFT/IDFT的次数是两个串的$a,b$数组DFT预处理$6+6=12$次，$30$个字符对的IDFT$30$次，一共只有$42$次。时间复杂度一下子就砍了一半。

最后求解的时候每一个位置遍历所有字符对$36$个，每一个字符对的操作是常数级（因为字符集很小），所以这个做法整体的时间复杂度可以认为是$O(42n\log n+36n)$，经过计算，是可以通过本题的。

代码采用NTT。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod = (1ll << 47) * 7 * 4451 + 1;
const LL g = 3;
const LL N = 300005;string ss, tt;
LL s[125005], t[125005];
LL f[N][8][8];
LL n, m, c;
LL sa[8][N], ta[8][N];
LL tmp[N];
LL rev[N], fa[8], rk[8];LL mul(LL x, LL y) {return (x * y - (LL)(x / (long double)mod * y + 1e-3) * mod + mod) % mod;
}template<class T> T power(T a, LL b) {T res = 1;for (; b; b >>= 1) {if (b & 1) res = mul(res, a);a = mul(a, a);}return res;
}void change(LL *y, int len) {for (int i = 0; i < len; ++i) {if (i < rev[i]) swap(y[i], y[rev[i]]);}
}void DFT(LL *y, int len, int on) {change(y, len);for (LL h = 2; h <= len; h <<= 1) {LL wn = power(g, (mod - 1) / h);if (on == -1) wn = power(wn, mod - 2); for (int j = 0; j < len; j += h) {LL w = 1ll;for (int k = j; k < j + h / 2; ++k) {LL u = y[k] % mod;LL t = mul(w, y[k + h / 2]);y[k] = (u + t) % mod;y[k + h / 2] = (u - t + mod) % mod;w = mul(w, wn);}}}
}int find(int x) {return ((fa[x] == x) ? x : (fa[x] = find(fa[x])));
}void merge(int i, int j) {int x = find(i), y = find(j);if (rk[x] <= rk[y]) fa[x] = y;else fa[y] = x;if (rk[x] == rk[y] && x != y) ++rk[y];
}void main2() {cin >> ss >> tt;n = ss.length(); m = tt.length();for (int i = 0; i < n; ++i) {s[i] = ss[i] - 'a';sa[s[i]][i] = 1;}for (int i = 0; i < m; ++i) {t[i] = tt[i] - 'a';ta[t[i]][m - i - 1] = 1;}LL len = 1;while (len <= (n + m)) len <<= 1;for (int i = 0; i < len; ++i) {rev[i] = rev[i >> 1] >> 1;if (i & 1) rev[i] |= len >> 1;}for (int i = 0; i < 6; ++i) {DFT(sa[i], len, 1);DFT(ta[i], len, 1);}for (int i = 0; i < 6; ++i) {for (int j = 0; j < 6; ++j) {if (i == j) continue;for (int k = 0; k < len; ++k) {tmp[k] = mul(sa[i][k], ta[j][k]); 	}DFT(tmp, len, -1);LL inv = power(len, mod - 2);for (int k = 0; k < len; ++k) {tmp[k] = mul(tmp[k], inv);if (k + 1 - m >= 0 and k + 1 - m < n) f[k + 1 - m][i][j] = tmp[k];}}}for (int i = 0; i < n - m + 1; ++i) {for (int j = 0; j < 6; ++j) {fa[j] = j; rk[j] = 1;}int sum = 0;for (int j = 0; j < 6; ++j) {for (int k = 0; k < 6; ++k) {if (j == k or !f[i][j][k]) continue;int fx = find(j), fy = find(k);if (fx != fy) {++sum;merge(fx, fy);}}}cout << sum << ' ';}
}int main() {ios::sync_with_stdio(false);cin.tie(0); cout.tie(0);LL _ = 1;
//	cin >> _;while (_--) main2();return 0;
}