Processing math: 100%

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

题目描述


多项式求逆

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

已知

A(x)B(x)1(modxn)

要求

A(x)C(x)1(modx2n)

两式相减可得

A(x)(B(x)C(x))1(modxn)

B(x)C(x)1(modxn)

平方一下

B2(x)+C2(x)2B(x)C(x)0(modx2n),因为x0~x2n拆分后一定有一个在n以内的

两边同乘A(x)

A(x)B2(x)+C(x)2B(x)0(modx2n)

C(x)2B(x)A(x)B2(x)(modx2n)

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

分治FFT/NTT

一般形式(本题与此无关):

f(i+j)=f(i)g(j)

cdq即可,log^2

题解

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

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

那么可得

f[i]=Lj=0[xj]Pi(x)f[i+j]

解方程可得

f[i]=11[x0]Pi(x)Lj=1Pi(x)f[i+j]

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

f[i]=11[x0]Pi(x)Lj=1Pi(x)f[ij]

接下来开始操♂作,定义

P[i]=pix+(1pi)(分清楚小括号和中括号),若i<1则P[i]=1

Pi(x)=iiL+1P[i]

Fi(x)=j<=ixjf(j)

Si(x)=j<=iP[j]

G[x][y]=xi=yL+1P[i]=Sx(x)/SyL(x)

(这里拓展了一下G的定义,为了方便之后的计算)

F[x][y]=Fx1(x)G[x][y]

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

F[1][n]=0G[1][n]=S1(x)/SnL+1(x),直接乘会爆掉所以这里也要分治NTT

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

左半边:

G[x][m]=G[x][y]yLi=mL+1P[i]=G[x][y]yi=m+1P[iL]

F[x][m]=Fx1(x)G[x][m]=Fx1(x)G[x][y]yi=m+1P[iL]

=F[x][y]yi=m+1P[iL]

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

右半边:

G[m+1][y]=G[x][y]m+1i=x+1P[i]=G[x][y]mi=xP[i+1]

F[m+1][y]=Fm(x)G[m+1][y]=Fm(x)G[x][y]mi=xP[i+1]

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

=Fm(x)G[x][y]mi=xP[i+1]Fx1(x)G[x][y]mi=xP[i+1]+

=(Fm(x)Fx1(x))G[x][y]mi=xP[i+1]+F[x][y]mi=xP[i+1]

=F[x][y]mi=xP[i+1]+(Fm(x)Fx1(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\n",(ans[n]+mod)%mod);
	
	fclose(stdin);
	fclose(stdout);
	
	return 0;
}
posted @   gmh77  阅读(308)  评论(0编辑  收藏  举报
编辑推荐:
· 深入理解 Mybatis 分库分表执行原理
· 如何打造一个高并发系统?
· .NET Core GC压缩(compact_phase)底层原理浅谈
· 现代计算机视觉入门之:什么是图片特征编码
· .NET 9 new features-C#13新的锁类型和语义
阅读排行:
· Spring AI + Ollama 实现 deepseek-r1 的API服务和调用
· 《HelloGitHub》第 106 期
· 数据库服务器 SQL Server 版本升级公告
· 深入理解Mybatis分库分表执行原理
· 使用 Dify + LLM 构建精确任务处理应用
点击右上角即可分享
微信分享提示