题面
题解
妈呀这辣鸡题目调了我整整三天……最后发现竟然是因为分治\(NTT\)之后的多项式长度不是\(2\)的幂导致把多项式的值存下来的时候发生了一些玄学错误……玄学到了我\(WA\)的点全都是\(WA\)在\(2\)的幂次行里……
看到这种题目二话不说先推倒
\[ \begin{aligned} [x^k]Ans &={1\over nm}\sum_{i=1}^n\sum_{j=1}^m\left(a_i+b_j\right)^k\\ &={1\over nm}\sum_{i=1}^n\sum_{j=1}^m\sum_{p=0}^k{k\choose p}{a_i}^p{b_j}^{k-p}\\ &={k!\over nm}\sum_{p=0}^k{\sum_{i=1}^n{a_i}^p\over p!}{\sum_{j=1}^m{b_j}^{k-p}\over (k-p)!}\\ \end{aligned} \]
然后这就被画成了一个卷积的形式
定义两个多项式\(A(x)=\sum_{i=0}^\infty x^i\sum_{j=1}^n{a_j}^i\),和\(B(x)=\sum_{i=0}^\infty x^i\sum_{j=1}^m{b_j}^i\),只要我们能求出这两个多项式的系数,然后一通乱搞之后就能求出\(Ans\)了
然后继续推倒
\[ \begin{aligned} A(x) &=\sum_{i=0}^\infty x^i\sum_{j=1}^n{a_j}^i\\ &=\sum_{j=1}^n\sum_{i=0}^\infty {a_j}^ix^i\\ &=\sum_{i=1}^n{1\over 1-a_ix}\\ &=\sum_{i=1}^n {a_i}^0+{a_i}^1x^1+{a_i}^2x^2+... \end{aligned} \]
所以……这玩意儿该咋算啊……
我们设
\[ \begin{aligned} G(x) &=\sum_{i=1}^n{-a_i\over 1-a_ix}\\ &=\sum_{i=1}^n-{a_i}^1-{a_i}^2x-{a_i}^3x^2-...\\ \end{aligned} \]
那么就有\(A(x)=-xG(x)+n\)
然而我还是不会算\(G\)啊……
那就继续推倒
\[ \begin{aligned} G(x) &=\sum_{i=1}^n{-a_i\over 1-a_ix}\\ &=\sum_{i=1}^n\ln'\left(1-a_ix\right)\\ &=\ln'\left(\prod_{i=1}^n (1-a_ix)\right) \end{aligned} \]
分治\(NTT\)就行啦
然后没有然后了
我错了多项式比计算几何难调多了
//minamoto#include#define R register#define fp(i,a,b) for(R int i=(a),I=(b)+1;i I;--i)#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)using namespace std;char buf[1<<21],*p1=buf,*p2=buf;inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}int read(){ R int res,f=1;R char ch; while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1); for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0'); return res*f;}char sr[1<<21],z[20];int C=-1,Z=0;inline void Ot(){fwrite(sr,1,C+1,stdout),C=-1;}void print(R int x){ if(C>1<<20)Ot();if(x<0)sr[++C]='-',x=-x; while(z[++Z]=x%10+48,x/=10); while(sr[++C]=z[Z],--Z);sr[++C]='\n';}const int N=(1<<18)+5,P=998244353,Gi=332748118;inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}inline int dec(R int x,R int y){return x-y<0?x-y+P:x-y;}inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}int ksm(R int x,R int y){ R int res=1; for(;y;y>>=1,x=mul(x,x))if(y&1)res=mul(res,x); return res;}vector r[21];int rt[2][N<<1],inv[N],fac[N],ifac[N],lim,d;inline void init(R int len){lim=1,d=0;while(lim <<=1,++d;}void Pre(){ fp(d,1,18){ r[d].resize(1< >1]>>1)|((i&1)<<(d-1)); } inv[0]=inv[1]=fac[0]=fac[1]=ifac[0]=ifac[1]=1; fp(i,2,262144){ fac[i]=mul(fac[i-1],i), inv[i]=mul(P-P/i,inv[P%i]), ifac[i]=mul(ifac[i-1],inv[i]); } for(R int t=(P-1)>>1,i=1,x,y;i<=262144;i<<=1,t>>=1){ x=ksm(3,t),y=ksm(Gi,t); rt[1][i]=rt[0][i]=1; fp(k,1,i-1){ rt[1][k+i]=mul(rt[1][k+i-1],x), rt[0][k+i]=mul(rt[0][k+i-1],y); } }}int rev[N];void NTT(int *A,int ty){ fp(i,0,lim-1)if(i >1),init(len<<1); static int A[N],B[N]; fp(i,0,len-1)A[i]=a[i],B[i]=b[i]; fp(i,len,lim-1)A[i]=B[i]=0; NTT(A,1),NTT(B,1); fp(i,0,lim-1)A[i]=mul(A[i],mul(B[i],B[i])); NTT(A,0); fp(i,0,len-1)b[i]=dec(add(b[i],b[i]),A[i]); fp(i,len,lim-1)b[i]=0;}void Ln(int *a,int *b,int len){ static int A[N],B[N]; fp(i,1,len-1)A[i-1]=mul(a[i],i);A[len-1]=0; Inv(a,B,len),init(len<<1);fp(i,len,lim-1)A[i]=B[i]=0; NTT(A,1),NTT(B,1); fp(i,0,lim-1)A[i]=mul(A[i],B[i]); NTT(A,0); fp(i,1,len-1)b[i]=mul(A[i-1],inv[i]);b[0]=0; fp(i,len,lim-1)b[i]=0;}int D[25][N];void solve(int *a,int d,int l,int r){ if(l==r)return D[d][0]=1,D[d][1]=P-a[l],void(); int mid=(l+r)>>1; solve(a,d,l,mid),solve(a,d+1,mid+1,r),init(r-l+1+1); static int A[N],B[N]; fp(i,0,mid-l+1)A[i]=D[d][i];fp(i,mid-l+2,lim-1)A[i]=0; fp(i,0,r-mid)B[i]=D[d+1][i];fp(i,r-mid+1,lim-1)B[i]=0; NTT(A,1),NTT(B,1); fp(i,0,lim-1)A[i]=mul(A[i],B[i]); NTT(A,0); fp(i,0,r-l+1)D[d][i]=A[i]; fp(i,r-l+2,lim-1)D[d][i]=0;}int a[N],b[N],ak[N],bk[N];int n,m,t;void calc(int *a,int *b,int n,int t){ static int A[N],B[N]; solve(a,1,1,n); init(t+1);int len=lim; fp(i,0,n)A[i]=D[1][i]; fp(i,n+1,len-1)A[i]=0; Ln(A,B,len); fp(i,1,len-1)B[i-1]=mul(B[i],i);B[len-1]=0; b[0]=n; fp(i,1,t)b[i]=mul(P-B[i-1],ifac[i]);}void Mul(int *a,int *b){ init(t<<1); NTT(a,1),NTT(b,1); fp(i,0,lim-1)a[i]=mul(a[i],b[i]); NTT(a,0); int invm=ksm(mul(n,m),P-2); fp(i,1,t)print(1ll*a[i]*fac[i]%P*invm%P);}int main(){// freopen("testdata.in","r",stdin); Pre(); n=read(),m=read(); fp(i,1,n)a[i]=read(); fp(i,1,m)b[i]=read(); t=read(); calc(a,ak,n,t),calc(b,bk,m,t); Mul(ak,bk); return Ot(),0;}