1. 程式人生 > >bzoj 3527: [Zjoi2014]力【FFT】

bzoj 3527: [Zjoi2014]力【FFT】

註意 math logs oid cpp 開始 ios 下標 ++

大力推公式,目標是轉成卷積形式:\( C_i=\sum_{j=1}^{i}a_jb_{i-j} \)
首先下標從0開始存,n--
\[ F_i=\frac{\sum_{j<i}\frac{q_jq_i}{(j-i)^2}-\sum_{j>i}\frac{q_jq_i}{(j-i)^2}}{q_i} \]
\[ F_i=\sum_{j<i}\frac{q_j}{(j-i)^2}-\sum_{j>i}\frac{q_j}{(j-i)^2} \]

\[ a_i=\sum_{j<i}\frac{q_j}{(j-i)^2} \]
\[ b_i=\sum_{j>i}\frac{q_j}{(j-i)^2} \]


\[ f_j=F_j \]
\[ g_i=\frac{1}{i^2} \]

\[ F_i=a_i-b_i \]
先推ai
!註意j==i的情況下g函數為0不影響結果,所以方便起見把大於小於都換成了大於等於小於等於
\[ a_i=\sum_{j<i}\frac{q_j}{(j-i)^2} \]
\[ =\sum_{j=0}^{i}\frac{q_j}{(i-j)^2} \]
\[ =\sum_{j=0}^{i}f_i*g_{i-j} \]
於是卷積*1get
再推bi
\[ b_i=\sum_{j>i}\frac{q_j}{(j-i)^2} \]
\[ =\sum_{j=i}^{n}\frac{q_j}{(j-i)^2} \]

然後用j+i替換j
\[ =\sum_{j+i=i}^{n}\frac{q_{j+i}}{(j+i-i)^2} \]
\[ =\sum_{j=0}^{n-i}\frac{q_{j+i}}{j^2} \]
\[ =\sum_{j=0}^{n-i}f_{j+i}\*g_j \]
然後用n-i-j替換j
\[ =\sum_{n-i-j=0}^{n-i}f_{n-i-j+i}*g_{n-i-j} \]
\[ =\sum_{0 \leq n-i-j \leq n-i}f_{n-i-j+i}*g_{n-i-j} \]
然後發現這樣轉化一下
\[ 0 \leq n-i-j\Rightarrow j \leq n-i \]

\[ n-i-j \leq n-i\Rightarrow 0 \leq j \]
\[ =\sum_{j=0}^{n-i}f_{n-j}*g_{n-i-j} \]

\[ f1_i=f_{n-i} \]
\[ =\sum_{j=0}^{n-i}f1_j*g_{n-i-j} \]

\[ t=n-i \]
\[ =\sum_{j=0}^{t}f1_j*g_{t-j} \]
於是卷積*2get
然後直接上FFT即可
!!!對於g函數,應該是g[i]=1.0/i/i !!g[i]=1.0/(i*i)會炸精度!

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int N=1000005;
int n,lm,bt,re[N];
struct cd
{
    double r,i;
    cd(double R=0,double I=0)
    {
        r=R,i=I;
    }
    cd operator + (cd &a) const
    {
        return cd(r+a.r,i+a.i);
    }
    cd operator - (cd &a) const
    {
        return cd(r-a.r,i-a.i);
    }
    cd operator * (cd &a) const
    {
        return cd(r*a.r-i*a.i,r*a.i+i*a.r);
    }
}f[N],f1[N],g[N],a[N],b[N];
void dft(cd a[],int f)
{
    for(int i=1;i<lm;i++)
        if(i<re[i])
            swap(a[i],a[re[i]]);
    for(int i=2;i<=lm;i<<=1)
    {
        cd wi=cd(cos(2.0*M_PI/i),f*sin(2.0*M_PI/i));
        for(int k=0;k<lm;k+=i)
        {
            cd w=cd(1,0),x,y;
            for(int j=0;j<(i>>1);j++)
            {
                x=a[j+k];
                y=a[j+k+(i>>1)]*w;
                a[j+k]=x+y;
                a[j+k+(i>>1)]=x-y;
                w=w*wi;
            }
        }
    }
    if(f==-1)
        for(int i=0;i<lm;i++)
            a[i].r/=lm;
}
int main()
{
    scanf("%d",&n);
    n--;
    for(int i=0;i<=n;i++)
    {
        scanf("%lf",&f[i].r);
        f1[n-i].r=f[i].r;
    }
    for(int i=1;i<=n;i++)
        g[i].r=1.0/i/i;
    for(int i=0;;i++)
        if((1<<i)>2*n)
        {
            bt=i;
            lm=(1<<i);
            break;
        }
    for(int i=0;i<lm;i++)
        re[i]=(re[i>>1]>>1)|((i&1)<<(bt-1));
    dft(f,1);
    dft(f1,1);
    dft(g,1);
    for(int i=0;i<lm;i++)
    {
        a[i]=f[i]*g[i];
        b[i]=f1[i]*g[i];
    }
    dft(a,-1);
    dft(b,-1);
    for(int i=0;i<=n;i++)
        printf("%.3f\n",a[i].r-b[n-i].r);
    return 0;
}

bzoj 3527: [Zjoi2014]力【FFT】