FFT模板使用

一、FFT模板（来自邝斌模板）

struct Complex{
    double x,y;
    Complex(double _x = 0.0, double _y = 0.0){
        x = _x; y = _y;
    }
    Complex operator -(const Complex &b)const{
        return Complex(x - b.x, y - b.y);
    }
    Complex operator +(const Complex &b)const{
        return Complex(x + b.x, y + b.y);
    }
    Complex operator *(const Complex &b)const{
        return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
    }
};

void change(Complex y[],int len){
    int i,j,k;
    for(i = 1, j = len / 2; i < len - 1; i++){
        if(i < j) swap(y[i],y[j]);
        k = len / 2;
        while(j >= k){
            j -= k;
            k /= 2;
        }
        if(j < k) j += k;
    }
}
//len必须是2^k
//on==1进行DFT；on==-1进行IDFT
void fft(Complex y[],int len,int on){
    change(y,len);
    for(int h = 2; h <= len; h <<= 1){
        Complex wn(cos(-on*2*pi/h),sin(-on*2*pi/h));
        for(int j = 0; j < len; j += h){
            Complex w(1,0);
            for(int k = j; k < j + h / 2; k++){
                Complex u = y[k];
                Complex t = w * y[k + h / 2];
                y[k] = u + t;
                y[k + h / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if(on == -1){
        for(int i = 0; i < len; i++) y[i].x /= len;
    }
}

Complex x[maxn<<2];
int len,len1;
len = 1;
while(len < 2 * len1) len <<= 1;
for(int i = 0; i < len1; i++)
   x[i] = a[i];
for(int i = len1; i < len; i++)
   x[i] = Complex(0,0);
fft(x,len,1);
fft(x,len,-1);

二、FFT详解请看大佬的blogGGN_2015

三、个人理解

1、FFT（Fast Fourier Transformation）是离散傅氏变换（DFT）的快速算法。

2、傅里叶变换：能将满足一定条件的某个函数表示成三角函数或他们的积分的线性组合。

3、卷积定理指出：傅里叶变换可以化复杂的卷积运算为简单的乘积运算，从而提供了计算卷积的一种简单手段。

4、卷积是两个变量在某范围内相乘后求和的结果。如果卷积的变量是序列x(n)和h(n)，则卷积的结果

5、可以看出FFT可以用于计算f(n)*g(n)的结果。

6、f(n)*g(n)的结果：

设：f(x) = x+2x^3+x^4，g(x) = x+2x^3+x^4，则

h(x) = f(n)*g(n)

       = (x+2x^3+x^4)*(x+2x^3+x^4)

       =(x*x+x*2x^3+x*x^4)+(2x^3*x+2x^3*2x^3+2x^3*x^4)+(x^4*x+x^4*2x^3+x^4*x^4)

       =x^2+(2+2)x^4+(1+1)x^5+4x^6+(2+2)x^7+x^8（请记住这一步，可以看做是不同位置两两相加）

       = x^2+4x^4+2x^5+4x^6+4x^7+x^8

7、把f(x)换个形式以符合计算机的使用。令x[i]等于幂为i的系数，如f(x) = 0*x^0+1*x^1+0*x^2+2*x^3+1*x^4，则x[] = {0,1,0,2,1}

四、例题

1、HDU4609 给出 n 条线段长度，问任取 3 根，组成三角形的概率。（详细看邝斌大神）

num[]表示长度为i的线段有多少根，现在求线段长度两两相加后长度为i的组合有多少种。如果暴力求线段相加的长度则复杂度为O(n^2)，会超时。这时考虑到卷积可以求不同位置两两相加的结果，而且可以用FFT以O(n)的复杂度求出卷积的结果。

现在我们构造f(x)和g(x)函数，令k*x^i表示长度为i的线段（组合）有k根（种），f(x) = sum(num[i] * x^i)。因为现在求线段两两组合所以g(x) = f(x)。

fft(x,len,1);//DFT使f(x)函数转化为三角函数或他们的积分的线性组合，便于计算

for(int i = 0; i < len; i++) x[i] = x[i] * x[i];//计算f(x)*g(x)的卷积结果，此时还是三角函数或他们的积分的线性组合的形式

fft(x,len,-1);//IDFT将三角函数或他们的积分的线性组合的形式转化为离散点的形式，得到结果

后面就是组合数学的内容了，详细请看邝斌大神的题解。

2、Gym-101667H 求B在A中选一个起点开始匹配，能匹配到的位置的最大数量。

简单想到，先将A的字符转化S->R,R->P,P->S，表示成哪些字符能得分。接下来也是两两匹配，所以我们可以对每种字符分别处理。针对字符S，x1[i]表示A位置i的字符为S则为1，否则为0。此时，如果位置i，j的字符相同则卷积位置i+j的结果为1，否则为0。当卷积结果计算完以后，我们只需查看sum(x[i+j]，x[i+j+2]，...，x[i+j+2*len2])(0<=i<len1,0<=j<len2)，sum的最大值。

但是这样复杂度不能满足。所以，反转B，则结果为max(x[i+j+len2]，x[i+1+j+len2-1]，...，x[i+len2+j+len2-len2])，即max(x[i+j+len2])。

#include<bits/stdc++.h>
#define ms(a, b) memset(a, b, sizeof(a));
#define lowbit(x) (x&(-x))
using namespace std;
typedef long long LL;
typedef unsigned long long uLL;
typedef pair<int, int> pii;
const LL Linf = 3e12 + 10;
const int inf = 1e9 + 7;
const int maxn = 1e5 + 10;
const int maxm = 2e3 + 10;
const int maxb = 26;
const int sigma_size = 62;
const double eps = 1e-8;
const int mod =  20071027;
const double pi = acos(-1.0);
const int Ha = 2333;

struct Complex{
    double x,y;
    Complex(double _x = 0.0, double _y = 0.0){
        x = _x; y = _y;
    }
    Complex operator -(const Complex &b)const{
        return Complex(x - b.x, y - b.y);
    }
    Complex operator +(const Complex &b)const{
        return Complex(x + b.x, y + b.y);
    }
    Complex operator *(const Complex &b)const{
        return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
    }
};

void change(Complex y[],int len){
    int i,j,k;
    for(i = 1, j = len / 2; i < len - 1; i++){
        if(i < j) swap(y[i],y[j]);
        k = len / 2;
        while(j >= k){
            j -= k;
            k /= 2;
        }
        if(j < k) j += k;
    }
}

void fft(Complex y[],int len,int on){
    change(y,len);
    for(int h = 2; h <= len; h <<= 1){
        Complex wn(cos(-on*2*pi/h),sin(-on*2*pi/h));
        for(int j = 0; j < len; j += h){
            Complex w(1,0);
            for(int k = j; k < j + h / 2; k++){
                Complex u = y[k];
                Complex t = w * y[k + h / 2];
                y[k] = u + t;
                y[k + h / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if(on == -1){
        for(int i = 0; i < len; i++) y[i].x /= len;
    }
}

Complex x1[maxn<<2],x2[maxn<<2];
char str1[maxn],str2[maxn];
int sum[maxn<<2];
int len,len1,len2;

void f(char c){
    for(int i = 0; i < len1; i++)
        x1[i] = Complex((str1[i] == c) ? 1 : 0, 0);
    for(int i = len1; i < len; i++)
        x1[i] = Complex(0,0);
    for(int i = 0; i < len2; i++)
        x2[i] = Complex((str2[i] == c) ? 1 : 0, 0);
    for(int i = len2; i < len; i++)
        x2[i] = Complex(0,0);
    fft(x1,len,1);
    fft(x2,len,1);
    for(int i = 0; i < len; i++)
        x1[i] = x1[i] * x2[i];
    fft(x1,len,-1);
    for(int i = 0; i < len; i++)
        sum[i] += (int)(x1[i].x + 0.5);
}

int main()
{
#ifdef LOCAL
    freopen("input.in","r",stdin);
//    freopen("out.in","w",stdout);
#endif
    while(~scanf("%d%d%s%s",&len1,&len2,str1,str2)){
        ms(sum,0);
        for(int i = 0; i < len2 / 2; i++)
            swap(str2[i],str2[len2 - i - 1]);

        for(int i = 0; i < len1; i++){
            if(str1[i] == 'R') str1[i] = 'P';
            else if(str1[i] == 'P') str1[i] = 'S';
            else if(str1[i] == 'S') str1[i] = 'R';
        }
        len = 1;
        while(len < 2 * len1 || len < 2 * len2) len <<= 1;

        f('R');
        f('P');
        f('S');

        int ans = 0;
        for(int i = len2 - 1; i < len; i++)
            ans = max(ans,sum[i]);
        printf("%d\n",ans);
    }
    return 0;
}