hdu_4656 Evaluation

题目:

  
  
  Given , calculate .

题目:

  由题意得:
  
  把 强行二项式展开:
  
  稍微整理一下:
  
  设 ,观察发现 可以通过一次多项式乘法得到,所以可以先处理出 ,接下来的推导中都认为 是常量。
  再把常量全部提出来,设 ,那么:
  
  接下来就是整个推导最神的地方了:
  
  
  这样的话,就可以分成两个次数分别为 的多项式,相乘即可。
  题目中给定的模数不满足NTT,所以要用任意模数的NTT。我用的是拆系数的写法,大致原理是:还是像普通的FFT一样相乘,通过把系数拆成两位来提高精度,但是板子有学长的大量优化,所以不可读。背下来当个黑箱算法.....复杂度和普通的NTT一样:

代码:

#include <bits/stdc++.h>
#define debug(x) 0
#define pw(x) ((ll)(x)*(x))
#define mset(a,b) memset(a,b,sizeof a)
#define rep(i,a,b) for(int i(a),i##_END_(b);i<=i##_END_;i++)
#define drep(i,a,b) for(int i(a),i##_END_(b);i>=i##_END_;i--)
template<class T> inline bool tomax(T &a,T b){return a<b?a=b,1:0;}
template<class T> inline bool tomin(T &a,T b){return b<a?a=b,1:0;}
typedef long long ll;
const int N=524288+5;
const double PI=acos(-1);

namespace FFT{
    struct Comp{
        double r,i;
        inline Comp operator+ (const Comp _) const {
            return (Comp){r+_.r,i+_.i};
        }
        inline Comp operator- (const Comp _) const {
            return (Comp){r-_.r,i-_.i};
        }
        inline Comp operator* (const Comp _) const {
            return (Comp){r*_.r-i*_.i,r*_.i+i*_.r};
        }
    };
    int rev[N];
    void DFT(Comp a[],Comp w[],int n){
        rep(i,0,n-1)
            if(i<rev[i])
                std::swap(a[i],a[rev[i]]);
        for(int l=2,m=1;l<=n;l<<=1,m<<=1)
            for(Comp *p=a;p!=a+n;p+=l)
                rep(i,0,m-1){
                    Comp t=w[n/l*i]*p[i+m];
                    p[i+m]=p[i]-t;
                    p[i]=p[i]+t;
                }
    }
    int log2(int x){
        int r=0;
        while((1<<r)<x)r++;
        return r;
    }
    Comp a[N],b[N],c[N],d[N],omg[N];
    void calc(int _a[],int _b[],int res[],int n,int m,int mod){
        int lim=log2(n+m),r=1<<lim;
        rep(i,0,r-1)a[i]=b[i]=c[i]=d[i]=(Comp){0,0};
        rep(i,0,n-1)a[i]=(Comp){_a[i]&32767,_a[i]>>15};
        rep(i,0,m-1)b[i]=(Comp){_b[i]&32767,_b[i]>>15};
        rep(i,0,r-1){
            omg[i]=(Comp){cos(PI*2/r*i),sin(PI*2/r*i)};
            rev[i]=(rev[i>>1]>>1)|((i&1)<<(lim-1));
        }
        DFT(a,omg,r);DFT(b,omg,r);
        rep(i,0,r-1){
            int j=(r-i)&(r-1);
            c[j]=(Comp){(a[i].r+a[j].r)/2,(a[i].i-a[j].i)/2}*b[i];
            d[j]=(Comp){(a[i].i+a[j].i)/2,(a[j].r-a[i].r)/2}*b[i];
        }
        DFT(c,omg,r);DFT(d,omg,r);
        rep(i,0,n+m-2){
            ll u=c[i].r/r+0.5,v=c[i].i/r+0.5;
            ll x=d[i].r/r+0.5,y=d[i].i/r+0.5;
            u%=mod;v%=mod;x%=mod;y%=mod;
            res[i]=(u+(v+x<<15)+(y<<30))%mod;
        }
        mset(a,0);mset(b,0);mset(c,0);mset(d,0);
    }
}

int fac[N],ifac[N];
const int MOD=1e6+3;
int qkpow(int x,int k){
    int ans=1;
    for(;k;k>>=1,x=(ll)x*x%MOD)
        if(k&1)ans=(ll)ans*x%MOD;
    return ans;
}
inline int inv(int x){return qkpow(x,MOD-2);}
void init(){
    int n=N-5;
    fac[0]=1;
    rep(i,1,n)fac[i]=(ll)fac[i-1]*i%MOD;
    ifac[n]=inv(fac[n]);
    drep(i,n-1,0)ifac[i]=(ll)ifac[i+1]*(i+1)%MOD;
}

int a[N],g[N],p[N],q[N],res[N],ans[N];
int main(){
    init();
    int n,b,c,d;
    scanf("%d%d%d%d",&n,&b,&c,&d);
    rep(i,0,n-1)scanf("%d",a+i);
    rep(i,0,n-1)p[i]=(ll)a[i]*fac[i]%MOD;
    rep(i,0,n-1)q[n-i]=(ll)ifac[i]*qkpow(d,i)%MOD;
    FFT::calc(p,q,res,n,n+1,MOD);
    rep(i,0,n-1)g[i]=res[i+n];
    p[0]=0;
    rep(i,0,n-1){
        p[n-i]=(ll)g[i]*qkpow(b,i)%MOD;
        p[n-i]=(ll)p[n-i]*ifac[i]%MOD*inv(qkpow(c,pw(i)%(MOD-1)))%MOD;
    }
    rep(i,0,2*n-1)
        q[i]=qkpow(c,pw(i)%(MOD-1));
    FFT::calc(p,q,res,n+1,2*n,MOD);
    rep(i,0,n-1)
        ans[i]=(ll)res[i+n]*inv(qkpow(c,pw(i)%(MOD-1)))%MOD;
    rep(i,0,n-1)
        printf("%d\n",ans[i]);
    return 0;
}