文档章节

【模板篇】NTT和三模数NTT

o
 osc_x4h57ch8
发布于 2018/04/24 08:29
字数 1393
阅读 0
收藏 0
ntt

之前写过FFT的笔记. 我们知道FFT是在复数域上进行的变换. 而且经过数学家的证明, DFT是复数域上唯一满足循环卷积性质的变换.

而我们在OI中, 经常遇到对xxxx取模的题目, 这就启发我们可不可以在模运算的意义下找一个这样的变换. 然后我们发现有个神奇的东西, 原根$g$, 这东西在模意义下相当于单位复根$-e^{\frac{2\pi i}{n}}$.

所以我们预处理一下$g$的幂和逆元, 然后改一下fft的代码就出现了快速数论变换ntt 懒得写了 直接上代码:

void getwn(){ //预处理原根的幂和逆元
	int x=qpow(3,p-2);
	for(int i=0;i<20;++i){
		wn[i]=qpow(3,(p-1)/(1<<i));
		inv[i]=qpow(x,(p-1)/(1<<i));
	}
}
void ntt(int *y,bool f){ rev(y); //翻转代码和fft无异
	for(int m=2,id=1;m<=n;m<<=1,++id){ //id用来记录转到第几下了
		for(int k=0;k<n;k+=m){
			int w=1,wm=f?wn[id]:inv[id]; //如果是dft就用幂, idft就用幂的逆元
			for(int j=0;j<m/2;++j){
                                //这里跟fft一样, 不过要对p取模
				int u=y[k+j]%p,t=1ll*w*y[k+j+m/2]%p;
				y[k+j]=u+t; if(y[k+j]>p) y[k+j]-=p;
				y[k+j+m/2]=u-t; if(y[k+j+m/2]<0) y[k+j+m/2]+=p;
				w=1ll*w*wm%p;
			}
		}
	}
	if(!f){
		int x=qpow(n,p-2);
		for(int i=0;i<n;++i)
			y[i]=1ll*y[i]*x%p;
	}
}

好像差不多呢~ 不过这样就要求我们找一个原根好求的数. 比如著名的uoj数: 998244353 还有1004535809和469762049等, 这三个数原根都是3~ 好像因为当时一看到模数不是1e9+7一般就会想到ntt, vfk为了防止这一点, 模数统一采用998244353, 现在看看收效不错.

不过 有些丧心病狂的人就是要用1e9+7作为ntt的模数, 甚至还出现了可以不模质数的情况! 那我们怎么解决任意模数ntt呢? 我们可以采用拆系数ntt或者三模数ntt. 这里介绍一下三模数ntt. 对于一般的数据范围, $n\leq10^5, a_i\leq10^9$, 这样可能会到$10^510^{9^2}=10^{23}$级别. 所以我们可以找三个乘积$>10^{23}$的ntt-friendly的数, 然后分别ntt再想办法合并. 我们假如答案是ans, 那我们做三次ntt后就能得到如下三个柿子. $$ \left{\begin{matrix} ans\equiv a_1(\mod m_1)\ ans\equiv a_2(\mod m_2)\ ans\equiv a_3(\mod m_3) \end{matrix}\right. $$ 我们把前两个柿子通过中国剩余定理合并, 就可以得到 $$ \left{\begin{matrix} ans\equiv A(\mod M)\
ans\equiv a_3(\mod m_3) \end{matrix}\right. $$ 其中, $M=m_1
m_2$ 这样我们设$ans=kM+A$, $$ kM+A\equiv a_3(\mod m_3) \ k=(a_3-A)*M^{-1} (\mod m_3) $$ 这样我们求出$k$然后代回到$ans=kM+A$就可以求对任意模数取模的结果了.

中国剩余定理合并的时候直接乘是可以爆long long的, 所以我们要用到$O(1)$快速乘~

下面上一波代码: luogu4245 【模板】MTT <span class="heimu" title="你自己还知道啊~">哎呀觉得自己码风有点丑啊qwq</span>

#include <cstdio>
#include <cstring>
#include <algorithm>
typedef long long LL;
const int N=600020,p0=469762049,p1=998244353,p2=1004535809;
const LL M=1ll*p0*p1;
int wn[20],nw[20],rev[N],n,lg,p;
int qpow(int a,int b,int p,int s=1){
    for(;b;b>>=1,a=1ll*a*a%p)
        if(b&1) s=1ll*s*a%p;
    return s;
}
LL mul(LL a,LL b,LL p){ a%=p; b%=p;
    return (a*b-(LL)((long double)a*b/p)*p+p)%p;
}
void calcw(int p){
    int x=qpow(3,p-2,p);
    for(int i=0;i<20;++i){
        wn[i]=qpow(3,(p-1)/(1<<i),p);
        nw[i]=qpow(x,(p-1)/(1<<i),p);
    }
}
void init(){
    for(int i=0;i<n;++i)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<lg);
}
void ntt(int *y,bool f,int p){ calcw(p);
    for(int i=0;i<n;++i) if(i<rev[i]) std::swap(y[i],y[rev[i]]);
    for(int m=2,id=1;m<=n;m<<=1,++id){
        for(int k=0;k<n;k+=m){
            int w=1,wm=f?wn[id]:nw[id];
            for(int j=0;j<m>>1;++j){
                int &a=y[k+j]; int &b=y[k+j+m/2];
                int u=a%p,t=1ll*w*b%p;
                a=u+t; if(a>p) a-=p;
                b=u-t; if(b<0) b+=p;
                w=1ll*w*wm%p;
            }
        }
    } int x=qpow(n,p-2,p);
    if(!f) for(int i=0;i<n;++i) y[i]=1ll*y[i]*x%p;
}
char c1[N],c2[N]; int a[N],b[N],c[N],d[N],ans[3][N];
int main(){
    int l1,l2; scanf("%d%d%d",&l1,&l2,&p);
    for(int i=0;i<=l1;++i) scanf("%d",&a[i]),a[i]%=p;
    for(int i=0;i<=l2;++i) scanf("%d",&b[i]),b[i]%=p;
    for(n=1;n<l1||n<l2;n<<=1,++lg); n<<=1; init();
    std::copy(a,a+n,c); std::copy(b,b+n,d);
    ntt(c,1,p0); ntt(d,1,p0);
    for(int i=0;i<n;++i) ans[0][i]=1ll*c[i]*d[i]%p0;
    std::copy(a,a+n,c); std::copy(b,b+n,d);
    ntt(c,1,p1); ntt(d,1,p1);
    for(int i=0;i<n;++i) ans[1][i]=1ll*c[i]*d[i]%p1;
    std::copy(a,a+n,c); std::copy(b,b+n,d);
    ntt(c,1,p2); ntt(d,1,p2);	
    for(int i=0;i<n;++i) ans[2][i]=1ll*c[i]*d[i]%p2;
    ntt(ans[0],0,p0); ntt(ans[1],0,p1); ntt(ans[2],0,p2);
    for(int i=0;i<n;++i){
        LL A=mul(1ll*ans[0][i]*p1%M,qpow(p1%p0,p0-2,p0),M)
            +mul(1ll*ans[1][i]*p0%M,qpow(p0%p1,p1-2,p1),M);
        if(A>M) A-=M;
        LL k=((ans[2][i]-A)%p2+p2)%p2*qpow(M%p2,p2-2,p2)%p2;
        a[i]=1ll*(k%p)*(M%p)%p+A%p;
        if(a[i]>p) a[i]-=p;
    }
    for(int i=0;i<=l1+l2;++i) printf("%d ",a[i]);
}
o
粉丝 0
博文 500
码字总数 0
作品 0
私信 提问
加载中
请先登录后再评论。

暂无文章

pycurl libcurl link-time ssl backend (nss)

pip uninstall pycurlecho 'pycurl==7.19.5.1 --global-option="--with-nss"' > requires.pypip install -r requires.py...

小红手
19分钟前
17
0
计算机网络性能衡量

1、速率 单位时间(s)内传输信息(bit)量 单位:KB/s, MB/s, Gb/s K = 10^3 ,M = 10^6, G=10^9 一般表示的是理想的传输速率 2、带宽 计算机网络中的带宽和通信等领域的带宽概念不一样,计算机网...

osc_np3y0rbq
19分钟前
3
0
互联网掀起农家乐,巨头上演AI掘金战

配图来自Canva **前有网易、阿里AI养猪,后有腾讯AI养鹅,互联网大佬们纷纷玩起了“农家乐”,互联网的生意在尖端技术的引领之下频频跨界,巨头之间的较量也从线上延伸至线下。**自古“民以食...

osc_5cok9i01
20分钟前
5
0
原来!我在4年前就开始体验雾游戏了!

前有云游戏后有雾游戏,游戏的方式看来起来越来越多种多样。那么“震撼业界”的雾游戏到底是什么来头?它依靠什么改变游戏界?它的原理又是什么? 本月月初,著名的日本游戏杂志《Fami通》表...

osc_j34n26zn
21分钟前
5
0
活动预告|田溯宁与你相约GSMA Thrive·万物生晖,分享5G风口下的创新与投资洞察

在万物互联的时代背景下,5G+AI+IoT的技术变革与融合,正在引发一场深刻的全产业创新与变革。5G技术创新、行业应用及投资机遇已成为科技行业所瞩目的焦点。 6月30日,宽带资本董事长田溯宁将...

osc_0qnrwmy3
23分钟前
7
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部