标签：re int omega 单位根 MTT 反演 Complex choose Sigma

[HNOI2019]白兔之舞
\[ 注:W是一个矩阵,表示题中w[i][j],下列式子中的W做加减法时表示W_{xy}\\ ans_t=\Sigma_{i\ mod\ k=t}^L{L \choose i}W^i\\ =\Sigma_{i=0}^L[k|(i-t)]{L \choose i}W^i\\ =\Sigma_{i=0}^L\frac1k\Sigma_{j=0}^{k-1}(\omega_k^{i-t})^j{L \choose i}W^i\\ =\Sigma_i^L\frac1k\Sigma_{j=0}^{k-1}\omega_k^{ij-tj}{L \choose i}W^i\\ =\frac1k\Sigma_{j=0}^{k-1}\omega_k^{{t \choose 2}+{j \choose 2}-{t+j \choose 2}}\Sigma_{i=0}^L\omega_k^{ij}{L \choose i}W^i\\ 设C(j)=\Sigma_{i=0}^L\omega_k^{ij}{L \choose i}W^i\\ 把它变成二项式的形式:\\ A(j)=\Sigma_{i=0}^L{L \choose i}(W\omega_k^j)^i*1^{L-i}\\ \therefore A(j)=(W\omega_k^j+1)^L\\ \therefore ans_t=\frac{\omega_k^{t \choose 2}}k\Sigma_{j=0}^{k-1}\omega_k^{j \choose 2}A(j)\omega_k^{-{t+j \choose 2}}\\ 这好像是还不是卷积,但是我们可以把它变成卷积\\ 令f(x)=\Sigma_{i=0}^n\omega_k^{i \choose 2}\\ 令g(x)=\Sigma_{i=0}^nA(i)\omega_k^{-{t+i \choose 2}}\\ 那么按照常规操作,将f倍增,再将g翻转.\\ \therefore ans_t=\frac{\omega_k^{t \choose 2}}k(f*g)(k+t)\\ 然后卷积用MTT算,就没了; \]

#pragma GCC optimize(3)
#include<bits/stdc++.h>
#define N (65536*4)
using namespace std;
const double pi=acos(-1);
int n,k,l,x,y,p;
struct Matrix{
    int a[3][3];
    void format(){
        memset(a,0,sizeof(a));
    }
    friend Matrix operator + (const Matrix &a,const Matrix &b){
        Matrix re;
        for(int i=0;i<3;i++)
        for(int j=0;j<3;j++){
            re.a[i][j]=(a.a[i][j]+b.a[i][j])%p;
        }
        return re;
    }
    friend Matrix operator * (const Matrix &a,const Matrix &b){
        Matrix re;
        re.format();
        for(int i=0;i<3;i++)
        for(int j=0;j<3;j++)
        for(int k=0;k<3;k++){
            re.a[i][j]=(re.a[i][j]+1ll*a.a[i][k]*b.a[k][j])%p;
        }
        return re;
    }
    friend Matrix operator * (const int &a,const Matrix &b){
        Matrix re;
        for(int i=0;i<3;i++)
        for(int j=0;j<3;j++){
            re.a[i][j]=(1ll*a*b.a[i][j])%p;
        }
        return re;
    }
    friend Matrix operator * (const Matrix &a,const int &b){return b*a;}
}I,W;
Matrix Pow(Matrix x,int y){
    Matrix re=I;
    while(y){
        if(y&1)re=re*x;
        y>>=1;
        x=x*x;
    }
    return re;
}
struct Complex{double x,y;};
Complex operator + (Complex a,Complex b){return (Complex){a.x+b.x,a.y+b.y};}
Complex operator - (Complex a,Complex b){return (Complex){a.x-b.x,a.y-b.y};}
Complex operator * (Complex a,Complex b){return (Complex){a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y};}
Complex a[N],b[N],c[N];

int Pow(int x,int y){
    int re=1;
    while(y){
        if(y&1)re=(1ll*re*x)%p;
        y>>=1;
        x=(1ll*x*x)%p;
    }
    return re;
}

int rev[N];
int pre(int n,int m){
    int high=0,cnt=1;;
    while(cnt<n+m)cnt<<=1,high++;
    for(int i=0;i<cnt;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(high-1));
    return cnt;
}
int root(int num){
    int div[128],sum=0;
    int phi=num-1;
    for(int i=2;i<=phi;i++){
        if(phi%i==0){
            div[++sum]=i;
            while(phi%i==0)phi/=i;
        }
    }
    if(phi!=1)div[++sum]=phi;
    int re=2;
    while(1){
        int flag=1;
        for(int i=1;i<=sum;i++){
            int pw=Pow(re,(num-1)/div[i]);
            if(pw==1){flag=0;break;}
        }
        if(flag)return re;
        re++;
    }
}
Complex w[N];
void fft(int lim,Complex *buf,int dft=1){
    for(int i=0;i<lim;i++)if(i<rev[i])swap(buf[i],buf[rev[i]]);
    for(int len=2;len<=lim;len<<=1){
        for(int s=0;s<lim;s+=len){
            for(int k=s;k<s+len/2;k++){
                Complex x,y;
                x=buf[k],y=w[lim/len*(k-s)]*buf[k+len/2];
                buf[k]=x+y;
                buf[k+len/2]=x-y;
            }
        }
    }
    if(dft==-1)for(int i=0;i<lim;i++)buf[i].x/=lim;
}
Complex A1[N],B1[N],A2[N],B2[N];
void Mul(int lim,int *A,int *B){
    for(int i=0;i<lim;i++)A1[i].x=A[i]>>15,B1[i].x=A[i]&32767;
    for(int i=0;i<lim;i++)A2[i].x=B[i]>>15,B2[i].x=B[i]&32767;
    for(int i=0;i<N;++i)w[i]=(Complex){cos(i*pi*2/lim),sin(i*pi*2/lim)};
    fft(lim,A1);
    fft(lim,B1);
    fft(lim,A2);
    fft(lim,B2);
    for(int i=0;i<lim;i++){
        Complex A1_=A1[i],A2_=A2[i];
        Complex B1_=B1[i],B2_=B2[i];
        A1[i]=A1_*A2_;
        A2[i]=B1_*B2_;
        B1[i]=A1_*B2_;
        B2[i]=A2_*B1_;
    }
    for(int i=0;i<lim;++i)w[i]=(Complex){cos(i*pi*2/lim),-sin(i*pi*2/lim)};
    fft(lim,A1,-1);
    fft(lim,B1,-1);
    fft(lim,A2,-1);
    fft(lim,B2,-1);
    int m=32768;
    for(int i=0;i<lim;i++){
        double A1_=A1[i].x,A2_=A2[i].x;
        double B1_=B1[i].x,B2_=B2[i].x;
        A[i]=((long long)(A1_+0.5)*m%p*m%p+(long long)(B1_+B2_+0.5)*m%p+(long long)(A2_+0.5))%p;
    }
}
int gen[N],f[N],g[N];
int omg(long long x){
    x=((x%k)+k)%k;
    return gen[x];
}
signed main(){
#ifdef KOG_juruo
    freopen("data.in","r",stdin);
    freopen("my.out","w",stdout);
#endif
    scanf("%d%d%d%d%d%d",&n,&k,&l,&x,&y,&p);
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            scanf("%d",&W.a[i][j]);
        }
    }
    for(int i=0;i<3;i++)I.a[i][i]=1;
    gen[0]=1,gen[1]=root(p);
    gen[1]=Pow(gen[1],(p-1)/k);
    for(int i=2;i<k;i++)gen[i]=(1ll*gen[i-1]*gen[1])%p;
    for(int i=0;i<=2*k;i++)f[i]=omg(1ll*-i*(i-1)/2);
    for(int i=0;i<k;i++)g[i]=1ll*omg(1ll*i*(i-1)/2)*(Pow(W*omg(i)+I,l).a[x-1][y-1])%p;
    reverse(g,g+k+1);
    int lim=pre(2*k,k);
    Mul(lim,f,g);
    int inv=Pow(k,p-2);
    for(int t=0;t<k;t++){
        printf("%lld\n",1ll*omg(1ll*t*(t-1)/2)*inv%p*f[t+k]%p);
    }
    return 0;
}
/*
3 4 100 3 3 998244353
1 1 1
1 1 1
1 1 1

*/

标签：re,int,omega,单位根,MTT,反演,Complex,choose,Sigma
来源： https://www.cnblogs.com/nlKOG/p/10754952.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。