两道FFT的应用题

BZOJ 快速傅里叶之二

题意

计算：

C[k]=∑k≤i<n(a[i]∗b[i−k])

思路

正好看到《具体数学》上处理和式的Tricks，虽然热身题也不会做，但碾OI题还是很稳的(Orz神犇高教授)…

对于：

∑k≤i<n(a[i]∗b[i−k])

将b数组倒置，即b[i]=b′[n−1−i]，原式变为：

c[k]=∑0≤i<n,0≤n+k−1−i<n(a[i]∗b′[n+k−1−i])

化简得到

c[k]=∑k≤i≤n+k−1(a[i]∗b′[n+k−1−i])

我们想得到形如

c′[k]=∑0≤i≤ka[i]∗b′[k−i]

的形式。只需要将后式中的k用k+n−1替换：

c′[k+n−1]=∑ia[i]∗b′[k+n−1−i]

条件是

0≤k+n−1−i<n⟺k≤i≤n+k−1

所以

c[k]=c′[k+n−1]=∑k≤i≤n+k−1a[i]∗b′[k+n−1−i]

然后用卷积做就好了。

Code

#include <bits/stdc++.h>
using namespace std;

typedef complex<double> Complex;
const int MAXN = 300005;
int rev[MAXN];
Complex A[MAXN];
Complex a[MAXN], b[MAXN], c[MAXN];
int n;

void fft(Complex a[], int n, int flag)
{
    int lgn = int(log2(n)+0.01);
    rev[0] = 0;
    for (int i = 1; i < n; i++)
        rev[i] = (rev[i>>1]>>1)|((i&1)<<(lgn-1));
    for (int i = 0; i < n; i++)
        A[rev[i]] = a[i];
    Complex u, t;
    for (int k = 2; k <= n; k <<= 1) {
        Complex dw = Complex(cos(2*M_PI/k), flag*sin(2*M_PI/k));
        for (int i = 0; i < n; i += k) {
            Complex w = 1;
            for (int j = 0; j < k>>1; j++) {
                u = A[i+j], t = w*A[i+j+(k>>1)];
                A[i+j] = u+t, A[i+j+(k>>1)] = u-t;
                w *= dw;
            }
        }
    }
    for (int i = 0; i < n; i++)
        a[i] = A[i]/Complex(flag==1?1:n,0);
}

int main()
{
    scanf("%d", &n);
    for (int i = 0; i < n; i++) {
        double x, y; scanf("%lf%lf", &x, &y);
        a[i] = Complex(x, 0), b[n-i-1] = Complex(y, 0);
    }
    int nn = 1;
    while (nn < n*2) nn <<= 1;
    fft(a, nn, 1), fft(b, nn, 1);
    for (int i = 0; i < nn; i++) c[i] = a[i]*b[i];
    fft(c, nn, -1);
    for (int i = 0; i < n; i++) printf("%d\n", int(c[i+n-1].real()+0.01));
    return 0;
}

BZOJ3160 万境人踪灭

题意

给定一个ab串，求所有不连续回文子序列的数量和。

思路

由于卷积可以处理关于一个对称轴两侧对称的字符总数，因此将串中的a、b变为0、1，自己和自己做卷积，就可以求出所有对称轴下b所对应数对称的字符总数。然后把a、b变为1、0，求出a的结果，相加即为关于对称轴两边对称的总字符数Ai，则2Ai就是关于其对称的字符串总数。

那么如何去除连续的呢？跑一遍Manacher就好了。复杂度O(nlgn)

Code

#include <bits/stdc++.h>
using namespace std;

const int MAXN = 300005;
typedef complex<double> Complex;
int rev[MAXN];
Complex A[MAXN];
const long long mod = 1000000007ll;

void fft(Complex a[], int n, int flag)
{
    rev[0] = 0;
    int lgn = int(log2(n)+0.01);
    for (int i = 1; i < n; i++)
        rev[i] = (rev[i>>1]>>1)|((i&1)<<(lgn-1));
    for (int i = 0; i < n; i++)
        A[rev[i]] = a[i];
    Complex u, t;
    for (int k = 2; k <= n; k <<= 1) {
        Complex dw = Complex(cos(2*M_PI/k), sin(flag*2*M_PI/k));
        for (int i = 0; i < n; i += k) {
            Complex w = 1;
            for (int j = 0; j < k>>1; j++) {
                u = A[i+j], t = A[i+j+(k>>1)]*w;
                A[i+j] = u+t, A[i+j+(k>>1)] = u-t;
                w *= dw;
            }
        }
    }
    for (int i = 0; i < n; i++)
        a[i] = A[i]/(flag == 1?1:Complex(n,0));
}

int p[MAXN];
long long manacher(char str[], int n)
{
    str[0] = '^', str[++n] = '#';
    int id = 0, mx = 0;
    long long ans = 0;
    for (int i = 1; i <= n; i++) {
        if (mx >= i) p[i] = min(mx-i, p[id-(i-id)]);
        else p[i] = 0;
        while (str[i+p[i]+1] == str[i-p[i]-1]) p[i]++;
        if (i+p[i] > mx) id = i, mx = i+p[i];
        (ans += p[i]+1) %= mod;
    }
    ans--;
    id = mx = 0;
    for (int i = 1; i < n; i++) {
        if (mx >= i) p[i] = min(mx-i, p[id-(i-id)]);
        else p[i] = 0;
        while (str[i+p[i]+1] == str[i-p[i]]) p[i]++;
        if (i+p[i] > mx) id = i, mx = i+p[i];
        (ans += p[i]) %= mod;
    }
    return ans;
}

long long power(int a, int n)
{
    if (n == 0) return 1;
    long long p = power(a, n>>1);
    (p *= p) %= mod;
    if (n&1) (p *= a) %= mod;
    return p;
}

char str[MAXN];
Complex a[MAXN], c[MAXN];
int cp[MAXN];

int main()
{
    scanf("%s", str+1);
    int len = strlen(str+1), n;
    for (n = 1; n < (len+1)*2; n<<=1);
    long long sub = manacher(str, len), ans = 0;
    for (int i = 1; i <= len; i++) a[i] = str[i] == 'a';
    fft(a, n, 1);
    for (int i = 0; i < n; i++) a[i] *= a[i];
    fft(a, n, -1);
    for (int i = 0; i < n; i++) cp[i] = int(a[i].real()+0.1), a[i] = 0;
    for (int i = 1; i <= len; i++) a[i] = str[i] == 'b';
    fft(a, n, 1);
    for (int i = 0; i < n; i++) a[i] *= a[i];
    fft(a, n, -1);
    for (int i = 0; i < n; i++) cp[i] += int(a[i].real()+0.1);
    for (int i = 0; i < n; i++) (ans += power(2, (cp[i]+1)/2)-1) %= mod;
    cout << ((ans-sub)%mod+mod)%mod << endl;
    return 0;
}