NTT
单位根的定义就不说了。
显然有:
带入这个可以直接证得:
我们用图像理性理解可得(就是绕着原点把向量转了180度):
依然根据图像可得(就是绕着原点旋转了360度):
由欧拉公式得:
所以有:
我们先带入\(\omega_n^1,\omega_n^2,...,\omega_n^n\)到多项式\(A\)中,求出\(A(\omega_n^1),A(\omega_n^2),...,A(\omega_n^n)\)。
为了方便,假设长度\(n\)为\(2^k\)(高位不够的添0即可)
把A下标奇偶分类:
显然有:
这两个式子只有后面一项是相反的,可以递归求解。
于是给出代码:
inline void FFT(complex<double> *a, int len) {
if (!len) return ; complex<double> a1[len], a2[len];
for (int i = 0; i < len; ++i) a1[i] = a[i << 1], a2[i] = a[i << 1 | 1];
FFT(a1, len >> 1); FFT(a2, len >> 1);
complex<double> w(cos(PI / len), sin(PI / len)), wk(1, 0);
for (int i = 0; i < len; ++i, wk *= w)
a[i] = a1[i] + wk * a2[i], a[i + len] = a1[i] - wk * a2[i];
}
考虑怎么从点值多项式转换到系数多项式。
我们钦定\(y_i=A(\omega_n^i)\),在有一多项式\(C\),满足:
则我们带入\(\omega_n^{-k}\),得到:
设:
当\(k\neq 0\)时为0,否则为\(n\)
则:
即当\(i=k\)时为\(n\),所以:
我们惊讶的发现这样对\(C\)做一次FFT之后点值除以n就是多项式的系数了。
代码结合一下:
inline void FFT(complex<double> *a, int len, int flag) {
if (!len) return ; complex<double> a1[len], a2[len];
for (int i = 0; i < len; ++i) a1[i] = a[i << 1], a2[i] = a[i << 1 | 1];
FFT(a1, len >> 1, flag); FFT(a2, len >> 1, flag);
complex<double> w(cos(PI / len), flag * sin(PI / len)), wk(1, 0);
for (int i = 0; i < len; ++i, wk *= w)
a[i] = a1[i] + wk * a2[i], a[i + len] = a1[i] - wk * a2[i];
}
发现递归版的码会T,手玩一下发现实际上奇偶变换后下标的操作相当于二进制反过来,可以改成非递归来模拟,自己对着码手玩看看就明白。
inline void FFT(complex *a, int type) {
for (int i = 0; i < lim; ++i)
if (i < rev[i]) swap(a[i], a[rev[i]]);
for (int mid = 1; mid < lim; mid <<= 1) {
complex wn; wn = complex(cos(pi / mid), type * sin(pi / mid));
for (int j = 0; j < lim; j += mid << 1) {
complex bas; bas = complex(1, 0);
for (int k = 0; k < mid; ++k, bas = bas * wn) {
complex x = a[j + k], y = bas * a[j + mid + k];
a[j + k] = x + y;
a[j + mid + k] = x - y;
}
}
}
}
解释一下,第一层循环枚举的是递归的层数,即当前合并的两个多项式的长度。第二层就是枚举当前要合并多项式的起点,第三层就是枚举的具体的那一个系数。这么说不是很清楚,还是自己造样例跟着代码手玩一下就明白了。
而NTT呢?设模数为p,g是p的原根,则不需要证明的给出,\(\omega_n^1\)等价于\(g^{p-1\over n} \bmod p\)。把上面的码代码里的wn换成这个就行了。一般p=998244353,此时g=3。
给个板子:
struct poly {
int n;
vector<ll> x;
inline void NTT(int flag) {
for (int i = 0; i < n; ++i)
if (i < rev[i]) swap(x[i], x[rev[i]]);
for (int mid = 1; mid < n; mid <<= 1) {
ll wn = power(flag == 1 ? G : Gi, (mod - 1) / (mid << 1));
for (int j = 0; j < n; j += mid << 1) {
ll bas = 1;
for (int k = 0; k < mid; ++k, bas = (bas * wn) % mod) {
ll xx = x[j + k], y = (bas * x[j + mid + k]) % mod;
x[j + k] = (xx + y) % mod;
x[j + mid + k] = ((xx - y) % mod + mod) % mod;
}
} cerr << endl;
}
}
};
inline int max_(int a, int b) {
return a > b ? a : b;
}
inline poly mul(poly A, poly B) {
poly a, b; a = A; b = B;
int tmp = a.n + b.n; a.n = 1;
int L = 0;
while (a.n <= tmp) a.n <<= 1, ++L;
for (int i = 0; i <= a.n; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L - 1);
b.n = a.n;
a.NTT(1); b.NTT(1);
for (int i = 0; i < a.n; ++i) a.x[i] = (a.x[i] * b.x[i]) % mod;
a.NTT(-1);
const ll inv = power(a.n, mod - 2);
a.n = tmp;
for (int i = 0; i <= a.n; ++i) a.x[i] = (a.x[i] * inv) % mod;
return a;
}
原文地址:https://www.cnblogs.com/wwlwakioi/p/15004166.html
- Enterprise Library Policy Injection Application Block 之三:PIAB的扩展—创建自定义CallHandler(提供Source Code下载)
- CodeSmith 创建Ado.Net自定义模版(四)
- TensorFlow图像分类教程
- Enterprise Library Policy Injection Application Block 之一: PIAB Overview
- Python教学——第七天
- 大数据将带来电视媒体生态式变革!大数据如何深度融合电视媒体?
- Silverlight SEO优化
- Silverlight性能优化
- WCF后续之旅(6): 通过WCF Extension实现Context信息的传递
- WCF后续之旅(6): 通过WCF Extension实现Context信息的传递
- 理性的相亲方法!精品课:《决策树》
- Asp.Net无刷新分页( jquery.pagination.js)
- 为什么网站需要用CDN来加速?
- Jmeter常用获取数据的几种方式
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- codeforces 1256C (贪心+构造)
- codeforces 722C(带权并查集+反向思维)
- codeforces 1144D(思维)
- 经典的SparkSQL/Hive-SQL/MySQL面试-练习题
- codeforces 1249E(dp)
- Redis-KV数据库Java连接以及Jedis包的使用
- codeforces 1203D1(暴力)
- codeforces 1366B(线段相交)
- 一文搞懂Python自动化测试框架
- codeforces 1005D(数学)
- JSP开发简单实例演示
- Linux笔记【003】| Linux系统目录结构与基本命令
- codeforces1322A(括号匹配)
- codeforces 1296D(贪心)
- codeforces 1399D