模板:多项式乘法(FFTNTT)

文章目录

  • 前言
  • 系数表示法和点值表示法
  • 单位根
  • 离散傅立叶变换(DFT)
  • 位逆序置换(蝴蝶变换)
  • 离散傅立叶逆变换(IDFT)
  • 代码
  • NTT
  • 代码

所谓快速傅立叶变换,就是傅立叶发明的一种快速的变换

(逃)

前言

FFT,用于在nlognnlognnlogn的复杂度内求两个多项式相乘
算是多项式的入门吧。
十一被rabbit的数学打蒙了所以开始碰这些东西。
不知道以后会不会直接鸽掉
(update on 2022.1.4:三个月后终于回来填坑了)。
很妙的算法。
代码写成迭代,精炼之后很优美。

系数表示法和点值表示法

系数表示法就是我们常用的函数的表示形式。
对于一个n次的函数,点值表示法则是给出了n+1个互异的点,通过这n+1个互异的点就可以唯一确定的描述出这个函数(可以想想高斯消元)。

这东西有啥用?
考虑如果两个多项式都表示成了点值形式,且选取的点的横坐标相同,我把它们的纵坐标对应乘起来,就能得到它们乘积的多项式的点值表示
乘起来的复杂度降到了O(n)O(n)O(n)

但是这个点值表示法求起来暴力似乎还是n方啊…
所以快速傅立叶变换其实解决的就是系数表示法和点值表示法快速互相转换的问题

单位根

规定:

若复数ω\omegaω满足ωn=1\omega^n=1ωn=1,就称ω\omegaωn次单位根

对于一个n项的多项式,我们让它的点值表示取的横坐标为ωnk(0≤k<n)\omega_n^k(0\leq k <n)ωnk0k<n)
然后就会有很多很妙的性质:

  1. ωnmkm=ωnk\omega_{nm}^{km}=\omega_n^kωnmkm=ωnk
  2. ωnk+n=ωnk\omega_{n}^{k+n}=\omega_n^kωnk+n=ωnk
  3. ωnk+n2=−wnk\omega_n^{k+\frac{n}{2}}=-w_n^kωnk+2n=wnk

都可以结合几何意义较为显然的证明。

离散傅立叶变换(DFT)

那么回到问题。
首先,我们要把多项式填零直到 n=2kn=2^kn=2k 的长度。
接下来,我们要求:
f(x)=a0+a1x+a2x2+a3x3+...+an−1xn−1f(x)=a_0+a_1x+a_2x^2+a_3x^3+...+a_{n-1}x^{n-1}f(x)=a0+a1x+a2x2+a3x3+...+an1xn1
让我们把这个函数分奇偶拆成两个函数:
A(x)=a0+a2x+a4x2+a6x3+...+an−1xn2A(x)=a_0+a_2x+a_4x^2+a_6x^3+...+a_{n-1}x^{\frac{n}{2}}A(x)=a0+a2x+a4x2+a6x3+...+an1x2n
B(x)=a1+a3x+a5x2+a7x3+...+an−2xn2B(x)=a_1+a_3x+a_5x^2+a_7x^3+...+a_{n-2}x^{\frac{n}{2}}B(x)=a1+a3x+a5x2+a7x3+...+an2x2n
那么就有:
f(x)=A(x2)+xB(x2)f(x)=A(x^2)+xB(x^2)f(x)=A(x2)+xB(x2)
我们把 x=ωnkx=\omega_n^kx=ωnk 带入。
分情况讨论:
(设:0≤k<n20\le k<\dfrac{n}{2}0k<2n
当指数 <n2< \dfrac{n}{2}<2n 时:
f(ωnk)=A((ωnk)2)+ωnkB((ωnk)2)f(\omega_n^k)=A((\omega_n^k)^2)+\omega_n^kB((\omega_n^k)^2)f(ωnk)=A((ωnk)2)+ωnkB((ωnk)2)
=A(ωn2k)+ωnkB(ωn2k)=A(\omega_n^{2k})+\omega_n^kB(\omega_n^{2k})=A(ωn2k)+ωnkB(ωn2k)
=A(ωn2k)+ωnkB(ωn2k)=A(\omega_{\frac{n}{2}}^{k})+\omega_n^kB(\omega_{\frac{n}{2}}^{k})=A(ω2nk)+ωnkB(ω2nk)
类似的,当指数 ≥n2\ge\dfrac{n}{2}2n 时:
f(ωnk+n2)=A(ωn2k+n)+ωnk+n2B(ωn2k+n)f(\omega_n^{k+\frac{n}{2}})=A(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}B(\omega_n^{2k+n})f(ωnk+2n)=A(ωn2k+n)+ωnk+2nB(ωn2k+n)
=A(ωn2k)−ωnkB(ωn2k)=A(\omega_n^{2k})-\omega_n^{k}B(\omega_n^{2k})=A(ωn2k)ωnkB(ωn2k)
=A(ωn2k)−ωnkB(ωn2k)=A(\omega_{\frac{n}{2}}^{k})-\omega_n^{k}B(\omega_{\frac{n}{2}}^{k})=A(ω2nk)ωnkB(ω2nk)
这样,我们就可以在 O(len)O(len)O(len) 的复杂度内合并两个长度为 lenlenlen 的多项式,就可以利用分治把复杂度降到 O(nlog⁡n)O(n\log n)O(nlogn)

位逆序置换(蝴蝶变换)

这个也很妙
有点线性求逆元那味了
通过观察发现,系数 i 调换位置之后会去往它的二进制表示反过来的位置(记为rir_iri
而上面那个东西可以通过:
r[i]=r[i>>1]>>1+[i&1]∗(n>>1)r[i]=r[i>>1]>>1+[i\&1]*(n>>1) r[i]=r[i>>1]>>1+[i&1](n>>1)
来线性的求
这个还是挺好理解的画一画就行了

离散傅立叶逆变换(IDFT)

我们变完之后当然还要变回来。
我会 n3n^3n3 高斯消元!
设之前DFT得到的点值表示 bk=f(ωnk)=∑i=0n−1ai(ωnk)ib_k=f(\omega_n^k)=\sum_{i=0}^{n-1}a_i(\omega_{n}^k)^ibk=f(ωnk)=i=0n1ai(ωnk)i
构造一个函数:
g(x)=b0+b1x+b2x2+b3x3+...+bn−1xn−1g(x)=b_0+b_1x+b_2x^2+b_3x^3+...+b_{n-1}x_{n-1}g(x)=b0+b1x+b2x2+b3x3+...+bn1xn1
然后我们用 ωn0,ωn−1,ωn−2,...,ωn−n\omega_n^{0},\omega_n^{-1},\omega_n^{-2},...,\omega_n^{-n}ωn0,ωn1,ωn2,...,ωnn 作为横坐标带入。(不难发现刚才 DFT 需要的性质依然成立,所以可以用同样的方法求解)
那么我们在第 kkk 位得到的新的点值就是:
g(ωn−k)=∑i=0n−1bi(ωn−k)ig(\omega_n^{-k})=\sum_{i=0}^{n-1}b_i(\omega_n^{-k})^ig(ωnk)=i=0n1bi(ωnk)i
=∑i=0n−1∑j=0n−1aj(ωni)j(ωn−k)i=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_{n}^i)^j(\omega_n^{-k})^i=i=0n1j=0n1aj(ωni)j(ωnk)i
=∑j=0j−1aj∑i=0n−1(ωnj−k)i=\sum_{j=0}^{j-1}a_j\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i=j=0j1aji=0n1(ωnjk)i
考虑 ∑i=0n−1(ωnj−k)i\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^ii=0n1(ωnjk)i,这一项,当 j=kj=kj=k 时,显然等于0;当 j≠kj\ne kj=k 时,根据等比公式,有:
∑i=0n−1(ωnj−k)i=(wnj−k)n−1wnj−k−1\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i=\frac{(w_n^{j-k})^n-1}{w_n^{j-k}-1}i=0n1(ωnjk)i=wnjk1(wnjk)n1
=(wnn)j−k−1wnj−k−1=1j−k−1wnj−k−1=0=\frac{(w_n^{n})^{j-k}-1}{w_n^{j-k}-1}=\frac{1^{j-k}-1}{w_n^{j-k}-1}=0=wnjk1(wnn)jk1=wnjk11jk1=0
所以原来的式子只有当 j=kj=kj=k 的时候有值,等于 nakna_knak
所以我们直接 NTT 完除 nnn 即可。

代码

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define il inline
const int N=5e6+100;
const int M=150;
const int mod=998244353;
const double pi=acos(-1.0);
inline ll read(){ll x=0,f=1;char c=getchar();while(!isdigit(c)){if(c=='-') f=-1;c=getchar();}while(isdigit(c)){x=x*10+c-'0';c=getchar();}return x*f;
}
int n,m,len,k;
struct node{double x,y;node(double a=0,double b=0){x=a;y=b;}
}A[N],B[N];
il node operator * (node a,node b){return (node){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
}
il node operator + (node a,node b){return (node){a.x+b.x,a.y+b.y};
}
il node operator - (node a,node b){return (node){a.x-b.x,a.y-b.y};
}
int r[N];
il void fft(node *x,int flag){for(int i=0;i<=len;i++){if(i<r[i]) swap(x[i],x[r[i]]);}for(int l=1;l<len;l<<=1){node o(cos(pi/l),flag*sin(pi/l));for(int st=0;st<len;st+=l<<1){node t(1,0);for(int j=0;j<l;j++,t=t*o){node u=x[st+j],v=t*x[st+j+l];x[st+j]=u+v;x[st+j+l]=u-v;}}}return;
}
int main(){n=read();m=read();for(int i=0;i<=n;i++) A[i].x=read();for(int i=0;i<=m;i++) B[i].x=read();len=n+m;while((1<<k)<=len) k++;len=1<<k;for(int i=1;i<=len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));fft(A,1);fft(B,1);for(int i=0;i<=len;i++) A[i]=A[i]*B[i];fft(A,-1);for(int i=0;i<=n+m;i++) printf("%d ",(int)(A[i].x/(len)+0.5));return 0;
}
/**/

NTT

NTT说简单也简单说难也难。
简单是因为板子几乎和FFT一模一样。
难是因为那个群的概念完全看不懂qwq。
然后我就选择调过了证明部分。

背下来背下来!——宋队

总的来说,写一写对背板子有帮助的我可以理解的。
对于一个质数P,若其存在原根(基本全是3)且可以表示为:
P=k∗2n+1P=k*2^n+1P=k2n+1
它就是一个ntt模数
ggg为什么叫原根?
因为它有一些关键性质:

1.gk∗n=1(modP)1. g^{k*n}=1(mod P)1.gkn=1(modP)
2.gi(modP)g^i(mod P)gi(modP)0≤i≤p0\leq i\leq p0ip 的范围内取值两两不同。

设:gkg^kgk单位根,记为gng_ngn
那么它就有很多熟悉的性质:

  1. gnn=1(modP)g_n^n=1(mod P)gnn=1(modP)
  2. g2k2k=gkkg_{2k}^{2k}=g_k^kg2k2k=gkk
  3. gnk+n/2=−gnkg_n^{k+n/2}=-g_n^kgnk+n/2=gnk

没错,推FFT时 ω\omegaω 需要有的性质它全有!
所以就可以直接上FFT的板子啦!
然后烦人的复数、浮点运算和精度问题全都拜拜了。
针不戳。
宋队诚不欺我。
逆变换还是倒数,从FFT的求共轭改成求逆元。
就像喝水一样。

代码

il void ntt(ll *x,int flag){for(int i=0;i<lim;i++){if(i<r[i]) swap(x[i],x[r[i]]);}for(register int l=1;l<lim;l<<=1){register ll g=ksm(3,(mod-1)/(l<<1));if(flag==-1) g=ksm(g,mod-2);for(register int st=0;st<lim;st+=l<<1){register ll now=1;for(register int i=0;i<l;i++,now=now*g%mod){ll u=x[st+i],v=x[st+i+l]*now%mod;x[st+i]=(u+v)%mod;x[st+i+l]=(u-v+mod)%mod;}}}if(flag==-1){register ll ni=ksm(lim,mod-2);for(register int i=0;i<lim;i++) x[i]=x[i]*ni%mod;}
}

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/317905.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

Loj #149. 01 分数规划(01分数规划模板题)

链接 题意&#xff1a; 题解&#xff1a; 详细解法看这里 这里说个点&#xff0c;eps一定要开足够小&#xff0c;我一开始开的1e-5&#xff0c;结果就过了90%的数据&#xff0c;开到1e-7就足够了 代码&#xff1a; #include<bits/stdc.h> typedef long long ll; usin…

CF986C AND Graph(图论+二进制连边)

CF986C AND Graph \(\color{yellow}{\bigstar\texttt{Hint}}\)&#xff1a;和每个点连接的点是这个数取反后的子集&#xff0c;考虑将这个点和它的反连边&#xff0c;那么所有对应的数的子集都是同一个连通块内的。 之后的一种简单的寻找连通块就是直接对所有没有访问过的点暴力…

AT2376-[AGC014D]Black and White Tree【结论,博弈论】

正题 题目链接:https://www.luogu.com.cn/problem/AT2376 题目大意 给出nnn个点的一棵树&#xff0c;先后手轮流选择一个未染色的点染上白色(先手)/黑色(后手)&#xff0c;如果最后有一个白色的点连接的都是白色的点则先手获胜&#xff0c;否则后手获胜。求是否先手必胜。 1≤…

如何测试 ASP.NET Core Web API

在本文中&#xff0c;我们将研究如何测试你的 ASP .NET Core 2.0 Web API 解决方案。我们将了解使用单元测试进行内部测试&#xff0c;使用全新的 ASP .NET Core 的集成测试框架来进行外部测试。本文要点正确理解并使用单元测试和你的 ASP .NET Core Web API 解决方案一样重要。…

[帝皇杯day 1] [NOIP2018模拟赛]小P的loI(暴力+素筛),【NOIP模拟赛】创世纪(贪心),无聊的数对(线段树)

文章目录T1&#xff1a;小P的loltitlesolutioncodeT2&#xff1a;创世纪titlesolutioncodeT3&#xff1a;无聊的数对titlesolutioncodeT1&#xff1a;小P的lol title solution 此题非常水… 先用素数筛&#xff0c;筛出[1,n][1,n][1,n]中的质数 质数越小&#xff0c;倍数的分…

P6144 [USACO20FEB]Help Yourself P(DP+线段树)

P6144 [USACO20FEB]Help Yourself P 将线段按照了 \(r\) 排序&#xff0c;设右端点为 \(r\) 的答案为 \(f_r\)&#xff0c;发现这样转移非常困难。 \(\color{yellow}{\bigstar\texttt{Trick}}\)&#xff1a;区间覆盖的题要按照左端点排序&#xff0c;记右端点为 \(r\) 时的答案…

Desert King POJ - 2728

题意&#xff1a; 给定N个平面上的点的坐标和它们的权值&#xff0c;任意两点之间的边的价值是它们的距离&#xff0c;费用是两点权值之差的绝对值&#xff0c;求该图的一棵生成树&#xff0c;使得该树所有边的费用之和与价值之和的比值最小&#xff08;只需求这个比值即可&am…

洛谷P3338:力(FFT)

传送门 解析 算是比较适合的FFT入门题了吧 一个重要的trick&#xff1a; 当函数无法表示成卷积时&#xff0c;可以把函数翻转过来 然后调一调就又是卷积了 一个重要的注意事项是FFT的lim一定是两多项式相乘结果多项式的项数&#xff01; 即使后面的项根本没有用也一样 其他的…

CF650E-Clockwork Bomb【并查集】

正题 题目链接:https://www.luogu.com.cn/problem/CF650E 题目大意 给出nnn个点的两棵树&#xff0c;你每次可以选择第一棵树上的一条边改成另一条边&#xff0c;但是改完之后必须还是一棵树&#xff0c;求最少的步数把第一棵树改成第二棵树。 1≤n≤51051\leq n\leq 5\times…

[2-sat专练]poj 3683,hdu 1814,hdu 1824,hdu 3622,hdu 4115,hdu 4421

文章目录Priest Johns Busiest DaycodePeaceful CommissioncodeLets go homecodeBomb GamecodeEliminate the ConflictcodeBit MagiccodePriest John’s Busiest Day 题目 司仪必须在婚礼开始或结束时出现&#xff0c;考虑把第iii场婚礼拆成两个点 iii&#xff1a;表示司仪在婚…

.NET Core中的一个接口多种实现的依赖注入与动态选择

最近有个需求就是一个抽象仓储层接口方法需要SqlServer以及Oracle两种实现方式&#xff0c;为了灵活我在依赖注入的时候把这两种实现都给注入进了依赖注入容器中&#xff0c;但是在服务调用的时候总是获取到最后注入的那个方法的实现&#xff0c;这时候就在想能不能实现动态的选…

Sightseeing Cows POJ - 3621

题意&#xff1a; L个点&#xff0c;P边的点边带权的有向图&#xff0c;求一个环点权和与边权和比值的最大值。 题解&#xff1a; 01分数规划判负环 详细看这里 还是套用01分数规划模型&#xff0c;点权为value[i],边权为cost[u],一个环为C&#xff0c;问题要求最大化 最…

CF559E Gerald and Path(DP)

CF559E Gerald and Path 设 \(dp(i,p)\) 表示完成前 \(i\) 条线段的覆盖&#xff0c;最右端位于 \(p\) 点的最大收益。 转移&#xff1f;向下一条线段转移时加上他们中间的距离&#xff1f;发现这样没有办法统计 \(p\) 点以前的空位了&#xff01; \(\color{yellow}{\bigstar\t…

多项式乘法:练习总结

文章目录前言力代码礼物代码差分和前缀和代码开拓者的知识代码总结前言 这两天由于国庆集训全是阴间的生成函数&#xff0c;所以就学了一点点相关的内容 其实就学了个FFT和NTT 也算是点开了一个小小的技能点吧 进入多项式才发现里面世界的广阔 然而由于这玩意没有一个是NOIP考…

CF1444C-Team-Building【可撤销并查集】

正题 题目链接:https://www.luogu.com.cn/problem/CF1444C 题目大意 给出nnn个点mmm条边的一张图&#xff0c;总共kkk个颜色&#xff0c;每个点有一个颜色。 询问有多少无序颜色对(x,y)(x,y)(x,y)满足x≠yx\neq yx​y且颜色为xxx或yyy的点构成的生成子图是一个二分图。 1≤…

C# Memory Cache 踩坑记录

背景前些天公司服务器数据库访问量偏高,运维人员收到告警推送,安排我团队小伙伴排查原因.我们发现原来系统定期会跑一个回归测试,该测运行的任务较多,每处理一条任务都会到数据库中取相关数据,高速地回归测试也带来了高频率的数据库读取.解决方案1我们认为每个任务要取的数据大…

模板:无旋treap

文章目录前言操作合并分裂插入删除查找第k大查询x的排名前驱后继完整代码所谓无旋treap&#xff0c;就是不带旋转的treap 前言 现在“理论上”我会四种平衡树了 之前说无旋treap功能弱&#xff0c;是我狗眼看银低了qwq 这玩意区间上该搞的也能搞 而且确实是非常的 好写&#x…

[贪心专题]CF549G,CF351E,CF226D,CF1276C,CF1148E,CF798D

文章目录T1&#xff1a;CF1276C Beautiful RectangletitlesolutioncodeT2&#xff1a;CF226D The tabletitlesolutioncodeT3&#xff1a;CF549G Happy LinetitlesolutioncodeT4&#xff1a;CF798D Mike and distributiontitlesolutioncodeT5&#xff1a;CF351E Jeff and Permut…

51nod1766-树上的最远点对【结论,线段树】

正题 题目链接:http://www.51nod.com/Challenge/Problem.html#problemId1766 题目大意 给出nnn个点的一棵树&#xff0c;mmm次询问给出两个区间&#xff0c;要求在两个区间中各选一个点使得他们之间距离最大。 1≤n,m≤1051\leq n,m\leq 10^51≤n,m≤105 解题思路 结论就是两…

P3385 【模板】负环

P3385 【模板】负环 题意&#xff1a; 给定一个 n 个点的有向图&#xff0c;请求出图中是否存在从顶点 1 出发能到达的负环。 负环的定义是&#xff1a;一条边权之和为负数的回路。 题解&#xff1a; 先说结论&#xff1a; 判断给定的有向图中是否存在负环。 利用 spfa 算…