「Luogu P4705」玩游戏

Luogu P4705 玩游戏

题意

有长度为 \(n\) 的序列 \(a\) 和长度为 \(m\) 的序列 \(b\)

对于 \(k=1,2,\dotsc,t\) 求随机两个元素 \(a_i\) 和 \(b_j\)，\((a_i+b_j)^k\) 的期望

模 \(998244353\)

\(n,m,t\le 10^5\)

做法

有点神

首先二项式展开

\[ \begin{align} & \sum_{i=1}^n \sum_{j=1}^m (a_i+b_j)^k \\ = & \sum_{i=1}^n \sum_{j=1}^m \sum_{x=0}^k a_i^xb_j^{k-x} \binom{k}{x} \\ = & k!\sum_{i=1}^n \sum_{j=1}^m \sum_{x=0}^k \frac{a_i^x}{x!} \times \frac{b_j^{k-x}}{(k-x)!} \\ = & k!\sum_{x=0}^k \left(\sum_{i=1}^n \frac{a_i^x}{x!} \right) \left(\sum_{i=1}^m \frac{b_i^{k-x}}{(k-x)!}\right) \end{align} \]

可以发现右边是卷积的形式，我们只要对于每一个 \(x\) 求出 \(\sum\limits_{i=1}^n a_i^x\) 和 \(\sum\limits_{i=1}^m b_i^x\) 即可

构造多项式 \(f(x)=\prod\limits_{i=1}^n (a_ix+1)\)，这可以用分治 NTT 在 \(\mathcal O(n\log^2n)\) 的复杂度内求出来

由于 \(\ln (ax+1) = -\sum\limits_{i=1}^\infty \frac{(-1)^i}{i}a^ix^i\)

证明如下

首先有

\[ \frac{1}{1-x}=\sum_{i=0}^\infty x^i \]

所以

\[ \begin{align} (\ln(1-x))' &= -\frac{1}{1-x} = -\sum_{i=0}^\infty x^i \\ \ln(1-x) &= -\sum_{i=1}^\infty \frac{x^i}{i} \end{align} \]

用 \(-ax\) 代入即可

所以 \(\ln(f(x)) = -\sum\limits_{j=1}^n \sum\limits_{i=1}^\infty \frac{(-1)^i}{i}a_j^ix^i\)

多余的系数去掉可以得到我们需要的每个 \(\sum\limits_{i=1}^n a_i^x\)

求这个的另一种方法是从 牛顿恒等式 推，本质似乎是一样的

总复杂度 \(\mathcal O(n\log^2n)\)

---

代码

#include<cstdio>
#include<algorithm>
#include<cctype>
#include<string.h>
#include<cmath>
#include<vector>

using namespace std;
#define ull unsigned long long

inline char read() {
	static const int IN_LEN = 1000000;
	static char buf[IN_LEN], *s, *t;
	return (s==t?t=(s=buf)+fread(buf,1,IN_LEN,stdin),(s==t?-1:*s++):*s++);
}
template<class T>
inline void read(T &x) {
	static bool iosig;
	static char c;
	for (iosig=false, c=read(); !isdigit(c); c=read()) {
		if (c == '-') iosig=true;
		if (c == -1) return;
	}
	for (x=0; isdigit(c); c=read()) x=((x+(x<<2))<<1)+(c^'0');
	if (iosig) x=-x;
}
const int OUT_LEN = 10000000;
char obuf[OUT_LEN], *ooh=obuf;
inline void print(char c) {
	if (ooh==obuf+OUT_LEN) fwrite(obuf, 1, OUT_LEN, stdout), ooh=obuf;
	*ooh++=c;
}
template<class T>
inline void print(T x) {
	static int buf[30], cnt;
	if (x==0) print('0');
	else {
		if (x<0) print('-'), x=-x;
		for (cnt=0; x; x/=10) buf[++cnt]=x%10+48;
		while(cnt) print((char)buf[cnt--]);
	}
}
inline void flush() { fwrite(obuf, 1, ooh - obuf, stdout); }


const unsigned P = 998244353;
struct Z{
	unsigned x;
	Z(const unsigned _x=0):x(_x){}
	inline Z operator +(const Z &rhs)const{ return x+rhs.x<P?x+rhs.x:x+rhs.x-P;}
	inline Z operator -(const Z &rhs)const{ return x<rhs.x?x-rhs.x+P:x-rhs.x;}
	inline Z operator -()const{ return x?P-x:0;}
	inline Z operator *(const Z &rhs)const{ return static_cast<ull>(x)*rhs.x%P;}
	inline Z operator +=(const Z &rhs){ return x=x+rhs.x<P?x+rhs.x:x+rhs.x-P, *this;}
	inline Z operator -=(const Z &rhs){ return x=x<rhs.x?x-rhs.x+P:x-rhs.x, *this;}
	inline Z operator *=(const Z &rhs){ return x=static_cast<ull>(x)*rhs.x%P, *this;}
};
int n, m, t;
vector<Z> a, b;

namespace Poly{
	const int MAX_LEN = 1<<18;

	Z w[MAX_LEN], Inv[MAX_LEN];// for DFT

	inline Z Pow(Z x, int y=P-2){
		Z ans=1;
		for(; y; y>>=1, x=x*x) if(y&1) ans=ans*x;
		return ans;
	}
	inline void Init(){
		for(int i=1; i<MAX_LEN; i<<=1){
			w[i]=1;
			Z t=Pow(3, (P-1)/i/2);
			for(int j=1; j<i; ++j) w[i+j]=w[i+j-1]*t;
		}
		Inv[1]=1;
		for(int i=2; i<MAX_LEN; ++i) Inv[i]=Inv[P%i]*(P-P/i);
	}
	inline int Get(int x){ int n=1; while(n<=x) n<<=1; return n;}
	inline void DFT(vector<Z> &f, int n){
		static ull F[MAX_LEN];
		if((int)f.size()!=n) f.resize(n);
		for(int i=0, j=0; i<n; ++i){
			F[i]=f[j].x;
			for(int k=n>>1; (j^=k)<k; k>>=1);
		}
		for(int i=1; i<n; i<<=1) for(int j=0; j<n; j+=i<<1){
			Z *W=w+i;
			ull *F0=F+j, *F1=F+j+i;
			for(int k=j; k<j+i; ++k, ++W, ++F0, ++F1){
				ull t=(*F1)*(W->x)%P;
				(*F1)=*F0+P-t, (*F0)+=t;
			}
		}
		for(int i=0; i<n; ++i) f[i]=F[i]%P;
	}
	inline void IDFT(vector<Z> &f, int n){
		f.resize(n), reverse(f.begin()+1, f.end());
		DFT(f, n);
		Z I=Pow(n);
		for(int i=0; i<n; ++i) f[i]=f[i]*I;
	}
	inline vector<Z> Add(const vector<Z> &f, const vector<Z> &g){
		vector<Z> ans=f;
		for(unsigned i=0; i<f.size(); ++i) ans[i]+=g[i];
		return ans;
	}
	inline vector<Z> Mul(const vector<Z> &f, const vector<Z> &g){
		if(f.size()*g.size()<=1000){
			vector<Z> ans;
			ans.resize(f.size()+g.size()-1);
			for(unsigned i=0; i<f.size(); ++i) for(unsigned j=0; j<g.size(); ++j)
				ans[i+j]+=f[i]*g[j];
			return ans;
		}
		static vector<Z> F, G;
		F=f, G=g;
		int p=Get(f.size()+g.size()-2);
		DFT(F, p), DFT(G, p);
		for(int i=0; i<p; ++i) F[i]*=G[i];
		IDFT(F, p);
		return F.resize(f.size()+g.size()-1), F;
	}
	vector<Z> &PolyInv(const vector<Z> &f, int n=-1){
		if(n==-1) n=f.size();
		if(n==1){
			static vector<Z> ans;
			return ans.clear(), ans.push_back(Pow(f[0])), ans;
		}
		vector<Z> &ans=PolyInv(f, (n+1)/2);
		vector<Z> tmp(&f[0], &f[0]+n);
		int p=Get(n*2-2);
		DFT(tmp, p), DFT(ans, p);
		for(int i=0; i<p; ++i) ans[i]=((Z)2-ans[i]*tmp[i])*ans[i];
		IDFT(ans, p);
		return ans.resize(n), ans;
	}
	inline vector<Z> Derivative(const vector<Z> &a){
		vector<Z> ans(a.size()-1);
		for(unsigned i=1; i<a.size(); ++i) ans[i-1]=a[i]*i;
		return ans;
	}
	inline vector<Z> Integral(const vector<Z> &a){
		vector<Z> ans(a.size()+1);
		for(unsigned i=0; i<a.size(); ++i) ans[i+1]=a[i]*Inv[i+1];
		return ans;
	}
	inline vector<Z> PolyLn(const vector<Z> &f){
		vector<Z> ans=Mul(Derivative(f), PolyInv(f));
		ans.resize(f.size()-1);
		return Integral(ans);
	}
	vector<Z> divide(int l, int r, const vector<Z> &f){
		if(l==r) return vector<Z>{1, f[l]};
		int mid=(l+r)>>1;
		return Mul(divide(l, mid, f), divide(mid+1, r, f));
	}
	inline vector<Z> solve(const vector<Z> &f, int t){
		vector<Z> ans=divide(0, f.size()-1, f);
		ans.resize(t+1), ans=PolyLn(ans);
		for(int i=1; i<=t; ++i) ans[i]*=(i&1?i:P-i);
		return ans[0]=f.size(), ans;
	}
}
int main() {
	Poly::Init();
	read(n), read(m), a.resize(n), b.resize(m);
	for(int i=0; i<n; ++i) read(a[i].x);
	for(int i=0; i<m; ++i) read(b[i].x);
	read(t);
	a=Poly::solve(a, t);
	b=Poly::solve(b, t);
	for(int i=1, k=1; i<=t; k=(ull)k*Poly::Inv[++i].x%P) a[i]=a[i]*k;
	for(int i=1, k=1; i<=t; k=(ull)k*Poly::Inv[++i].x%P) b[i]=b[i]*k;
	a=Poly::Mul(a, b), a.resize(t+1);
	for(int i=1, k=Poly::Pow((ull)n*m%P).x; i<=t; k=(ull)k*++i%P) print((a[i]*k).x), print('\n');
	return flush(), 0;
}