NTT学习笔记

粗略的学习笔记~

NTT(快速数论变换)

  • 解决多项式卷积在取模意义下的值,并且不使用复数(浮点数)而是使用整数(在取模的意义下)来减少精度误差。

其实FFT预处理单位根和将较大数据拆成两部分也可以减少精度误差,但是还是没有NTT的精度高,而且不方便取模。


这里假设你会了FFT。
不会的话看这里FFT

FFT之所以精度误差大,就是在于单位根为无理数 e e e的多少次方,所以必须用浮点数,那么由于计算机的特性,误差会很大。

而这里是在取模的一个前提下,那么如果我们找到了一个新的且为整数的单位根(这里是在取模意义下),我们记其为 G G G,那么用这个根去替代原来的 e e e,进行运算。由于所有运算就会变成整数运算,所以精度误差大大降低。

如何找到这么一个神奇 G G G满足原来的单位根的性质呢?


原根

原根:对于一个数字 m m m,当 g c d ( a , m ) = 1 gcd(a,m)=1 gcd(a,m)=1 a a a m m m的最大公约数为1,也就是互质)时,我们记满足 a x ≡ 1 ( m o d   m ) a^x\equiv 1(\rm mod\ m) ax1(mod m)式子的最小 x x x o r d m ( a ) ord_m(a) ordm(a),然后当 o r d m ( a ) = φ ( m ) ord_m(a)=\varphi(m) ordm(a)=φ(m)时称 a a a m m m的原根。(不一定所有的数字都有原根)。

扩展:
费马小定理: a p − 1 ≡ 1 ( m o d   p ) a^{p-1}\equiv 1(\rm mod\ p) ap11(mod p) p p p为素数。(这个可以用来求逆元,也可写为 a φ ( p ) ≡ 1 ( m o d   p ) a^{\varphi(p)}\equiv 1(\rm mod\ p) aφ(p)1(mod p)因为 φ ( p ) = p − 1 \varphi(p)=p-1 φ(p)=p1 p p p为质数)
费马小定理还有个用途就是当幂的指数特别大时,需要对 p p p取模,我们可以将指数对 p p p取模然后求(因为根据费马小定理,将模出来等于一的可以先模掉)。
离散对数:满足 a x ≡ k ( m o d   m ) a^x\equiv k(\rm mod\ m) axk(mod m)我们将 x x x称为 a a a在模 m m m的意义下的离散对数,我们可以简记为 x = i n d a ( k ) x=ind_a(k) x=inda(k)

性质(我们将原根一般用 g g g表示):

  1. g i m o d   m ( i ∈ [ 0 , m − 1 ] ) g^i\rm mod\ m(i\in[0,m-1]) gimod m(i[0,m1])可以将 [ 0 , m − 1 ] [0,m-1] [0,m1]内的数字全部取到。
  2. ( g m − 1 p i m o d   m ) ≠ 1 (g^{\frac{m-1}{p_i}}\rm mod\ m)≠1 (gpim1mod m)̸=1,其中 p i p_i pi m m m的素数因子。

  • 原根的求法:
    因为一般原根很小,我们采用暴力枚举的方法。
    对于一个数字 m m m我们求它的原根,根据性质2,我们先线性筛出 m m m的所有素数因子,然后从 2 2 2开始枚举原根 g g g,如果对于一个枚举的 g g g,满足对于 m m m的所有质因数 p i p_i pi,满足性质2,那么这个数就是原根,然后break掉就好了。

要求原根的一个题目BZOJ3992

原根的详细讲解-hdxrie的友链


转换

通过离散对数的定义,我们可以将原来FFT单位根的方程 x n = 1 x^n=1 xn=1写为 x n ≡ 1 ( m o d   m ) x^n\equiv 1(\rm mod\ m) xn1(mod m),那么我们求的新的整数 x x x不就是原根了吗?

所以我们用 g g g来代替了单位根,那么 我 们 仍 旧 学 F F T 一 样 记 ω n k = g k ( m − 1 ) n ( k ∈ Z , ) 我们仍旧学FFT一样记\omega_n^k=g^{\frac{k(m-1)}{n}}(k\in Z,) FFTωnk=gnk(m1)(kZ,) m m m为模数,这样就代替了单位复数根,但是只有当 m m m为素数时,原根才一定存在,所以我们取模数一般取素数,而且这种质数一般取费马质数(NTT质数)。

一些常用的取模素数以及原根:IN There


但是模数并不总是素数,而是一些普通的数字,那么我们怎么办呢?

那么对于模数,我们可以用几个乘积大于最大卷积的系数的模数,分别求取取模后的卷积,然后用中国剩余定理合并即可。

不会中国剩余定理的->百度

但是有时卷积的系数会非常之大,导致你不能直接合并,否则爆long long,所以我们要使用一个小技巧:

假如我们有三个模数 m 1 , m 2 , m 3 m_1,m_2,m_3 m1,m2,m3,其中 m 1 × m 2 m_1\times m_2 m1×m2不会爆long long,但是 m 1 × m 2 × m 3 m_1\times m_2\times m_3 m1×m2×m3会爆,那么我们先直接合并前两个得到 M = m 1 × m 2 M=m_1\times m_2 M=m1×m2的卷积,然后我们就会剩下以下这个方程式:

{ x ≡ a 1   ( m o d   M ) x ≡ a 2   ( m o d   m 3 ) \begin{cases} x \equiv a_1\ ({mod}\ M) \\ x\equiv a_2\ ({mod}\ m_3) \end{cases} {xa1 (mod M)xa2 (mod m3)

那么转换一下就可以写成:

{ x = a 1 + k 1 × M x = a 3 + k 3 × m 3 \begin{cases} x = a_1+k_1\times M \\ x= a_3+k_3\times m_3 \end{cases} {x=a1+k1×Mx=a3+k3×m3

然后可以写成:

x = k 1 × M + a 1 = k 3 × m 3 + a 3 x=k_1\times M+a_1=k_3\times m_3+a_3 x=k1×M+a1=k3×m3+a3

那么:

k 1 × M ≡ a 3 − x ( m o d   m 3 ) k_1\times M\equiv a_3-x(\rm mod\ m_3) k1×Ma3x(mod m3)

k 3 × m 3 k_3\times m_3 k3×m3 m 3 m_3 m3模掉了。

所以我们解出 k 1 k_1 k1,然后带入原式即可算出答案 x x x,那么这个问题就解决了。

最后我们求出的是在 m o d ( m 1 × m 2 × m 3 ) \rm mod(m_1\times m_2\times m_3) mod(m1×m2×m3)的意义下的卷积,如果原来让你求 a n s   m o d   P ( a n s &lt; ( m 1 × m 2 × m 3 ) ) ans\ \rm mod\ P(ans&lt;(m_1\times m_2\times m_3)) ans mod P(ans<(m1×m2×m3)),那么你再将 a n s ans ans P P P取模即可。


模板代码:

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;

const int M=3e6+10;
const ll p1=469762049ll,p2=998244353ll,p3=1004535809ll,g=3;
const ll Mod=468937312667959297ll;
//三个乘积大于模数且大于最大卷积的系数P*(n-1)^2,原根均为3
int n,m,p;
ll t1[M],t2[M],t3[M],t4[M],ans[3][M];
int R[M];
int fpow(int a,int b,int mod){
	int ans=1;
	for(;b;b>>=1,a=(1ll*a*a)%mod){
		if(b&1) ans=(1ll*ans*a)%mod;
	}
	return ans;
}
ll mul(ll a,ll b,ll mod){
    a%=mod,b%=mod;
    return ((a*b-(ll)((ll)((long double)a/mod*b+1e-3)*mod))%mod+mod)%mod;
}

void DNT(ll *a,int n,int f,int mod){
	for(int i=0;i<n;i++) if(i<R[i]) swap(a[i],a[R[i]]);
	for(int i=2;i<=n;i<<=1){
		int now=i>>1;
		int wn=fpow(g,(mod-1)/i,mod);
		if(f==-1) wn=fpow(wn,mod-2,mod);//逆变换时,逆元,或者采用如下的交换法
		for(int j=0;j<n;j+=i){
			int w=1,x,y;
			for(int k=j;k<j+now;k++,w=1ll*w*wn%mod){
				x=a[k],y=1ll*w*a[k+now]%mod;
				a[k]=(x+y)%mod;a[k+now]=(x-y+mod)%mod;
			}
		}
	}
	if(f==-1){
		int inv=fpow(n,mod-2,mod);
		for(int i=0;i<=n;i++){
			a[i]=1ll*a[i]*inv%mod;
		}
//		for(int i=1;i<=n/2;i++) swap(a[i],a[n-i]);//交换法,因为g^(-1)=g^(n-1),所以可以先正变换,翻转就为逆变换。
	}
}

void NTT(ll *a,ll *b,int k,int mod){
	DNT(a,k,1,mod),DNT(b,k,1,mod);
	for(int i=0;i<=k;i++) a[i]=1ll*a[i]*b[i]%mod;
//	DNT(a,k,-1,mod);
}
int k;
void mcpy(int d){
	for(int i=0;i<=k;i++) ans[d][i]=t3[i];
	if(d==2) return;
	memset(t3,0,sizeof(t3));memset(t4,0,sizeof(t4));
	for(int i=0;i<=n;i++) t3[i]=t1[i];
	for(int i=0;i<=m;i++) t4[i]=t2[i];
}

void REV(){
	DNT(ans[0],k,-1,p1);
	DNT(ans[1],k,-1,p2);
	DNT(ans[2],k,-1,p3);
}

int lg;
int main(){
	scanf("%d%d%d",&n,&m,&p);
	for(int i=0;i<=n;i++) scanf("%lld",&t1[i]),t3[i]=t1[i];
	for(int i=0;i<=m;i++) scanf("%lld",&t2[i]),t4[i]=t2[i];
	k=1;while(k<=n+m)k<<=1,++lg;
	for(int i=0;i<k;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(lg-1));
	NTT(t3,t4,k,p1);
	mcpy(0);
	NTT(t3,t4,k,p2);
	mcpy(1);
	NTT(t3,t4,k,p3);
	mcpy(2);
	REV();
	for(int i=0;i<=n+m;i++){//中国剩余定理合并+后面推导出的公式
		ll x=((mul(1ll*ans[0][i]*p2%Mod,fpow(p2%p1,p1-2,p1),Mod))+
			  (mul(1ll*ans[1][i]*p1%Mod,fpow(p1%p2,p2-2,p2),Mod)))%Mod;
		ll y=((((ans[2][i]-x)%p3+p3)%p3)*fpow(Mod%p3,p3-2,p3))%p3;
//		printf("%lld %lld %lld ",x,y,ans[2][i]);
		printf("%lld ",(1ll*(y%p)*(Mod%p)%p+x%p)%p);
	}
	return 0;
}

一些例题

【模板】任意模数NTT

HDU5829 Rikka with Subset - 题解

[SDOI2015]序列统计


End

目前初学,还有很多不懂,望大佬指点错误!


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

VictoryCzt

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值