FFT可以用来计算多项式乘法,但是复数的运算中含有大量的浮点数,精度较低。对于只有整数参与运算的多项式,有时,\(\text{NTT(Number-Theoretic Transform)}\)会是更好的选择。

若\(a,p\)互素,且\(p>1\),对于\(a^k \equiv 1 (\mod p)\)的最小的\(k\),称为\(a\)模\(p\)的,记做\(\sigma_p(a)\)。

\(E.g.\) \(\sigma_7(2)=3\)

\(2^1\equiv 2(\mod 7)\)

\(2^2\equiv 4(\mod 7)\)

\(2^3\equiv 1(\mod 7)\)

对于一个数\(g\),\(g\)的阶一定是\(p-1\)的约数。

证明:

假设最小的\(k\)不是\(p-1\)的约数,找到\(x\)满足\(xk>p-1>(x-1)k\),由费马小定理可知

\[g^{xk}\equiv g^{p-1}\equiv 1 \equiv g^{xk-(p-1)} (\mod p)
\]

\(xk-(p-1)<k\),与假设矛盾。

原根

定义

\(FFT\)中,我们使用单位复根\(\omega_n^k=\cos k\frac{2\pi}{n}+i\sin k\frac{2\pi}{n}\),那有没有什么能够代替单位复根且解决精度问题呢?这就是原根。

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

若\(p\)为质数,设\(g\)为\(p\)的原根,那么\(g^i \mod p(1<j<p,1\le i\le p-1)\)的结果两两不同。且其等价于\(g^{p-1}\equiv 1(\mod p)\)当且仅当指数为\(p-1\)的时候成立。(这里\(p\)是素数)

简单证明一下:

显然\(g^0 \equiv 1(\mod p)\)

由原根的定义可知满足\(g^{i} \equiv 1(\mod p)\)的最小正整数为\(\varphi(p)=p-1\)

故由指数循环节可知,\(g^i \mod p(1<j<p,1\le i\le p-1)\)的结果两两不同。

性质

考虑在FFT当中我们需要单位复根的以下性质:

  1. \(\omega_n^t\)互不相同,保证点值的合法性;

  2. \(\omega_{2n}^{2k} = \omega_n^k\),用于分治;

  3. \(\omega_n^{k+\frac{n}{2}} = -\omega_n^k\),用于分治;

  4. 当\(k\neq 0\)时,\(1+\omega_n^k+(\omega_n^k)^2+\dots +(\omega_n^k)^{n-1}=0\),用于逆变换。

性质一

令\(\omega_n=g^q\),\(1,g^q,g^{2q},\dots,g^{(n-1)q}\)互不相同,满足性质一

性质二

由\(\omega_n = g^q\)可知,\(\omega_{2n}=g^{\frac{q}{2}}(p=\frac{q}{2} \times 2n + 1)\),故\(\omega_{2n}^{2k} = g^{2k \frac{q}{2}}=g^{kq}=g^q\),满足性质二

性质三

根据费马小定理得

\[\omega_n^n=g^{nq}=g^{p-1}\equiv 1(\mod p)
\]

又因为\((\omega_n^{\frac{n}{2}})^2=\omega_n^n\),所以\(\omega_n^{\frac{n}{2}}\equiv \pm 1 (\mod p)\),而根据性质一可得\(\omega_n^{\frac{n}{2}}\neq \omega_n^0\),即\(\omega_n^{\frac{n}{2}}\equiv -1(\mod p)\)。可推出\(\omega_n^{k+\frac{n}{2}}=\omega_n^k \times \omega_n^{\frac{n}{2}}=-\omega_n^k (\mod p)\),满足性质三

性质四

当\(k\neq 0\)时

\[S(\omega_n^k)=1+\omega_n^k+(\omega_n^k)^2+\dots +(\omega_n^k)^{n-1}
\]

\[\omega_n^k S(\omega_n^k)=\omega_n^k+(\omega_n^k)^2+(\omega_n^k)^3+\dots +(\omega_n^k)^n
\]

\[\omega_n^k S(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^n-1
\]

\[S(\omega_n^k)=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}
\]

性质三中的推论可知,\((\omega_n^k)^n-1=(\omega_n^n)^k-1\equiv \omega_n^n-1\equiv 0 (\mod p)\),故\(S(\omega_n^k)=0\),性质四成立。

求原根

求一个质数的原根,可以用枚举法——枚举\(g\),检验\(g\)是否为\(p\)的原根。

根据前面的关于阶知识可知,检验时,只需枚举\(p-1\)的所有约数\(q\),检验\(g^q\equiv 1(\mod p)\)即可。

代码实现

将\(FFT\)里所有关于\(\omega_n\)的运算替换成\(g^q\)在模意义下的运算即可,注意\(\div n\)要改为\(\times n^{-1}\)

#include <bits/stdc++.h>
using namespace std; typedef long long ll;
inline ll ty() {
char ch = getchar(); ll x = 0, f = 1;
while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
return x * f;
} const int _ = 4e6 + 10;
const ll P = 998244353, G = 3, Gx = 332748118;
int N, M, r[_];
ll A[_], B[_]; ll ksm(ll a, ll b) {
ll ret = 1;
for ( ; b; b >>= 1) {
if (b & 1) ret = ret * a % P;
a = a * a % P;
}
return ret;
} void ntt(int lim, ll *a, int op) {
for (int i = 0; i < lim; ++i)
if (i < r[i]) swap(a[i], a[r[i]]);
for (int len = 2; len <= lim; len <<= 1) {
int mid = len >> 1;
ll Wn = ksm(op == 1 ? G : Gx, (P - 1) / len);
for (int i = 0; i < lim; i += len) {
ll w = 1;
for (int j = 0; j < mid; ++j, w = (w * Wn) % P) {
ll x = a[i + j], y = w * a[i + j + mid] % P;
a[i + j] = (x + y) % P;
a[i + j + mid] = (x - y + P) % P;
}
}
}
} int main() {
#ifndef ONLINE_JUDGE
freopen("ntt.in", "r", stdin);
freopen("ntt.out", "w", stdout);
#endif
N = ty(), M = ty();
for (int i = 0; i <= N; ++i) A[i] = (ty() + P) % P;
for (int i = 0; i <= M; ++i) B[i] = (ty() + P) % P;
int lim = 1, k = 0;
while (lim <= N + M) lim <<= 1, ++k;
for (int i = 0; i < lim ; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (k - 1));
ntt(lim, A, 1);
ntt(lim, B, 1);
for (int i = 0; i < lim; ++i) A[i] = (A[i] * B[i]) % P;
ntt(lim, A, -1);
ll invx = ksm(lim, P - 2);
for (int i = 0; i <= N + M; ++i)
printf("%lld ", (A[i] * invx) % P);
return 0;
}

参考资料

从傅里叶变换到数论变换 | Menci's Blog

快速数论变换(NTT)小结 - 自为风月马前卒 - 博客园

最新文章

  1. cmake 编译 c++ dll 的一个例子(更新1)
  2. webform组合查询和分页
  3. jQuery Ajax 简单的实现跨域请求
  4. Html5中的跨页面消息传输
  5. Count the Colors(线段树染色)
  6. poj 1273 Drainage Ditches 最大流入门题
  7. 上传Test Result和attachment到ALM
  8. Qt 显示透明flash和编写QtWebkit插件
  9. ios 中的UI控件学习总结(1)
  10. Hadoop Hive sql语法详解
  11. 工作线程基类TaskSvc
  12. java关键字中文对比
  13. hive的简单理解--笔记
  14. HashMap 你真的了解吗?
  15. java自动化-数据驱动juint演示,上篇
  16. 使用Keepalived配置主从热备实现Nginx高可用(HA)
  17. Android jni中回调java的方法
  18. 使用GO开发ChainCode
  19. Object与Class的区别
  20. 读取 Excel 之 Epplus

热门文章

  1. MySQL的多表联查
  2. Redis Python(二)
  3. vue应用调试工具 vue-devtools安装
  4. Linux(Centos7)下redis5集群搭建和使用
  5. Nginx基础知识点总结和优化项
  6. Nginx + FastCGI + Django在windows上部署及nginx常用命令
  7. SpringMVC 简单限流方案设计
  8. 阿里面试官:字符串在JVM中如何存放?90%的人就真的只回答在哪里存放
  9. jQuery-点击返回顶部
  10. c++实现通讯录管理系统(控制台版)