CodeChef BINOMSUM

题目大意

假定有$T$天,每天有$K$个小时,你要安排每一天的日程。
第$i$天有$D+i−1$道菜,第一个小时你必须选择$L$道菜吃,接下来每个小时你可以选择吃一道菜或者选择$A$个活动中的一个参加,不能连续两个小时吃菜(但是一天的最后一个小时和下一天的第一个小时连着吃是可以的)。
$K,A$是事先给定的值。有$Q$次询问,每次给出$L,D,T$,问这$T$天每一天安排的方案数之和
答案对$P$(事先给定)取模。

$2\leq K\leq 10^5,1\leq Q\leq 500,1\leq L\leq D\leq D+T-1\leq10^7,10^8+7\leq P\leq10^9+7,P\in\mathbb P$


题目分析

首先根据题意很快就可以写出方案数的式子:
$$\sum_{i=1}^T{D+i+1\choose L}\sum_{j=0}^{K-1}\left(D+i-1\right)^jA^{K-j-1}{K-j-1\choose j}$$可以发现第二个$\Sigma$后面的是一个度数为$K/2$的多项式,但是我们后面懒得管这么多直接当$K$。
然后是一波操作:考虑将后面的多项式写成阶乘多项式的形式:
$$F(x)=\sum_{i=0}^{K-1}a_i\prod_{j=1}^i(x+j)$$也就是答案是
$$\sum_{i=D}^{D+T-1}{i\choose L}F(i)$$然后推一推
$$
\begin{align}
\sum_{i=D}^{D+T-1}{i\choose L}F(i)&=\sum_{i=D}^{D+T-1}{i\choose L}\sum_{j=0}^{K-1}a_j\prod_{k=1}^j(i+k)\\
&=\sum_{i=D}^{D+T-1}{i\choose L}\sum_{j=0}^{K-1}a_j\frac{(i+j)!}{i!}\\
&=\sum_{i=0}^{K-1}a_i\sum_{j=D}^{D+T-1}{j\choose L}\frac{(i+j)!}{j!}\\
&=\frac{\sum_{i=0}^{K-1}a_i\sum_{j=D}^{D+T-1}\frac{(i+j)!}{(j-L)!}}{L!}\\
&=\frac{\sum_{i=0}^{K-1}a_i(i+L)!\sum_{j=D}^{D+T-1}{i+j\choose i+L}}{L!}\\
&=\frac{\sum_{i=0}^{K-1}a_i(i+L)!\left({D+T+i\choose L+i+1}-{D+i\choose L+i+1}\right)}{L!}
\end{align}
$$如果我们知道了$\{a_n\}$,单次询问就可以直接$O(K)$计算。
然后我们来算这个$\{a_n\}$。显然有$F(x)=\sum_{i=0}^{K-1}a_i\cdot i!{x+i\choose i}$。
我们考虑求出$F(x)$在$-1,-2,…,-K$处的点值,目的是利用里面的负号凑出二项式反演的形式:
$$
\begin{align}
F(-n-1)&=\sum_{i=0}^{K-1}a_i\cdot i!{-n-1+i\choose i}\\
&=\sum_{i=0}^na_i\cdot i!(-1)^i{n\choose i}
\end{align}
$$
反演得到$$a_n\cdot n!=\sum_{i=0}^{K-1}(-1)^i{n\choose i}F(-i-1)$$也就是说$$a_n=\sum_{i=0}^n\frac{(-1)^iF(-i-1)}{i!}\frac1{(n-i)!}$$直接FFT做卷积就好了。
至于这个$F(x)$怎么快速地求值,观察$F(x)$原本对应的式子可以发现,只要将组合数用杨辉三角拆开,就能得到递推式:$F_k(x)=AF_{k-1}(x)+xAF_{k-2}(x)$。用矩阵快速幂就好了。
时间复杂度$O(K\log K+QK+(D+T))$。


代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cctype>
#include <cmath>
using namespace std;
typedef long double db;
typedef long long LL;
inline int read()
{
int x=0,f=1;
char ch=getchar();
while (!isdigit(ch)) f=ch=='-'?-1:f,ch=getchar();
while (isdigit(ch)) x=x*10+ch-'0',ch=getchar();
return x*f;
}
int buf[30];
inline void write(int x)
{
if (x<0) putchar('-'),x=-x;
for (;x;x/=10) buf[++buf[0]]=x%10;
if (!buf[0]) buf[++buf[0]]=0;
for (;buf[0];putchar('0'+buf[buf[0]--]));
}
const db pi=acos(-1);
const int N=15000000;
const int L=262144;
int fact[N+5],invf[N+5];
int K,A_,P,T,q,len;
inline int quick_power(int x,int y)
{
int ret=1;
for (;y;y>>=1,x=1ll*x*x%P) if (y&1) ret=1ll*ret*x%P;
return ret;
}
inline int sig(int x){return x&1?P-1:1;}
inline int C(int n,int m){return n>=m?1ll*fact[n]*invf[m]%P*invf[n-m]%P:0;}
void pre()
{
fact[0]=1;
for (int i=1;i<=N;++i) fact[i]=1ll*fact[i-1]*i%P;
invf[N]=quick_power(fact[N],P-2);
for (int i=N;i>=1;--i) invf[i-1]=1ll*invf[i]*i%P;
}
struct matrix
{
int num[2][2];
int r,c;
inline matrix(){r=c=2,num[0][0]=num[1][1]=1,num[0][1]=num[1][0]=0;}
inline matrix operator*(matrix const mat)
{
matrix ret;memset(ret.num,0,sizeof ret.num);
ret.r=r,ret.c=mat.c;
for (int i=0;i<ret.r;++i)
for (int j=0;j<ret.c;++j)
for (int k=0;k<c;++k)
(ret.num[i][j]+=1ll*num[i][k]*mat.num[k][j]%P)%=P;
return ret;
}
}trans,fk;
inline matrix operator^(matrix x,int y)
{
matrix ret=matrix();
for (;y;y>>=1,x=x*x) if (y&1) ret=ret*x;
return ret;
}
inline int Fk(int x)
{
fk.r=1,fk.c=2,fk.num[0][0]=0,fk.num[0][1]=1;
trans.r=trans.c=2,trans.num[0][0]=0,trans.num[0][1]=1ll*x*A_%P,trans.num[1][0]=1,trans.num[1][1]=A_;
return (fk*(trans^K)).num[0][0];
}
struct Z
{
db x,y;
inline Z(db x_=0.,db y_=0.){x=x_,y=y_;}
inline Z operator+(Z const z)const{return Z(x+z.x,y+z.y);}
inline Z operator-(Z const z)const{return Z(x-z.x,y-z.y);}
inline Z operator*(Z const z)const{return Z(x*z.x-y*z.y,x*z.y+y*z.x);}
inline Z operator*(db const k)const{return Z(x*k,y*k);}
}omega[L+5],A[2][L],B[2][L],cnt[4][L],t[L];
int f[L],g[L],a[L];
int trs[L];
void FFT_pre()
{
for (len=1;len<(K<<1)-1;len<<=1);
for (int i=0,ret;i<len;++i)
{
ret=0;
for (int x=i,j=1;j<len;j<<=1,x>>=1) ret=(ret<<1)|(x&1);
trs[i]=ret;
}
for (int i=0;i<=len;++i) omega[i]=Z(cos(2.*pi*i/len),sin(2.*pi*i/len));
}
void DFT(Z *a,int sig)
{
for (register int i=0;i<len;++i) t[trs[i]]=a[i];
for (register int l=2;l<=len;l<<=1)
for (register int j=0,h=l>>1;j<len;j+=l)
for (register int i=0;i<h;++i)
{
Z w=omega[sig>0?len/l*i:len-len/l*i],u=t[i+j],v=t[i+j+h]*w;
t[i+j]=u+v,t[i+j+h]=u-v;
}
for (register int i=0;i<len;++i) a[i]=t[i];
if (sig<0) for (register int i=0;i<len;++i) a[i]=a[i]*(1./len);
}
void mult(int *a,int *b,int *c)
{
for (register int i=0;i<len;++i) A[0][i]=Z(a[i]/T,0),A[1][i]=Z(a[i]%T,0);
for (register int i=0;i<len;++i) B[0][i]=Z(b[i]/T,0),B[1][i]=Z(b[i]%T,0);
DFT(A[0],1),DFT(A[1],1),DFT(B[0],1),DFT(B[1],1);
for (register int i=0;i<len;++i) cnt[0][i]=A[0][i]*B[0][i];
for (register int i=0;i<len;++i) cnt[1][i]=A[0][i]*B[1][i]+A[1][i]*B[0][i];
for (register int i=0;i<len;++i) cnt[2][i]=A[1][i]*B[1][i];
DFT(cnt[0],-1),DFT(cnt[1],-1),DFT(cnt[2],-1);
for (register int i=0;i<len;++i) c[i]=((LL)(cnt[0][i].x+.5)%P*T%P*T%P+(LL)(cnt[1][i].x+.5)%P*T%P+(LL)(cnt[2][i].x+.5)%P)%P;
}
void calc()
{
pre();
for (int i=0;i<K;++i) f[i]=1ll*Fk(P-i-1)*invf[i]%P*sig(i)%P,g[i]=invf[i];
FFT_pre(),mult(f,g,a);
}
int main()
{
freopen("binomsum.in","r",stdin),freopen("binomsum.out","w",stdout);
for (K=read(),A_=read(),P=read(),T=trunc(sqrt(P)),calc(),q=read();q--;)
{
int l=read(),d=read(),t_=read(),ans=0;
for (int i=0;i<K;++i) (ans+=1ll*a[i]*((C(d+i+t_,l+i+1)-C(d+i,l+i+1)+P)%P)%P*fact[i+l]%P)%=P;
ans=1ll*ans*invf[l]%P,write(ans),putchar('\n');
}
fclose(stdin),fclose(stdout);
return 0;
}