快速数论变换(NTT)小结

NTT

在FFT中,我们需要用到复数,复数虽然很神奇,但是它也有自己的局限性——需要用double类型计算,精度太低

那有没有什么东西能够代替复数且解决精度问题呢?

这个东西,叫原根

原根

原根的定义

m是正整数,a是整数,若am的阶等于\phi(m),则称a为模m的一个原根

定义中用到了群论的一些知识,不过不会也没关系,不影响接下来的学习

我们定义P为素数,gP的原根

接下来不加证明的扔出一个很重要定理

  • P为素数,假设一个数gP的原根,那么g^i \mod P (1<g<P,0<i<P)的结果两两不同

不要问我为什么,因为我也不知道。。

考虑原根为什么能代替单位根进行运算,(这部分可以跳过)

原因很简单,因为它具有和单位根相同的性质

在FFT中,我们用到了单位根的四条性质,而原根也满足这四条性质

1 . 对于所有\omega_n ^ t (0 \leq t \leq n - 1)均不相同

这一条可以由上面的定理得到

2 .\omega_{2n} ^ {2k} = \omega_n ^ k

通过代换可以得到

3 .\omega_n ^ { k + \frac{n}{2} } = -\omega_n ^ k

根据费马小定理和性质1可以得到

4 .1 + \omega_n ^ k + (\omega_n ^ k) ^ 2 + \dots + (\omega_n ^ k) ^ {n - 1} = 0

由性质3和FFT中傅里叶逆变换的定理可以得到

这样我们最终可以得到一个结论

\omega_n \equiv g^\frac{p-1}{n} \mod p

然后把FFT中的\omega\_n都替换掉就好了

p建议取998244353,它的原根为3

如何求任意一个质数的原根呢?

可以证明满足g^r \equiv 1(\mod p)的最小的r一定是p-1的约数

对于质数p,质因子分解p−1,若g^{\frac{p-1}{p_i}} \neq 1 \pmod p恒成立,gp的原根

实现

NTT求卷积代码:

确实比FFT快了不少

#include<cstdio>
#define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1<<21, stdin), p1 == p2) ? EOF : *p1++)
#define swap(x,y) x ^= y, y ^= x, x ^= y
#define LL long long 
const int MAXN = 3 * 1e6 + 10, P = 998244353, G = 3, Gi = 332748118; 
char buf[1<<21], *p1 = buf, *p2 = buf;
inline int read() { 
    char c = getchar(); int x = 0, f = 1;
    while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
int N, M, limit = 1, L, r[MAXN];
LL a[MAXN], b[MAXN];
inline LL fastpow(LL a, LL k) {
	LL base = 1;
	while(k) {
		if(k & 1) base = (base * a ) % P;
		a = (a * a) % P;
		k >>= 1;
	}
	return base % P;
}
inline void NTT(LL *A, int type) {
	for(int i = 0; i < limit; i++) 
		if(i < r[i]) swap(A[i], A[r[i]]);
	for(int mid = 1; mid < limit; mid <<= 1) {	
		LL Wn = fastpow( type == 1 ? G : Gi , (P - 1) / (mid << 1));
		for(int j = 0; j < limit; j += (mid << 1)) {
			LL w = 1;
			for(int k = 0; k < mid; k++, w = (w * Wn) % P) {
				 int x = A[j + k], y = w * A[j + k + mid] % P;
				 A[j + k] = (x + y) % P,
				 A[j + k + mid] = (x - y + P) % P;
			}
		}
	}
}
int main() {
	N = read(); M = read();
	for(int i = 0; i <= N; i++) a[i] = (read() + P) % P;
	for(int i = 0; i <= M; i++) b[i] = (read() + P) % P;
	while(limit <= N + M) limit <<= 1, L++;
	for(int i = 0; i < limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));	
	NTT(a, 1);NTT(b, 1);	
	for(int i = 0; i < limit; i++) a[i] = (a[i] * b[i]) % P;
	NTT(a, -1);	
	LL inv = fastpow(limit, P - 2);
	for(int i = 0; i <= N + M; i++)
		printf("%d ", (a[i] * inv) % P);
	return 0;
}

原创声明,本文系作者授权云+社区发表,未经许可,不得转载。

如有侵权,请联系 yunjia_community@tencent.com 删除。

编辑于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏ml

2015编程之美(资格赛)--基站选址

题目3 : 基站选址 时间限制:2000ms 单点时限:1000ms 内存限制:256MB 描述 需要在一个N × M的网格中建立一个通讯基站,通讯基站仅必须建...

33450
来自专栏斑斓

Spark 1.4为DataFrame新增的统计与数学函数

Spark一直都在快速地更新中,性能越来越快,功能越来越强大。我们既可以参与其中,也可以乐享其成。 目前,Spark 1.4版本在社区已经进入投票阶段,在Gi...

37370
来自专栏web前端教室

javascript 算法初识

最近有空,想学习下算法。一直感觉它很高深的样子,尤其我数学又不好。 但我还是想学学看,万一能学到点东西呢,,, 先来了解下算法的定义:是指解题方案的准确而完整的...

21060
来自专栏杨建荣的学习笔记

对于随机数的一些分析

多年前我朋友圈的一个朋友公司年会抽奖出现了下面的这样一幕:CTO现场review代码。本来带着一丝娱乐精神,结果被无限放大了。所以年会中大家都会很自然想revi...

36480
来自专栏小樱的经验随笔

BZOJ 2038: [2009国家集训队]小Z的袜子(hose)【莫队算法裸题&&学习笔记】

2038: [2009国家集训队]小Z的袜子(hose) Time Limit: 20 Sec  Memory Limit: 259 MB Submit: 9...

31360
来自专栏落影的专栏

OpenGLES进阶教程7-天空盒效果

教程 OpenGLES入门教程1-Tutorial01-GLKit OpenGLES入门教程2-Tutorial02-shader入门 OpenGLES入门...

40460
来自专栏数据小魔方

动态地理信息可视化——leaflet在线地图简介

最近稍微涉猎了一下leaflet这个包,突然感到发现了动态可视化的新大门,这个包所提供的地图类型、动态效果、图层展示方式都大大扩展了ggplot作图系统的在数据...

56340
来自专栏技术专栏

二维数组的DP问题

问题:平面上有N*M个格子,每个格子中放着一定数量的苹果。你从左上角的格子开始,每一步只能向下走或是向右走,每次走到一个格子上就把格子里的苹果收集起来,这样下去...

16230
来自专栏利炳根的专栏

学习笔记CB003:分块、标记、关系抽取、文法特征结构

分块,根据句子的词和词性,按照规则组织合分块,分块代表实体。常见实体,组织、人员、地点、日期、时间。名词短语分块(NP-chunking),通过词性标记、规则识...

23390
来自专栏大数据

有向无环图检测

01 — Spark背景介绍 Apache Spark 是专为大规模数据处理而设计的快速通用的计算引擎。Spark 是一种与 Hadoop 相似的开源集群计算环...

49970

扫码关注云+社区

领取腾讯云代金券