上一章:浅谈 FFT (1){:target=”_blank”}
3.4 多项式乘法的迭代实现
递归实现 FFT 的代码,看着很优美,也的确能过掉这道板子题{:target=”_blank”}
但是,常数太大
面对 105 以上的数据,有很大几率会 TLE,需要优化
观察下我们 FFT 中分治得到的结果,下面以 n=8 为例
[a _ 0, a _ 1, a _ 2, a _ 3, a _ 4, a _ 5, a _ 6, a _ 7] \left\\{ \begin {aligned} [a _ 0, a _ 2, a _ 4, a _ 6] [a _ 0, a _ 2, a _ 4, a _ 6] \left\\{ \begin {aligned} [a _ 0, a _ 4] [a _ 0, a _ 4] \left\\{ \begin {aligned} a _ 0 \\\ a _ 4 \end {aligned} \right. \\\ [a _ 2, a _ 6] \left\\{ \begin {aligned} a _ 2 \\\ a _ 6 \end {aligned} \right. \end {aligned} \right. \\\ [a _ 1, a _ 3, a _ 5, a _ 7] \left\\{ \begin {aligned} [a _ 1, a _ 5] [a _ 1, a _ 5] \left\\{ \begin {aligned} a _ 1 \\\ a _ 5 \end {aligned} \right. \\\ [a _ 3, a _ 7] \left\\{ \begin {aligned} a _ 3 \\\ a _ 7 \end {aligned} \right. \end {aligned} \right. \end {aligned} \right.
(看下面这张图可能会更好懂一点?)
[\color{red}{a _ 0}, \color{blue}{a _ 1}, \color{red}{a _ 2}, \color{blue}{a _ 3}, \color{red}{a _ 4}, \color{blue}{a _ 5}, \color{red}{a _ 6}, \color{blue}{a _ 7}]
\left\\{
\begin {aligned}
[a _ 0, a _ 2, a _ 4, a _ 6]
[\color{red}{a _ 0}, \color{blue}{a _ 2}, \color{red}{a _ 4}, \color{blue}{a _ 6}]
\left\\{
\begin {aligned}
[a _ 0, a _ 4]
[\color{red}{a _ 0}, \color{blue}{a _ 4}]
\left\\{
\begin {aligned}
a _ 0
\\\
a _ 4
\end {aligned}
\right.
\\\
[\color{red}{a _ 2}, \color{blue}{a _ 6}]
\left\\{
\begin {aligned}
a _ 2
\\\
a _ 6
\end {aligned}
\right.
\end {aligned}
\right.
\\\
[\color{red}{a _ 1}, \color{blue}{a _ 3}, \color{red}{a _ 5}, \color{blue}{a _ 7}]
\left\\{
\begin {aligned}
[\color{red}{a _ 0}, \color{blue}{a _ 4}]
[\color{red}{a _ 1}, \color{blue}{a _ 5}]
\left\\{
\begin {aligned}
a _ 1
\\\
a _ 5
\end {aligned}
\right.
\\\
[\color{red}{a _ 3}, \color{blue}{a _ 7}]
\left\\{
\begin {aligned}
a _ 3
\\\
a _ 7
\end {aligned}
\right.
\end {aligned}
\right.
\end {aligned}
\right.
对比我们的原序列和操作之后得到的序列:(这里只写下标)
原序列:0,1,2,3,4,5,6,7
后序列:0,4,2,6,1,5,3,7
啥也看不出来对吧,我们把所有数转成二进制之后再看:
原序列:000,001,010,011,100,101,110,111
后序列:000,100,010,110,001,101,011,111
这个就不难找出规律了吧?
每个位置分治后的最终位置,为其二进制位翻转后得到的位置
我只能看出这个规律,不会证。。wtcl
我们发现了这个规律,但是要 O(n) 处理出每个数的二进制表示翻转之后的结果,还需要对位运算有一定了解
关于位运算的一些基本知识{:target=”_blank”}
我们这么考虑这个问题:
对于一个数 i,我们将其二进制表示记做 f[i],考虑如何求 f[i]
假设我们已经算出来了 f[0∼i−1]
结论:
若 i 是偶数,那么 f[i]=f[i/2]+
'0'
若 i 是奇数,那么 f[i]=f[i/2]+
'1'
这个也不难理解,i/2 即为 i 去掉末位之后的结果,然后如果 i 是奇数,则末位为 1,否则末位为 0
然后再考虑如何求每个数的二进制表示翻转之后的结果
由于它是翻转的,所以我们直接把结论得到的东西反过来即可
即
若 i 是偶数,那么 f[i]=
'0'
+f[i/2]若 i 是奇数,那么 f[i]=
'1'
+f[i/2]
然后就不难写出 O(n) 求二进制翻转的代码了:
int bit = log_2(n) - 1;
for (int i = 0; i < n; i ++ ) f[i] = f[i >> 1] >> 1 | (i & 1) << bit;
O(loglogn) 求 log2(n) 的代码在这里{:target=”_blank”}
这样,我们就可以把原序列中的每个数,先放在最终的位置上,再一步一步向上合并
C ++ 代码
#include <cmath>
#include <cstdio>
// complex 板子
struct cp
{
double a, b;
cp(double x = 0, double y = 0){a = x, b = y;}
cp operator + (const cp &t) const {return cp(a + t.a, b + t.b);}
cp operator - (const cp &t) const {return cp(a - t.a, b - t.b);}
cp operator * (const cp &t) const {return cp(a * t.a - b * t.b, a * t.b + b * t.a);}
};
// loglogn 求 logn 板子
int log_2(int x)
{
int res = 0;
if (x & 0xffff0000) res += 16, x >>= 16;
if (x & 0xff00) res += 8, x >>= 8;
if (x & 0xf0) res += 4, x >>= 4;
if (x & 0xc) res += 2, x >>= 2;
if (x & 2) res ++ ;
return res;
}
// 交换两个虚数
void swap(cp &a, cp &b)
{
static cp t;
t = a, a = b, b = t;
}
const int N = 400005;
const double pi = 2 * acos(-1);
int n, m, k;
int r[N]; // 存变换之后的下标,即上面的 f[i]
cp a[N], b[N];
void FFT(cp a[], int f)
{
// 将 a 变换成分治后的序列
for (int i = 0; i < k; i ++ ) if (i < r[i]) swap(a[i], a[r[i]]);
// 要用到的四个变量,t 存单位根,w 存单位根的 i 次方,x 和 y 做临时变量
cp t, w, x, y;
// 枚举区间长度
for (int len = 2; len <= k; len <<= 1)
{
// 找到单位根
t = cp(cos(pi / len), f * sin(pi / len));
// 枚举所有区间左端点
for (int l = 0; l < k; l += len)
{
// 单位根的 0 次方为 1
w = cp(1, 0);
// 和递归 FFT 一样分治处理
for (int i = 0, j = len >> 1; j < len; i ++ , j ++ , w = w * t)
{
// 先用 x 和 y 表示 a[l + i] 和 w * a[l + j],然后再计算,不然会覆盖掉前面的
x = a[l + i], y = w * a[l + j];
a[l + i] = x + y, a[l + j] = x - y;
}
}
}
}
int main()
{
scanf("%d %d", &n, &m);
for (int i = 0; i <= n; i ++ ) scanf("%lf", &a[i].a);
for (int i = 0; i <= m; i ++ ) scanf("%lf", &b[i].a);
int bit = log_2(n + m);
// k 的意义和递归 FFT 中的 k 一样
k = 1 << bit + 1;
// 先将 r 数组处理好
for (int i = 0; i < k; i ++ ) r[i] = r[i >> 1] >> 1 | (i & 1) << bit;
// 然后一顿 FFT 就可以了
FFT(a, 1), FFT(b, 1);
for (int i = 0; i < k; i ++ ) a[i] = a[i] * b[i];
FFT(a, -1);
for (int i = 0; i <= n + m; i ++ ) printf("%d ", (int)(a[i].a / k + 0.5));
return 0;
}
3.5 三变二优化 这个不知道学名叫啥,瞎起个名字
设 ck=ak+ibk
那么 c2k=a2k−b2k+2iakbk
那么我们就可以把两个多项式存到一个数组里面,然后只用做两次 FFT 即可
另外,注意到上面重复调用了很多次三角函数,我们可以直接预处理出 ωin,可以减少很多调用三角函数的次数
#include <cmath>
#include <cstdio>
struct cp
{
double a, b;
cp(double x = 0, double y = 0){a = x, b = y;}
cp operator + (const cp &t) const {return cp(a + t.a, b + t.b);}
cp operator - (const cp &t) const {return cp(a - t.a, b - t.b);}
cp operator * (const cp &t) const {return cp(a * t.a - b * t.b, a * t.b + b * t.a);}
};
int log_2(int x)
{
int res = 0;
if (x & 0xffff0000) res += 16, x >>= 16;
if (x & 0xff00) res += 8, x >>= 8;
if (x & 0xf0) res += 4, x >>= 4;
if (x & 0xc) res += 2, x >>= 2;
if (x & 2) res ++ ;
return res;
}
void swap(cp &a, cp &b)
{
static cp t;
t = a, a = b, b = t;
}
const int N = 400005;
const double pi = 2 * acos(-1);
int n, m, k, r[N];
cp w[N], a[N]; // w 预处理 omega ^ k,a 存两个多项式
void FFT(bool type)
{
for (int i = 0; i < k; i ++ ) if (i < r[i]) swap(a[i], a[r[i]]);
static cp x, y;
for (int len = 2; len <= k; len <<= 1)
for (int l = 0; l < k; l += len)
for (int i = 0, j = len >> 1; j < len; i ++ , j ++ )
{
// w[k / len * i] = omega _ len ^ i
x = a[l + i], y = w[k / len * i];
if (type) y = y * a[l + j];
else y = cp(y.a, -y.b) * a[l + j];
a[l + i] = x + y, a[l + j] = x - y;
}
}
int main()
{
scanf("%d %d", &n, &m);
for (int i = 0; i <= n; i ++ ) scanf("%lf", &a[i].a);
for (int i = 0; i <= m; i ++ ) scanf("%lf", &a[i].b);
int bit = log_2(n + m);
k = 1 << bit + 1;
for (int i = 0; i < k; i ++ ) r[i] = r[i >> 1] >> 1 | (i & 1) << bit;
// 预处理 omega ^ k
w[0].a = 1;
cp t(cos(pi / k), sin(pi / k));
for (int i = 1; i < k; i ++ ) w[i] = w[i - 1] * t;
// 只做两次 FFT
FFT(1);
for (int i = 0; i < k; i ++ ) a[i] = a[i] * a[i];
FFT(0);
for (int i = 0; i <= n + m; i ++ ) printf("%d ", (int)(a[i].b / k / 2.0 + 0.5));
return 0;
}
是的这就是 FFT 的终极形式了(应该没有别的神仙优化了8?)
练习题:【模板】A*B Problem升级版(FFT快速傅里叶){:target=”_blank”}
这名字怎么这么长。。
板子题,直接贴代码
C ++ 代码
#include <cmath>
#include <cstdio>
#include <cstring>
struct cp
{
double a, b;
cp(double x = 0, double y = 0){a = x, b = y;}
cp operator + (const cp &t) const {return cp(a + t.a, b + t.b);}
cp operator - (const cp &t) const {return cp(a - t.a, b - t.b);}
cp operator * (const cp &t) const {return cp(a * t.a - b * t.b, a * t.b + b * t.a);}
};
int log_2(int x)
{
int res = 0;
if (x & 0xffff0000) res += 16, x >>= 16;
if (x & 0xff00) res += 8, x >>= 8;
if (x & 0xf0) res += 4, x >>= 4;
if (x & 0xc) res += 2, x >>= 2;
if (x & 2) res ++ ;
return res;
}
void swap(cp &a, cp &b)
{
static cp t;
t = a, a = b, b = t;
}
const int N = 4000005;
const double pi = 2 * acos(-1);
int n, m, r[N];
char f[N], g[N], c[N];
cp w[N], a[N];
void init()
{
m += n;
int bit = log_2(m);
n = 1 << bit + 1;
for (int i = 0; i < n; i ++ ) r[i] = r[i >> 1] >> 1 | (i & 1) << bit;
w[0].a = 1;
cp t(cos(pi / n), sin(pi / n));
for (int i = 1; i < n; i ++ ) w[i] = w[i - 1] * t;
}
void FFT(bool type)
{
for (int i = 0; i < n; i ++ ) if (i < r[i]) swap(a[i], a[r[i]]);
static cp x, y;
for (int len = 2; len <= n; len <<= 1)
for (int l = 0; l < n; l += len)
for (int i = 0, j = len >> 1; j < len; i ++ , j ++ )
{
x = a[l + i], y = w[n / len * i];
if (type) y = y * a[l + j];
else y = cp(y.a, -y.b) * a[l + j];
a[l + i] = x + y, a[l + j] = x - y;
}
}
int main()
{
scanf("%s\n%s", f, g);
n = strlen(f) - 1, m = strlen(g) - 1;
for (int i = n; ~i; i -- ) a[i].a = (double) (f[n - i] ^ 48);
for (int i = m; ~i; i -- ) a[i].b = (double) (g[m - i] ^ 48);
init();
FFT(1);
for (int i = 0; i < n; i ++ ) a[i] = a[i] * a[i];
FFT(0);
int t = 0;
for (int i = 0; i <= m; i ++ )
{
t += (int) (a[i].b / n / 2.0 + 0.5);
c[i] = t % 10 ^ 48, t /= 10;
}
if (t) printf("%d", t);
for (int i = m; ~i; i -- ) putchar(c[i]);
return 0;
}
在 TLE 的边缘疯狂试探
3.6 FFT 在字符串匹配中的应用
例题:KMP 字符串{:target=”_blank”}
是的这是KMP的板子题,我们现在要用 FFT 解决它
对于给定的两个字符串 s[m],p[n],
我们可以定义匹配函数 c(i,j)=s[i]−p[j]
再定义完全匹配函数 C(x)=n−1∑i=0c(x+i,i)
那么如果 s[x∼x+n−i] 与 p 完全匹配,则 C(x)=0
但是我们可以只通过 C(x) 是否为 0,直接判断 s[x x+n−i] 与 p 是否完全匹配吗?
用完全匹配函数判断一下 abc
和 cba
试试?
似乎在我们的定义下,abc
和 cba
是匹配的?
匹配函数没错,如果两个字符不同,一定会返回一个非零值
所以问题出在了我们的完全匹配函数上
在我们的定义下,只要两个字符串所包含的字符集合是相同的,不管排列如何,完全匹配函数都会返回 0
(当然,两个字符串所包含的字符集合不同,完全匹配函数也是有可能返回 0 的,比如 add
和 ccc
也是匹配的)
我们的完全匹配函数,直接将所有匹配函数的结果相加,而没有考虑正负号
对于判定正负号,我们可以给匹配函数,直接暴力地加一个绝对值上去
但是这样就只能 O(nm) 暴力算了
不难想到,我们可以将匹配函数,变为 c(i,j)=(s[i]−p[j])2
这样就不会出问题了,但是会难算一些
先将其代入完全匹配函数
C(x)=n−1∑i=0(s[x+i]−p[i])2
将平方展开
C(x)=n−1∑i=0(s[x+i]2+p[i]2−2s[x+i]p[i])
再将求和展开
C(x)=n−1∑i=0s[x+i]2+n−1∑i=0p[i]2−2n−1∑i=0s[x+i]p[i]
对于第一项,我们可以用前缀和处理;
对于第二项,是一个常数,可以预处理出来;
所以我们剩下要考虑的,就是如何快速求出第三项
n−1∑i=0s[x+i]p[i]
这么看确实没什么思路,我们令 j=n−i−1 并将其代入,得
∑i+j=n−1,0≤i<ns[x+i]p[n−j−1]
然后我们可以将 p 翻转,那么就有
∑i+j=n−1,0≤j<ns[x+i]p[j]
也就是
∑i+j=n+x−1,0≤j<ns[i]p[j]
这不就是我们熟悉的卷积形式了吗
那么我们就可以用 FFT 解决掉这项了
我们也就可以用 O(nlogn) 的复杂度,求出所有的 C 函数的值了
C ++ 代码
#include <cmath>
#include <cstdio>
// 前缀和有可能爆 int
typedef long long ll;
struct cp
{
double a, b;
cp(double x = 0, double y = 0) {a = x, b = y;}
cp operator + (const cp &t) const {return cp(a + t.a, b + t.b);}
cp operator - (const cp &t) const {return cp(a - t.a, b - t.b);}
cp operator * (const cp &t) const {return cp(a * t.a - b * t.b, a * t.b + b * t.a);}
};
int log_2(int x)
{
int res = 0;
if (x & 0xffff0000) res += 16, x >>= 16;
if (x & 0xff00) res += 8, x >>= 8;
if (x & 0xf0) res += 4, x >>= 4;
if (x & 0xc) res += 2, x >>= 2;
if (x & 2) res ++ ;
return res;
}
void swap(cp &a, cp &b)
{
static cp t;
t = a, a = b, b = t;
}
const int N = 100005;
const int M = 4000005;
const double pi = 2 * acos(-1);
int n, m, k;
int r[M], c[M]; // c 即最后求出来的 C 函数
char s[N], p[M];
ll sum, f[M]; // sum 用于预处理上面 C 函数中第二项的值。f 用于预处理 C 函数第一项的值。
cp a[M], w[M];
void init()
{
k = n + m - 2;
int bit = log_2(k);
k = 1 << bit + 1;
for (int i = 0; i < k; i ++ ) r[i] = r[i >> 1] >> 1 | (i & 1) << bit;
const cp t(cos(pi / k), sin(pi / k)); (*w).a = 1;
for (int i = 1; i < k; i ++ ) w[i] = w[i - 1] * t;
}
void FFT(bool type)
{
for (int i = 0; i < k; i ++ ) if (i < r[i]) swap(a[i], a[r[i]]);
static cp x, y;
for (int len = 2; len <= k; len <<= 1)
for (int l = 0; l < k; l += len)
for (int i = 0, j = len >> 1; j < len; i ++ , j ++ )
{
x = a[l + i], y = w[k / len * i];
if (type) y = y * a[l + j];
else y = cp(y.a, -y.b) * a[l + j];
a[l + i] = x + y, a[l + j] = x - y;
}
}
int main()
{
scanf("%d\n%s\n%d\n%s", &n, p, &m, s);
// 预处理 sum 和 f
for (int i = 0; i < n; i ++ ) sum += p[i] * 1ll * p[i];
for (int i = 1; i <= m; i ++ ) f[i] = f[i - 1] + s[i - 1] * 1ll * s[i - 1];
// 将 s 和 p 放到 a 中做 FFT,其中 p 要反着放
for (int i = 0; i < n; i ++ ) a[i].a = p[n - i - 1];
for (int i = 0; i < m; i ++ ) a[i].b = s[i];
// 一顿 FFT
init();
FFT(1);
for (int i = 0; i < k; i ++ ) a[i] = a[i] * a[i];
FFT(0);
// 求出 c
for (int i = 1; i < n + m; i ++ ) c[i] = (int) (a[i - 1].b / k / 2.0 + 0.5);
for (int i = n; i <= m; i ++ )
// 如果 C(i) 为 0,则输出 i - n
if (f[i] - f[i - n] + sum - 2 * c[i] == 0)
printf("%d ", i - n);
return 0;
}
很不幸,最后一个点 TLE 了。。
但是如果你是卡常数带师,应该也可以卡常过掉这题
卡常后的 AC 代码
(不保证每次提交都能过,AC 看运气)
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC target("avx,sse2,sse3,sse4,mmx")
#pragma GCC optimize("Ofast")
#pragma GCC optimize("inline")
#pragma GCC optimize("-fgcse")
#pragma GCC optimize("-fgcse-lm")
#pragma GCC optimize("-fipa-sra")
#pragma GCC optimize("-ftree-pre")
#pragma GCC optimize("-ftree-vrp")
#pragma GCC optimize("-fpeephole2")
#pragma GCC optimize("-ffast-math")
#pragma GCC optimize("-fsched-spec")
#pragma GCC optimize("unroll-loops")
#pragma GCC optimize("-falign-jumps")
#pragma GCC optimize("-falign-loops")
#pragma GCC optimize("-falign-labels")
#pragma GCC optimize("-fdevirtualize")
#pragma GCC optimize("-fcaller-saves")
#pragma GCC optimize("-fcrossjumping")
#pragma GCC optimize("-fthread-jumps")
#pragma GCC optimize("-funroll-loops")
#pragma GCC optimize("-fwhole-program")
#pragma GCC optimize("-freorder-blocks")
#pragma GCC optimize("-fschedule-insns")
#pragma GCC optimize("inline-functions")
#pragma GCC optimize("-ftree-tail-merge")
#pragma GCC optimize("-fschedule-insns2")
#pragma GCC optimize("-fstrict-aliasing")
#pragma GCC optimize("-fstrict-overflow")
#pragma GCC optimize("-falign-functions")
#pragma GCC optimize("-fcse-skip-blocks")
#pragma GCC optimize("-fcse-follow-jumps")
#pragma GCC optimize("-fsched-interblock")
#pragma GCC optimize("-fpartial-inlining")
#pragma GCC optimize("no-stack-protector")
#pragma GCC optimize("-freorder-functions")
#pragma GCC optimize("-findirect-inlining")
#pragma GCC optimize("-fhoist-adjacent-loads")
#pragma GCC optimize("-frerun-cse-after-loop")
#pragma GCC optimize("inline-small-functions")
#pragma GCC optimize("-finline-small-functions")
#pragma GCC optimize("-ftree-switch-conversion")
#pragma GCC optimize("-foptimize-sibling-calls")
#pragma GCC optimize("-fexpensive-optimizations")
#pragma GCC optimize("-funsafe-loop-optimizations")
#pragma GCC optimize("inline-functions-called-once")
#pragma GCC optimize("-fdelete-null-pointer-checks")
#include <cmath>
#include <cstdio>
typedef long long ll;
struct cp
{
double a, b;
inline cp(double x = 0, double y = 0) {a = x, b = y;}
inline cp operator + (const cp &t) const {return cp(a + t.a, b + t.b);}
inline cp operator - (const cp &t) const {return cp(a - t.a, b - t.b);}
inline cp operator * (const cp &t) const {return cp(a * t.a - b * t.b, a * t.b + b * t.a);}
};
inline int log_2(register int x)
{
int res = 0;
if (x & 0xffff0000) res += 16, x >>= 16;
if (x & 0xff00) res += 8, x >>= 8;
if (x & 0xf0) res += 4, x >>= 4;
if (x & 0xc) res += 2, x >>= 2;
if (x & 2) res ++ ;
return res;
}
inline void swap(cp &a, cp &b)
{
static cp t;
t = a, a = b, b = t;
}
const int N = 100005, M = 4000005;
const double pi = 2 * acos(-1);
int n, m, k;
int r[M];
char s[M >> 2], p[N];
ll sum, f[M >> 2];
cp a[M], w[M];
inline void init()
{
k = n + m - 2;
register int bit = log_2(k), *i, *j;
for (i = j = r, k = 1 << bit + 1, bit = k >> 1; i < r + k; i += 2, j ++ )
*(i + 1) = (*i = *j >> 1) | bit;
cp t(cos(pi / k), sin(pi / k)), *pt; (*w).a = 1;
for (pt = w + 1; pt < w + k; pt ++ ) *pt = *(pt - 1) * t;
}
inline void FFT(bool type)
{
register int i, len, l;
for (i = 0; i < k; i ++ ) if (i > r[i]) swap(a[i], a[r[i]]);
static cp *x, *y, v;
for (len = 2; len <= k; len <<= 1)
for (l = 0; l < k; l += len)
for (i = len >> 1, x = a + l, y = x + i; i < len; x ++ , y ++ , i ++ )
{
v = w[k / len * (i - (len >> 1))];
if (type) v.b = -v.b;
v = v * *y;
*y = *x - v, *x = *x + v;
}
}
int main()
{
scanf("%d\n%s\n%d\n%s", &n, p, &m, s);
for (int i = 0; i < n; i ++ ) p[i] -= 96;
for (int i = 0; i < m; i ++ ) s[i] -= 96;
register char *c; register ll *it, *jt; cp *pt;
for (c = p; *c; c ++ ) sum += *c * 1ll * *c;
for (it = f, c = s; *c; c ++ , it ++ ) *(it + 1) = *it + *c * 1ll * *c;
for (pt = a, c = p + n - 1; c >= p; c -- , pt ++ ) pt->a = *c;
for (pt = a, c = s; *c; c ++ , pt ++ ) pt->b = *c;
init();
FFT(0);
for (pt = a; pt < a + k; pt ++ ) *pt = *pt * *pt;
FFT(1);
k <<= 1;
for (pt = a + n - 1, it = f + n, jt = f; it <= f + m; pt ++ , it ++ , jt ++ )
if (((ll) (pt->b / k + 0.5)) << 1 == *it - *jt + sum)
printf("%d ", it - f - n);
return 0;
}
是的我卡了半天常并加了火车头才过
练习题:残缺的字符串{:target=”_blank”}
这题可以使用 FFT 做字符串匹配
刚才那题中当且仅当两个字符相同才算匹配,
我们定义了匹配函数 c(i,j)=si−pj
在这题中,如果 s[i]= *
,或者 p[j]= ‘*’,也算这两个字符是匹配的
那么不难想到,若 s[i]= *
,则令 s[i]=0,若 p[i]= *
,则令 p[i]=0
然后定义匹配函数 c(i,j)=(s[i]−p[j])2s[i]p[j] 即可
完全匹配函数 C(x)=n−1∑i=0c(x+i,i) 不变,但是要重新推下式子
C(x)= n−1∑i=0c(x+i,i)2 = n−1∑i=0(s[x+i]−p[i])2s[x+i]p[i] = n−1∑i=0(s[x+i]3p[i]+s[x+i]p[i]3−2s[x+i]2p[i]2) = n−1∑i=0s[x+i]3p[i]+n−1∑i=0s[x+i]p[i]3−2n−1∑i=0s[x+i]2p[i]2
同样,代入 j=n−i−1,翻转 p,得
C(x)=(∑i+j=n−1s[x+i]3p[j])+(∑i+j=n−1s[x+i]p[j]3)−2(∑i+j=n−1s[x+i]2p[j]2)
(脑补一下 0≤j<n,求和符号下面写不下了。。)也就是
C(x)=(∑i+j=n+x−1s[i]3p[j])+(∑i+j=n+x−1s[i]p[j]3)−2(∑i+j=n+x−1s[i]2p[j]2)
然后我们发现我们并不能预处理任何东西。。要做三遍 FFT。。
C ++ 代码
FFT 好像无 O2
不可过,这题应该是故意把 FFT 卡掉了的。
#include <cmath>
#include <cstdio>
using namespace std;
#define double long double
struct cp
{
double a, b;
cp operator + (const cp &t) const {return (cp){a + t.a, b + t.b};}
cp operator - (const cp &t) const {return (cp){a - t.a, b - t.b};}
cp operator * (const cp &t) const {return (cp){a * t.a - b * t.b, a * t.b + b * t.a};}
};
int lg2(int x)
{
int res = 0;
if (x & 0xffff0000) res += 16, x >>= 16;
if (x & 0xff00) res += 8, x >>= 8;
if (x & 0xf0) res += 4, x >>= 4;
if (x & 0xc) res += 2, x >>= 2;
if (x & 2) res ++ ;
return res;
}
void swap(cp &a, cp &b)
{
static cp t;
t = a, a = b, b = t;
}
const int N = 300005, M = N << 2;
const double pi = 2 * acos(-1);
int n, m, k;
int r[M];
int res[N], cnt;
char s[N], p[N];
cp a[M], c[M];
void init()
{
k = n + m - 2;
int bit = lg2(k); k = 1 << bit + 1;
for (int i = 0; i < k; i ++ ) r[i] = r[i >> 1] >> 1 | (i & 1) << bit;
}
void FFT(bool type)
{
for (int i = 0; i < k; i ++ ) if (i > r[i]) swap(a[i], a[r[i]]);
static cp *x, *y, v, w1, w;
for (register int len = 2, l, i; len <= k; len <<= 1)
{
w1 = (cp){cos(pi / len), type ? sin(pi / len) : -sin(pi / len)};
for (l = 0; l < k; l += len)
{
w = (cp){1, 0};
for (i = len >> 1, x = a + l, y = x + i; i < len; x ++ , y ++ , i ++ , w = w * w1)
v = w * *y, *y = *x - v, *x = *x + v;
}
}
}
int main()
{
scanf("%d%d\n%s\n%s", &n, &m, p, s);
for (int i = 0; i < n; i ++ ) p[i] = p[i] == '*' ? 0 : p[i] ^ 96;
for (int i = 0; i < m; i ++ ) s[i] = s[i] == '*' ? 0 : s[i] ^ 96;
for (int i = 0, j = n, t; -- j > i; i ++ ) t = p[j], p[j] = p[i], p[i] = t;
init();
for (int i = 0; i < n; i ++ ) a[i].a = p[i];
for (int i = 0; i < m; i ++ ) a[i].b = s[i] * 1.0 * s[i] * s[i];
FFT(1);
for (int i = 0; i < k; i ++ ) c[i] = c[i] + a[i] * a[i], a[i].a = a[i].b = 0;
for (int i = 0; i < n; i ++ ) a[i].a = p[i] * 1.0 * p[i] * p[i];
for (int i = 0; i < m; i ++ ) a[i].b = s[i];
FFT(1);
for (int i = 0; i < k; i ++ ) c[i] = c[i] + a[i] * a[i], a[i].a = a[i].b = 0;
for (int i = 0; i < n; i ++ ) a[i].a = p[i] * 1.0 * p[i];
for (int i = 0; i < m; i ++ ) a[i].b = s[i] * 1.0 * s[i];
FFT(1);
for (int i = 0; i < k; i ++ ) a[i] = a[i] * a[i], a[i] = c[i] - a[i] - a[i];
FFT(0);
for (int i = n - 1; i < m; i ++ )
if (a[i].b < k)
res[cnt ++ ] = i - n + 2;
printf("%d\n", cnt);
for (int i = 0; i < cnt; i ++ ) printf("%d ", res[i]);
return 0;
}
一些经典题:
[AH2017/HNOI2017]礼物{:target=”_blank”}
万径人踪灭{:target=”_blank”}
[ZJOI2014]力{:target=”blank”}
公式挂了
写的很清楚,但有个问题,在luogu加了那个三变二优化后跑的更慢了是为什么。。。。
同时催更NTT,FWT
坐等NTT、FWT
可能要半年后了。。抱歉QAQ
没事
反正半年后我也不一定就会了(doge)卡常宗师orz
Orz您来了
orz
把O2改成Ofast
个人认为
Ofast
作用不是很大,没有O2
效果好,开了Ofast
之后还是会TLE
开O2Wa说明你代码有Ub的地方
FFT
会有精度丢失,NTT
就不会了哦哦,好吧
卡常变量加register
卡常的那段代码加了
register
了,还是过不去%%%%%
建议再讲一下蝴蝶变换吧……感觉好晕
蝴蝶变换是啥?
emm我百度了下,那个蝴蝶变换就是在递归改迭代的时候,把我写的递归 FFT 中的那个
t
数组去掉。。但实际上改成迭代的时候,我代码里面已经把那个去掉了。。
so,相当于是已经做了蝴蝶变换,只是没提出来这个是蝴蝶变换鹅已
而且我感觉,要是提出来那个是什么蝴蝶变换,然后专门去写一下,可能好多人一看到蝴蝶变换这个名字就被吓跑了
哦哦原来如此
我谔谔
完全看不懂
哪没看懂啊QAQ?可能是我写的不好
建议先看一下上一章QAQ
看不懂