zoukankan      html  css  js  c++  java
  • 【数学】【多项式】快速数论变换(NTT)

    引入

    求两个多项式的卷积

    Description

    给定两个多项式 (Fleft(x ight), Gleft(x ight)) 的系数表示法,求两个多项式的卷积。

    如:

    [Fleft(x ight) = 2x + 1 ]

    [Gleft(x ight) = x^2 + 2x + 1 ]

    得到的结果即为:

    [Fleft(x ight)Gleft(x ight) = left(2x + 1 ight)left(x^2 + 2x + 1 ight) = 2x^3 + 5x^2 + 4x + 1 ]

    直接由系数计算的复杂度是 (mathcal Oleft(n ^ 2 ight)) 级别的,而我们发现在点值表示下,计算卷积的复杂度是 (mathcal Oleft(n ight)) 的,因此考虑把系数表示转化为点值表示然后再做乘法。

    这样为了实现这个转化,我们引入了快速傅里叶变换(FFT)。通过把多项式按奇偶次分成两个关于 (x^2) 的多项式,并不断递归下去实现这个过程,并通过求得的一对 (x) 两个子多项式分别的值,一相加,一相减得到原多项式中 (x) 满足正负配对的一对点值。这个过程中为了满足 (x) 取值的正负配对,我们把 (x) 的取值扩充到复数域,取 (n) 次单位根。这样实现了由系数向点值的转化 DFT。

    在由点值转化回系数的时候,我们从矩阵运算的角度,意识到 IDFT 矩阵应为 DFT 矩阵的逆矩阵,而 DFT 矩阵的逆矩阵又只需要对每一项取其倒数,并乘上 (n^{-1})。因此我们在 DFT 上稍加修改就实现了 IDFT。

    FFT 存在的一点不优美之处

    根据 FFT 的过程,我们可知因为复数运算,FFT 对运算的精度要求很高,这也就影响了算法的时空开销。

    而这个瓶颈似乎就是出在使用了复数运算上。

    考虑是否能在 (mathbf R) 甚至是 (mathbf Z) 范围内找到一种能够替代 (n) 次单位根的某种神奇的元素。

    它存在吗?当然存在。这就是快速数论变换(NTT)用到的原根

    一点数论知识和概念

    若对于 (a, p),存在正整数 (l) 使得 (a^l equiv 1 left(mod p ight)),则把满足这个条件的最小的正整数 (l) 称为 (a) 在模 (p) 下的

    举例:

    (a = 2, p = 5) 时,因为 (2 ^ 4 = 16 equiv 1 left(mod 5 ight)),可以证明,(4) 是满足该条件的最小正整数。因此 (4) 就是 (2) 在模 (5) 意义下的阶。

    给定 (a, p)(l) 的方法是一个典型的高次同余方程,可以使用 BSGS 算法求解。

    原根

    (g) 在模 (p) 意义下的阶是 (varphileft(p ight)),其中 (varphi) 是欧拉函数,则 (g) 称为模 (p) 意义下的原根

    如何求一个原根呢?

    关于原根,由如下的结论:

    (gcdleft(g, p ight) = 1),设 (p_1, p_2,cdots,p_k)(varphileft(p ight)) 所有的质因数,则 (g) 是模 (p) 的一个原根,当且仅当对于任意的 (p_i),都有 (g^{frac{varphileft(p ight)}{p_i}} otequiv 1 left(mod p ight))

    可以证明,(p) 的最小原根在 (p ^ {0.25}) 级别。因此枚举原根的时间一般都是可以接受的。

    快速数论变换(NTT)

    原根的性质

    根据阶和原根的定义,我们可以发现原根一个很优美的性质,那就是 ({g, g ^ 2, g ^ 3, cdots, g ^ {varphileft(p ight)}}) 在模 (p) 下互不相同。次数若继续增大,则在模 (p) 下就会出现循环。这似乎与单位根的性质非常相似。

    特殊的,若 (p) 是质数,(varphileft(p ight) = p - 1)。若 (p - 1) 的分解中含有较多的 (2) 这个因子,就可以把它的原根作拿来代替单位根了!

    这个模数 (p) 一般取 (998244353left(2^{23} imes 7 imes 17 + 1 ight)) 即可。此时原根 (g = 3)

    只要多项式的长度不要求在 (8 imes 10^6) 以上,这个模数和原根基本就够用。

    实现

    为了实现 (n) 次单位根乘 (n) 次方会出现循环的性质,联系上文提到的原根的性质,直接用 (g^{frac{varphileft(p ight)}{n}})(p) 为质数时为 (g^{frac{p - 1}{n}}))来替换 (n) 次单位根即可。

    这样,我们在 FFT 上稍作修改,就实现了 NTT。

    NTT 实现 IDFT 的时候,存在两种效果等价的方式,一种是直接对转换后的序列 reverse 一下,另一种是取原根的逆元。我选择了后者,在 (p = 998244353) 时这个数为 (332748118)

    代码实现

    void NTT(LL a[], int len, int type) {
    	for(int i = 0; i < len; ++i) if(i < rev[i]) swap(a[i], a[rev[i]]);
    	for(int h = 2; h <= len; h <<= 1) {
    		LL gn = qpow(type == 1 ? g : invg , ((Mod - 1) / (LL)h));
    		for(int j = 0; j < len; j += h) {
    			LL gk = 1;
    			for(int k = j; k < j + h / 2; ++k) {
    				LL e = a[k], o = gk * a[k + h / 2] % Mod;
    				a[k] = (e + o) % Mod; a[k + h / 2] = ((e - o) % Mod + Mod) % Mod;
    				gk = gk * gn % Mod;
    			}
    		}
    	}
    	if(type == -1) {
    		LL inv = qpow(len, Mod - 2);
    		for(int i = 0; i < len; ++i) a[i] = a[i] * inv % Mod;
    	}
    }
    

    其他

    完整代码

    洛谷 P3803 多项式乘法(NTT)

    #include <bits/stdc++.h>
    #define LL long long
    
    template <typename Temp> inline void read(Temp & res) {
    	Temp fh = 1; res = 0; char ch = getchar();
    	for(; !isdigit(ch); ch = getchar()) if(ch == '-') fh = -1;
    	for(; isdigit(ch); ch = getchar()) res = (res << 3) + (res << 1) + (ch ^ '0');
    	res = res * fh;
    }
    
    using namespace std;
    
    const int Maxn = 2097200;
    const LL Mod = 998244353, g = 3, invg = 332748118;
    
    inline LL qpow(LL A, LL P) {
    	LL res = 1;
    	while(P) {
    		if(P & 1) res = res * A % Mod;
    		A = A * A % Mod;
    		P >>= 1;
    	}
    	return res;
    }
    
    int rev[Maxn];
    void Init(int len) {
    	for(int i = 1; i < len; ++i) {
    		rev[i] = rev[i >> 1] >> 1;
    		if(i & 1) rev[i] |= len >> 1;
    	}
    }
    
    void NTT(LL a[], int len, int type) {
    	for(int i = 0; i < len; ++i) if(i < rev[i]) swap(a[i], a[rev[i]]);
    	for(int h = 2; h <= len; h <<= 1) {
    		LL gn = qpow(type == 1 ? g : invg , ((Mod - 1) / (LL)h));
    		for(int j = 0; j < len; j += h) {
    			LL gk = 1;
    			for(int k = j; k < j + h / 2; ++k) {
    				LL e = a[k], o = gk * a[k + h / 2] % Mod;
    				a[k] = (e + o) % Mod; a[k + h / 2] = ((e - o) % Mod + Mod) % Mod;
    				gk = gk * gn % Mod;
    			}
    		}
    	}
    	if(type == -1) {
    		LL inv = qpow(len, Mod - 2);
    		for(int i = 0; i < len; ++i) a[i] = a[i] * inv % Mod;
    	}
    }
    
    void polymul(LL a[], LL b[], int lenA, int lenB) {
    	int L = lenA + lenB, len = 1; while(L) {len <<= 1; L >>= 1;}
    	Init(len); NTT(a, len, 1); NTT(b, len, 1);
    	for(int i = 0; i < len; ++i) a[i] = a[i] * b[i];
    	NTT(a, len, -1);
    }
    
    int n, m;
    LL A[Maxn], B[Maxn];
    
    signed main() {
    	read(n); read(m);
    	for(int i = 0; i <= n; ++i) read(A[i]);
    	for(int i = 0; i <= m; ++i) read(B[i]);
    	polymul(A, B, n, m);
    	for(int i = 0; i <= n + m; ++i) printf("%lld ", A[i]);
    	return 0;
    }
    

    鸣谢

    (排名不分先后)

    OI-Wiki

  • 相关阅读:
    WCF 第十三章 可编程站点 为站点创建操作
    WCF 第十三章 可编程站点 所有都与URI相关
    WCF 第十二章 对等网 使用自定义绑定实现消息定向
    WCF 第十三章 可编程站点 使用WebOperationContext
    Using App.Config for user defined runtime parameters
    WCF 第十三章 可编程站点
    WCF 第十三章 可编程站点 使用AJAX和JSON进行网页编程
    WCF 第十二章 总结
    WCF 第十三章 可编程站点 使用WebGet和WebInvoke
    WCF 第十三章 可编程站点 URI和UriTemplates
  • 原文地址:https://www.cnblogs.com/zimujun/p/14472471.html
Copyright © 2011-2022 走看看