浅谈 FFT (2)

作者： 作者的头像   垫底抽风 ,  2020-08-24 09:27:39 ,  所有人可见 ,  阅读 2283

32


42

上一章：浅谈 FFT (1){:target=”_blank”}

3.4 多项式乘法的迭代实现

递归实现 $\text{FFT}$ 的代码，看着很优美，也的确能过掉这道板子题{:target=”_blank”}

但是，常数太大

面对 $10^5$ 以上的数据，有很大几率会 $\text{TLE}$，需要优化

观察下我们 $\text{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 \sim 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(\log \log n)$ 求 $\log _ 2(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 三变二优化 这个不知道学名叫啥，瞎起个名字

设 $c _ k = a _ k + i b _ k$

那么 $c _ k ^ 2 = a _ k ^ 2 - b _ k ^ 2 + 2 i a _ k b _ k$

那么我们就可以把两个多项式存到一个数组里面，然后只用做两次 $\text{FFT}$ 即可

另外，注意到上面重复调用了很多次三角函数，我们可以直接预处理出 $\omega _ n ^ i$，可以减少很多调用三角函数的次数

#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;
}

是的这就是 $\text{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的板子题，我们现在要用 $\text{FFT}$ 解决它

对于给定的两个字符串 $s[m], p[n]$，

我们可以定义匹配函数 $c(i, j) = s[i] - p[j]$

再定义完全匹配函数 $ C(x) = \begin {aligned} \sum _ {i = 0} ^ {n - 1} c(x + i, i) \end {aligned} $

那么如果 $s[x \sim 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) = \begin {aligned} \sum _ {i = 0} ^ {n - 1} (s[x + i] - p[i]) ^ 2 \end {aligned} $$

将平方展开

$$ C(x) = \begin {aligned} \sum _ {i = 0} ^ {n - 1} (s[x + i] ^ 2 + p[i] ^ 2 - 2 s[x + i] p[i]) \end {aligned} $$

再将求和展开

$$ C(x) = \begin {aligned} \sum _ {i = 0} ^ {n - 1} s[x + i] ^ 2 + \sum _ {i = 0} ^ {n - 1} p[i] ^ 2 - 2 \sum _ {i = 0} ^ {n - 1} s[x + i] p[i] \end {aligned} $$

对于第一项，我们可以用前缀和处理；

对于第二项，是一个常数，可以预处理出来；

所以我们剩下要考虑的，就是如何快速求出第三项

$$ \begin {aligned} \sum _ {i = 0} ^ {n - 1} s[x + i] p[i] \end {aligned} $$

这么看确实没什么思路，我们令 $j = n - i - 1$ 并将其代入，得

$$ \begin {aligned} \sum _ {i + j = n - 1, 0 \leq i < n} s[x + i] p[n - j - 1] \end {aligned} $$

然后我们可以将 $p$ 翻转，那么就有

$$ \begin {aligned} \sum _ {i + j = n - 1, 0 \leq j < n} s[x + i] p[j] \end {aligned} $$

也就是

$$ \begin {aligned} \sum _ {i + j = n + x - 1, 0 \leq j < n} s[i] p[j] \end {aligned} $$

这不就是我们熟悉的卷积形式了吗

那么我们就可以用 $\text{FFT}$ 解决掉这项了

我们也就可以用 $O(n \log n)$ 的复杂度，求出所有的 $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;
}

很不幸，最后一个点 $\text{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”}

这题可以使用 $\text{FFT}$ 做字符串匹配

刚才那题中当且仅当两个字符相同才算匹配，

我们定义了匹配函数 $c(i, j) = s _ i - p _ j$

在这题中，如果 $s[i] = $ *，或者 $p[j] = $ ‘*’，也算这两个字符是匹配的

那么不难想到，若 $s[i] = $ *，则令 $s[i] = 0$，若 $p[i] = $ *，则令 $p[i] = 0$

然后定义匹配函数 $c(i, j) = (s[i] - p[j]) ^ 2 s[i] p[j]$ 即可

完全匹配函数 $ C(x) = \begin {aligned} \sum _ {i = 0} ^ {n - 1} c(x + i, i) \end {aligned} $ 不变，但是要重新推下式子

$$ \begin {aligned} C(x) = &\ \sum _ {i = 0} ^ {n - 1} c (x + i, i) ^ 2 \\\ \\\ = &\ \sum _ {i = 0} ^ {n - 1} (s[x + i] - p[i]) ^ 2 s[x + i] p[i] \\\ \\\ = &\ \sum _ {i = 0} ^ {n - 1} (s[x + i] ^ 3 p[i] + s[x + i] p[i] ^ 3 - 2 s[x + i] ^ 2 p[i] ^ 2) \\\ \\\ = &\ \sum _ {i = 0} ^ {n - 1} s[x + i] ^ 3 p[i] + \sum _ {i = 0} ^ {n - 1} s[x + i] p[i] ^ 3 - 2 \sum _ {i = 0} ^ {n - 1} s[x + i] ^ 2 p[i] ^ 2 \\\ \\\ \end {aligned} $$

同样，代入 $j = n - i - 1$，翻转 $p$，得

$$ C(x) = \begin {aligned} (\sum _ {i + j = n - 1} s[x + i] ^ 3 p[j]) + (\sum _ {i + j = n - 1} s[x + i] p[j] ^ 3) - 2 (\sum _ {i + j = n - 1} s[x + i] ^ 2 p[j] ^ 2) \\\ \end {aligned} $$

（脑补一下 $0 \leq j < n$，求和符号下面写不下了。。）也就是

$$ C(x) = \begin {aligned} (\sum _ {i + j = n + x - 1} s[i] ^ 3 p[j]) + (\sum _ {i + j = n + x - 1} s[i] p[j] ^ 3) - 2 (\sum _ {i + j = n + x - 1} s[i] ^ 2 p[j] ^ 2) \\\ \end {aligned} $$

然后我们发现我们并不能预处理任何东西。。要做三遍 $\text{FFT}$。。

C ++ 代码

$\text{FFT}$ 好像无 O2 不可过，这题应该是故意把 $\text{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”}

$$ $$

25 评论

用户头像
Conan15   4个月前         踩      回复

公式挂了

用户头像
_Undertaker_   2022-07-26 17:12         踩      回复

写的很清楚，但有个问题，在luogu加了那个三变二优化后跑的更慢了是为什么。。。。

同时催更NTT，FWT

用户头像
RingweEH   2020-09-17 19:19         踩      回复

坐等NTT、FWT

用户头像
垫底抽风   2020-09-18 19:39         踩      回复

可能要半年后了。。抱歉QAQ

用户头像
RingweEH   2020-09-19 07:12    回复了 垫底抽风 的评论         踩      回复

没事反正半年后我也不一定就会了（doge）

用户头像
vcrlwpx   2020-08-26 14:31         踩      回复

卡常宗师orz

用户头像
垫底抽风   2020-08-26 14:44         踩      回复

Orz您来了

用户头像
itdef   2020-08-24 19:52         踩      回复

orz

用户头像
Hicode002   2020-08-24 17:19         踩      回复

把O2改成Ofast

用户头像
垫底抽风   2020-08-24 17:31         踩      回复

个人认为 Ofast 作用不是很大，没有 O2 效果好，开了 Ofast 之后还是会 TLE

用户头像
Hicode002   2020-08-24 17:17         踩      回复

开O2Wa说明你代码有Ub的地方

用户头像
垫底抽风   2020-08-24 17:30         踩      回复

FFT 会有精度丢失，NTT 就不会了

用户头像
Hicode002   2020-08-24 17:34    回复了 垫底抽风 的评论         踩      回复

哦哦，好吧

用户头像
Hicode002   2020-08-24 17:17         踩      回复

卡常变量加register

用户头像
垫底抽风   2020-08-24 17:30         踩      回复

卡常的那段代码加了 register 了，还是过不去

用户头像
cht   2020-08-24 14:28         踩      回复

%%%%%
建议再讲一下蝴蝶变换吧……感觉好晕

用户头像
垫底抽风   2020-08-24 16:17         踩      回复

蝴蝶变换是啥？

用户头像
垫底抽风   2020-08-24 16:27         踩      回复

emm我百度了下，那个蝴蝶变换就是在递归改迭代的时候，把我写的递归 FFT 中的那个 t 数组去掉。。
但实际上改成迭代的时候，我代码里面已经把那个去掉了。。
so，相当于是已经做了蝴蝶变换，只是没提出来这个是蝴蝶变换鹅已

用户头像
垫底抽风   2020-08-24 16:28         踩      回复

而且我感觉，要是提出来那个是什么蝴蝶变换，然后专门去写一下，可能好多人一看到蝴蝶变换这个名字就被吓跑了

用户头像
cht   2020-08-24 18:16    回复了 垫底抽风 的评论         踩      回复

哦哦原来如此

用户头像
cht   2020-08-24 18:16    回复了 垫底抽风 的评论         踩      回复

我谔谔

用户头像
Kue   2020-08-24 09:51         踩      回复

完全看不懂

用户头像
垫底抽风   2020-08-24 10:02         踩      回复

哪没看懂啊QAQ？可能是我写的不好

用户头像
垫底抽风   2020-08-24 10:04         踩      回复

建议先看一下上一章QAQ

用户头像
dys   2020-09-02 15:55    回复了 垫底抽风 的评论         踩      回复

看不懂