【题目】
LG
给定两个序列
a
,
b
a,b
a,b,定义
k
k
k次价值为分别随机从
a
,
b
a,b
a,b中取出一个数
a
x
,
b
x
a_x,b_x
ax,bx,
(
a
x
+
b
x
)
k
(a_x+b_x)^k
(ax+bx)k的期望。
对于所有
k
∈
[
1
,
t
]
k\in[1,t]
k∈[1,t],求序列的
k
k
k次价值。
∣
a
∣
,
∣
b
∣
,
k
≤
1
0
5
|a|,|b|,k\leq 10^5
∣a∣,∣b∣,k≤105
【解题思路】
化柿子吧,这个形式显然就是二项式展开一下吧。
a
n
s
k
⋅
n
m
=
∑
i
=
1
n
∑
j
=
1
m
(
a
i
+
b
j
)
k
=
∑
x
=
0
k
∑
i
=
1
n
∑
j
=
1
m
(
x
k
)
a
i
x
b
j
k
−
x
=
k
!
×
∑
x
=
0
k
(
∑
i
=
1
n
a
i
x
x
!
)
(
∑
j
=
1
m
b
j
k
−
x
(
k
−
x
)
!
)
\begin{aligned} ans_k \cdot nm=& \sum_{i=1}^n\sum_{j=1}^m(a_i+b_j)^k\\ =& \sum_{x=0}^k\sum_{i=1}^n\sum_{j=1}^m{x\choose k}a_i^xb_j^{k-x}\\ =& k!\times \sum_{x=0}^k(\sum_{i=1}^n\frac {a_i^x}{x!})(\sum_{j=1}^m\frac {b_j^{k-x}} {(k-x)!}) \end{aligned}
ansk⋅nm===i=1∑nj=1∑m(ai+bj)kx=0∑ki=1∑nj=1∑m(kx)aixbjk−xk!×x=0∑k(i=1∑nx!aix)(j=1∑m(k−x)!bjk−x)
右边是一个卷积,那么我们只需要对于每个
k
k
k,求出一个数列每个数字的
k
k
k次方和。
这其实是一个比较经典的问题,写出关于
k
k
k的生成函数,则显然为:
∑
i
=
1
n
1
1
−
a
i
x
\sum_{i=1}^n\frac 1 {1-a_ix}
i=1∑n1−aix1
这个问题有一个比较好看的做法:
我们现在不妨求一下
∑
i
=
1
n
a
i
1
−
a
i
x
\sum_{i=1}^n\frac {a_i} {1-a_ix}
i=1∑n1−aixai
实际上这个东西求出来
0
0
0次项系数显然就是原来柿子的
1
1
1次项系数。
对它每个部分上下除以
a
i
a_i
ai,然后通分,令
r
i
=
1
a
i
r_i=\frac 1 {a_i}
ri=ai1可以得到:
∑
i
=
1
n
∏
j
≠
i
(
r
j
−
x
)
∏
i
=
1
n
(
r
i
−
x
)
\frac {\sum_{i=1}^n\prod _{j\neq i} (r_j-x)}{\prod_{i=1}^n (r_i-x)}
∏i=1n(ri−x)∑i=1n∏j̸=i(rj−x)
你会发现这个函数的分子是分母的导数。
然后分治
FFT
+
\text{FFT}+
FFT+多项式求逆,就可以得到整个东西了。只不过这个做法无法得到
0
0
0次项系数,不过这不重要,全部加一次就可以了。
复杂度
O
(
n
log
2
n
)
O(n\log ^2n)
O(nlog2n)
更直观的理解?令上面我们要求的和为
f
(
x
)
f_(x)
f(x)
可以发现
ln
′
(
1
−
a
i
x
)
=
−
a
i
1
−
a
i
x
\ln'(1-a_ix)=\frac {-a_i} {1-a_ix}
ln′(1−aix)=1−aix−ai,
则令
g
(
x
)
=
∑
i
=
1
n
ln
′
(
1
−
a
i
x
)
=
(
ln
(
∏
i
=
1
n
(
1
−
a
i
x
)
)
)
′
g(x)=\sum_{i=1}^n\ln' (1-a_ix)=(\ln(\prod_{i=1}^n(1-a_ix)))'
g(x)=∑i=1nln′(1−aix)=(ln(∏i=1n(1−aix)))′,有
f
(
x
)
=
−
x
×
g
(
x
)
+
n
f(x)=-x\times g(x)+n
f(x)=−x×g(x)+n
那么
g
g
g可以分治
FFT
\text{FFT}
FFT,就求完了。
对了其实原来那个柿子直接分治 FFT \text{FFT} FFT也可以做的,不过常数巨大而已。
我写了我提到的两种,其中一个注释掉了。
【参考代码】
#include<bits/stdc++.h>
using namespace std;
const int N=262333,mod=998244353,G=3;
namespace IO
{
int read()
{
int ret=0;char c=getchar();
while(!isdigit(c)) c=getchar();
while(isdigit(c)) ret=ret*10+(c^48),c=getchar();
return ret;
}
void write(int x){if(x>9)write(x/10);putchar(x%10^48);}
void writeln(int x){write(x);putchar('\n');}
}
using namespace IO;
namespace Math
{
int fac[N],ifac[N],inv[N];
int upm(int x){return x>=mod?x-mod:(x<0?x+mod:x);}
void up(int &x,int y){x=upm(x+y);}
int mul(int x,int y){return 1ll*x*y%mod;}
int qpow(int x,int y){int res=1;for(;y;y>>=1,x=mul(x,x))if(y&1)res=mul(res,x);return res;}
int getinv(int x){return qpow(x,mod-2);}
void initmath()
{
fac[0]=1;for(int i=1;i<N;++i)fac[i]=mul(fac[i-1],i);
ifac[N-1]=getinv(fac[N-1]);for(int i=N-2;~i;--i)ifac[i]=mul(ifac[i+1],i+1);
inv[1]=1;for(int i=2;i<N;++i)inv[i]=mod-mul(mod/i,inv[mod%i]);
}
}
using namespace Math;
namespace Poly
{
int m,L,rev[N];
void ntt(int *a,int n,int f)
{
for(int i=0;i<n;++i) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1)
{
int wn=qpow(G,(mod-1)/(i<<1));
if(!~f) wn=getinv(wn);
for(int j=0;j<n;j+=i<<1)
{
int w=1;
for(int k=0;k<i;++k,w=mul(w,wn))
{
int x=a[j+k],y=mul(w,a[i+j+k]);
a[j+k]=upm(x+y);a[i+j+k]=upm(x-y);
}
}
}
if(!~f) for(int i=0,iv=getinv(n);i<n;++i) a[i]=mul(a[i],iv);
}
void reget(int n)
{
for(m=1,L=0;m<=n;m<<=1,++L);
for(int i=0;i<m;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
}
void polymul(int *a,int *b,int *c)
{
static int A[N],B[N];
copy(a,a+m,A);copy(b,b+m,B);
ntt(A,m,1);ntt(B,m,1);
for(int i=0;i<m;++i) c[i]=mul(A[i],B[i]);
ntt(c,m,-1);
fill(A,A+m,0);fill(B,B+m,0);
}
void polymult(int *a,int *b,int *c,int dega,int degb)
{
static int A[N];
reget(dega+degb);polymul(a,b,A);
for(int i=0;i<m;++i) c[i]=A[i];
fill(A,A+m,0);
}
void polyder(int *a,int deg)
{
for(int i=0;i<deg-1;++i) a[i]=mul(a[i+1],i+1);
a[deg-1]=0;
}
void polyint(int *a,int deg)
{
for(int i=deg-1;i;--i) a[i]=mul(a[i-1],inv[i]);
a[0]=0;
}
void polyinv(int *a,int *b,int deg)
{
if(deg==1){b[0]=getinv(a[0]);return;}
polyinv(a,b,(deg+1)>>1);
static int A[N],B[N];
reget(deg*2);copy(a,a+deg,A);copy(b,b+deg,B);fill(A+deg,A+m,0);fill(B+deg,B+m,0);
ntt(A,m,1);ntt(B,m,1);
for(int i=0;i<m;++i) A[i]=mul(A[i],mul(B[i],B[i]));
ntt(A,m,-1);
for(int i=0;i<deg;++i) b[i]=upm(mul(2,b[i])-A[i]);
fill(A,A+m,0);fill(B,B+m,0);
}
void polyln(int *a,int *b,int deg)
{
static int A[N],B[N];
polyinv(a,A,deg);copy(a,a+deg,B);polyder(B,deg);
reget(deg*2);polymul(A,B,A);
for(int i=0;i<deg;++i) b[i]=A[i];
polyint(b,deg);
fill(A,A+m,0);fill(B,B+m,0);
}
}
using namespace Poly;
namespace DreamLolita
{
int n,m,T;
int f[25][N],coe[N],A[N],B[N],C[N],a[N],b[N];
/*void solve(int l,int r,int dep)
{
if(l==r){f[dep][0]=1;f[dep][1]=mod-coe[l];return;}
int mid=(l+r)>>1;
solve(l,mid,dep);solve(mid+1,r,dep+1);
polymult(f[dep],f[dep+1],f[dep],r-l+2,0);
fill(f[dep+1],f[dep+1]+Poly::m,0);//wrong because write m
}
void GETK(int *a,int *b,int n)
{
memset(f,0,sizeof(f));memset(coe,0,sizeof(coe));
for(int i=1;i<=n;++i) coe[i]=a[i];
solve(1,n,0);
polyln(f[0],b,T+1);
polyder(b,T+1);
for(int i=T;i;--i) b[i]=mod-b[i-1];b[0]=n;
for(int i=0;i<=T;++i) b[i]=mul(b[i],ifac[i]);
}*/
void solve(int l,int r,int dep)
{
if(l==r){f[dep][0]=coe[l];f[dep][1]=mod-1;return;}
int mid=(l+r)>>1;
solve(l,mid,dep);solve(mid+1,r,dep+1);
polymult(f[dep],f[dep+1],f[dep],r-l+2,0);
fill(f[dep+1],f[dep+1]+Poly::m,0);
}
int Q[N],R[N],P[N];
void GETK(int *a,int *b,int n)
{
memset(f,0,sizeof(f));memset(coe,0,sizeof(coe));
memset(P,0,sizeof(P));memset(Q,0,sizeof(Q));memset(R,0,sizeof(R));
for(int i=1;i<=n;++i) coe[i]=getinv(a[i]);
solve(1,n,0);
for(int i=0;i<=n;++i) Q[i]=f[0][i];
if(T+1<=n) fill(Q+T+1,Q+n+1,0);//wrong because no fill
polyinv(Q,R,T+1);polyder(Q,T+1);
polymult(Q,R,P,T+1,T+1);
for(int i=T;i;--i) b[i]=upm(-P[i-1]); b[0]=n;
for(int i=0;i<=T;++i) b[i]=mul(b[i],ifac[i]);
}
void solution()
{
initmath();
n=read();m=read();
for(int i=1;i<=n;++i) a[i]=read();
for(int i=1;i<=m;++i) b[i]=read();
T=read();
GETK(a,A,n);GETK(b,B,m);
polymult(A,B,C,T,T);
for(int i=1,iv=mul(inv[n],inv[m]);i<=T;++i) writeln(mul(iv,mul(fac[i],C[i])));
}
}
int main()
{
#ifdef Durant_Lee
freopen("LGP4705.in","r",stdin);
freopen("LGP4705.out","w",stdout);
#endif
DreamLolita::solution();
return 0;
}