luogu5075 [JSOI2012]分零食
题目链接：luogu
首先有一个很显然的dp，记(f_{i,j})为分给前(i)个人(j)块糖的答案，那么有(f_{i,j}=sum_{k=0}^jf_{i-1,k} imes g_{j-k}).其中(g_i)表示分给一个人(i)块糖时的快乐度（特别的,(g_0=0)）
这个dp直接做很没有前途，但是发现它的转移是类似卷积的形式，于是考虑往生成函数的方向思考。
记(g_i)的生成函数为(G(x))，(f_{i,j})的生成函数为(F_i(x))，那么每次转移其实就是(F_i(x)=F_{i-1}(x)G(x)),初始条件为(F_0(x)=1), 最后的答案为(sum_{i=1}^m[x^m]F_i(x)).那么有(F_i(x)=G^i(x)).
事实上答案可以看成是等比数列求和的形式，于是可以用等比数列求和转化后使用多项式求逆+快速幂，可以做到一个(log)，但是模数是一个比较麻烦的问题。
考虑一个不依赖模数的做法：记(F_n(x)=sum_{i=1}^n G^i(x)). 使用类似快速幂的方法求(F_n(x)).
- (n)为偶数时
已知(F_{frac{n}{2}}(x),G^{frac{n}{2}}(x)).
[egin{aligned}
F_n(x)=&sum_{i=1}^{frac{n}{2}}G^i(x)+sum_{i=frac{n}{2}+1}^nG^i(x)\
=&F_{frac{n}{2}}(x)+G^{frac{n}{2}}(x)sum_{i=1}^{frac{n}{2}}F^i(x)\
=&F_{frac{n}{2}}(x)+G^{frac{n}{2}}(x)F_{frac{n}{2}}(x)\
G^n(x)=&G^{frac{n}{2}}(x)G^{frac{n}{2}}(x)
end{aligned}
]
- (n)为奇数时
发现我们上面求的其实是(n-1)的情况，把(n)的贡献单独加上去即可。
[G^n(x)=G^{n-1}(x)G^1(x)\
F_n(x)=F_{n-1}(x)+G^n(x)
]
#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<vector>
#include<bitset>
#include<math.h>
#include<stack>
#include<queue>
#include<set>
#include<map>
using namespace std;
typedef long long ll;
typedef long double db;
typedef vector<int> vi;
typedef pair<int,int> pii;
const int N=100000+100;
const db pi=acos(-1.0);
#define lowbit(x) (x)&(-x)
#define sqr(x) (x)*(x)
#define rep(i,a,b) for (register int i=a;i<=b;i++)
#define per(i,a,b) for (register int i=a;i>=b;i--)
#define go(u,i) for (register int i=head[u];i;i=sq[i].nxt)
#define fir first
#define sec second
#define mp make_pair
#define pb push_back
#define eps 1e-8
int maxd;
inline int read()
{
int x=0,f=1;char ch=getchar();
while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
return x*f;
}
namespace My_Math{
#define N 100000
int fac[N+100],invfac[N+100];
int add(int x,int y) {return x+y>=maxd?x+y-maxd:x+y;}
int dec(int x,int y) {return x<y?x-y+maxd:x-y;}
int mul(int x,int y) {return 1ll*x*y%maxd;}
ll qpow(ll x,int y)
{
ll ans=1;
while (y)
{
if (y&1) ans=mul(ans,x);
x=mul(x,x);y>>=1;
}
return ans;
}
int getinv(int x) {return qpow(x,maxd-2);}
int C(int n,int m)
{
if ((n<m) || (n<0) || (m<0)) return 0;
return mul(mul(fac[n],invfac[m]),invfac[n-m]);
}
void math_init()
{
fac[0]=invfac[0]=1;
rep(i,1,N) fac[i]=mul(fac[i-1],i);
invfac[N]=getinv(fac[N]);
per(i,N-1,1) invfac[i]=mul(invfac[i+1],i+1);
}
#undef N
}
using namespace My_Math;
int n,m,O,S,U,f1[N<<2],f[N<<2],g[N<<2],h[N<<2];
namespace Polynomial{
struct complex{
double x,y;
complex(double _x=0.0,double _y=0.0) {x=_x;y=_y;}
};
complex operator +(complex a,complex b)
{
return complex(a.x+b.x,a.y+b.y);
}
complex operator -(complex a,complex b)
{
return complex(a.x-b.x,a.y-b.y);
}
complex operator *(complex a,complex b)
{
return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
complex A[N<<2],B[N<<2];
int r[N<<2];
void fft(int lim,complex *a,int typ)
{
rep(i,0,lim-1)
if (i<r[i]) swap(a[i],a[r[i]]);
for (int mid=1;mid<lim;mid<<=1)
{
int len=(mid<<1);
complex wn=complex(cos(pi/mid),sin(pi/mid)*typ);
for (int sta=0;sta<lim;sta+=len)
{
complex w=complex(1,0);
for (int j=0;j<mid;j++,w=w*wn)
{
complex x=a[sta+j],y=a[sta+j+mid]*w;
a[sta+j]=x+y;a[sta+j+mid]=x-y;
}
}
}
}
int calcr(int len)
{
int lim=1,cnt=0;
while (lim<len) {lim<<=1;cnt++;}
rep(i,0,lim-1)
r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
return lim;
//rep(i,0,lim-1) cout << r[i] << " ";cout << endl;
}
void mul(int *a,int *b,int *c,int lim)
{
//rep(i,0,lim-1) cout << a[i] << " ";cout << endl;
//rep(i,0,lim-1) cout << b[i] << " ";cout << endl;
rep(i,0,lim-1)
{
A[i]=complex(a[i],0);
B[i]=complex(b[i],0);
}
fft(lim,A,1);fft(lim,B,1);
rep(i,0,lim-1) A[i]=A[i]*B[i];
fft(lim,A,-1);
rep(i,1,m) c[i]=(int)(A[i].x/lim+0.5)%maxd;
//rep(i,1,m) cout << c[i] << " ";cout << endl << endl;
}
void solve(int k,int lim)
{
if (k==1)
{
rep(i,0,lim-1) f[i]=g[i]=f1[i];
return;
}
solve(k>>1,lim);
mul(f,g,h,lim);
rep(i,0,lim-1) g[i]=add(g[i],h[i]);
mul(f,f,f,lim);
if (k&1)
{
mul(f,f1,f,lim);
rep(i,0,lim-1)
g[i]=add(g[i],f[i]);
}
}
}
using namespace Polynomial;
int calc(int x) {return (1ll*O*sqr(x)+S*x+U)%maxd;}
int main()
{
m=read();maxd=read();n=read();O=read();S=read();U=read();
rep(i,1,m) f1[i]=calc(i);
int lim=calcr(m<<1);
solve(min(n,m),lim);
printf("%d",g[m]);
return 0;
}