FFT+NTT

前置补充：

复数(Complex Number)
定义 $i=\sqrt{-1}$ 所以 $i^2=-1$ 。
那么复数的基本形式可以表示为 $z=a*i+b$，其中 $a ,b$ 为实数.
对于其加减乘,其实就是两个部分分别考虑.
$x_1=a_1*i+b_1,x_2=a_2*i+b_2$
加法:
$x_1+x_2=(a_1+a_2)*i+(b_1+b_2)$
乘法:
$x_1*x_2=(a_1*i+b_1)*(a_2*i) +(a_1*i+b_1)*b_2$
$x_1*x_2=a_1a_2*i^2+a_2b_1*i+a_1b_2*i+b_1b_2$
$\because i^2=-1$
$\therefore x_1*x_2=(a_2b_1+a_1b_2)*i+(b_1b_2-a_1a_2)$
差不多代码里重载一下就是这个样子:
多项式的定义
一个环R上的关于 $x$ 的多项式可以写作
$A(x)=\sum^n_{i=0}a_ix^i$
其中$a_i\in R$,$ x $ 被称为该多项式的自由元.
多项式的运算
加法(对应的系数相加减):
$A(x)\pm B(x)=\sum^n_{i=0}{(a_i\pm b_i)x^i}$
数列卷积(交叉得值)(非多项式)
$c_k=\sum_{i+j=k}{a_ib_j}$
$c_0=a_0b_0$
$c_1=a_0b_1+a_1b_0$
$c_2=a_0b_2+a_1b_1+a_2b_0$
多项式乘法
$A(x)B(x)=\sum^n_{i=0}\sum^n_{j=0}a_ib_jx^{i+j}=\sum^{2n}_{i=0}c_kx^{k}$
其中 $c$ 是 $a$ 和 $b$ 的卷积.
很明显,如果你直接算就是 $O(n^2)$ 的时间复杂度.
多项式的点值和系数表示法
我们知道,你可以用三个点确定一条二次函数曲线,那么 $n$ 个点能不能确定一个最高项为 $n-1$ 的多项式呢?
答案是肯定的.这两个玩意是可以互相转化的.高斯消元可以从点值转成系数
如果给出 $A(x)$ 和 $B(x)$ 分别在 $x_0,x_1,...,x_n$ 下的点值,则可以在$O(n)$的时间内得到 $A(x)\pm B(x)(点值直接加再转回去),A(x)B(x)(乘起来再转回去)$
单位根
满足 $x^n-1=0$ 的 $x$ 被称为 $n$ 次单位根.
实数的单位根只会有-1,1
放到复数上,就会有其他答案
复数模长:$|c|=\sqrt{a^2+b^2}$
复数辐角:$\theta$

$c=|c|\cos \theta+i|c|\sin \theta$
复数乘法有一个重要性质：模长相乘，辐角相加，自己推推即可证明
记一些特殊的模长为1的复数为$ω^k_n$
其中当 $n=8$ 时,可以变成如下图

这八个都叫做8次单位根
对于下面的应用,我们的$ω_n=\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}$

fft

离散傅里叶变换(DFT)
设 $a$ 是长度为 $n$ 的数列,对于 $0\le k<n$ 令
$b_k=\sum^{n-1}_{i=0}a_i.ω_n^{k_i}$
则称 $b$ 为 $a$ 的离散傅里叶变换,记作 $b=F(a)$
就是 $a$ 的每一项乘上 $n$ 次单位根的第 $k_i$ 的位置
可以发现,将 $a_i$ 分离出来,$b_k$就是$A(x)$在$ω_n^k$位置的点值
因此计算 $a$ 的 DFT 与计算 $A(x)$ 在$ω^0_n,ω^1_n,...,ω^{n-1}_n$处的点值表示是等价的.
离散傅里叶逆变换(IDFT)

$a_k=\frac{1}{n}\sum^{n-1}_{i=0}b_i\cdot ω_n^{-k_i}(0\le k < n)$

将一式中的定义带进二式,完成推导

DFT:
$b_k=\sum^{n-1}_{i=0}a_i\cdotω_n^{ki}(0\le k<n)$
IDFT:
$a_k=\frac{1}{n}\sum^{n-1}_{i=0}b_i\cdot ω_n^{-ki}(0\le k< n)$

我们把DFT的定义式带进IDFT，就会有以下的式子

$\sum^{n-1}_{i=0}b_i\cdotω_n^{-k_i}=\sum^{n-1}_{i=0}ω^{-ki}_n\sum^{n-1}_{j=0}ω_n^{ij}a_j$

$\qquad\qquad\qquad\;=\sum_{j=0}^{n-1}a_j\sum^{n-1}_{i=0}ω_n^{-ki}\cdotω_n^{ij}$

$\qquad\qquad\qquad\;=\sum^{n-1}_{j=0}a_j\sum^{n-1}_{i=0}ω_n^{i(j-k)}$

考虑 $\sum_{i=0}^{n-1}ω_n^{i(j-k)}$ 这一部分

当 $j=k$ 时 $\sum_{i=0}^{n-1}ω_n^{i(j-k)}=\sum_{i=0}^{n-1}1=n$

$otherwise:$ $\because 0\le j,k\le n,|j-k|\le{n},ω_n^{j-k}\neq1$

$\therefore \sum^{n-1}_{i=0}ω_n^{i(j-k)}=\sum^{n-1}_{i=0}(ω_n^{j-k})^i$

$\qquad\qquad\qquad\;\;=\frac{1-(ω_n^{j-k})^n}{1-ω_n^{j-k}}$

$\qquad\qquad\qquad\;\;=\frac{1-(ω_n^n)^{j -k}}{1-ω_n^{j-k}}=0$

因此会有以下推导:

$\frac{1}{n}\sum_{i=0}^{n-1}b_i\cdotω_n^{-ki}=\sum_{i=0}^{n-1}ω_n^{-ki}\sum^{n-1}_{j=0}ω^{ij}_{a_j}$

$\qquad\qquad\qquad\quad\;=\frac{1}{n}\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}ω_n^{i(j-k)}$

$\qquad\qquad\qquad\;\quad=\frac{1}{n}\cdot na_k=a_k$

两式是可以互推的,是等价的

所以这就完成了点值和差值的互相转换.
3. 循环卷积
对于两个长度为 $n$ 的序列 $a,b$,定义$$c_k=\sum_{(i+j) mod n=k}a_ib_j$$
则称 $c$ 为 $a$ 与 $b$ 的循环卷积,记作 $c=a*b$
对于循环卷积与DFT有定理:
$F(a*b)=F(a)\cdot F(b)$其中$\cdot$表示逐点相乘
4. 两个关于单位根的推广:
$ω_n\quad ω_{2n}$
$(ω_{2n}^k)^2=ω_n^k$
$ω_{2n}^{n+k}=-ω_{2n}^k$
把图画出来推导即可

蝴蝶操作
按照 $n=2m$ 将 $A(x)$ 的项按次数的奇偶性分类

$A(x)=\sum_{0\le i<n}a_ix^i=\sum_{0\le i<m}a_{2i}x^{2i}+\sum_{0\le i<m}a_{2i+1}+x^{2i+1}$

$\qquad\;=\sum_{0\le i<m}a_{2i}x^{2i}+x\sum_{0\le i<m}a_{2i+1}x^{2i}$

定义:
$A_0(x)=\sum_{0\le i<m}a_{2i}x^i$

$A_1(x)=\sum_{0\le i<m}a_{2i+1}+x^i$

那么 $A_0(x),A_1(x)$ 都是次数不超过 $m$ 的多项式,并且有:

$A(x)=A_0(x^2)+xA_1(x^2)$

考虑单位根的性质:

$A(ω_n^k)=A_0((ω_n^k)^2)+w_n^kA_1((ω_n^k)^2)$

$\qquad\quad=A_0(ω_m^k)+w_n^kA_1(ω_m^k)$

$A(ω_n^{m+k})=A_0((ω_n^{m+k})^2)+ω_n^{m+k}A_1((ω_n^{m+k})^2)$

$\qquad\qquad=A_0(ω_m^k)-w_n^kA_1(ω_m^k)$

从上面可以看出,当我们取得了 $A_0(x),A_1(x)$ 在 $ω_m^0,ω_m^1,...,ω_m^{m-1}$ 的点权,那么就可以在 $O(n)$ 的时间内计算出 $A(x)$ 在 $ω_m^0,ω_m^1,...,ω_m^{m-1}$ 的点权.

计算 $A_0,A_1$ 的点值可以递归分治.

主定理可证复杂度为 $O(n\log n)$

IDFT的转化
蝴蝶操作是针对DFT的，如何在转成点值后变回系数？
只需要将FFT中的 $ω_n$ 全部变成 $ω_n^{-1}$ ,最后乘上 $\frac{1}{n}$ 即可,观察式子变换得到

例题

P3803 【模板】多项式乘法（FFT）

/*
雲璃猫猫が好きです
すべての生命よ,歌のように輝いています
截剣式、斬、断、破です!
*/
#include<bits/stdc++.h>
#include<bits/extc++.h>
#define INF 1e18
#define lb double
#define ls (id<<1)
#define rs (id<<1|1)
#define rep(i,l,r,k) for(int i=(l);i<=(r);i+=(k))
#define dep(i,r,l,k) for(int i=(r);i>=(l);i-=(k))
#define tep(ch,cr) for(auto ch:cr)
#define mk(a,b) make_pair(a,b)
#define me(a,b) memset(a,b,sizeof(a))
#define pb(x) push_back(x)
#define pr putchar
#define fi first
#define se second
#define wl while
#define max(a,b)((a)>(b)?(a):(b))
#define min(a,b)((a)<(b)?(a):(b))
using namespace std;
random_device rd;
unsigned int seed=rd();
mt19937 Rand(seed);
typedef pair<int,int> pii;
const int M=5e6+110,mod=1e9+7,Mod=998244353;
__gnu_pbds::gp_hash_table<string,int>ml;
inline int read(){int sum=0,k=1;char c=getchar();
	while(c>'9'||c<'0'){if(c=='-')k=-1;c=getchar();
	}while(c>='0'&&c<='9'){sum=sum*10+c-48;c=getchar();
	}return sum*k;
}inline void wr(int x){if(x<0) putchar('-'),x=-x;
	if(x>9) wr(x/10);return void(putchar(x%10+'0'));}
const lb pi=acos(-1.0);
struct Node{
	lb x,y;
	//A=x+yi
	//复数加减乘
	inline Node operator+(const Node &X)const{
		Node C;
		C.x=X.x+x;C.y=y+X.y;
		return C;
	}
	inline Node operator-(const Node &X)const{
		Node C;
		C.x=x-X.x;C.y=y-X.y;
		return C;
	}
	inline Node operator*(const Node &X)const{
		Node C;
		C.x=(x*X.x)-(y*X.y);
		C.y=(x*X.y)+(y*X.x);
		return C;
	}
}a[M],b[M];
inline void FFT(Node *res,int n,int id){
	if(n==1) return void();//边界
	Node a1[n],b1[n];//分割
	rep(i,0,n-1,2)
		a1[i>>1]=res[i],//奇偶分
		b1[i>>1]=res[i+1];
	FFT(a1,n>>1,id);FFT(b1,n>>1,id);//递归
	Node w=(Node){//增加量
		cos(2.0*pi/(lb)n),
		id*sin(2.0*pi/(lb)n)
	},x=(Node){1,0};
	rep(i,0,(n>>1)-1,1){//回推
		res[i]=a1[i]+x*b1[i];
		res[i+(n>>1)]=a1[i]-x*b1[i];
		x=x*w;
	}
}
int n,m,mx=1;
signed main(){
	n=read(),m=read();
	rep(i,0,n,1) cin>>a[i].x;rep(i,0,m,1) cin>>b[i].x;
	wl(mx<=m+n) mx<<=1;
	FFT(a,mx,1);FFT(b,mx,1);//变点值
	rep(i,0,mx,1) a[i]=a[i]*b[i];//转出新点值
	FFT(a,mx,-1);//转回系数
	rep(i,0,m+n,1)
		wr((int)(a[i].x/mx+0.5)),pr(' ');
	return 0;
}

P1919 【模板】高精度乘法 / A*B Problem 升级版
版子,考虑倒序存储后处理进位即可
因为可以看作$A(x)=a_1\cdot10^0+a_2\cdot10^1+a_3\cdot10^2+a_4\cdot10^3+...+a_n\cdot10^{n-1}$
两个数卷起来即可

/*
雲璃猫猫が好きです
すべての生命よ,歌のように輝いています
截剣式、斬、断、破です!
*/
#include<bits/stdc++.h>
#include<bits/extc++.h>
#define int long long
#define INF 1e18
#define lb   double
#define ls (id<<1)
#define rs (id<<1|1)
#define rep(i,l,r,k) for(int i=(l);i<=(r);i+=(k))
#define dep(i,r,l,k) for(int i=(r);i>=(l);i-=(k))
#define tep(ch,cr) for(auto ch:cr)
#define mk(a,b) make_pair(a,b)
#define me(a,b) memset(a,b,sizeof(a))
#define pb(x) push_back(x)
#define pr putchar
#define fi first
#define se second
#define wl while
#define max(a,b)((a)>(b)?(a):(b))
#define min(a,b)((a)<(b)?(a):(b))
using namespace std;
random_device rd;
unsigned int seed=rd();
mt19937 Rand(seed);
typedef pair<int,int> pii;
const int M=3e6+110,mod=1e9+7,Mod=998244353;
const lb pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679;
__gnu_pbds::gp_hash_table<string,int>ml;
inline int read(){int sum=0,k=1;char c=getchar();
	while(c>'9'||c<'0'){if(c=='-')k=-1;c=getchar();
	}while(c>='0'&&c<='9'){sum=sum*10+c-48;c=getchar();
	}return sum*k;
}inline void wr(int x){if(x<0) putchar('-'),x=-x;
	if(x>9) wr(x/10);return void(putchar(x%10+'0'));}
struct Node{
	lb x,y;//z=x+yi
	Node operator+(const Node &X)const{
		return (Node){x+X.x,y+X.y};
	}
	Node operator-(const Node &X)const{
		return (Node){x-X.x,y-X.y};
	}
	Node operator*(const Node &X)const{
		return (Node){x*X.x-y*X.y,x*X.y+y*X.x};
	}
}a[M],b[M];
inline void fft(Node *res,int n,int id){
	if(n==1) return void();
	Node a1[n],a2[n];
	rep(i,0,n-1,2) a1[i>>1]=res[i],a2[i>>1]=res[i+1];
	fft(a1,n>>1,id);fft(a2,n>>1,id);
	Node w=(Node){
		cos(2.0*pi/(lb)n),sin(2.0*pi/(lb)n)*id
	},x=(Node){1,0};
	rep(i,0,(n>>1)-1,1) res[i]=a1[i]+x*a2[i],
	res[i+(n>>1)]=a1[i]-x*a2[i],x=x*w;
}
int n,m,mx=1,ans[M],tot;
string sa,sb;
signed main(){
	cin>>sa>>sb;
	n=(int)sa.size();m=(int)sb.size();
	rep(i,0,n-1,1) a[i].x=sa[n-i-1]-'0';
	rep(i,0,m-1,1) b[i].x=sb[m-i-1]-'0';
	wl(mx<=m+n) mx<<=1;
	fft(a,mx,1),fft(b,mx,1);
	rep(i,0,mx,1) a[i]=a[i]*b[i];
	fft(a,mx,-1);
	tot=0;
	rep(i,0,mx,1){
		ans[i]+=(int)(a[i].x/mx+0.5);
		if(ans[i]>=10)
			ans[i+1]+=ans[i]/10,ans[i]%=10,mx+=(i==mx);
	}
	wl(!ans[mx]&&mx>=1) mx--;mx++;
	wl(--mx>=0) wr(ans[mx]);
	return 0;
}
/*

*/

NTT（快速数论变化）

原根
由费马小定理:
$若p\in prime 且 p|a 则有: a^{p-1}\equiv 1(mod p)$
对于 $m \in \mathbf{N}_{+}$ ，如果存在 $g \in \mathbf{Z}$ 且 $g \perp m$ 使得 $\delta_{m}(g)=\left|\mathbf{Z}_{m}^{*}\right|=\varphi(m)$ ，就称 $g$ 为模 $m$ 的原根（primitive root modulo m ）。其中， $\varphi(m)$ 是欧拉函数。
简单来说,就是就是一个最小整数 $g$ 使得其 $1~\varphi (m)$ 次幂在取模 $m$ 时互不相同.
优势
FFT使用了复数,要使用三角函数,所以会有精度的误差,ntt在常数大,精度误差小,特定模数才方便.
素数表(g为原根)：

/*
雲璃猫猫が好きです
すべての生命よ,歌のように輝いています
截剣式、斬、断、破です!
*/
#include<bits/stdc++.h>
#include<bits/extc++.h>
#define int long long
#define INF 1e18
#define lb long double
#define ls (id<<1)
#define rs (id<<1|1)
#define rep(i,l,r,k) for(int i=(l);i<=(r);i+=(k))
#define dep(i,r,l,k) for(int i=(r);i>=(l);i-=(k))
#define tep(ch,cr) for(auto ch:cr)
#define mk(a,b) make_pair(a,b)
#define me(a,b) memset(a,b,sizeof(a))
#define pb(x) push_back(x)
#define pr putchar
#define fi first
#define se second
#define wl while
#define max(a,b)((a)>(b)?(a):(b))
#define min(a,b)((a)<(b)?(a):(b))
using namespace std;
random_device rd;
unsigned int seed=rd();
mt19937 Rand(seed);
typedef pair<int,int> pii;
const int M=4e6+110,mod=998244353,Mod=998244353;
__gnu_pbds::gp_hash_table<string,int>ml;
inline int read(){int sum=0,k=1;char c=getchar();
	while(c>'9'||c<'0'){if(c=='-')k=-1;c=getchar();
	}while(c>='0'&&c<='9'){sum=sum*10+c-48;c=getchar();
	}return sum*k;
}inline void wr(int x){if(x<0) putchar('-'),x=-x;
	if(x>9) wr(x/10);return void(putchar(x%10+'0'));}
inline int qp(int x,int y,int res=1){
	wl(y){
		if(y&1) res=res*x%Mod;
		x=x*x%Mod,y>>=1;
	}return res;
}
int g=3,gi=332748118;
inline void ntt(int *res,int n,int id){
	if(n==1) return void();
	int a[n],b[n];
	rep(i,0,n,2) a[i>>1]=res[i],b[i>>1]=res[i+1];
	ntt(a,n>>1,id);ntt(b,n>>1,id);
	int w=1,
	x=((id==1)?qp(g,(Mod-1)/n):qp(gi,(Mod-1)/n));
	rep(i,0,(n>>1)-1,1)
		res[i]=(a[i]+w*b[i]+Mod)%Mod,
		res[i+(n>>1)]=(a[i]-w*b[i]+Mod)%Mod,
		w=w*x%Mod;
}
int mx=1,a[M],b[M];
signed main(){
	int n=read(),m=read(); 
	rep(i,0,n,1) a[i]=(read()+Mod)%Mod;
	rep(i,0,m,1) b[i]=(read()+Mod)%Mod;
	wl(mx<=m+n) mx<<=1;
	ntt(a,mx,1);ntt(b,mx,1);
	rep(i,0,mx,1) a[i]=a[i]*b[i]%Mod;
	ntt(a,mx,-1);
	int inv=qp(mx,Mod-2);
	rep(i,0,m+n,1) wr((((a[i]*inv)+Mod)%Mod+Mod)%Mod),pr(32);
	return 0;
}
/*

*/

posted @ 2025-11-14 14:49 rerecloud 阅读(26) 评论(14) 收藏举报

刷新页面返回顶部

入梦云端,可兮明虹

recloud-重返青云关注国家大事喵,注意身体健康喵,学好喵,云璃猫猫好看喵

FFT+NTT

前置补充：

fft

例题

NTT（快速数论变化）

公告

入梦云端,可兮明虹

recloud-重返青云 关注国家大事喵,注意身体健康喵,学好喵,云璃猫猫好看喵

FFT+NTT

前置补充：

fft

例题

NTT（快速数论变化）

公告

recloud-重返青云关注国家大事喵,注意身体健康喵,学好喵,云璃猫猫好看喵