ZJOI2014 力

题目:

  给出 个数 ,给出 的定义如下:
  
  令 ,求 .

思路:

  由题意得: ,因为 ,所以可以考虑列出两个多项式相乘 。因为两次对答案的贡献系数不同(一个是 ,另一个是 ),所以分别做两次FFT即可。

代码:

#include <bits/stdc++.h>
#define pw(x) ((ll)(x)*(x))
#define mset(a,b) memset(a,b,sizeof a)
#define _(...) (void)(__VA_ARGS__)
#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 double PI=acos(-1);
const int N=(1<<18)+5;

struct Complex{
    double r,i;
    Complex operator + (const Complex _) const {
        return (Complex){r+_.r,i+_.i};
    }
    Complex operator - (const Complex _) const {
        return (Complex){r-_.r,i-_.i};
    }
    Complex operator * (const Complex _) const {
        return (Complex){r*_.r-i*_.i,r*_.i+i*_.r};
    }
};

int rev[N];
void fft(Complex a[],int n,bool inv){
    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){
        Complex base=(Complex){cos(PI*2/l),(inv?-1:1)*sin(PI*2/l)};
        for(Complex *p=a;p!=a+n;p+=l){
            Complex w=(Complex){1,0};
            rep(i,0,m-1){
                Complex t=w*p[i+m];
                p[i+m]=p[i]-t;
                p[i]=p[i]+t;
                w=w*base;
            }
        }
    }
    if(inv)
        rep(i,0,n-1)
            a[i].r/=n;
}

int log2(int x){
    int lim=0;
    while((1<<lim)<x)lim++;
    return lim;
}

Complex a[N],b[N];
double p[N],res[N];
int main(){
    int n;
    scanf("%d",&n);
    rep(i,1,n){
        scanf("%lf",p+i);
        a[i].r=p[i];
    }

    int lim=log2((n+1)*2),d=1<<lim;
    rep(i,0,d-1)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(lim-1));

    fft(a,d,false);
    rep(i,1,n-1)
        b[n-i].r=1.0/pw(i);
    fft(b,d,false);
    rep(i,0,d-1)
        b[i]=a[i]*b[i];
    fft(b,d,true);
    rep(i,n+1,n*2)
        res[i-n]-=b[i].r;

    mset(b,0);
    rep(i,1,n-1)
        b[i].r=1.0/pw(i);
    fft(b,d,false);
    rep(i,0,d-1)
        b[i]=a[i]*b[i];
    fft(b,d,true);
    rep(i,1,n)
        res[i]+=b[i].r;

    rep(i,1,n){
        if(fabs(res[i])<1e-8)puts("0.0000");
        else printf("%.4lf\n",res[i]);
    }
    return 0;
}