粗略的学习笔记~
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) ax≡1(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) ap−1≡1(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)=p−1当 p p p为质数)
费马小定理还有个用途就是当幂的指数特别大时,需要对 p p p取模,我们可以将指数对 p p p取模然后求(因为根据费马小定理,将模出来等于一的可以先模掉)。
离散对数:满足 a x ≡ k ( m o d m ) a^x\equiv k(\rm mod\ m) ax≡k(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表示):
- g i m o d m ( i ∈ [ 0 , m − 1 ] ) g^i\rm mod\ m(i\in[0,m-1]) gimod m(i∈[0,m−1])可以将 [ 0 , m − 1 ] [0,m-1] [0,m−1]内的数字全部取到。
- ( g m − 1 p i m o d m ) ≠ 1 (g^{\frac{m-1}{p_i}}\rm mod\ m)≠1 (gpim−1mod 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
转换
通过离散对数的定义,我们可以将原来FFT单位根的方程 x n = 1 x^n=1 xn=1写为 x n ≡ 1 ( m o d m ) x^n\equiv 1(\rm mod\ m) xn≡1(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(m−1)(k∈Z,), 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} {x≡a1 (mod M)x≡a2 (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×M≡a3−x(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 < ( m 1 × m 2 × m 3 ) ) ans\ \rm mod\ P(ans<(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;
}
一些例题
HDU5829 Rikka with Subset - 题解
End
目前初学,还有很多不懂,望大佬指点错误!