题面在这里
description
对于(k=1,2,...,t)，求$$frac{1}{nm}sum_{i=1}{n}sum_{j=1}{m}(a_i+b_j)^k$$
对(998244353)取模。
data range
[1le n,m,kle 10^5,0le a_i,b_ile 998244352
]
solution
由于要求多项式
[egin{aligned}
Ans(x)&=sum_{k=1}^{t}sum_{i=1}^{n}sum_{j=1}^{m}(a_i+b_j)^kx^k\
&=sum_{k=1}^{t}sum_{i=1}^{n}sum_{j=1}^{m}sum_{s=1}^{k}inom{k}{s}a_i^sb_j^{k-s}x^k\
&=sum_{k=1}^{t}k!sum_{i=1}^{n}sum_{j=1}^{m}sum_{s=1}^{k}frac{a_i^s}{s!}frac{b_j^{k-s}}{(k-s)!}x^k\
&=sum_{k=1}^{t}k!sum_{s=1}^{k}(sum_{i=1}^{n}frac{a_i^s}{s!})(sum_{j=1}^{m}frac{b_j^{k-s}}{(k-s)!})x^k\
&=sum_{k=1}^{t}k!sum_{s=1}^{k}frac{sum_{i=1}^{n}a_i^s}{s!}frac{sum_{j=1}^{m}b_j^{k-s}}{(k-s)!}x^k\
end{aligned}]
我们记
[f_k(x)=sum_{s=1}^{k}frac{sum_{i=1}^{n}a_i^s}{s!}frac{sum_{j=1}^{m}b_j^{k-s}}{(k-s)!}x^s
]
则$$Ans(x)=sum_{k=1}{t}f_k(1)xk$$
即其系数的前缀和
于是现在我们要求出$$g_a(x)=sum_{i=1}{t}sum_{j=1}{n}a_jixi$$
构造函数(h(x)=prod_{i=1}^{n}(a_ix+1))，因为
[ln[h(x)]=sum_{i=1}^{n}ln(a_ix+1)
]
而
[ln'(a_ix+1)=frac{a_i}{a_ix+1}=sum_{j=1}^{infty}(-1)^ja_i^{j+1}x^j
]
对上面这个式子求积分，我们有
[ln(a_ix+1)=sum_{j=1}^{infty}(-1)^{j-1}frac{a_i^j}{j}x^j
]
于是我们有
[egin{aligned}
ln[h(x)]&=sum_{i=1}^{n}sum_{j=1}^{infty}(-1)^{j-1}frac{a_i^j}{j}x^j\
&=sum_{j=1}^{infty}(-1)^{j-1}frac{sum_{i=1}^{n}a_i^j}{j}x^j\
end{aligned}]
(h(x))我们可以使用分治(FFT)在(O(nlog^2n))的时间内求出；
(ln[h(x)])可以使用多项式求(ln)在(O(nlogn))的时间内求出，
系数稍加处理即可得到(g_a(x))和(g_b(x))；
最后将(g_a(x))和(g_b(x))做一遍(NTT)即可得到(f(x))。
总时间复杂度为(O(nlog^2n))
常数巨大
调试(3+day)后终于过了第二个样例
code
#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cmath>
#include<cstring>
#define RG register
#define il inline
using namespace std;
typedef long long ll;
typedef double dd;
const int N=400010;
const int mod=998244353;
const dd pi=acos(-1);
il int read(){
RG int d=0,w=0;char ch=getchar();
while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
if(ch=='-')w=1,ch=getchar();
while(ch>='0'&&ch<='9')d=(d<<3)+(d<<1)+(ch^48),ch=getchar();
return w?-d:d;
}
il int poww(int a,int b){
RG int ret=1;
for(;b;b>>=1,a=1ll*a*a%mod)
if(b&1)ret=1ll*ret*a%mod;
return ret;
}
int r[N];
il void NTT(int *a,int n,int opt){
for(RG int i=0;i<n;i++)if(i<r[i])swap(a[i],a[r[i]]);
for(RG int i=2;i<=n;i<<=1){
RG int wn=poww((opt==1)?3:((mod+1)/3),(mod-1)/i);
for(RG int j=0;j<n;j+=i){
RG int w=1;
for(RG int k=j;k<j+(i>>1);k++,w=1ll*w*wn%mod){
RG int x=1ll*a[k+(i>>1)]*w%mod;
a[k+(i>>1)]=(a[k]-x+mod)%mod;
a[k]=(a[k]+x)%mod;
}
}
}
if(opt==-1)
for(RG int i=0,rev=poww(n,mod-2);i<n;i++)
a[i]=1ll*a[i]*rev%mod;
}
int inv[N];
il void initinv(int n){
inv[0]=inv[1]=1;
for(RG int i=2;i<n;i++)inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
}
int x[N],y[N];
#define M ((L+R)>>1)
void solve(int *a,int L,int R){//分治FFT
if(L==R)return;solve(a,L,M);solve(a,M+1,R);
RG int n=M-L+1,m=R-M,len=1,l=0;
for(;len<=(n+m);len<<=1,l++);
for(RG int i=0;i<len;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
for(RG int i=0;i<len;i++)x[i]=y[i]=0;x[0]=y[0]=1;
for(RG int i=1;i<=n;i++)x[i]=a[L+i-1];
for(RG int i=1;i<=m;i++)y[i]=a[M+i];
NTT(x,len,1);NTT(y,len,1);
for(RG int i=0;i<len;i++)x[i]=1ll*x[i]*y[i]%mod;
NTT(x,len,-1);
for(RG int i=1;i<=n+m;i++)a[L+i-1]=x[i];
}
il void getdao(int *a,int *x,int n){//多项式求导
for(RG int i=0;i<n;i++)x[i]=1ll*a[i+1]*(i+1)%mod;x[n-1]=0;
}
il void getjifen(int *a,int *x,int n){//多项式求积分
for(RG int i=n-1;i;i--)x[i]=1ll*a[i-1]*inv[i]%mod;x[0]=0;
}
int xi[N],yi[N];
void getinv(int *f,int *g,int n,int l){//多项式求逆
if(n==1){g[0]=poww(f[0],mod-2);return;}getinv(f,g,n>>1,l-1);
for(RG int i=0;i<(n<<1);i++)r[i]=(r[i>>1]>>1)|((i&1)<<l);
for(RG int i=0;i<(n<<1);i++)xi[i]=yi[i]=0;
for(RG int i=0;i<n;i++)xi[i]=f[i];
for(RG int i=0;i<(n>>1);i++)yi[i]=g[i];
NTT(xi,n<<1,1);NTT(yi,n<<1,1);
for(RG int i=0;i<(n<<1);i++)
xi[i]=1ll*(2-1ll*xi[i]*yi[i]%mod+mod)%mod*yi[i]%mod;
NTT(xi,n<<1,-1);
for(RG int i=0;i<n;i++)g[i]=xi[i];
}
void getln(int *a,int *f,int n,int l){//多项式求ln
memset(x,0,sizeof(x));memset(y,0,sizeof(y));
getdao(a,x,n);getinv(a,y,n,l);
for(RG int i=0;i<(n<<1);i++)r[i]=(r[i>>1]>>1)|((i&1)<<l);
NTT(x,n<<1,1);NTT(y,n<<1,1);
for(RG int i=0;i<(n<<1);i++)x[i]=1ll*x[i]*y[i]%mod;
NTT(x,n<<1,-1);
getjifen(x,f,n);
}
int n,m,t,a[N],b[N],f[N],g[N],l,rv;
int main()
{
n=read();m=read();rv=1ll*poww(n,mod-2)*poww(m,mod-2)%mod;
for(RG int i=1;i<=n;i++)a[i]=read();
for(RG int i=1;i<=m;i++)b[i]=read();
solve(a,1,n);solve(b,1,m);//分治FFT得到h(x);
t=read();a[0]=b[0]=1;
for(l=0;(1<<l)<=t;l++);initinv(1<<l);
getln(a,f,(1<<l),l);getln(b,g,(1<<l),l);
for(RG int i=1;i<=t;i++){
f[i]=(i&1)?(1ll*f[i]*i%mod):((mod-1ll*f[i]*i%mod)%mod);
g[i]=(i&1)?(1ll*g[i]*i%mod):((mod-1ll*g[i]*i%mod)%mod);
}
//多项式求ln得到ga(x)和gb(x);
for(RG int i=1;i<=t;i++){
inv[i]=1ll*inv[i-1]*inv[i]%mod;
f[i]=1ll*f[i]*inv[i]%mod;
g[i]=1ll*g[i]*inv[i]%mod;
}
f[0]=n;g[0]=m;m=n=t;
for(m+=n,n=1,l=0;n<=m;l++,n<<=1);
for(RG int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
NTT(f,n,1);NTT(g,n,1);
for(RG int i=0;i<n;i++)f[i]=1ll*f[i]*g[i]%mod;
NTT(f,n,-1);
for(RG int i=1,fac=1;i<=t;i++,fac=1ll*fac*i%mod){
f[i]=1ll*f[i]*fac%mod;
printf("%lld
",1ll*f[i]*rv%mod);
}
return 0;
}