有A位小朋友，有三个共同的欢乐系数O，S和U。如果有一位小朋友得到了x个糖果，那么她的欢乐程度就是f(x)=Ox2+Sx+U。
现在校长开始分糖果了，一共有M个糖果。有些小朋友可能得不到糖果，对于那些得不到糖果的小朋友来说，欢乐程度就是1。如果一位小朋友得不到糖果，那么在她身后的小朋友们也都得不到糖果。（即这一列得不到糖果的小朋友一定是最后的连续若干位）
所有分糖果的方案都是等概率的。现在问题是：期望情况下，所有小朋友的欢乐程度的乘积是多少？呆呆同学很快就有了一个思路，只要知道总的方案个数T和所有方案下欢乐程度乘积的总和S，就可以得到答案Ans=S/T。现在他已经求出来了T的答案，但是S怎么求呢？他就不知道了。你能告诉他么？
因为答案很大，你只需要告诉他S对P取模后的结果。
M<=10000，P<=255，A<=108，O<=4，S<=300，U<=100。
Solution
将一个人的选择写成生成函数
[G(x)=1+(O+S+U)x+cdots+(Om^2+Sm+U)x^m
]
于是答案为
[sum_{i=1}^A[x^m]G^i(x)=[x^m]sum_{i=1}^AG^i(x)
]
考虑到这里模数不是质数，如果用多项式快速幂和多项式求逆处理会很麻烦，我们考虑一种类似倍增的思路
令 (n=a+1)，设三个辅助多项式，(R_i,P_i,Q_i)，有递推关系
[R_i=f^{2^i}quad P_i =(R_{i-1}+1)P_{i-1}quad Q_i=Q_{i-1}R_{i-1}^{n_i}
]
其中 (n_i) 表示 (n) 的第 (i) 个二进制位的值，统计答案时则有
[ans=sum_{n_i=1} P_iQ_i
]
为了方便理解，给出一个例子，(a=4,m=4) 的情况
| i | P | Q | R | ans |
|---|---|---|---|---|
| (0) | (1) | (1) | (f) | (1) |
| (1) | (1+f) | (f) | (f^2) | 1 |
| (2) | (1+f+f^2+f^3) | (f) | (f^4) | (1+f+f^2+f^3+f^4) |
| (3) | (1+f+f^2+f^3) | (f) | (f^4) | (1+f+f^2+f^3+f^4) |
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 100005;
namespace po {
const int mod = 998244353;
int qpow(int p,int q) {return (q&1?p:1)*(q?qpow(p*p%mod,q/2):1)%mod;}
int rev[N],inv[N],w[N],sz;
void presolve(int l) {
int len=1;
sz=0;
while(len<l) len<<=1, ++sz;
for(int i=1;i<len;i++) {
inv[i]=(i==1?1:inv[mod%i]*(mod-mod/i)%mod);
rev[i]=(rev[i>>1]>>1)|((i&1)<<(sz-1));
}
int wn=qpow(3,(mod-1)/len);
w[len/2]=1;
for(int i=len/2+1;i<len;i++) w[i]=w[i-1]*wn%mod;
for(int i=len/2-1;i;i--) w[i]=w[i<<1];
}
int pre(int l) {int g; for(g=1;g<l;g<<=1); return g;}
void ntt(int *a,int o,int n) {
static unsigned long long s[N];
int t=sz-__builtin_ctz(n),x;
for(int i=0;i<n;i++) s[rev[i]>>t]=a[i];
for(int l=1;l<n;l<<=1) for(int i=0;i<n;i+=l<<1) for(int j=0;j<l;j++) {
x=s[i+j+l]*w[j+l]%mod;
s[i+j+l]=s[i+j]+mod-x;
s[i+j]+=x;
}
for(int i=0;i<n;i++) a[i]=s[i]%mod;
if(o) {
x=qpow(n,mod-2);
for(int i=0;i<n;i++) a[i]=a[i]*x%mod;
reverse(a+1,a+n);
}
}
void mult(int n,int *x,int *y,int *z) {
static int a[N],b[N];
int l=pre(n<<1);
for(int i=0;i<l;i++) {
a[i]=(i<n?x[i]:0);
b[i]=(i<n?y[i]:0);
}
ntt(a,0,l); ntt(b,0,l);
for(int i=0;i<l;i++) z[i]=a[i]*b[i]%mod;
ntt(z,1,l);
for(int i=n;i<l;i++) z[i]=0;
}
}
int f[N],p[N],q[N],r[N],g[N],t[N],m,a,mod,o,s,u;
signed main() {
cin>>m>>mod>>a>>o>>s>>u;
a=min(a,m);
++a;
po::presolve(2*(m+1));
for(int i=1;i<=m;i++) r[i]=(o*i*i+s*i+u)%mod;
p[0]=1;
q[0]=1;
r[0]=0;
for(int i=0;i<20;i++) {
if((a>>i)&1) {
po::mult(m+1,p,q,t);
for(int i=0;i<=m;i++) g[i]=(g[i]+t[i])%mod;
po::mult(m+1,q,r,t);
for(int i=0;i<=m;i++) q[i]=t[i]%mod;
}
po::mult(m+1,r,p,t);
for(int i=0;i<=m;i++) p[i]=(p[i]+t[i])%mod;
po::mult(m+1,r,r,t);
for(int i=0;i<=m;i++) r[i]=t[i]%mod;
}
cout<<g[m];
}