[CTSC2010]性能优化

[CTSC2010]性能优化 

循环卷积快速幂


两个注意点:
n+1不是2^k*P+1形式,任意模数又太慢?n=2^k1*3^k2*5^k3*7^k4

多路分治!深刻理解FFT运算本质:分治,推式子得到从下往上的迭代公式

最后求的是w_n^i的点值

 

快速幂:

循环卷积快速幂比较特殊,就是G*F,>=n的项的系数加到-n位置上

所以,由于w(n,p+n)=w(n,p),点值相乘直接得到G*F的点值表达

F,B点值,快速幂相乘即可

不用把n扩充系数等。

 

// luogu-judger-enable-o2
// luogu-judger-enable-o2
#include<bits/stdc++.h>
#define reg register int
#define il inline
#define fi first
#define se second
#define mk(a,b) make_pair(a,b)
#define numb (ch^'0')
using namespace std;
typedef long long ll;
template<class T>il void rd(T &x){
    char ch;x=0;bool fl=false;
    while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);
    for(x=numb;isdigit(ch=getchar());x=x*10+numb);
    (fl==true)&&(x=-x);
}
template<class T>il void output(T x){if(x/10)output(x/10);putchar(x%10+'0');}
template<class T>il void ot(T x){if(x<0) putchar('-'),x=-x;output(x);putchar(' ');}
template<class T>il void prt(T a[],int st,int nd){for(reg i=st;i<=nd;++i) ot(a[i]);putchar('\n');}

namespace Miracle{
const int N=5e5+5;
int pri[4]={2,3,5,7};
int n,mod,C;
int G,GI;
int ci[10];
int ad(int x,int y){
    return x+y>=mod?x+y-mod:x+y;
}
int qm(int x,int y){
    int ret=1;
    while(y){
        if(y&1) ret=(ll)ret*x%mod;
        x=(ll)x*x%mod;
        y>>=1;
    }
    return ret;
}
int dep[33],cnt;
void divi(){
    int tmp=n;
    while(tmp%2==0) tmp/=2,dep[++cnt]=2,++ci[2];
    while(tmp%3==0) tmp/=3,dep[++cnt]=3,++ci[3];
    while(tmp%5==0) tmp/=5,dep[++cnt]=5,++ci[5];
    while(tmp%7==0) tmp/=7,dep[++cnt]=7,++ci[7];
}
bool che(int x){
    for(reg i=0;i<4;++i){
        if(ci[pri[i]]){
            if(qm(x,n/pri[i])==1) return false;
        }
    }
    return true;
}
void fin(){
    G=2;
    while(!che(G)) ++G;
}

int f[N];
int b[N];
int pos[N];
int getpos(int x){
//    cout<<" getpos "<<x<<endl;
    int len=n,l=0;
    for(reg i=1;i<=cnt;++i){
//        cout<<" i "<<i<<" x "<<x<<" ll "<<l<<endl;
        int be=(x-l)%dep[i],th=(x-l)/dep[i];
        len=len/dep[i];
        x=l+len*be+th;
        l=x-th;
    }
//    cout<<" bac "<<x<<endl;
    return x;
}
int g[N];
int pw[7*N][2];
int mem[10];
void FFT(int *f,int n,int c){
    for(reg i=0;i<n;++i) g[i]=f[i];
    for(reg i=0;i<n;++i) f[pos[i]]=g[i];
    int las=1;
    for(reg i=cnt;i>=1;--i){
        int p=dep[i];
        int st=(mod-1)/(las*p);
        for(reg l=0;l<n;l+=las*p){
            for(reg b=0;b<las;++b){
                for(reg j=0;j<p;++j){
                    mem[j]=0;
                    for(reg i=0;i<p;++i){
                        mem[j]=ad(mem[j],(ll)pw[st*(b+j*las)*i][c]*f[l+b+i*las]%mod);
                    }
                }
                for(reg j=0;j<p;++j){
                    f[l+b+j*las]=mem[j];
                }
            }
        }
        las*=p;
    }
}
int main(){
    rd(n);rd(C);mod=n+1;
    divi();fin();GI=qm(G,mod-2);
//    cout<<" G "<<G<<" GI "<<GI<<endl;//
    pw[0][0]=pw[0][1]=1;
    for(reg i=1;i<=7*n;++i) pw[i][0]=(ll)pw[i-1][0]*GI%mod,pw[i][1]=(ll)pw[i-1][1]*G%mod;
    for(reg i=0;i<n;++i) rd(f[i]),f[i]%=mod;
    for(reg i=0;i<n;++i) rd(b[i]),b[i]%=mod;
    for(reg i=0;i<n;++i) pos[i]=getpos(i);
//    cout<<" pos "<<endl;
//    prt(pos,0,n-1);
    FFT(f,n,1);
//    cout<<" ff "<<endl;
//    prt(f,0,n-1);
    
    FFT(b,n,1);
//    cout<<" bb "<<endl;
//    prt(b,0,n-1);
    
    for(reg i=0;i<n;++i) f[i]=(ll)f[i]*qm(b[i],C)%mod;
    FFT(f,n,0);
    for(reg i=0;i<n;++i){
        f[i]=(ll)f[i]*qm(n,mod-2)%mod;
        printf("%d\n",f[i]);
    }
    return 0;
}

}
signed main(){
    Miracle::main();
    return 0;
}

/*
   Author: *Miracle*
   Date: 2019/3/19 18:42:47
*/

算是对FFT本质的理解吧。

 

posted @ 2019-03-19 20:04  *Miracle*  阅读(441)  评论(0编辑  收藏  举报