bzoj 4332: JSOI2012 分零食 快速傅立葉變換


題目:

Description
同學們依次排成了一列,其中有A位小朋友,有三個共同的歡樂系數O,S和U。如果有一位小朋友得到了x個糖果,那么她的歡樂程度就是\(f(x)=O*x^2+S*x+U\)
現在校長開始分糖果了,一共有M個糖果。有些小朋友可能得不到糖果,對於那些得不到糖果的小朋友來說,歡樂程度就是1。如果一位小朋友得不到糖果,那么在她身后的小朋友們也都得不到糖果。(即這一列得不到糖果的小朋友一定是最后的連續若干位)
所有分糖果的方案都是等概率的。現在問題是:期望情況下,所有小朋友的歡樂程度的乘積是多少?呆呆同學很快就有了一個思路,只要知道總的方案個數T和所有方案下歡樂程度乘積的總和S,就可以得到答案Ans=S/T。現在他已經求出來了T的答案,但是S怎么求呢?他就不知道了。你能告訴他么?
因為答案很大,你只需要告訴他S對P取模后的結果。

題解:

首先這道題我們可以考慮枚舉一下有多少人得到了零食
\(g[i][j]\)表示\(i\)個人里分下去了\(j\)個零食得到的值,n為人數,m為零食數
這樣我們有\(ans = \sum_{i=1}^ng[i][m]\)
\(g[i][j]\)的遞推我們有
\[g[i][j] = \sum_{k=1}^{j-1}g[i-1][j-k]*F(k)\]
其中\(F(k)\)表示將\(k\)個零食分給一個人得到的權
然后我們驚奇地發現后面的式子是一個卷積的形式
所以我們可以得到\(g_i = g_{i-1}*F\)
由於卷積滿足結合律,所以我們有\(g_i = g_0*F^i\)
這樣的話我們能夠完成\(nlogn\)的單點求值,但是我們要求的是\(g_i\)的一個和
我們把答案\(ans\)拓展為一個多項式\(f(x)\),答案儲存在第f(n)的第\(m\)
則有\(f_n = \sum_{i=1}^{n}g_i\)
然后我們發現:
\[f_n = f_{\frac{n}{2}} + \sum_{i=\frac{n}{2}+1}^{n}g_i\]
\[f_n = f_{\frac{n}{2}} + \sum_{i=1}^{\frac{n}{2}}g_{\frac{n}{2}+i}\]

在繼續推導之前首先我們需要證明: \(g_{i+j} = g_i*g_j\)
\(g_i = g_0*F^i\)可得:\(g_i*g_j = g_0*g_0*F^i*F^j\)
因為:\(g_0*g_0 = g_0\)所以有\(g_0*g_0*F^i*F^j = g_0*F^i*F^j = g_{i+j}\)得證

所以繼續上式的推導我們有
\[f_n = f_{\frac{n}{2}} + \sum_{i=1}^{\frac{n}{2}}g_{\frac{n}{2}}*g_i\]
我們將卷積的形式再拆解開來:
\[f_n = f_{\frac{n}{2}} + \sum_{i=1}^{\frac{n}{2}}\sum_{j=1}^{m-1}g_{\frac{n}{2},m-j}*g_{i,j}\]
\[f_n = f_{\frac{n}{2}} + \sum_{j=1}^{m-1}g_{\frac{n}{2},m-j}*\sum_{i=1}^{\frac{n}{2}}g_{i,j}\]
\[f_n = f_{\frac{n}{2}} + \sum_{j=1}^{m-1}g_{\frac{n}{2},m-j}*f_{\frac{n}{2},j}\]
\[f_n = f_{\frac{n}{2}} + g_{\frac{n}{2}}*f_{\frac{n}{2}}\]
我們又知道:
\[g_n = g_{\frac{n}{2}}*g_{\frac{n}{2}}\]
所以迭代倍增即可求解
復雜度\(O(nlog^2n)\)

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
inline void read(int &x){
    x=0;char ch;bool flag = false;
    while(ch=getchar(),ch<'!');if(ch == '-') ch=getchar(),flag = true;
    while(x=10*x+ch-'0',ch=getchar(),ch>'!');if(flag) x=-x;
}
const int maxn = 41000;
const double pi = acos(-1);
int mod;
struct complex{
    long double x,y;
    complex(){x=y=0;}
    complex(long double a,long double b){x=a;y=b;}
    complex operator + (const complex &r){return complex(x+r.x,y+r.y);}
    complex operator - (const complex &r){return complex(x-r.x,y-r.y);}
    complex operator * (const complex &r){return complex(x*r.x-y*r.y,x*r.y+y*r.x);}
    complex operator / (const long double &r){return complex(x/r,y/r);}
};
inline void FFT(complex *x,int n,int p){
    for(int i=0,t=0;i<n;++i){
        if(i > t) swap(x[i],x[t]);
        for(int j=n>>1;(t^=j) < j;j>>=1);
    }
    for(int m=2;m<=n;m<<=1){
        int k = m>>1;
        complex wn(cos(p*2*pi/m),sin(p*2*pi/m));
        for(int i=0;i<n;i+=m){
            complex w(1,0),u;
            for(int j=0;j<k;++j,w=w*wn){
                u = x[i+j+k]*w;
                x[i+j+k] = x[i+j] - u;
                x[i+j] = x[i+j] + u;
            }
        }
    }
    if(p == -1 ) for(int i=0;i<n;++i) x[i] = x[i]/n;
}
complex ca[maxn],cb[maxn],cc[maxn];int len,m;
inline int mul(int *a,int *b,int *c){
    for(int i=0;i<len;++i){
        ca[i] = complex((long double)a[i],0);
        cb[i] = complex((long double)b[i],0);
    }
    FFT(ca,len,1);FFT(cb,len,1);
    for(int i=0;i<len;++i) cc[i] = ca[i]*cb[i];
    FFT(cc,len,-1);
    for(int i=0;i<=m;++i){
        c[i] = ((int)floor(cc[i].x + 0.5)) % mod;
    }
}
int f[maxn],g[maxn],arr[maxn],tmp[maxn];
inline void qpow(int k){
    if(k == 1){
        for(int i=0;i<=m;++i) f[i] = g[i] = arr[i];
        return ;
    }qpow(k>>1);
    mul(f,g,tmp);mul(g,g,g);
    for(int i=0;i<=m;++i){
        f[i] += tmp[i];
        if(f[i] >= mod) f[i] -= mod;
    }
    if(k&1){
        mul(g,arr,g);
        for(int i=0;i<=m;++i){
            f[i] += g[i];
            if(f[i] >= mod) f[i] -= mod;
        }
    }
}
int main(){
    read(m);read(mod);
    for(len = 1;(len) <= (m<<1);len<<=1);
    int n,a,b,c;read(n);read(a);read(b);read(c);
    a %= mod;b %= mod;c %= mod;
    for(int i=1;i<=m;++i) arr[i] = ((a*i*i % mod) + (b*i % mod) + c) % mod;
    qpow(n);printf("%d\n",f[m]);
    getchar();getchar();
    return 0;
}

注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
  © 2014-2022 ITdaan.com 联系我们: