FFT 【JSOI2012】bzoj4332 分零食


題目大意:

有n個小朋友,m塊糖。
給小朋友分糖,如果一個小朋友分不到糖,那他后面的小朋友也分不到糖。
每個小朋友有一個喜悅值,有三個參數,O,S,U,設一個小朋友分到糖數為x,則這個小朋友的喜悅值為O*x x+ S x +U,分不到糖的小朋友的喜悅值為1。
求所有分糖方案下 所有小朋友喜悅值乘積 的和。

題目分析:

如果沒有小朋友必須是前面的連續的一段這個要求,那就是一個FFT模板,把喜悅值數組自乘n次即可。
但是有這個要求之后就不能這么做了。

這個時候我們考慮DP。
設數組g[i][j]為給前i個小朋友分j塊糖(每個小朋友都能分到糖)的答案。
設數組f[i]為一個小朋友分到i塊糖的喜悅值。

那么轉移方程就可以寫成g[i][j]=∑(1<=k<=m)g[i-1][j-k]*f[k]。
這樣就是標准的卷積的形式,就可以用快速傅里葉變換來優化,可以用快速冪迅速求出g[min(n,m)][m],括弧但是我們並沒有用到快速冪括回笑。

我們求出g[min(n,m)][m]並沒有什么用,因為不一定會分到第n個小朋友。

所以我們要求的是∑(1<=i<=min(n,m))g[i][m]。
設數組p[i][j]=∑(1<=k<=i)g[k][j]。
我們只要求出p[min(n,m)][m]即可。

p可以遞推來求
當 i mod 2 為0時 ,p[i]=p[i/2]+p[i/2]*g[i/2];
當 i mod 2 為1時 ,p[i]=p[i/2]+p[i/2]*g[i/2]+g[i];

這個可以感性的證明一下:
首先g[a+b]=g[a]*g[b]。
給a+b個小朋友分糖,就相當於給前a個小朋友分糖和給a+1到a+b個小朋友分糖的答案求一個卷積,在點值表達式的意義下計算是合法的。

那么p[i/2]為g[1]到g[i/2]的和,把p[i/2]乘以g[i/2],那我們就得到了g[i/2+1]到g[i]的和,再加上p[i/2],就是所求的p[i]。

當i為奇數時,我們就先求p[i-1],然后再加上g[i]就可以了,括弧還是挺好理解的吧括回笑。

設top=min(n,m).
我們要求p[top]和g[top],那就要求出p[top/2] 和 g[top/2]。
那就可以遞歸處理啦。

時間復雜度為O(mlogmlogmin(n,m))。

此篇題解部分借鑒自sengxian’s blog.

注意:

  1. 亘古不變的精度問題!!!!
  2. 在乘的時候不要連乘啊,要不然會變成循環卷積啊什么的雲雲,會錯,就是說每次在做乘法的時候先DFT過去,乘完之后立刻IDFT回來,然后把超過m的部分全部去掉,再進行下一次乘法。

代碼:

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <iostream>
#define N 1200000
using namespace std;
const double pi=acos(-1);
const double DFT=2.0;
const double IDFT=-2.0;
long long n,m,O,S,U,mod;
int len;
int pos[N];
struct complex{
    double a,b;
    complex(double a=0,double b=0):a(a),b(b){}
    complex operator + (const complex &c) { return complex(a+c.a,b+c.b); }
    complex operator - (const complex &c) { return complex(a-c.a,b-c.b); }
    complex operator * (const complex &c) { return complex(a*c.a-b*c.b,a*c.b+b*c.a); }
}g[N],p[N],f[N],tmp[N];
void init(int len)
{
    for(int i=1;i<len;i++)
    {
        pos[i]=pos[i>>1]>>1;
        if(i&1) pos[i]|=len>>1;
    }
    return;
}
void Fast_Fourier_Transform(complex x[],int len,double mode)
{
    for(int i=1;i<len;i++)
    if(i<pos[i]) swap(x[i],x[pos[i]]);
    for(int i=2;i<=len;i<<=1)
    {
        int step=i>>1;
        complex wm=complex(cos(2*pi/i),sin(mode*pi/i));
        for(int j=0;j<len;j+=i)
        {
            complex w=complex(1,0);
            for(int k=j;k<j+step;k++)
            {
                complex a=x[k],b=w*x[k+step];
                x[k]=a+b;x[k+step]=a-b;
                w=w*wm;
            }
        }
    }
    if(mode==IDFT)
        for(int i=0;i<len;i++) x[i].a=(x[i].a)/len;
    return;
}
#define FFT Fast_Fourier_Transform
void modify(complex x[])
{
    for(int i=0;i<len;i++)
    {
        if(i<=m) x[i].a=((long long)(x[i].a+0.5)%mod+mod)%mod,x[i].b=0;
        else x[i]=complex(0,0);
    }
}
void solve(long long top)
{
    if(top==1)
    {
        for(int i=0;i<len;i++) p[i]=f[i],g[i]=f[i];
        FFT(p,len,IDFT); FFT(g,len,IDFT);
        modify(p); modify(g);
        return;
    }
    solve(top/2);
    FFT(p,len,DFT); FFT(g,len,DFT);
    for(int i=0;i<len;i++) p[i]=p[i]+p[i]*g[i];
    for(int i=0;i<len;i++) g[i]=g[i]*g[i];
    FFT(p,len,IDFT); FFT(g,len,IDFT);
    modify(p); modify(g);
    if(top&1)
    {
        FFT(g,len,DFT); FFT(p,len,DFT);
        for(int i=0;i<len;i++) g[i]=g[i]*f[i];
        for(int i=0;i<len;i++) p[i]=p[i]+g[i];
        FFT(g,len,IDFT); FFT(p,len,IDFT);
        modify(p); modify(g);
    }
    return;
}
int main()
{
    scanf("%lld%lld",&m,&mod);
    scanf("%lld",&n);
    scanf("%lld%lld%lld",&O,&S,&U);
    for(int i=1;i<=m;i++) f[i].a=(O*i*i+S*i+U)%mod;
    long long top=min(n,m);
    len=1;
    while(len<=(m<<1)) len<<=1;
    init(len);
    FFT(f,len,DFT);
    solve(top);
    printf("%lld\n",((long long)(p[m].a))%mod);
    return 0;
}

注意!

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



 
粤ICP备14056181号  © 2014-2020 ITdaan.com