[#LOJ6388] [THUPC2018] 赛艇 快速傅里叶变换 FFT

题目链接
今天$ShichengXiao$讲了卷积相关..这是$PPT$上的题

首先这个题 大意就是给你一个大的和一个小的$0/1$矩阵，然后你把小的对齐到大的上面，统计有多少种方案不会出现有两个$1$重叠

这个题感觉很套路(巧妙)啊 如果两项是1那这两项乘起来就是$1$

我们就可以考虑把它转换为乘法形式，就定义一个点的权值就是把小矩阵的右下角固定在这个点，每一项对应相乘求和得到的值，显然只有这个值是$0$的时候这个点才属于合法的方案，二维不好做我们把它拆成一维的，就是把第二行接在第一行后面，第三行接在第二行后面…那么$f_x = \sum_{i = 0}^{n - 1}a_ib_i$不难发现这是一个卷积的形式，翻转$b$序列，最后统计答案就好了

还是讲详细一点，小矩阵要补成和大矩阵一样大，空的地方补$0$就好了，这样子模拟一下运算过程，发现每个点移动的时候大矩阵的对应项就乘上了小矩阵补的$0$，所以是不会影响答案的

复杂度$O(nmlog(nm))$

hint: $FFT$复数类还是手写吧…我用自带的$TLE$了

其实$NTT$要比$FFT$快的多，可是还要有$LSTete$才更快啊

Codes

#include<bits/stdc++.h>
#include<tr1/unordered_map>

#define map tr1::unordered_map
#define id(i, j) ((i) * m + (j))

using namespace std;

template<class T>inline bool chkmin(T &_, T __) {return _ > __ ? _ = __, 1 : 0;}
template<class T>inline bool chkmax(T &_, T __) {return _ < __ ? _ = __, 1 : 0;}

typedef long double db;

struct Complex {
    db x, y;
    Complex operator + (const Complex &T) {
        return Complex{x + T.x, y + T.y};
    }
    Complex operator - (const Complex &T) {
        return Complex{x - T.x, y - T.y};
    }
    Complex operator * (const Complex &T) {
        return Complex{x * T.x - y * T.y, x * T.y + y * T.x};
    }
};

const int N = 1500 + 10, M = 1e7 + 10;
const db Pi = acos(-1);
int a[N][N], mp[N][N];
int n, m, k, hb, lb;

map<int, map<int, int> >b;
Complex A[M], B[M];

void Init() {
    int hmin, hmax, lmin, lmax; char c; 
    hb = lb = hmin = hmax = lmin = lmax = b[1][1] = 1;

    scanf("%d%d%d", &n, &m, &k);

    for(int i = 1; i <= n; ++ i) 
        for(int j = 1; j <= m; ++ j)
            scanf("%1d", a[i] + j);

    for(int i = 1; i <= k; ++ i) {
        cin >> c;
        if(c == 'w') {b[-- hb][lb] = 1; chkmin(hmin, hb);}
        if(c == 'a') {b[hb][-- lb] = 1; chkmin(lmin, lb);}
        if(c == 's') {b[++ hb][lb] = 1; chkmax(hmax, hb);}
        if(c == 'd') {b[hb][++ lb] = 1; chkmax(lmax, lb);}
    }


    for(int i = hmin; i <= hmax; ++ i)
        for(int j = lmin; j <= lmax; ++ j) 
            mp[i - hmin + 1][j - lmin + 1] = b[i][j];

    hb = hmax - hmin + 1, lb = lmax - lmin + 1;

    for(int i = 0; i < n; ++ i)
        for(int j = 0; j < m; ++ j) 
            A[id(i, j)].x = a[i + 1][j + 1], B[id(i, j)].x = mp[i + 1][j + 1];

    reverse(B, B + n * m);
}

void FFT_Init(int limit, Complex *a) {
    int k = log2(limit);
    for(int j = 0; j < limit; ++ j) {
        int t = 0;
        for(int i = 0; i < k; ++ i)
            if((1 << i) & j)
                t |= 1 << (k - i - 1);
        if(j < t) swap(a[t], a[j]);
    }
}

void FFT(int n, Complex *a, int fh) {
    FFT_Init(n, a);
    for(int limit = 2; limit <= n; limit <<= 1) {
        double theta = 2.0 * Pi / limit;
        Complex Wn = Complex{cos(theta), sin(theta) * fh}, W = Complex{1, 0};
        for(int j = 0; j < n; j += limit, W = Complex{1, 0}) 
            for(int i = j; i < j + (limit >> 1); ++ i, W = W * Wn) {
                Complex a1 = a[i], a2 = a[i + (limit >> 1)] * W;
                a[i] = a1 + a2, a[i + (limit >> 1)] = a1 - a2;
            }
    }
    if(fh == -1) 
        for(int i = 0; i < n; ++ i)
            a[i].x /= n;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("6388.in", "r", stdin);
    freopen("6388.out", "w", stdout);
#endif
    int res = 0, limit = 1; Init();
    while(limit <= (n * m * 2)) limit <<= 1;
    //for(int i = 0; i < limit; ++ i)
    //  printf("%d ", int(A[i].real()));
    //puts("");
    //for(int i = 0; i < limit; ++ i)
    //  printf("%d ", int(B[i].real()));
    //return 0;
    FFT(limit, A, 1), FFT(limit, B, 1); 
    for(int i = 0; i < limit; ++ i) A[i] = B[i] * A[i];
    FFT(limit, A, -1);
    //for(int i = 0; i < limit; ++ i)
    //  printf("%d ", int(A[i].real() + 0.5));
    for(int i = hb - 1; i < n; ++ i)
        for(int j = lb - 1; j < m; ++ j) {
            if(A[id(i, j) - id(hb - 1, lb - 1) + n * m - 1].x < 0.5)
                ++ res;
        }
    printf("%d\n", res);
    cerr << 1.0 * clock() / CLOCKS_PER_SEC << endl;
    return 0;
}