luogu P4239 【模板】多项式求逆(加强版)

18 篇文章 0 订阅

背景:

多项式全家桶 eating... \text{eating...} eating...
复制模板 10 min 10\text{min} 10min,开始卡常 1h \text{1h} 1h

题目传送门:

https://www.luogu.org/problemnew/show/P4239

题意:

给定一个多项式 F ( x ) F(x) F(x),求一个多项式 G ( x ) G(x) G(x)使得 F ( x ) ∗ G ( x ) ≡ 1 ( m o d    x n ) F(x)*G(x)≡1(\mod x^n) F(x)G(x)1(modxn)。最后对系数取模(任意模数)。

思路:

任意模数?
2 k + 1 ( k ∈ N + ) 2^k+1(k∈N_{+}) 2k+1(kN+)模数请看 luogu P4238 \text{luogu P4238} luogu P4238【模板】多项式求逆
那不就是倍增多项式求逆 + \text{+} +任意模数 NTT \text{NTT} NTT吗?

又是三个模数 NTT \text{NTT} NTT啊。

可是直接做中国剩余定理时可能有问题啊(有负数 − A 2 F -A^2F A2F出现)?
怎么解决?
我们当前求逆中需要计算的式子是 2 A − A 2 F = A ( 2 − A F ) 2A-A^2F=A(2-AF) 2AA2F=A(2AF)
那么先把 A ∗ F A*F AF用三个模数算出来然后中国剩余定理合并出模 1 0 9 + 7 10^9+7 109+7下的值,取负,然后常数项取 1 1 1,然后转换成模三个模数下的值再和 A A A卷一次,然后出合并成模 1 0 9 + 7 10^9+7 109+7的值即可。

理论时间复杂度: Θ ( n log ⁡ n ) \Theta(n\log n) Θ(nlogn),但常数大得飞起,无奈之下对于所用能用 int \text{int} int代替 long long \text{long long} long long的都改了过来,竟然卡过了,无奈。
大力用卡常技巧,相信奇迹。

代码:

#include<cstdio>
#include<cstring>
#include<algorithm>
#define R register
#define I inline
using namespace std;
typedef long long LL;
typedef long double LD;
const int Mod[3]={469762049,998244353,1004535809},MOD=1000000007,G=3;
LL M=(LL)Mod[0]*Mod[1];
	int a[300010],b[300010],f[300010],g[300010],ans[3][300010];
	int limit,n,l,r[300010];
	int mod,inv_G;
I int read()
{
	int x=0,f=1;
	char ch=getchar();
	for(;ch<'0'||ch>'9';ch=getchar())
		if(ch=='-') f=-1;
	for(;ch>='0'&&ch<='9';x=(x<<3)+(x<<1)+(ch^48),ch=getchar());
	return x*f;
}
LL dg(int x,int k,LL mod)
{
	if(!k) return 1;
	LL op=dg(x,(k>>1),mod);
	if(k&1) return op*op%mod*x%mod; else return op*op%mod;
}
I LL times(LL a,LL b,LL mod)
{
	a%=mod,b%=mod;
	return ((a*b-(LL)((LL)((LD)a/mod*b+1e-3)*mod))%mod+mod)%mod;
}
I int inv(int x,int mod)
{
	return dg(x,mod-2,mod);
}
I void NTT(int *now,int limit,int op)
{
	for(R int i=0;i<limit;i++)
		if(i<r[i]) swap(now[i],now[r[i]]);
	for(R int mid=1;mid<limit;mid<<=1)
	{
		int wn=(int)dg(op==1?G:inv_G,(mod-1)/(mid<<1),mod);
		for(int j=0;j<limit;j+=(mid<<1))
		{
			int w=1;
			for(int k=0;k<mid;k++,w=(int)(((LL)w*wn)%mod))
			{
				LL x=now[j+k],y=(int)((LL)w*now[j+k+mid]%mod);
				now[j+k]=(int)(((LL)x+y)%mod);
				now[j+k+mid]=(int)(((LL)x-y+mod)%mod);
			}
		}
	}
}
I int merge(LL x,LL y,LL z)
{
	LL A=(times(x*Mod[1]%M,dg(Mod[1]%Mod[0],Mod[0]-2,Mod[0]),M)+
		  times(y*Mod[0]%M,dg(Mod[0]%Mod[1],Mod[1]-2,Mod[1]),M))%M;
	LL K=((z-A)%Mod[2]+Mod[2])%Mod[2]*dg(M%Mod[2],Mod[2]-2,Mod[2])%Mod[2];
	return (int)(((K%MOD)*(M%MOD)%MOD+A%MOD)%MOD);
}
I void mul(int *f,int *g,int *h,int n,int limit,int md)
{
	for(int i=n;i<limit;i++)
		a[i]=b[i]=h[i]=0;
	for(R int i=0;i<n;i++)
		a[i]=f[i],b[i]=g[i];
	mod=md;
	inv_G=inv(G,mod);
	NTT(a,limit,1);
	NTT(b,limit,1);
	for(R int i=0;i<limit;i++)
		a[i]=((LL)a[i]*b[i])%mod;
	NTT(a,limit,-1);
	int INV=(int)inv(limit,mod);
	for(R int i=0;i<n;i++)
		h[i]=(LL)a[i]*INV%mod;
}
void work(int n,int *f,int *g)
{
	if(n==1)
	{
		g[0]=inv(f[0],MOD);
		return;
	}
	work((n+1)>>1,f,g);
	limit=1,l=0;
	while(limit<(n<<1))
		limit<<=1,l++;
	for(R int i=1;i<limit;i++)
		r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
	mul(f,g,ans[0],n,limit,Mod[0]),mul(f,g,ans[1],n,limit,Mod[1]),mul(f,g,ans[2],n,limit,Mod[2]);
	ans[0][0]=ans[1][0]=ans[2][0]=1;
	for(R int i=1;i<=n;i++)
		ans[0][i]=ans[1][i]=ans[2][i]=(-merge(ans[0][i],ans[1][i],ans[2][i])+MOD)%MOD;
	mul(ans[0],g,ans[0],n,limit,Mod[0]),mul(ans[1],g,ans[1],n,limit,Mod[1]),mul(ans[2],g,ans[2],n,limit,Mod[2]);
	for(R int i=0;i<n;i++)
		g[i]=merge(ans[0][i],ans[1][i],ans[2][i]);
}
int main()
{
	n=read();
	for(R int i=0;i<n;i++)
		f[i]=read();
	work(n,f,g);
	for(R int i=0;i<n;i++)
		printf("%d ",g[i]);
}
  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值