[多维FFT Bluestein′s Algorithm] Codechef October Challenge 2017 .Chef and Horcrux

题目里那个 img

其实不重要,只要能算出 pi 就行了

pi 的话发现就是个多维FFT的转移形式

Hillan大佬博客里拷了个代码改一改就好了…

Bluestein算法里的FFT可以用暴力卷积代替,因为数据范围小FFT常数大

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;

const int P=330301441;
int rev[3200005];
int w[2][3200005];

typedef long long ll;

inline int Pow(int x,ll y){
    int ret=1;
    for(;y;y>>=1,x=1LL*x*x%P) if(y&1) ret=1LL*ret*x%P;
    return ret;
}

int c[3200005];
int C,K,G,InvG;
int tmp[3200005],tmp1[3200005],B1[3200005],B2[3200005],nxt[1110];

int N,inv;

inline void DFT(int* a,int d,int n,int fl)
{
    int Nf=2*n,*A=tmp,*B=tmp1,*BB=B1;
    if(fl==-1) swap(A,B),BB=B2;
    for(int i=0,j=0;i<n;i++,j=nxt[i]) 
            c[n-i]=1LL*a[i*d]*B[j]%P;
    for(int i=n;i<2*n;i++){
        int *pos=a+(i-n)*d;
        *pos=0;
        for(int j=0;j<=i;j++)
            *pos=(*pos+1LL*BB[i-j]*c[j])%P;
    }
    for(int i=0,j=0;i<n;i++,j=nxt[i]) 
            a[i*d]=1LL*a[i*d]*B[j]%P;
    if(fl==-1)
      for(int i=0;i<n;i++)a[i*d]=1LL*inv*a[i*d]%P;
    for(int i=0;i<(1<<N);i++) c[i]=0;
}

void FWT(int* a,int n,int m,int fl)
{
    int len=1;
    for(int i=0;i<m;len*=K,i++)
        for(int j=0;j<n;j+=len*K)
            for(int k=0;k<len;k++)
                DFT(a+j+k,len,K,fl);
}
int Da[3200005];
int no[3200005];

int num[2200010];
int ans[3200005];

inline char nc(){
    static char buf[100000],*p1=buf,*p2=buf;
    return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++; 
}

inline void rea(int &x){
    char c=nc(); x=0;
    for(;c>'9'||c<'0';c=nc());for(;c>='0'&&c<='9';x=x*10+c-'0',c=nc());
}

inline void rea(long long &x){
    char c=nc(); x=0;
    for(;c>'9'||c<'0';c=nc());for(;c>='0'&&c<='9';x=x*10+c-'0',c=nc());
}

int main()
{
    int T; rea(T);
    while(T--)
    {
        int n,m=0; long long x;
        rea(n); rea(K); rea(x);
        for(int i=1,x;i<=n;i++)
            rea(x),num[x]++;
        ll pro=1,INV=Pow(Pow(2,n),P-2);
        for(n=1;n<=100000;n*=K) m++;
        for(int i=0;i<=n;i++) pro=1LL*pro*(num[i]=Pow(2,num[i]))%P;
        for(int i=0,*j=Da;i<=n;i++,j++){
            pro=pro*Pow(num[i],P-2)%P;
            *j=pro*INV%P;
            pro=pro*(num[i]-1)%P;
            }
            C=n*2;
            G=Pow(22,(P-1)/2/K),InvG=Pow(G,P-2);

        N=1; inv=Pow(K,P-2);
        for(;(1<<N)<(3*K);N++);
        int Nf=2*K,G=Pow(22,(P-1)/Nf),InvG=Pow(G,P-2);

        tmp[0]=tmp1[0]=1;
        for(int i=1;i<Nf;i++) tmp[i]=1LL*tmp[i-1]*G%P,tmp1[i]=1LL*tmp1[i-1]*InvG%P;
            for(int i=0,j=0;i<2*K;i++,j=(j+(i<<1)+1)%Nf) B1[i]=tmp[j],B2[i]=tmp1[j];
        for(int i=0,j=0;i<K;i++,j=nxt[i]) nxt[i+1]=((j+(i+1<<1)+1)%Nf);

            FWT(Da,n,m,1);
      for(int i=0;i<=n;i++) Da[i]=Pow(Da[i],x%(P-1));

      FWT(Da,n,m,-1);
      int ANS=0;
      for(int i=0;i<=n;i++)
        ANS=(ANS+1LL*Pow(i,2*i)*Pow(Da[i],3*i))%P;
      cout<<ANS<<endl;
      for(int i=0;i<=n;i++) num[i]=0;
    }
    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值