6508. 【GDOI2020模拟03.11】我的朋友们（多项式求逆、分治NTT）

题目描述

多项式求逆

第一次写就是这么毒瘤的题目

已知

(A(x)B(x)equiv 1(mod;x^n))

要求

(A(x)C(x)equiv 1(mod;x^{2n}))

两式相减可得

(A(x)(B(x)-C(x))equiv 1(mod;x^{n}))

(B(x)-C(x)equiv 1(mod;x^{n}))

平方一下

(B^2(x)+C^2(x)-2B(x)C(x)equiv 0(mod;x^{2n}))，因为x^0~x^2n拆分后一定有一个在n以内的

两边同乘A(x)

(A(x)B^2(x)+C(x)-2B(x)equiv 0(mod;x^{2n}))

(C(x)equiv 2B(x)-A(x)B^2(x)(mod;x^{2n}))

虽然做了log n次乘法，但n的总和还是n级别的，所以还是(O(n log n))

分治FFT/NTT

一般形式（本题与此无关）：

(f(i+j)=sum{f(i)*g(j)})

cdq即可，log^2

题解

期望的一般算法：计算以某个状态为起点到终点的期望步数

设f[i]表示从i开始的答案，Pi(x)表示在滑块[i,i+L-1]中的期望多项式，其中x^k表示有k个不合法

那么可得

(f[i]=sum_{j=0}^{L}{[x^j]P_i(x)f[i+j]})

解方程可得

(f[i]=frac{1}{1-[x^0]P_i(x)}sum_{j=1}^{L}{P_i(x)f[i+j]})

不好算，所以把数组p翻转，变成（Pi(x)变成[i-L+1,i]的期望）

(f[i]=frac{1}{1-[x^0]P_i(x)}sum_{j=1}^{L}{P_i(x)f[i-j]})

接下来开始操♂作，定义

(P[i]=p_ix+(1-p_i))（分清楚小括号和中括号），若i<1则P[i]=1

则(P_i(x)=prod_{i-L+1}^{i}{P[i]})

(F_i(x)=sum_{j<=i}{x^j*f(j)})

(S_i(x)=prod_{j<=i}{P[j]})

(G[x][y]=prod_{i=y-L+1}^{x}{P[i]}=S_x(x)/S_{y-L}(x))

（这里拓展了一下G的定义，为了方便之后的计算）

(F[x][y]=F_{x-1}(x)*G[x][y])

如果能求出每个F[x][x]，那么([x^x]F[x][x])的系数就是(f[x])，字母不够用意会一下

(F[1][n]=0)，(G[1][n]=S_1(x)/S_{n-L+1}(x))，直接乘会爆掉所以这里也要分治NTT

向下分治，设m=(x+y)/2

左半边：

(G[x][m]=G[x][y]*prod_{i=m-L+1}^{y-L}{P[i]}=G[x][y]*prod_{i=m+1}^{y}{P[i-L]})

(F[x][m]=F_{x-1}(x)*G[x][m]=F_{x-1}(x)*G[x][y]*prod_{i=m+1}^{y}{P[i-L]})

(=F[x][y]*prod_{i=m+1}^{y}{P[i-L]})

有可能i-L<1，但是可以发现这样f必然为0，所以可以把P设为1

右半边：

(G[m+1][y]=G[x][y]*prod_{i=x+1}^{m+1}{P[i]}=G[x][y]*prod_{i=x}^{m}{P[i+1]})

(F[m+1][y]=F_{m}(x)*G[m+1][y]=F_{m}(x)*G[x][y]*prod_{i=x}^{m}{P[i+1]})

为了方便接下来的操作，需要把Fm(x)的系数限制到x~y范围内

(=F_{m}(x)*G[x][y]*prod_{i=x}^{m}{P[i+1]}-F_{x-1}(x)*G[x][y]*prod_{i=x}^{m}{P[i+1]}+前面那一坨)

(=(F_{m}(x)-F_{x-1}(x))*G[x][y]*prod_{i=x}^{m}{P[i+1]}+F[x][y]*prod_{i=x}^{m}{P[i+1]})

(=F[x][y]*prod_{i=x}^{m}{P[i+1]}+(F_{m}(x)-F_{x-1}(x))*G[m+1][y])

看似时间复杂度不变，但实际上F只需要保留以x^y结尾的2(y-x+1)项，G只用保留前(y-x+1)项即可，剩余的不会对新的区间产生影响

因为最终只需要x^x

还有一个要注意的，乘到F里的G[m+1][y]要保留(y-x+1)项

关键在于怎么实现，只可意会不可言传

一个比较方便的写法是写一个 把A(x)的一段乘上B(x)的一段后把某一段（删掉开头若干项）丢到C(x)上 的函数，要仔细考虑具体删掉多少项

极其难写+卡时间，写之前做好心理准备

code

我相信比expln之类的要难写的多（虽然还没学）

#include <bits/stdc++.h>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define max(a,b) (a>b?a:b)
#define min(a,b) (a<b?a:b)
#define mod 998244353
#define Mod 998244351
#define ll long long
#define G 3
#define file
using namespace std;

ll aa[262144],bb[262144],bb2[262144],bb3[262144],cc[262144],A[262144],p[262144],p1[4194304],p2[4194304]; //p1=i-l p2=i+1
ll f[8388608],g[4194304],a[8388608],b[8388608],c[262144],ans[262144],P[262144]; //f=F g=G ans=f
int bg[262144],ed[262144],bg1[262144],ed1[262144],bg2[262144],ed2[262144],A2[19][262144],len2[262145],n2[262145],n3[262145],W1[19],W2[19],n,L,I,J,i,j,k,l,s,s2,N,N2,len,tot,tot1,tot2;

void look(ll *a,int x,int y)
{
	int i;
	fd(i,y,x) cout<<(a[i]+mod)%mod<<" ";cout<<endl;
}

ll qpower(ll a,int b)
{
	ll ans=1;
	
	while (b)
	{
		if (b&1)
		ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	
	return ans;
}

void dft(ll *a,int type,int N,int len)
{
	int i,j,k,l,S=N,s1=2,s2=1;
	ll w,W,u,v;
	
	fo(i,0,N-1) A[A2[len][i]]=a[i];
	memcpy(a,A,N*8);
	
	fo(i,1,len)
	{
		S>>=1;
		w=(type==1)?W1[i]:W2[i];
		
		fo(j,0,S-1)
		{
			W=1;
			fo(k,0,s2-1)
			{
				u=a[j*s1+k];
				v=a[j*s1+k+s2]*W;
				
				a[j*s1+k]=(u+v)%mod;
				a[j*s1+k+s2]=(u-v)%mod;
				
				W=W*w%mod;
			}
		}
		
		s1<<=1,s2<<=1;
	}
}

void mul(ll *x,ll *y,int n,int m)
{
	int i,j,k,l,N,N2,len;
	
	len=len2[n+m];N=n2[n+m];N2=n3[n+m];
	
	if (n<m)
	fo(i,n,m) x[i]=0;
	else
	fo(i,m,n) y[i]=0;
	fo(i,max(n,m),N-1) x[i]=y[i]=0;
	memcpy(bb3,y,N*8);
	
	dft(x,1,N,len);
	dft(bb3,1,N,len);
	fo(i,0,N-1) x[i]=x[i]*bb3[i]%mod;
	dft(x,-1,N,len);
	fo(i,0,N-1) x[i]=x[i]*N2%mod;
}
void Mul(ll *x,ll *y,ll *z,int X1,int X2,int Y1,int Y2,int Z1,int Z2,int st)
{
	int i;
	
	fo(i,X1,X2) aa[i-X1]=x[i];
	fo(i,Y1,Y2) bb[i-Y1]=y[i];
	mul(aa,bb,X2-X1+1,Y2-Y1+1);
	fo(i,Z1,Z2) z[i]=aa[i-Z1+st];
}

void ny(ll *a,int n)
{
	int i,j,k,l;
	
	cc[0]=qpower(a[0],Mod);
	for (i=2; i<=n; i*=2)
	{
		memcpy(bb2,a,i*8);
		
		mul(bb2,cc,i,i);mul(bb2,cc,i,i);
		fo(j,0,i+i-1)
		bb2[j]=(2*cc[j]-bb2[j])%mod;
		
		memcpy(cc,bb2,i*8);
	}
	memcpy(a,cc,n*8);
}

void work(int t,int x,int y)
{
	int i,mid=(x+y)/2;
	
	if (x==y)
	{
		bg[t]=tot+1;ed[t]=tot+2;
		
		if (x>L)
		p1[tot+1]=1-p[x-L],p1[tot+2]=p[x-L];
		else
		p1[tot+1]=1,p1[tot+2]=0;
		if (x<n)
		p2[tot+1]=1-p[x+1],p2[tot+2]=p[x+1];
		else
		p2[tot+1]=1,p2[tot+2]=0;
		
		tot+=2;
		return;
	}
	
	work(t*2,x,mid);
	work(t*2+1,mid+1,y);
	
	bg[t]=tot+1;ed[t]=tot+(y-x+1)+1;
	Mul(p1,p1,p1,bg[t*2],ed[t*2],bg[t*2+1],ed[t*2+1],bg[t],ed[t],0);
	Mul(p2,p2,p2,bg[t*2],ed[t*2],bg[t*2+1],ed[t*2+1],bg[t],ed[t],0);
	tot+=(y-x+1)+1;
}

void Work(int t,int x,int y)
{
	int i,mid=(x+y)/2;
	
	if (x==y)
	{
		if (x>=L)
		ans[x]=(f[bg1[t]+1]+1)*qpower(1-P[x],Mod)%mod;
		return;
	}
	
	bg2[t*2]=tot2+1;ed2[t*2]=tot2+(mid-x+1);tot2+=mid-x+1;
	Mul(g,p1,g,bg2[t],ed2[t],bg[t*2+1],ed[t*2+1],bg2[t*2],ed2[t*2],0);
	
	bg1[t*2]=tot1+1;ed1[t*2]=tot1+min(n+1,2*(mid-x+1));tot1+=min(n+1,2*(mid-x+1));
	Mul(f,p1,f,bg1[t],ed1[t],bg[t*2+1],ed[t*2+1],bg1[t*2],ed1[t*2],max(0,2*x-mid-1)-max(0,2*x-y-1));
	Work(t*2,x,mid);
	
//	---
	
	bg2[t*2+1]=tot2+1;ed2[t*2+1]=tot2+(y-mid);tot2+=y-mid;
	Mul(g,p2,a,bg2[t],ed2[t],bg[t*2],ed[t*2],bg2[t],ed2[t],0); //Ïȱ£Áôy-x+1Ïî
	fo(i,0,y-mid-1) g[bg2[t*2+1]+i]=a[bg2[t]+i];
	
	bg1[t*2+1]=tot1+1;ed1[t*2+1]=tot1+min(n+1,2*(y-mid));tot1+=min(n+1,2*(y-mid));
	Mul(f,p2,b,bg1[t],ed1[t],bg[t*2],ed[t*2],bg1[t*2+1],ed1[t*2+1],max(0,2*mid-y+1)-max(0,2*x-y-1));
	fo(i,x,mid) c[i]=ans[i];
	Mul(c,a,a,x,mid,bg2[t],ed2[t],bg1[t*2+1],ed1[t*2+1],max(0,2*mid-y+1)-x);
	
	fo(i,bg1[t*2+1],ed1[t*2+1]) f[i]=(a[i]+b[i])%mod;
	Work(t*2+1,mid+1,y);
}

int main()
{
	freopen("friends.in","r",stdin);
	#ifdef file
	freopen("friends.out","w",stdout);
	#endif
	
	I=2;J=1;s=2;len2[1]=1;n2[1]=2;n3[1]=499122177;
	fo(len,0,18)
	{
		if (len<18)
		W1[len+1]=qpower(G,(mod-1)/s),W2[len+1]=qpower(W1[len+1],Mod);
		
		s2=qpower(s,Mod);
		fo(i,0,J-1)
		{
			if (I<=262144)
			len2[I]=len+1,n2[I]=s,n3[I]=s2,++I;
			
			j=i;k=0;
			fo(l,1,len)
			k=(k<<1)+(j&1),j>>=1;
			
			A2[len][i]=k;
		}
		J*=2;s=s*2;
	}
	
	scanf("%d%d",&n,&L);
	fd(i,n,1)
	{
		scanf("%d%d",&p[i],&j);
		p[i]=p[i]*qpower(j,Mod)%mod;
	}
	P[L]=1;
	fo(i,1,L) P[L]=P[L]*(1-p[i])%mod;
	fo(i,L+1,n) P[i]=P[i-1]*(1-p[i])%mod*qpower(1-p[i-L],Mod)%mod;
	
	work(1,1,n);
	
	tot1=n+1;bg1[1]=1;ed1[1]=n+1;
	tot2=n+1;bg2[1]=1;ed2[1]=n;
	a[1]=p[1],a[0]=1-p[1];
	fo(i,0,n)
	b[i]=p1[bg[1]+i];
	N=pow(2,ceil(log2(n+1)));
	ny(b,N);
	Mul(a,b,g,0,1,0,N-1,1,n,0);
	
	Work(1,1,n);
	
	printf("%lld
",(ans[n]+mod)%mod);
	
	fclose(stdin);
	fclose(stdout);
	
	return 0;
}