两道FFT的应用题

最新推荐文章于 2020-07-29 13:05:30 发布

OIljt12138

最新推荐文章于 2020-07-29 13:05:30 发布

阅读量394

点赞数

分类专栏： ---数学问题--- 高等数学文章标签： fft 数学

本文链接：https://blog.csdn.net/OIljt12138/article/details/65702989

版权

---数学问题--- 同时被 2 个专栏收录

15 篇文章 2 订阅

订阅专栏

高等数学

11 篇文章 0 订阅

订阅专栏

BZOJ 快速傅里叶之二

题意

计算：

C [k] = \sum k \leq i < n (a [i] * b [i - k])

$C[k]=\sum_{k\le i<n}(a[i]*b[i-k])$

思路

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

对于：

\sum k \leq i < n (a [i] * b [i - k])

$\sum_{k\le i<n}(a[i]*b[i-k])$

将b数组倒置，即 $b[i] = b'[n-1-i]$ ，原式变为：

c [k] = \sum 0 \leq i < n, 0 \leq n + k - 1 - i < n (a [i] * b' [n + k - 1 - i])

$c[k] = \sum_{0\le i<n, 0\le n+k-1-i<n}(a[i]*b'[n+k-1-i])$

化简得到

c [k] = \sum k \leq i \leq n + k - 1 (a [i] * b' [n + k - 1 - i])

$c[k] = \sum_{k\le i\le n+k-1}(a[i]*b'[n+k-1-i])$

我们想得到形如

c' [k] = \sum 0 \leq i \leq k a [i] * b' [k - i]

$c'[k] = \sum_{0\le i\le k}a[i]*b'[k-i]$

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

c' [k + n - 1] = \sum i a [i] * b' [k + n - 1 - i]

$c'[k+n-1] = \sum_i a[i]*b'[k+n-1-i]$

条件是

0 \leq k + n - 1 - i < n ⟺ k \leq i \leq n + k - 1

$0\le k+n-1-i<n\iff k\le i \le n+k-1$

所以

c [k] = c' [k + n - 1] = \sum k \leq i \leq n + k - 1 a [i] * b' [k + n - 1 - i]

$c[k] = c'[k+n-1] = \sum_{k\le i\le n+k-1} a[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的结果，相加即为关于对称轴两边对称的总字符数 $A_i$ ，则 $2^{A_i}$ 就是关于其对称的字符串总数。

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

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