题意
给出两个序列(a,b),长度分别为(n,m),请求出(x=1,2,ldots,t),下式的值
[dfrac{sum_{i=1}^nsum_{j=1}^m(a_i+b_j)^x}{nm}
]
题解
看到还能发就来水一篇
本人数学比较菜有错误请联系这个屑
我们记(x)的答案为(ans_x),那么就有:
[ans_x imes nm=sum_{i=1}^nsum_{j=1}^m(a_i+b_j)^x
]
好像也做不了什么,那就把((a_i+b_j)^x)用二项式定理展开再说:
[ans_x imes nm=sum_{i=1}^nsum_{j=1}^msum_{k=0}^x binom{x}{k}a_i^{k}b_j^{x-k}
]
交换顺序,再把万恶的组合数展开。
[ans_x imes nm=sum_{k=0}^xsum_{i=1}^nsum_{j=1}^m binom{x}{k}a_i^{k}b_j^{x-k}=sum_{k=0}^xsum_{i=1}^nsum_{j=1}^mfrac{x!}{k!(x-k)!}a_i^{k}b_j^{x-k}
]
[frac{ans_x imes nm}{x!}=sum_{k=0}^xfrac{sum_{i=1}^na_i^k}{k!}frac{sum_{j=1}^mb_j^{x-k}}{(x-k)!}
]
如果我们记(A(k)=sum_{i=1}^na_i^k),(B(k)=sum_{i=1}^mb_i^k),那么就有:
[frac{ans_x imes nm}{x!}=sum_{k=0}^xfrac{A(k)}{k!}frac{B(x-k)}{(x-k)!}
]
这样卷积就很明显了。
于是剩下的问题就变成求(A(0),A(1),A(2),ldots,A(t)),考虑构造下面一个多项式:
[F(x)=prod_{i=1}^n(1+a_ix)
]
这个可以分治得到,复杂度显然正确
两边同取(ln)得到:
[ln F(x)=sum_{i=1}^nln(1+a_ix)
]
两边求导:
[(ln F(x))^prime=sum_{i=1}^nfrac{a_i}{1+a_ix}
]
又知,(frac{1}{1+a_ix}=frac{1}{1-(-a_ix)}),即为一个首项为(1),公比为(-a_ix)的数列的和,那么(frac{1}{1+a_ix}=sum_{j=0}^{+infty}(-a_ix)^j)
带带带回去:
[(ln F(x))^prime=sum_{i=1}^na_isum_{j=0}^{+infty}(-a_ix)^j
]
再积回去:
[ln F(x)=int sum_{i=1}^na_isum_{j=0}^{+infty}(-a_ix)^jdx=sum_{i=1}^na_isum_{j=1}^{+infty}frac{(-a_i)^{j-1}x^j}{j}
]
[=sum_{i=1}^nsum_{j=1}^{+infty}frac{(-1)^{j-1}a_i^{j}x^j}{j}=sum_{j=1}^{+infty}frac{(-1)^{j-1}x_j}{j}sum_{i=1}^na_i^j=sum_{j=1}^{+infty}frac{(-1)^{j-1}A(j)}{j}x^j
]
于是此时求出(A)也不难了,随便乘一下就好了。
(B)也同理。把(A,B)带回去,不难得到答案。
代码
#include<bits/stdc++.h>
namespace in{
#ifdef slow
inline int getc(){return getchar();}
#else
char buf[1<<21],*p1=buf,*p2=buf;
inline int getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
#endif
template <typename T>inline void read(T& t){
t=0;int f=0;char ch=getc();while (!isdigit(ch)){if(ch=='-')f = 1;ch=getc();}
while(isdigit(ch)){t=t*10+(ch-48);ch = getc();}if(f)t=-t;
}
template <typename T,typename... Args> inline void read(T& t, Args&... args){read(t);read(args...);}
}
namespace out{
char buffer[1<<21];int p1=-1;const int p2 = (1<<21)-1;
inline void flush(){fwrite(buffer,1,p1+1,stdout),p1=-1;}
inline void putc(const char &x) {if(p1==p2)flush();buffer[++p1]=x;}
template <typename T>void write(T x) {
static char buf[15];static int len=-1;if(x>=0){do{buf[++len]=x%10+48,x/=10;}while (x);}else{putc('-');do {buf[++len]=-(x%10)+48,x/=10;}while(x);}
while (len>=0)putc(buf[len]),--len;
}
}
using namespace std;
template<const int mod>
struct modint{
int x;
modint<mod>(int o=0){x=o;}
modint<mod> &operator = (int o){return x=o,*this;}
modint<mod> &operator +=(modint<mod> o){return x=x+o.x>=mod?x+o.x-mod:x+o.x,*this;}
modint<mod> &operator -=(modint<mod> o){return x=x-o.x<0?x-o.x+mod:x-o.x,*this;}
modint<mod> &operator *=(modint<mod> o){return x=1ll*x*o.x%mod,*this;}
modint<mod> &operator ^=(int b){
modint<mod> a=*this,c=1;
for(;b;b>>=1,a*=a)if(b&1)c*=a;
return x=c.x,*this;
}
modint<mod> &operator /=(modint<mod> o){return *this *=o^=mod-2;}
modint<mod> &operator +=(int o){return x=x+o>=mod?x+o-mod:x+o,*this;}
modint<mod> &operator -=(int o){return x=x-o<0?x-o+mod:x-o,*this;}
modint<mod> &operator *=(int o){return x=1ll*x*o%mod,*this;}
modint<mod> &operator /=(int o){return *this *= ((modint<mod>(o))^=mod-2);}
template<class I>friend modint<mod> operator +(modint<mod> a,I b){return a+=b;}
template<class I>friend modint<mod> operator -(modint<mod> a,I b){return a-=b;}
template<class I>friend modint<mod> operator *(modint<mod> a,I b){return a*=b;}
template<class I>friend modint<mod> operator /(modint<mod> a,I b){return a/=b;}
friend modint<mod> operator ^(modint<mod> a,int b){return a^=b;}
friend bool operator ==(modint<mod> a,int b){return a.x==b;}
friend bool operator !=(modint<mod> a,int b){return a.x!=b;}
bool operator ! () {return !x;}
modint<mod> operator - () {return x?mod-x:0;}
modint<mod> &operator++(int){return *this+=1;}
};
const int N=4e6+5;
const int mod=998244353;
const modint<mod> GG=3,Ginv=modint<mod>(1)/3,I=86583718;
struct poly{
vector<modint<mod>>a;
modint<mod>&operator[](int i){return a[i];}
int size(){return a.size();}
void resize(int n){a.resize(n);}
};
int rev[N];
inline int ext(int n){int k=0;while((1<<k)<n)k++;return k;}
inline void init(int k){int n=1<<k;for(int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));}
inline void ntt(poly&a,int k,int typ){
int n=1<<k;
for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int mid=1;mid<n;mid<<=1){
modint<mod> wn=(typ>0?GG:Ginv)^((mod-1)/(mid<<1));
for(int r=mid<<1,j=0;j<n;j+=r){
modint<mod> w=1;
for(int k=0;k<mid;k++,w=w*wn){
modint<mod> x=a[j+k],y=w*a[j+k+mid];
a[j+k]=x+y,a[j+k+mid]=x-y;
}
}
}
if(typ<0){
modint<mod> inv=modint<mod>(1)/n;
for(int i=0;i<n;i++)a[i]*=inv;
}
}
inline poly one(){poly a;a.a.push_back(1);return a;}
poly operator +(poly a,poly b){
int n=max(a.size(),b.size());a.resize(n),b.resize(n);
for(int i=0;i<n;i++)a[i]+=b[i];return a;
}
poly operator -(poly a,poly b){
int n=max(a.size(),b.size());a.resize(n),b.resize(n);
for(int i=0;i<n;i++)a[i]-=b[i];return a;
}
inline poly operator*(poly a,poly b){
int n=a.size()+b.size()-1,k=ext(n);
a.resize(1<<k),b.resize(1<<k),init(k);
ntt(a,k,1);ntt(b,k,1);for(int i=0;i<(1<<k);i++)a[i]*=b[i];
ntt(a,k,-1),a.resize(n);return a;
}
inline poly operator*(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]*=b;return a; }
inline poly operator/(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]/=b;return a; }
inline poly operator-(poly a){for(int i=0;i<a.size();i++)a[i]=-a[i];return a; }
poly inv(poly F,int k){
int n=1<<k;F.resize(n);
if(n==1){F[0]=modint<mod>(1)/F[0];return F;}
poly G,H=inv(F,k-1);
G.resize(n),H.resize(n<<1),F.resize(n<<1);
for(int i=0;i<n/2;i++)G[i]=H[i]*2;
init(k+1),ntt(H,k+1,1),ntt(F,k+1,1);
for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i]*F[i];
ntt(H,k+1,-1),H.resize(n);
for(int i=0;i<n;i++)G[i]-=H[i];return G;
}
inline poly inv(poly a){
int n=a.size();
a=inv(a,ext(n)),a.resize(n);return a;;
}
inline poly deriv(poly a){//求导
int n=a.size()-1;
for(int i=0;i<n;i++)a[i]=a[i+1]*(i+1);
a.resize(n);return a;
}
inline poly inter(poly a){//求原
int n=a.size()+1;a.resize(n);
for(int i=n;i>=1;i--)a[i]=a[i-1]/i;
a[0]=0;return a;
}
inline poly ln(poly a){
int n=a.size();
a=inter(deriv(a)*inv(a));
a.resize(n);return a;
}
#define mid (l+r>>1)
inline poly F(int l,int r,int *a){
poly G;
if(l==r){G.resize(2);G[0]=1;G[1]=a[l];}
else{
G=F(l,mid,a)*F(mid+1,r,a);
}
//printf("[%d,%d] ",l,r);for(int i=0;i<G.size();i++)printf("%d ",G[i].x);puts("");
return G;
}
int n,m,t,a[N],b[N];
poly A,B;
modint<mod>fac[N];
signed main(){
fac[0]=1;for(int i=1;i<N;i++)fac[i]=fac[i-1]*i;
in::read(n,m);
for(int i=1;i<=n;i++)in::read(a[i]);
for(int i=1;i<=m;i++)in::read(b[i]);
in::read(t);
A=F(1,n,a);A.resize(t+1);A=ln(A);
A[0]=n;for(int i=1;i<=t;i++)if(i&1)A[i]*=i;else A[i]*=i,A[i]=-A[i];
B=F(1,m,b);B.resize(t+1);B=ln(B);
B[0]=m;for(int i=1;i<=t;i++)if(i&1)B[i]*=i;else B[i]*=i,B[i]=-B[i];
//for(int i=0;i<=t;i++)cout<<A[i].x<<" ";cout<<endl;
//for(int i=0;i<=t;i++)cout<<B[i].x<<" ";cout<<endl;
for(int i=0;i<=t;i++)A[i]/=fac[i];
for(int i=0;i<=t;i++)B[i]/=fac[i];
poly ans=A*B;for(int i=0;i<=t;i++)ans[i]=(ans[i]*fac[i]/n)/m;
for(int i=1;i<=t;i++)cout<<ans[i].x<<endl;
}