zoukankan      html  css  js  c++  java
  • 浅谈$NTT$

    (NTT),快速数论变换,可以理解为带模数的FFT


    原根 & 阶

    先来补一点数论。(这里讲的应该很少,都是针对(ntt)胡的,具体的话可以去看《初等数论》那本小黄书)。

    阶(指数)

    如果(m > 1, (a,m) = 1),那么必有整数(d),使得下面这个柿子成立

    [a^d equiv 1 (mod m) ]

    我们称令这个柿子成立的最小的(d)(a mod m)(或指数),记作(delta_m(a))

    一些性质

    • ((1).) 因为(a,m)互质,所以根据欧拉定理(a^{varphi(m)} equiv 1 (mod m)),立即得到(delta_m(a) le varphi(m))

    • $(2). $ 若(d_0)也满足(a^{d_0} equiv 1 (mod m)),那么必有(delta_m(a) mid d_0),这个显然。

    • $(3). $ 由((1))((2))可以立即推出(delta_m(a) mid varphi(m))

    • (cdots)


    原根

    (m > 1, (a,m) = 1),且(delta_m(a) = varphi(m)),那么则称(a)(mod m)的一个原根

    特别(m)是质数(p)的时候若(delta_p(a) = p-1),则(a)是原根。

    性质

    (g)(m)的原根,那么(g, g^2,g^3,cdots,g^{varphi(m)})一定是在(mod m)意义下与(m)互质的(varphi(m))个数的一个排列。这个根据定义也显然。

    特别的当(m)是一个质数(p)的时候,有(g, g^2, g^3, cdots, g^{p-1})(mod p)意义下是(1...p-1)(p-1)个数的一个排列。

    偷偷告诉你,原根一般都很小,默认<=100就好了

    那对于一个素数,如何找原根呢?

    直接暴力,枚举(g = 1...100),然后(check),根据阶性质((1))我们在枚举(d)的时候可以枚举(p-1)的因子,注意不要枚举到(p-1)

    #include <bits/stdc++.h>
    using namespace std;
    int p;
    int fpow(int a,int b,int mod=p){
    	int ret=1; for (;b;b>>=1,a=1ll*a*a%mod) if (b&1) ret=1ll*ret*a%mod;
    	return ret;
    }
    int chk(int g){
    	for (int i=2;i*i<=p-1;i++) // !!从2开始
    		if ((p-1)%i==0&&(fpow(g,i)==1||fpow(g,(p-1)/i)==1)) return 0;
    	return 1;
    }
    int main(){
    	scanf("%d",&p);
    	for (int i=1;i<=100;i++)
    		if (chk(i)){printf("G: %d
    ",i); return 0;}
    	return 1;
    	return 0;
    }
    

    安利一个巨佬原根表


    NTT

    P.S. 以下的原根都是在一个形如(k2^t + 1)的质数(p)的同余系下定义的。((t)足够大,一般(t ge 20)

    (g)(mod p)的原根,那么(g, g^2, cdots, g^{p-1})(mod p)意义下是(p-1)个不同的数。


    回去看一下FFT,随便交一发就(T)了,你看这个FFT,他就是逊啦!

    发现(FFT)还有很多不足之处

    • 复数运算,浮点数乘法,自带无穷大常数

    • 多次调用三角函数,常数大,还掉精度

    • 复数乘法异常难背,还要手算一遍

    • (cdots)

    总之我们得想办法把频繁的复数运算和单位根去掉。

    在颓(FFT)的时候我们用到了单位根的几个性质

      1. (omega_{n}^{k} = (omega_{n}^{1})^k)
      1. (omega_{n}^{k} = omega_{n}^{k mod n})
      1. (omega_{n}^{k + n/2} = -omega_{n}^{k})
      1. (omega_{2n}^{2k} = omega_{n}^{k})

    当然根据单位根的定义,(omega_{n}^{0}, omega_{n}^{1}, cdots, omega_{n}^{n-1})必须是(n-1)个不同的数。

    结论就是我们在(mod p)意义下,我们用原根(g)((p-1)/n)次方代替(omega_{n}^{1})就对了,即(omega_{n}^1 = g^{(p-1)/n})

    我们来证一下性质(一下运算均是在(mod p)意义下的)。(1. 就是定义,即(omega_{n}^{k} = (omega_{n}^1)^k = g^{k(p-1)/n}),我们不用证)。

      1. (omega_{n}^{k} = omega_{n}^{k \% n})

        因为(p)是质数,有小费马(a^{p-1} equiv 1 (mod p)),立刻得到(g^{(p-1)k/n} = g^{(p-1)k/n \% (p-1)})

        这时候并不能直接取模,因为(k/n)不一定是整数。我们令({ x })表示(x)的小数部分。

        ((p-1)k/n \% (p-1) = (p-1){ k/n } = (p-1)(k \% n) /n),代上去即证。

      1. (omega_{n}^{k+n/2} = -omega_{n}^{k})

        根据定义,有(omega_{n}^{k+n/2} = omega_{n}^{k} cdot omega_{n}^{n/2})

        由于((omega_{n}^{n/2})^2 = omega_{n}^{n} = 1),所以(omega_{n}^{n/2})(mod p)意义下只可能是(1)(-1),又因为(omega_{n}^{n/2} ot= omega_{n}^{0}),所以(omega_{n}^{n/2} equiv -1 ( mod p))。代上去就好了。

      1. (omega_{2n}^{2k} = omega_{n}^{k})

        这个应该是最好证的。。。

        (g^{(p-1)/n})代上去,(g^{(p-1)2k/2n} = g^{(p-1)k/n})。没了。

    综上所述,在(mod p)意义下用(g^{(p-1)/n})代替(omega_{n}^{1})一点毛病都没有,在(IDFT)的时候用((omega_{n}^{1})^{-1})代替(omega_{n}^{-1}),算逆元就好了。

    所以再来看一下板子题 才发现这里系数都<10

    我们取一个信仰膜数(p = 998244353),然后在(mod p)意义下做(NTT),没了。

    #include <bits/stdc++.h>
    using namespace std;
    inline int read(){
    	char ch=getchar();int x=0,f=1;
    	while (!isdigit(ch)){f=ch=='-'?-f:f;ch=getchar();}
    	while (isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    	return x*f;
    }
    const int N=4e6+10,P=998244353,G=3,IG=332748118;
    int fpow(int a,int b,int mod=P){
    	int ret=1; for (;b;b>>=1,a=1ll*a*a%P) if (b&1) ret=1ll*a*ret%P;
    	return ret;
    }
    #define add(x,y) ((x)+(y)>=P?(x)+(y)-P:(x)+(y))
    #define sub(x,y) ((x)-(y)<0?(x)-(y)+P:(x)-(y))
    int rev[N];
    void ntt(int *f,int n,int flg){
    	for (int i=0;i<n;i++) if (i<rev[i]) swap(f[i],f[rev[i]]);
    	for (int len=2;len<=n;len<<=1){
    		int w1=fpow(flg==1?G:IG,(P-1)/len); // w^{-1} = (w^1)^{-1}
    		for (int i=0;i<n;i+=len){
    			int w=1;
    			for (int j=i;j<(len>>1)+i;j++){
    				int tmp=1ll*w*f[j+(len>>1)]%P;
    				f[j+(len>>1)]=sub(f[j],tmp);
    				f[j]=add(f[j],tmp);
    				w=1ll*w*w1%P;
    			}
    		}
    	}
    }
    int f[N],g[N];
    int main(){
    	int n=read(),m=read();
    	for (int i=0;i<=n;i++) f[i]=read();
    	for (int i=0;i<=m;i++) g[i]=read();
    	int limit=1; while (limit<=n+m)limit<<=1;
    	for (int i=0;i<limit;i++) rev[i]=rev[i>>1]>>1|((i&1)?limit>>1:0);
    	ntt(f,limit,1),ntt(g,limit,1);
    	for (int i=0;i<limit;i++) f[i]=1ll*f[i]*g[i]%P;
    	ntt(f,limit,-1); int ivn=fpow(limit,P-2);
    	for (int i=0;i<=m+n;i++) printf("%d ",1ll*f[i]*ivn%P);
    	return 0;
    }
    

    提交记录,总共(1.7s),nice!

  • 相关阅读:
    SqlServer:创建索引
    SqlServer:使用视图 View
    SqlServer:修改和删除数据
    网络管理:SNMPv1
    《剑指 Offer》学习记录:题 9:用两个栈实现队列
    《剑指 Offer》学习记录:题 27:二叉树的镜像
    团队冲刺9
    团队冲刺8
    团队冲刺7
    团队冲刺6
  • 原文地址:https://www.cnblogs.com/wxq1229/p/12235505.html
Copyright © 2011-2022 走看看