快速数论变换(NTT)
来源:互联网 发布:python web ui 编辑:程序博客网 时间:2024/04/30 01:54
FFT的复数精度问题让我很不爽,算法导论上确实提了一下可以用数论的方法实现傅立叶变换,可是我一直不知道怎么搞,现在终于找到了质料
以下内容转自ACdreamers (http://blog.csdn.net/acdreamers/article/details/39026505)
在上一篇文章中 http://blog.csdn.net/acdreamers/article/details/39005227 介绍了用快速傅里叶变
换来求多项式的乘法。可以发现它是利用了单位复根的特殊性质,大大减少了运算,但是这种做法是对复数系数的矩阵
加以处理,每个复数系数的实部和虚部是一个正弦及余弦函数,因此大部分系数都是浮点数,我们必须做复数及浮点数
的计算,计算量会比较大,而且浮点数的计算可能会导致误差增大。
今天,我将来介绍另一种计算多项式乘法的算法,叫做快速数论变换(NTT),在离散正交变换的理论中,已经证明在
复数域内,具有循环卷积特性的唯一变换是DFT,所以在复数域中不存在具有循环卷积性质的更简单的离散正交变换。
因此提出了以数论为基础的具有循环卷积性质的快速数论变换。
回忆复数向量,其离散傅里叶变换公式如下
离散傅里叶逆变换公式为
今天的快速数论变换(NTT)是在上进行的,在快速傅里叶变换(FFT)中,通过次单位复根来运算的,即满
足的,而对于快速数论变换来说,则是可以将看成是的等价,这里是模素数
的原根(由于是素数,那么原根一定存在)。即
所以综上,我们得到数论变换的公式如下
而数论变换的逆变换公式为
这样就把复数对应到一个整数,之后一切都是在系统内考虑。
上述数论变换(NTT)公式中,要求是素数且必须是的因子。由于经常是2的方幂,所以可以构造形
如的素数。通常来说可以选择为费马素数,这样的变换叫做费马数数论变换。
这里我们选择,,这样得到模的原根值为。
题目:http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1028
分析:题目意思就是大数相乘,此处用快速数论变换(NTT)实现。
代码:
- #include <iostream>
- #include <string.h>
- #include <stdio.h>
- using namespace std;
- typedef long long LL;
- const int N = 1 << 18;
- const int P = (479 << 21) + 1;
- const int G = 3;
- const int NUM = 20;
- LL wn[NUM];
- LL a[N], b[N];
- char A[N], B[N];
- LL quick_mod(LL a, LL b, LL m)
- {
- LL ans = 1;
- a %= m;
- while(b)
- {
- if(b & 1)
- {
- ans = ans * a % m;
- b--;
- }
- b >>= 1;
- a = a * a % m;
- }
- return ans;
- }
- void GetWn()
- {
- for(int i=0; i<NUM; i++)
- {
- int t = 1 << i;
- wn[i] = quick_mod(G, (P - 1) / t, P);
- }
- }
- void Prepare(char A[], char B[], LL a[], LL b[], int &len)
- {
- len = 1;
- int len_A = strlen(A);
- int len_B = strlen(B);
- while(len <= 2 * len_A || len <= 2 * len_B) len <<= 1;
- for(int i=0; i<len_A; i++)
- A[len - 1 - i] = A[len_A - 1 - i];
- for(int i=0; i<len - len_A; i++)
- A[i] = '0';
- for(int i=0; i<len_B; i++)
- B[len - 1 - i] = B[len_B - 1 - i];
- for(int i=0; i<len - len_B; i++)
- B[i] = '0';
- for(int i=0; i<len; i++)
- a[len - 1 - i] = A[i] - '0';
- for(int i=0; i<len; i++)
- b[len - 1 - i] = B[i] - '0';
- }
- void Rader(LL a[], int len)
- {
- int j = len >> 1;
- for(int i=1; i<len-1; i++)
- {
- if(i < j) swap(a[i], a[j]);
- int k = len >> 1;
- while(j >= k)
- {
- j -= k;
- k >>= 1;
- }
- if(j < k) j += k;
- }
- }
- void NTT(LL a[], int len, int on)
- {
- Rader(a, len);
- int id = 0;
- for(int h = 2; h <= len; h <<= 1)
- {
- id++;
- for(int j = 0; j < len; j += h)
- {
- LL w = 1;
- for(int k = j; k < j + h / 2; k++)
- {
- LL u = a[k] % P;
- LL t = w * (a[k + h / 2] % P) % P;
- a[k] = (u + t) % P;
- a[k + h / 2] = ((u - t) % P + P) % P;
- w = w * wn[id] % P;
- }
- }
- }
- if(on == -1)
- {
- for(int i = 1; i < len / 2; i++)
- swap(a[i], a[len - i]);
- LL Inv = quick_mod(len, P - 2, P);
- for(int i = 0; i < len; i++)
- a[i] = a[i] % P * Inv % P;
- }
- }
- void Conv(LL a[], LL b[], int n)
- {
- NTT(a, n, 1);
- NTT(b, n, 1);
- for(int i = 0; i < n; i++)
- a[i] = a[i] * b[i] % P;
- NTT(a, n, -1);
- }
- void Transfer(LL a[], int n)
- {
- int t = 0;
- for(int i = 0; i < n; i++)
- {
- a[i] += t;
- if(a[i] > 9)
- {
- t = a[i] / 10;
- a[i] %= 10;
- }
- else t = 0;
- }
- }
- void Print(LL a[], int n)
- {
- bool flag = 1;
- for(int i = n - 1; i >= 0; i--)
- {
- if(a[i] != 0 && flag)
- {
- printf("%d", a[i]);
- flag = 0;
- }
- else if(!flag)
- printf("%d", a[i]);
- }
- puts("");
- }
- int main()
- {
- GetWn();
- while(scanf("%s%s", A, B)!=EOF)
- {
- int len;
- Prepare(A, B, a, b, len);
- Conv(a, b, len);
- Transfer(a, len);
- Print(a, len);
- }
- return 0;
- }
- 快速数论变换(NTT)
- 快速数论变换(NTT)
- NTT(快速数论变换)
- 快速数论变换模板(NTT)
- 【模板】快速数论变换ntt
- NTT(快速数论变换)模板
- [SDOI2015]序列统计 (NTT:快速数论变换 + 快速幂)
- 快速傅里叶变换(FFT)和数论变换(NTT)模板
- Gym 100341C AVL TREE(NTT快速数论变换)
- 【快速数论变换NTT】AVL tree
- NTT FFT 数论变换 快速傅里叶变换 模板
- A*B NTT快速数论变换
- [NTT] 快速数论变换学习笔记
- 【快速傅立叶变换fft&数论变换ntt学习小记】
- 【模板】快速数论变换ntt(long long版)
- 【从0开始ACM】【LV3】【数学】【 FFT 快速傅里叶变换 && NTT 更高精度的 快速数论变换】
- NTT(快速数论变换)用到的各种素数及原根
- 多项式乘法运算 NTT(数论变换)实现
- sqlserver日期函数
- Cassandra 2.1.0 安装详解(单点与单个集群)
- 65446413123112
- 自己写的表格的动态加载,有些不足地方
- linux采用编译内核的方法增加系统功能调用
- 快速数论变换(NTT)
- WPF中控件弹出动画
- android 滚动条 相关属性
- 中序线索化二叉树
- C++数组名a和&a的区别
- 配置django运行环境(apache+wsgi)
- Monkey源码分析番外篇之Android注入事件的三种方法比较
- 屏幕方向android:screenOrientation
- 使用Fabric部署网站应用