luogu5075 [JSOI2012]分零食
题目链接：luogu
首先有一个很显然的dp，记\(f_{i,j}\)为分给前\(i\)个人\(j\)块糖的答案，那么有\(f_{i,j}=\sum_{k=0}^jf_{i-1,k}\times 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)\).
已知\(F_{\frac{n}{2}}(x),G^{\frac{n}{2}}(x)\).
发现我们上面求的其实是\(n-1\)的情况，把\(n\)的贡献单独加上去即可。
#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;
}
原文：https://www.cnblogs.com/encodetalker/p/12670689.html