BestCoder#44
NOI2015团灭记

BZOJ1919: [Ctsc2010]性能优化

SkyDec posted @ 2015年6月16日 12:57 in 杂乱无章 , 1542 阅读

练练tex

首先这道题是循环卷积,下标是取模的,我们先倍长,然后再模回来,设长度为n,设a与b的点值为GA,GB

[tex]2*n*F_k=\sum_{i=0}^{2*n-1}GA_i*GB_i*\omega_{2*n}^{-k*i}[/tex]

设实际上我们循环卷积搞出来的数组为T

[tex]2*n*T_k[/tex]

[tex]=2*n*(F_k+F_{n+k})[/tex]

[tex]=\sum_{i=0}^{2*n-1}GA_i*GB_i*(\omega_{2*n}^{-k*i}+\omega_{2*n}^{-(k+n)*i})[/tex]

[tex]=\sum_{i=0}^{2*n-1}GA_i*GB_i*(\omega_{2*n}^{-k*i}+\omega_{2*n}^{-k*i}*\omega_{2*n}^{-n*i})[/tex]

[tex]=\sum_{i=0}^{2*n-1}GA_i*GB_i*(\omega_{2*n}^{-k*i}+\omega_{2*n}^{-k*i}*\omega_{2}^{-i})[/tex]

当i为奇数时,[tex]\omega_{2}^{-i}[/tex]为-1,否则为1

[tex]=\sum_{i=0}^{2*n-1}GA_i*GB_i*(\omega_{2*n}^{-k*i}*2*[i\%2==0])[/tex]

[tex]=\sum_{i=0}^{n-1}GA_{2*i}*GB_{2*i}*(\omega_{2*n}^{-2*k*i}*2)[/tex]

[tex]=\sum_{i=0}^{n-1}GA_{2*i}*GB_{2*i}*(\omega_{n}^{-k*i}*2)[/tex]

所以:

[tex]n*T_k=\sum_{i=0}^{n-1}GA_{2*i}*GB_{2*i}*\omega_{n}^{-k*i}[/tex]

其中:

[tex]GA_{2*i}[/tex]

[tex]=\sum_{j=0}^{2*n-1}a_i*\omega_{2*n}^{2*i*j}[/tex]

[tex]=\sum_{j=0}^{2*n-1}a_i*\omega_{n}^{i*j}[/tex]

可以发现当[tex]n\leq j[/tex]时,[tex]a_j=0[/tex]

[tex]=\sum_{j=0}^{n-1}a_i*\omega_{n}^{i*j}[/tex]

我们可以发现这就是做了一遍长度为n的FFT

于是下标取模的卷积就是直接做长度为n的FFT就好啦o(^▽^)o

有了这个性质,我们算出b的点值形式GB后,直接令[tex]GB_{i}=GB_{i}^{k}[/tex]即可

然而长度是一堆质数的乘积,这个有点麻烦

设d是n的最小质因子

[tex]GA_{k}[/tex]

[tex]=\sum_{i=0}^{n-1}a_i*\omega_{n}^{i*k}[/tex]

[tex]=\sum_{i=0}^{d-1}\sum_{j=0}^{n/d}a_{j*d+i}*\omega_{n}^{(j*d+i)*k}[/tex]

[tex]=\sum_{i=0}^{d-1}(\omega_{n}^{k})^{i}\sum_{j=0}^{n/d}a_{j*d+i}*(\omega_{n}^{d*k})^{j}[/tex]

我们可以发现

[tex]\omega_{n}^{k*d}=\omega_{n}^{(k+\frac{n}{d})*d}[/tex]

于是像基于二进制的FFT一样,系数和点值都缩小了,递归下去做即可,稍微改改就可以变成非递归的了

/**************************************************************
    Problem: 1919
    User: SKYDEC
    Language: C++
    Result: Accepted
    Time:14024 ms
    Memory:18744 kb
****************************************************************/
 
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#define rep(i,j,k) for(int i=(int)j;i<=(int)k;i++)
#define per(i,j,k) for(int i=(int)j;i>=(int)k;i--)
using namespace std;
typedef long long LL;
typedef double db;
inline void read(int &x){
    x=0;char p=getchar();
    while(!(p<='9'&&p>='0'))p=getchar();
    while(p<='9'&&p>='0')x*=10,x+=p-48,p=getchar();
}
int P;
int G;
const int N=510000;
inline int Pow(int a,int b){
    int c=1;
    for(;b;b>>=1,a=a*1ll*a%P)if(b&1)c=c*1ll*a%P;
    return c;
}
int pr[N];
int mul[N];
int suf[N];
int w[2][N];
int m=0;
void init(int n){
    int x=n;
    for(int i=2;i*i<=x;i++){
        while(x%i==0){
            pr[++m]=i;
            x/=i;
        }
    }
    if(x!=1)pr[++m]=x;
    for(G=2;G;G++){
        bool flag=1;
        rep(i,1,m)if(Pow(G,(P-1)/pr[i])==1)flag=0;
        if(flag)break;
    }
    mul[0]=1;rep(i,1,m)mul[i]=mul[i-1]*pr[i];
    suf[m+1]=1;per(i,m,1)suf[i]=suf[i+1]*pr[i];
    w[0][0]=w[1][0]=1;
    int V=Pow(G,(P-1)/n);
    int VV=Pow(V,P-2);
    rep(i,1,n-1){
        w[0][i]=w[0][i-1]*1ll*V%P;
        w[1][i]=w[1][i-1]*1ll*VV%P;
    }
}
int d[N];
int tmp[N];
void Rev(int *a,int l,int r,int x){
    if(l==r)return;
    rep(i,l,r)tmp[i]=a[i];
    rep(v,0,pr[x]-1){
        rep(k,0,(r-l+1)/pr[x]-1)
        a[l+v*((r-l+1)/pr[x])+k]=tmp[l+v+k*pr[x]];
    }
    int len=(r-l+1)/pr[x];
    rep(v,0,pr[x]-1){
        Rev(a,l,l+len-1,x+1);
        l+=len;
    }
}
inline int fft(int *a,int n,int f){
    Rev(a,0,n-1,1);
    for(int i=m;i>=1;i--){
        int ww=Pow(G,(P-1)/pr[i]);
        if(f)ww=Pow(ww,P-2);
        for(int j=0,t=n/(suf[i]);j<n;j+=suf[i])
        for(int k=0,l=0;k<suf[i+1];k++,l+=t){
            int bt=1;
            rep(v,0,pr[i]-1){
                d[v]=a[j+k+(suf[i+1])*v]*1ll*bt%P;
                bt=bt*1ll*w[f][l]%P;
            }
            int base=1;
            rep(v,0,pr[i]-1){
                int ret=0;
                int vv=1;
                rep(l,0,pr[i]-1){
                    ret=(ret+d[l]*1ll*vv)%P;
                    vv=vv*1ll*base%P;
                }
                a[j+k+suf[i+1]*v]=ret;
                base=base*1ll*ww%P;
            }
        }
    }
    if(f){
        int vv=Pow(n,P-2);
        rep(i,0,n-1)a[i]=a[i]*1ll*vv%P;
    }
}
int n,C;
int a[N],b[N];
int main(){
    read(n);read(C);C%=n;
    P=n+1;
    init(n);
    rep(i,0,n-1)read(a[i]);
    fft(a,n,0);
    rep(i,0,n-1)read(b[i]);
    fft(b,n,0);
    rep(i,0,n-1)b[i]=Pow(b[i],C)%P;
    rep(i,0,n-1)a[i]=a[i]*1ll*b[i]%P;
    fft(a,n,1);
    rep(i,0,n-1)printf("%d\n",a[i]);
    return 0;
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter