分治FFT&&多项式求逆学习笔记

最新推荐文章于 2021-07-06 19:56:42 发布
最新推荐文章于 2021-07-06 19:56:42 发布
阅读量422 收藏 2
点赞数
CC 4.0 BY-SA版权
分类专栏： 分治fft 多项式求逆
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/ylsoi/article/details/86761042
本文介绍分治FFT算法解决特殊卷积问题及多项式求逆的高效方法，包括算法原理、优化技巧和具体实现。通过实例展示如何利用分治FFT优化多项式运算，以及如何快速求解多项式的逆。

摘要生成于 C知道 ，由 DeepSeek-R1 满血版支持， 前往体验 >

分治FFT：

解决的是形似以下的问题：

给定n次多项式g(x)g(x)g(x)，求多项式f(x)f(x)f(x)，其中fff的第iii项系数的表达式为∑fj×gi−j\sum f_{j}\times g_{i-j}∑fj​×gi−j​。

解法：

不难发现式子也是卷积的形式，但是与普通多项式乘法不一样的是，每一项的系数依赖前面的项的系数，使得普通的FFT无法起作用。

考虑分治，将区间[l,r][l,r][l,r]分为[l,mid],[mid+1,r][l,mid],[mid+1,r][l,mid],[mid+1,r]两个区间计算，计算完[l,mid][l,mid][l,mid]中的多项式的系数之后，可以很方便的将[l,mid][l,mid][l,mid]对[mid+1,r][mid+1,r][mid+1,r]的贡献给算出来，贡献的计算同上面的递推式。

由于贡献的计算也需要依赖前面的系数，所以分治时的进程是：
divide(l,mid);
calc([l,mid]->[mid+1,r]);
divide(mid+1,r);

不难发现每个点之前的部分会被拆成若干个区间累加到这个位置上，所以正确性显而易见。

/*=======================================
 * Author : ylsoi
 * Time : 2019.1.29
 * Problem : luogu4721
 * E-mail : ylsoi@foxmail.com
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
    freopen("luogu4721.in","r",stdin);
    freopen("luogu4721.out","w",stdout);
}

template<typename T>void read(T &_){
    _=0; T fl=1; char ch=getchar();
    for(;!isdigit(ch);ch=getchar())if(ch=='-')fl=-1;
    for(;isdigit(ch);ch=getchar())_=(_<<1)+(_<<3)+(ch^'0');
    _*=fl;
}

const int maxn=1e5+10;
const int mod=998244353;

ll qpow(ll x,ll y){
    ll ret=1; x%=mod;
    while(y){
        if(y&1)ret=ret*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return ret;
}

struct ksslbh{
    int lim,cnt,dn[maxn<<2];
    ll g[maxn<<2],ig[maxn<<2];
    void ntt(ll *A,int ty){
        REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
        for(int len=1;len<lim;len<<=1){
            ll w= ty==1 ? g[len<<1] : ig[len<<1];
            for(int L=0;L<lim;L+=len<<1){
                ll wk=1;
                REP(i,L,L+len-1){
                    ll u=A[i],v=A[i+len]*wk%mod;
                    A[i]=(u+v)%mod;
                    A[i+len]=(u-v)%mod;
                    wk=wk*w%mod;
                }
            }
        }
    }
    void work(int n,int m,ll *a,ll *b){
        lim=1,cnt=0;
        while(lim<=n+m)lim<<=1,++cnt;
        if(!cnt)cnt=1;
        REP(i,0,lim){
            dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
            if(i>n)a[i]=0;
            if(i>m)b[i]=0;
        }
        g[lim]=qpow(3,(mod-1)/lim);
        ig[lim]=qpow(g[lim],mod-2);
        for(int i=lim>>1;i;i>>=1){
            g[i]=g[i<<1]*g[i<<1]%mod;
            ig[i]=ig[i<<1]*ig[i<<1]%mod;
        }
        ntt(a,1),ntt(b,1);
        REP(i,0,lim-1)a[i]=a[i]*b[i]%mod;
        ntt(a,-1);
        ll inv=qpow(lim,mod-2);
        REP(i,0,lim-1)a[i]=a[i]*inv%mod;
    }
}T;

int n;
ll f[maxn<<2],g[maxn<<2],a[maxn<<2],b[maxn<<2];

void divide(int l,int r){
    if(l==r)return;
    int mid=(l+r)>>1;
    divide(l,mid);
    int nn=mid-l,mm=r-l-1;
    //nn means the pows of f/a
    //mm means the pows of g/b
    REP(i,0,nn)a[i]=f[l+i];
    REP(i,0,mm)b[i]=g[i+1];
    T.work(nn,mm,a,b);
    REP(i,mid+1,r)f[i]=(f[i]+a[i-1-l])%mod;
    divide(mid+1,r);
}

int main(){
    //File();
    read(n);
    REP(i,1,n-1)read(g[i]);
    f[0]=1;
    divide(0,n);
    REP(i,0,n-1)printf("%lld ",(f[i]+mod)%mod);
    return 0;
}

多项式求逆：

对于一个多项式A(x)A(x)A(x)，找到一个多项式A−1(x)A^{-1}(x)A−1(x)，使得A(x)×A−1(x)A(x)\times A^{-1}(x)A(x)×A−1(x)只有常数项的系数为1，其余项的系数全部为0，那么我们称A−1(x)A^{-1}(x)A−1(x)为多项式A(x)A(x)A(x)的逆。

当常数项存在逆元的时候，利用递推的方法去求解多项式的逆，不难发现所得的多项式是一个无限项的级数，所以我们一般只取前面n−1n-1n−1次，即对xnx^nxn取模。

暴力求法：

从低次依次向高次递推，每次求解nnn次项的系数ana_nan​时，将前nnn项带入，所得的和需为000，直接解出来即可。

优化：

观察上面的递推过程，不难发现后面的项求解需要用到前面的项的信息，利用分治FFT可以优化到nlog⁡2nn\log^2nnlog2n。

更加优秀的解法：

直接一项一项地递推显然太浪费时间，考虑利用倍增的思想来利用前面的项推出后面的项。

假设我们现在得到了前n2\frac{n}{2}2n​项的系数，那么我们就要想办法得出前nnn项的系数。

假设前n2\frac{n}{2}2n​项的多项式为B(x)B(x)B(x)，最终的多项式为A−1(x)A^{-1}(x)A−1(x)，那么我们可以得到
A−1(x)−B(x)≡0mod&ThinSpace;&ThinSpace;xn2 A^{-1}(x)-B(x)\equiv 0 \mod x^{\frac{n}{2}} A−1(x)−B(x)≡0modx2n​
为了使模数提升至xnx^{n}xn，可以将等式两边平方得：
(A−1(x)−B(x))2≡0mod&ThinSpace;&ThinSpace;xn (A^{-1}(x)-B(x))^2\equiv 0 \mod x^{n} (A−1(x)−B(x))2≡0modxn
考虑证明这一步，任意一个i(i&lt;n)i (i &lt; n)i(i<n)次项的系数的拆分必定有一个数&lt;n2&lt; \frac{n}{2}<2n​，而任意&lt;n2&lt; \frac{n}{2}<2n​次项的系数都为0，所以i(i&lt;n)i (i&lt;n)i(i<n)次项的系数都为000。

同时注意这种变化只有在等号右边等于000时才成立，这也是为什么要拿两式相减的理由。

然后将平方项拆开：
A−2(x)−2×A−1(x)×B(x)+B2(x)≡0mod&ThinSpace;&ThinSpace;xn A^{-2}(x)-2\times A^{-1}(x)\times B(x)+B^{2}(x)\equiv 0 \mod x^n\\ A−2(x)−2×A−1(x)×B(x)+B2(x)≡0modxn
两边同时乘以A(x)A(x)A(x)：
A−1(x)−2×B(x)+A(x)×B2(x)≡0mod&ThinSpace;&ThinSpace;xnA−1(x)≡2×B(x)−A(x)×B2(x)mod&ThinSpace;&ThinSpace;xn A^{-1}(x)-2\times B(x)+A(x)\times B^2(x)\equiv 0 \mod x^n\\ A^{-1}(x)\equiv 2\times B(x)-A(x)\times B^{2}(x) \mod x^n\\ A−1(x)−2×B(x)+A(x)×B2(x)≡0modxnA−1(x)≡2×B(x)−A(x)×B2(x)modxn
然后我们就成功地得到了在模xnx^nxn意义下的多项式A−1(x)A^{-1}(x)A−1(x)。

具体实现：

考虑递归实现上述过程。

实现过程中，理论上A(x)A(x)A(x)和B(x)B(x)B(x)项数都很多，但为了保证复杂度，A(x)A(x)A(x)每次取对计算有贡献的前nnn项，而B(x)B(x)B(x)原本是在模xn2x^\frac{n}{2}x2n​意义下的多项式，只取前n2\frac{n}{2}2n​项即可，所以最后乘出来的多项式为2n2n2n次的，时间复杂度满足T(n)=nlog⁡n+T(n2)T(n)=n\log n+T(\frac{n}{2})T(n)=nlogn+T(2n​)。

所以最终时间复杂度为nlog⁡nn\log nnlogn。


 /*=======================================
 * Author : ylsoi
 * Time : 2019.2.2
 * Problem : luogu4238
 * E-mail : ylsoi@foxmail.com
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
    freopen("luogu4238.in","r",stdin);
    freopen("luogu4238.out","w",stdout);
}

template<typename T>void read(T &_){
	_=0; T fl=1; char ch=getchar();
	for(;!isdigit(ch);ch=getchar())if(ch=='-')fl=-1;
	for(;isdigit(ch);ch=getchar())_=(_<<1)+(_<<3)+(ch^'0');
	_*=fl;
}

const int maxn=1e5+10;
const int mod=998244353;

int lim,cnt,dn[maxn<<2];
ll g[maxn<<2],ig[maxn<<2];

ll qpow(ll x,ll y){
	ll ret=1; x%=mod;
	while(y){
		if(y&1)ret=ret*x%mod;
		x=x*x%mod;
		y>>=1;
	}
	return ret;
}

void ntt(ll *A,int ty){
	REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
	for(int len=1;len<lim;len<<=1){
		ll w= ty==1 ? g[len<<1] : ig[len<<1];
		for(int L=0;L<lim;L+=len<<1){
			ll wk=1;
			REP(i,L,L+len-1){
				ll u=A[i],v=A[i+len]*wk%mod;
				A[i]=(u+v)%mod;
				A[i+len]=(u-v)%mod;
				wk=wk*w%mod;
			}
		}
	}
}

int n;
ll f[maxn<<2],a[maxn<<2],b[maxn<<2];

void init(){
	read(n);
	REP(i,0,n-1)read(f[i]);
	lim=1;
	while(lim<=n+n)lim<<=1;
	g[lim]=qpow(3,(mod-1)/lim);
	ig[lim]=qpow(g[lim],mod-2);
	for(int i=lim>>1;i;i>>=1){
		g[i]=g[i<<1]*g[i<<1]%mod;
		ig[i]=ig[i<<1]*ig[i<<1]%mod;
	}
}

void poly_inv(int k){
	if(k==1){
		b[0]=qpow(f[0],mod-2);
		return;
	}
	poly_inv((k+1)>>1);
	lim=1,cnt=0;
	while(lim<=k+k)lim<<=1,++cnt;
	if(!cnt)cnt=1;
	REP(i,0,lim-1)dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
	REP(i,0,k-1)a[i]=f[i];
	REP(i,k,lim-1)a[i]=0;
	ntt(a,1),ntt(b,1);
	REP(i,0,lim-1)b[i]=b[i]*(2ll-a[i]*b[i]%mod)%mod;
	ntt(b,-1);
	ll inv=qpow(lim,mod-2);
	REP(i,0,k-1)b[i]=b[i]*inv%mod;
	REP(i,k,lim)b[i]=0;
}

int main(){
    File();
	init();
	poly_inv(n);
	REP(i,0,n-1)printf("%lld ",(b[i]%mod+mod)%mod);
    return 0;
}
【模板】分治FFT -cdq分治/多项式求逆
远行客
11-21 386
分治 FFT
[学习笔记]FFT之多项式求逆与一类组合计数 DP 问题
女装的你，如此好看！
07-10 711
一、引入 有时候，我们会遇到这种形式的问题： 给定一个定义域为 [0,n)[0,n)[0,n) 且自变量为整数的函数（数组） fff ，有一个未知函数 ggg ， 给定 fff 与 ggg 的卷积 hhh ，求 g=hfg=hfg=\frac hf 。 这时候要用到一个技巧——多项式求逆。 二、概念 一个 nnn 项 n−1n−1n-1 次多项式 AAA 的逆元，定义为一个多项式 ...
多项式求逆/分治FFT 学习笔记
weixin_34080951的博客
01-11 320
一、多项式求逆 给定一个多项式\(F(x)\)，请求出一个多项式\(G(x)\)， 满足\(F(x) * G(x) \equiv 1 ( \mathrm{mod\:} x^n )\)。系数对\(998244353\)取模。 考虑递归求解，当\(F\)的最高次为\(0\)时，\(G_0=F_0^{-1}\) 假设我们知道了\(F(x)\)在模\(x^{\left \lceil \fra...
洛谷 P4721 [模板]分治FFT —— 分治FFT / 多项式求逆
aodan5477的博客
11-30 174
题目：https://www.luogu.org/problemnew/show/P4721 分治做法，考虑左边对右边的贡献即可； 注意最大用到的 a 的项也不过是 a[r-l] ，所以 NTT 可以只做到 2*(r-l)，能快一倍。 代码如下： #include<iostream> #include<cstdio> #include<cst...
多项式求逆(倍增或分治FFT)
jziwjxjd的博客
12-21 398
传送门 F(x)∗H(x)≡ 1 (mod  x⌈n2⌉)由题意得F(x)∗G(x)≡1(mod  x⌈n2⌉)那么F(x)∗[G(x)−H(x)]≡0(mod  x⌈n2⌉)就是G(x)−H(x)≡0(mod  x⌈n2⌉)两边同时平方得到G(x)2+H(x)2−2G(x)H(x)≡0(mod  xn)模数变成xn是有原因的,我们定义f(x)=G(x)−H(x),f(x)2为f(x)∗f(x)因为fi2=∑j=0ifj∗fi−j其中i和i−j不可能同时大于n2,所以一定有一项为0得证我
【学习笔记】超简单的多项式求逆（含全套证明）
繁凡さん的博客
12-30 3529
整理的算法模板合集： ACM模板 目录多项式求逆一、分治FFT二、倍增法及其证明三、多项式求逆例题[P4238 【模板】多项式乘法逆](https://www.luogu.com.cn/problem/P4238) 点我看多项式全家桶(●^◡_◡◡​^●) 多项式求逆 一、分治FFT 二、倍增法及其证明 三、多项式求逆例题 P4238 【模板】多项式乘法逆 直接按照上述的思路实现一下模板即可。注意要中间 h(x)h(x)h(x) 要取的是模 x⌈n⌉2x^{\frac{\lceil n \rceil}
【学习笔记】分治FFT
繁凡さん的博客
07-06 1506
分治FFT 1. Luogu P4721 【模板】分治 FFT Problem 给定序列 g1…n−1g_{1\dots n - 1}g1…n−1​ ，求序列 f0…n−1f_{0\dots n - 1}f0…n−1​ 。 其中 fi=∑j=1ifi−jgjf_i=\sum_{j=1}^if_{i-j}g_jfi​=∑j=1i​fi−j​gj​ ，边界为 f0=1f_0=1f0​=1。 答案对 998244353998244353998244353 取模。 Solution 显然直接暴力FFT是 O(n2l
[学习笔记][省选算法]多项式乘法之 FFT & NTT 及卷积
女装的你，如此好看！
06-29 1528
一、开头 （Place ：山东省神犇协会第 998244353 会议厅） SD 神犇 998244353 号（会长）：所有神犇，洛谷黑题都切完了吗？ 全体神犇：切完了！ 神犇 1004535809 号： 998244353 ，我昨天看到一个弱鸡，洛谷名叫 xyz32768 ，他一道黑题都没做，您认为应该怎样呢？ 神犇 998244353 号：这个问题您们考虑一下，应该用怎样的方式把 xy...
FFT多项式乘法学习笔记
过路人丶的博客
05-22 534
基础知识了解:点击打开链接公式证明:点击打开链接综合了这两个博客我才懂了什么是FFT，所以这两个博客讲的比较细的我就不说了，说一些综合的.为什么选择n个n次单位复根？公式1： w(n.k)^2 = w(n/2,k)应用了这个公式才可以把原来n大小的问题变成n/2的问题,这就是分治得到log(n)的由来.将n个数分成奇数偶数分开处理得到公式F(x) = G(x^2) + xH(x^2),其中G(x)...
【分治FFT+多项式求逆】LGP4705 玩游戏
MEET YOU FIRST TIME
04-12 239
【题目】 LG 给定两个序列a,ba,ba,b，定义kkk次价值为分别随机从a,ba,ba,b中取出一个数ax,bxa_x,b_xax​,bx​，(ax+bx)k(a_x+b_x)^k(ax​+bx​)k的期望。 对于所有k∈[1,t]k\in[1,t]k∈[1,t]，求序列的kkk次价值。 ∣a∣,∣b∣,k≤105|a|,|b|,k\leq 10^5∣a∣,∣b∣,k≤105 【解题思路】 化...
[n点无向图个数 分治FFT || 多项式求逆] BZOJ 3456 城市规划
雯舞
09-14 874
然后就分治FFT了 #include #include #include using namespace std; typedef long long ll; inline char nc(){ static char buf[100000],*p1=buf,*p2=buf; if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin)
洛谷 P4239 【模板】多项式求逆（加强版）任意模数fft
liangzihao1的博客
08-19 738
题目描述 给定一个多项式F(x)F(x)F(x)，请求出一个多项式G(x)G(x)G(x)，满足F(x)∗G(x)≡1(mod&amp;amp;nbsp;xn)F(x)∗G(x)≡1(mod&amp;amp;nbsp;xn)F(x)∗G(x)≡1(mod\ x^n )。系数对109+7109+710^9+7取模。 输入输出格式 输入格式： 首先输入一个整数nnn，表示输入多项式的次数。 接着输入nnn个整数，第iii个...
多项式求逆
qmmctzs的博客
02-18 316
前置知识 NTT 多项式取模: 多项式模 xnx^{n}xn 表示取多项式的前 nnn 位 多项式求逆 给定 f(x)=a0+a1x1+a2x2+⋅⋅⋅+anxnf(x)=a_{0}+a_{1}x^{1}+a_{2}x^{2}+···+a_{n}x^{n}f(x)=a0​+a1​x1+a2​x2+⋅⋅⋅+an​xn 求出 g(x)=b0+b1x1+b2x2+⋅⋅⋅+bkxk   (k≤n)g(x)=b_{0}+b_{1}x^{1}+b_{2}x^{2}+···+b_{k}x^{
FFT算法实现多项式乘法的程序解析
资源摘要信息:"快速傅里叶变换（FFT）是数字信号处理领域中一种非常重要的算法，主要用于实现多项式的快速相乘。本资源中的FFT算法相关文件提供了实现两个多项式系数序列通过FFT算法进行相乘的方法，并给出了相应的...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值