FFT模板 & 瞎做

参考资料

1. 快速傅里叶变换(FFT)详解
2. FFT 学习笔记
3. FFT模板以及简单的FFT入门题总结

模板

typedef cp complex<double>;
namespace fft
{
    cp o[4*MAXN],oi[4*MAXN];
    void init(int n)
    {
        for(int i=0;i<n;i++)
        {
            o[i] =cp(cos(2*PI/n*i),sin(2*PI/n*i));
            oi[i]=conj(o[i]);
        }
    }
    void transform(cp a[],const int n,cp o[])
    {
        //二进制翻转
        int k=0;
        while((1<<k)<n)k++;
        for(int i=0;i<n;i++)
        {
            int t=0;
            for(int j=0;j<k;j++)
                if(i&(1<<j))
                    t|=(1<<(k-j-1));
            if(i<t) swap(a[i],a[t]);
        }
        for(int l=2;l<=n;l*=2)
        {
            int m=l/2;
            //将两个长度为 m 的序列的答案合并为长度为 l 的序列的答案
            for(cp *p=a;p!=a+n;p+=l)
            {
                for(int i=0;i<m;i++)
                {
                    cp t=o[n/l*i]*p[m+i];
                    p[m+i]=p[i]-t,p[i]+=t;
                }
            }
        }
    }
    void dft (cp a[],int n){transform(a,n,o);}
    void idft(cp a[],int n){transform(a,n,oi);for(int i=0;i<n;i++) a[i]/=n;}
};

//n1 n2为多项式a1 a2系数的个数 不是最高次数(注意+1)

void mul(int a1[],int n1,int a2[],int n2,int res[])
{
    int n=1;
    while(n<n1+n2) n=n<<1;
    static cp c1[4*MAXN],c2[4*MAXN];
    for(int i=0;i< n;i++) c1[i]=c2[i]=0;
    for(int i=0;i<n1;i++) c1[i].real(a1[i]);
    for(int i=0;i<n2;i++) c2[i].real(a2[i]);
    fft::init(n);
    fft::dft(c1,n),fft::dft(c2,n);
    for(int i=0;i<n;i++) c1[i]*=c2[i];
    fft::idft(c1,n);
    for(int i=0;i<n1+n2-1;i++)
        res[i]=int(floor(c1[i].real()+0.5));
}


/*
void mul(double a1[],int n1,double a2[],int n2,double res[])
{
    int n=1;
    while(n<n1+n2) n=n<<1;
    static cp c1[4*MAXN],c2[4*MAXN];
    for(int i=0;i< n;i++) c1[i]=c2[i]=0;
    for(int i=0;i<n1;i++) c1[i].real(a1[i]);
    for(int i=0;i<n2;i++) c2[i].real(a2[i]);
    fft::init(n);
    fft::dft(c1,n),fft::dft(c2,n);
    for(int i=0;i<n;i++) c1[i]*=c2[i];
    fft::idft(c1,n);
    for(int i=0;i<n1+n2-1;i++)
        res[i]=c1[i].real();
}
*/

常数优化

struct cp
{
    double r,i;
    cp(double _r = 0,double _i = 0){r = _r; i = _i;}
    cp operator +(const cp  &b){return cp(r+b.r,i+b.i);}
    cp operator -(const cp  &b){return cp(r-b.r,i-b.i);}
    cp operator *(const cp  &b){return cp(r*b.r-i*b.i,r*b.i+i*b.r);}
    cp operator /(const int &n){return cp(r/n,i/n);}
};

namespace fft
{
    int R[MAXN*4];
    void init(int n)
    {
        int L=0;
        while((1<<L)<n) L++;
        for(int i=0;i<n;i++)
            R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    }
    void transform(cp a[],const int n,int f)
    {
        for(int i=0;i<n;i++)
            if(i<R[i])
                swap(a[i],a[R[i]]);
        for(int l=2;l<=n;l=l<<1)
        {
            int m=l/2;
            cp wn=cp(cos(2*PI/l),f*sin(2*PI/l));
            for(cp *p=a;p!=a+n;p+=l)
            {
                cp w=cp(1,0);
                for(int i=0;i<m;i++,w=w*wn)
                {
                    cp t=w*p[m+i];
                    p[m+i]=p[i]-t,p[i]=p[i]+t;
                }
            }
        }
    }
    void dft (cp a[],int n){transform(a,n, 1);}
    void idft(cp a[],int n){transform(a,n,-1);for(int i=0;i<n;i++) a[i]=a[i]/n;}
};

cp c1[MAXN*4],c2[MAXN*4];
void mul(LL a1[],int n1,LL a2[],int n2,LL res[])
{
    int n=1;
    while(n<n1+n2) n=n<<1;
    for(int i=0;i<=n;i++) c1[i]=c2[i]=cp(0,0);
    for(int i=0;i<n1;i++) c1[i]=cp(a1[i],0);
    for(int i=0;i<n2;i++) c2[i]=cp(a2[i],0);
    fft::init(n);
    fft::dft(c1,n),fft::dft(c2,n);
    for(int i=0;i<n;i++) c1[i]=c1[i]*c2[i];
    fft::idft(c1,n);
    for(int i=0;i<n1+n2-1;i++)
        res[i]=LL(floor(c1[i].r+0.5));
}

FFT瞎做 & 简要题解
1.luogu P3803 【模板】多项式乘法（FFT)
贴板子AC一气呵成，，，

2.CF Gym101002E E - K-Inversions
下标差的绝对值为K的一对字符A和B称为K逆序对，依次输出字符串中1~n-1逆序对的个数。

万恶之源，Wannafly Camp Day3开的比赛的E题，讲题时dls说这是FFT板子题一眼秒了，，，

构造两个函数
$$f(x)=\{\begin{array}{ll}1,& \text{str[x]==’A’}\\ 0,& \text{else}\end{array}g(x)=\{\begin{array}{ll}1,& \text{str[x]==’B’}\\ 0,& \text{else}\end{array}$$
问题等价于求函数
$$h(j)=\sum _{i=j+1}^{n}f(i)\ast g(i-j)$$
在1~n-1上的取值
设函数$g^{\prime }(x)$是由$g(x)$翻转得到的(字符串下标由1开始)，即：
$$g^{\prime }(n-x+1)=g(x)$$
原式变为
$$h(j)=\sum _{i=j+1}^{n}f(i)\ast g^{\prime }(n-i+j+1)$$
注意到$g(n-i+j+1)$和$f(i)$的自变量之和等于$n+j+1$，有点类似于我们熟悉的卷积，考虑扩展上下界至$[0,n+j+1]$：

因为$$f(i)=0,i\in [n+1,n+j+1]$$
$$g(n-i+j+1)=0,i\in [0,j]$$
所以
$$h(j)=\sum _{i=0}^{n+j+1}f(i)\ast g^{\prime }(n-i+j+1)$$
将$f(x)$和$g^{\prime }(x)$写成多项式后相乘，$x^{n+j+1}$项的系数就是$h(j)$的值，FFT板子贴上搞一搞就行。

int n;
int a[MAXN],b[MAXN],res[MAXN];
char s[MAXN];

int work()
{
    scanf("%s",s+1);
    n=strlen(s+1);
    for(int i=1;i<=n;i++) a[i]=(s[i]=='A'),b[n-i+1]=(s[i]=='B');
    mul(a,n+1,b,n+1,res);
    for(int j=1;j<n;j++)
        printf("%d\n",res[n+j+1]);
    return 0;
}

3.BZOJ2194 快速傅立叶之二

定义 $C_k={\sum }_{i=k}^{n-1}{A}_i{B}_{i-k}$ ，分别求 $C_0$ ~ $C_{n-1}$ 的值。

设数组$B^{\prime }$由数组$B$翻转所得，即：
$$B_{n-i-1}^{\prime }=B_i$$
所以
$$C_k=\sum _{i=k}^{n-1}{A}_i{B}_{n-i+k-1}^{\prime }$$

因为
$$B_{n-i+k+1}^{\prime }=0,i\in [0,k-1]$$
$$A_i=0,i\in [n,n+k-1]$$

所以可将上下界拓展至$[0,n+k-1]$，即：
$$C_k=D_{n+k-1}=\sum _{i=0}^{n+k-1}{A}_i{B}_{n-i+k-1}^{\prime }$$
将 $A_i$ 和 $B_i^{\prime }$ 作为两个多项式的系数，相乘后得到的$x^{n+k-1}$项的系数即为所求的 $C_k$

int n;
int a[MAXN],b[MAXN],res[MAXN];
 
int work()
{
    n=read();
    for(int i=0;i<n;i++)
        a[i]=read(),b[n-i-1]=read();
    mul(a,n,b,n,res);
    for(int i=0;i<n;i++)
        printf("%d\n",res[n+i-1]);
    return 0;
}

小结：形如 $C_j={\sum }_{i=0}^j{A}_i{B}_{j-i}$ 次数之和为定值的形式为标准的卷积，但次数之差为定值的形式可以通过翻转其中一个数组并扩展上下界来化为标准的卷积形式。

4.BZOJ3527 [Zjoi2014]力
定义
$$E_j=\sum _{i=1}^{j-1}\frac{q_i}{(i-j{)}^2}-\sum _{i=j+1}^n\frac{q_i}{(i-j{)}^2}$$
求$E{}_1$ ~ $E_n$的值。

将两个求和部分分开考虑，定义
$$f(x)=\{\begin{array}{ll}{q}_x& 1<=x<=n\\ 0& x=0\end{array}g(x)=\{\begin{array}{ll}\frac{1}{x^2}& x̸=0\\ 0& x=0\end{array}$$
$$\sum _{i=1}^{j-1}\frac{q_i}{(i-j{)}^2}$$就可以写作
$$\sum _{i=1}^{j-1}f(i)g(i-j)$$

因为 $f(0)=0$ 并且 $g(j-j)=g(0)=0$，所以可将上下界扩展至$[0,j]$，标准的卷积形式
$$\sum _{i=0}^{j}f(i)g(i-j)$$
第二部分可以通过第二/第三题的翻转序列并扩展上下界技巧转换为标准卷积形式，两部分分别FFT并相减即可。

int n;
double q[MAXN],p[MAXN];
double res1[MAXN],res2[MAXN];
 
int work()
{
    cin >> n;
    for(int i=1;i<=n;i++) cin >> q[i];
    for(int i=1;i<=n;i++) p[i]=(1.0/(double(i)*double(i)));
    mul(q,n+1,p,n+1,res1);
    for(int i=1;i<=(n/2);i++) swap(p[i],p[n-i+1]);
    mul(q,n+1,p,n+1,res2);
    for(int i=1;i<=n;i++)
        printf("%f\n",res1[i]-res2[n+i+1]);
    return 0;
}

5.HDU4609 3-idiots
给你一个包含n个数的数列 $a_1$ ~ $a_n$，求任选三个数构成三角形的概率(方案)。

构造一个多项式，$A_i$ 为 $x^i$ 的系数，并且$A_i$等于数列中 $i$ 出现的次数，将这个多项式和自己乘一下，所得新的多项式 $x^i$ 的系数 $B_i$ 即为数列中选两个数和为 $i$ 的方案数。但题目要求选的两个数不能相同，因此对于每一个 $a_i$ 都要在 $B_{2\ast a_i}$ 中减去选两个$a_i$ 的方案。此外 $A_i$ 所包含的方案中 $a_i+a_j$ 与 $a_j+a_i$ 被视为两种不同的方案(因为在多项式相乘中$x^{a_i+a_j}$的系数$B_{a_i+a_j}=A_{a_i}{A}_{a_j}+A_{a_j}{A}_{a_i})$，因此除以二后便得到了任选两个不相同的数构成 $i$ 的方案数 $B_i$ 。
考虑枚举三角形的最大边来计算所有不能构成三角形的方案数量，因为没有负数和0，不合法方案中的较小的两条边之和一定小于等于第三边，因此枚举最大边的长度 $i$ 并统计两边之和$<=i$的方案数即可(对系数$B_i$做前缀和)。
数组开小交了20多发

LL n,m,p,q;
int x[MAXN];
LL num[N],res[N];

void init()
{
    p=q=n*(n-1)*(n-2)/6,m=0;
    mst(num,0);
}

int work()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d",&n);
        init();
        for(int i=1;i<=n;i++)
            scanf("%d",&x[i]),num[x[i]]++,m=max(m,(LL)x[i]);
        mul(num,m+1,num,m+1,res);
        for(int i=1;i<=n  ;i++) res[x[i]*2]--;
        for(int i=1;i<=2*m;i++) res[i]/=2,res[i]+=res[i-1];
        for(int i=1;i<=2*m;i++) p-=res[i]*num[i];
        printf("%.7lf\n",(double)p/q);
    }
    return 0;
}

正向考虑似乎有点麻烦，可以参考kuangbin大大的博客，，，

6.UVA - 12298 Super Poker II
现在有一副只有合数(没有1)的扑克牌，写着某个数的扑克牌有四种(Spade Heart Club Diamond)且各有一张，现在有一些扑克牌丢了，求由四种扑克牌构成 $[a,b]$ 上每一个数的方案数。
和上一题的思想类似，构造四个多项式分别代表四种扑克牌，若某一种数字为 $x$ 的扑克牌存在就在对应的多项式系数$A_x$上让它等于 $1$，三次多项式乘法后输出次数为 $[a,b]$ 项的系数就行。
注意开long double 有点卡精度
上午交的下午才看到结果 神仙UVA评测姬

map<char,int> h;
int l,r,m;
LL  num[4][MAXN];
void pre()
{
    for(int i=1;i<=r;i++)
    {
        bool flag=0;
        for(int j=2;!flag&&j<i;j++)
            if(i%j==0)
                flag=1;
        for(int k=0;k<4;k++)
            num[k][i]=flag;
    }
}

int work()
{
    h['S']=0,h['H']=1,h['C']=2,h['D']=3;
    while(cin >> l >> r >> m && l && r)
    {
        pre();
        while(m--)
        {
            string card;
            cin >> card;
            int t=0;
            for(auto ch:card)
                if(ch>='0'&&ch<='9')
                    t=t*10+ch-'0';
            num[h[card.back()]][t]--;
        }
        for(int k=1;k<4;k++)
            mul(num[0],r+1,num[k],r+1,num[0]);
        for(int i=l;i<=r;i++)
            cout << num[0][i] << endl;
        cout << endl;
    }
    return 0;
}

然后翻到了学长的博客，又补了几题，，，

7.2018宁夏网络赛H题 Rock Paper Scissors Lizard Spock

题意抄下学长的ORZ，，，

给你类似剪刀石头布游戏的五种手势和十种克制关系。每种手势克制其他两种手势并被另外两种手势克制。
给你两个字符串分别表示A，B的手势顺序，A长B短，你可以随意挪动B相对于A的位置，求B最多能赢多少次？

首先可以考虑字符串中某一种手势赢的次数，以P为例，P可以赢K和R，所以从位置 $j$ 开始游戏能赢的次数等价于求函数:
$$h(j)=\sum _{i=j}^{j+m-1}f(i)g(i-j+1)$$
在 $[1,n]$ 的最大值(字符串下标从1开始)。

其中
$$f(x)=\{\begin{array}{ll}1,& \text{str1[x]=’K’ or str1[x]=’R’}\\ 0,& \text{else}\end{array}g(x)=\{\begin{array}{ll}1,& \text{str2[x]=’P’}\\ 0,& \text{else}\end{array}$$

老套路，翻转 $g(x)$ 扩展上下界：
$$g^{\prime }(m-x+1)=g(x)$$
$$h(j)=\sum _{i=j}^{j+m-1}f(i)g^{\prime }(m-i+j)=\sum _{i=0}^{j+m}f(i)g^{\prime }(m-i+j)$$

因为 $f(0)=0,g^{\prime }(m-j-m+j)=g^{\prime }(0)=g(m+1)=0$

其他四种手势也类似，最后把 $x^{m+j},j\in [1,n]$ 的系数累加取个max就是答案了。

int n,m;
char s1[MAXN],s2[MAXN];
int a[MAXN*2],b[MAXN*2],ans[MAXN];
map<char,vector<char>> h;
vector<char> type={'P','R','S','L','K'};

void pre()
{
    vector<char> t;
    t.clear(),t.pb('R'),t.pb('K'),h['P']=t; //P-R K
    t.clear(),t.pb('P'),t.pb('L'),h['S']=t; //S-P L
    t.clear(),t.pb('L'),t.pb('S'),h['R']=t; //R-L S
    t.clear(),t.pb('K'),t.pb('P'),h['L']=t; //L-K P
    t.clear(),t.pb('S'),t.pb('R'),h['K']=t; //K-S R
}

int work()
{
    pre();
    scanf("%s",s1+1);
    scanf("%s",s2+1);
    n=strlen(s1+1),m=strlen(s2+1);
    for(auto cur:type)
    {
        for(int i=1;i<=m;i++) b[m-i+1]=(s2[i]==cur);
        for(int i=1;i<=n;i++) a[i]=(find(all(h[cur]),s1[i])!=h[cur].end());
        mul(a,n+1,b,m+1,a);
        for(int j=1;j<=n;j++) ans[j]+=a[j+m];
    }
    int t=0;
    for(int j=1;j<=n;j++)
        t=max(t,ans[j]);
    return printf("%d\n",t);
}

8.BZOJ4259 残缺的字符串
题意中文的。。。大致就是带通配符的字符串匹配
定义
$$f(x)=\{\begin{array}{ll}0,& \text{str1[x]=’*’}\\ str1[x],& \text{else}\end{array}g(x)=\{\begin{array}{ll}0,& \text{str2[x]=’*’}\\ str2[x],& \text{else}\end{array}$$
字符串$str2$在字符串$str1$的 $j$ 位置开始能匹配上等价于
$$h(j)=\sum _{i=j}^{j+m-1}(f(i)-g(i-j+1){)}^2\ast f(i)\ast g(i-j+1)=0$$
展开
$$h(j)=\sum _{i=j}^{j+m-1}f(i{)}^3\ast g(i-j+1)-2\ast \sum _{i=j}^{j+m-1}f(i{)}^2\ast g(i-j+1{)}^2+\sum _{i=j}^{j+m-1}f(i)\ast g(i-j+1{)}^3$$
三项分别把 $g(x)$ 翻转，求卷积，合并，看下哪位是0就行。

这题卡常，，，抄板子请用下头那个

int n,m;
char s1[MAXN],s2[MAXN];
LL a[MAXN*4],b[MAXN*4],ans[MAXN];
LL f[MAXN],g[MAXN];
int work()
{
    scanf("%d%d",&m,&n);
    scanf("%s",s2+1);
    scanf("%s",s1+1);
    for(int i=1;i<=n;i++) f[i]=((s1[i]=='*')?0:(s1[i]-'a'+1));
    for(int i=1;i<=m;i++) g[i]=((s2[i]=='*')?0:(s2[i]-'a'+1));
 
    for(int i=1;i<=n;i++) a[i]=(LL)f[i]*f[i]*f[i];
    for(int i=1;i<=m;i++) b[m-i+1]=(LL)g[i];
    mul(a,n+1,b,m+1,a);
    for(int i=1;i<=n-m+1;i++) ans[i]+=a[i+m];
 
    for(int i=1;i<=n;i++) a[i]=f[i]*f[i];
    for(int i=1;i<=m;i++) b[m-i+1]=g[i]*g[i];
    mul(a,n+1,b,m+1,a);
    for(int i=1;i<=n-m+1;i++) ans[i]-=2*a[i+m];
 
    for(int i=1;i<=n;i++) a[i]=f[i];
    for(int i=1;i<=m;i++) b[m-i+1]=g[i]*g[i]*g[i];
    mul(a,n+1,b,m+1,a);
    for(int i=1;i<=n-m+1;i++) ans[i]+=a[i+m];
 
    int cnt=0;
    for(int i=1;i<=n-m+1;i++) cnt+=(ans[i]==0);
    printf("%d\n",cnt);
    for(int i=1;i<=n-m+1;i++)
        if(ans[i]==0)
            write(i),putchar(' ');
    if(cnt) putchar('\n');
    return 0;
}

9.A+B Problem Kattis - aplusb
给定n个数，求$i$ $j$ $k$两两不同的三元组$(i,j,k)$，使得$a_i+a_j=a_k$成立的三元组的个数。
由于又负数和0，首先把所有数加上一个$N$，并且把 $0$ 拿出单独考虑。
然后就和和5题的思想类似，构造多项式使得 $x^i$ 的系数 $A_i$ 为数列中 $i$ 出现的次数，做乘法后减去选两个 $a_i$ 的方案，由于没有0 所以 $a_i$ $a_j$ $a_k$ 一定两两不相同，枚举 $a_k$ 统计答案。
包含0的话分三种情况讨论，选三个0，选两个相同数和一个0，选两个相反数和一个0，加上这些情况就行。

int n,m;
LL ans=0,num[MAXN],res[MAXN*2];
LL zero=0;
vector<int> s;

LL A(LL n,LL m)
{
    LL ans=1;
    for(LL i=0;i<m;i++)
        ans*=(n-i);
    return ans;
}

int work()
{
    n=read();
    for(int i=1;i<=n;i++)
    {
        int t=read()+N;
        if(t==N) zero++;
        else num[t]++,m=max(m,t),s.pb(t);
    }
    mul(num,m+1,num,m+1,res);
    for(auto x:s) res[x*2]--;
    //for(int i=2*N;i<=2*m;i++) ans+=num[i-N]*res[i];
    for(auto x:s) ans+=res[x+N];

    for(int x=0;x<=2*N;x++)
        if(num[x]>=2)
            ans+=num[x]*(num[x]-1)*zero*2;
    for(int x=N;x<=2*N;x++)
        if(num[x]&&num[2*N-x])
            ans+=num[x]*num[2*N-x]*2*zero;
    if(zero>=3) ans+=A(zero,3);

    write(ans),putchar('\n');
    return 0;
}

10 P1919 【模板】A*B Problem升级版 凑够十题

10 BZOJ3160 万径人踪灭

我已经想到了一个绝妙的解释，但写不下了，贴个题解跑路，，，

回炉重造.jpg

int manacher(int len,char s[],int p[])
{
    static char a[MAXN];
    int n=1,id=0,mx=0;
    mst(a,0),a[0]='$',a[1]='#';
    for(int i=1;i<=len;i++)
        a[++n]=s[i],a[++n]='#';
    a[n+1]='\0';
    for(int i=1;i<=n;i++)
    {
        if(mx>i) p[i]=min(p[2*id-i],mx-i-1);
        else p[i]=1;
        while(a[i+p[i]]==a[i-p[i]]) p[i]++;
        if(i+p[i]-1>mx) mx=i+p[i]-1,id=i;
    }
    return n;
}

int n;
LL a[MAXN],b[MAXN],res[MAXN];
int p[MAXN],f[MAXN];
LL  pw[MAXN],ans;
char s[MAXN];

const LL MOD = 1e9+7;

void solve(char c)
{
    for(int i=1;i<=n;i++) a[i]=(s[i]==c);
    for(int i=1;i<=n;i++) b[i]=(s[i]==c);
    mul(a,n+1,b,n+1,res);
    for(int i=0;i<=2*n;i++) f[i]+=res[i];
}

void init()
{
    pw[0]=1;
    for(int i=1;i<=2*n;i++)
        pw[i]=pw[i-1]*2%MOD;
}

int work()
{
    scanf("%s",s+1);
    n=strlen(s+1);
    solve('a'),solve('b');
    int len=manacher(n,s,p);
    init();
    for(int i=1;i<=len;i++)
        ans+=pw[(f[i]+1)/2]-1-p[i]/2,
        ans=((ans%MOD)+MOD)%MOD;
    printf("%lld\n",ans);
    return 0;
}