首页 > 其他分享 >浅谈 DFT、IDFT、NTT

浅谈 DFT、IDFT、NTT

时间:2024-10-01 08:53:27浏览次数:8  
标签:begin end matrix DFT NTT cdots 浅谈 aligned vdots

DFT(离散傅里叶变换)

多项式分治
最早可能是由高斯发现的多项式可以分治,但他的手稿并未作为论文发表。

考虑多项式 \(F(x) = a_0 + a_1 x^{1} + a_2 x^{2} + \cdots + a_{n - 1} x^{n - 1}\) 其中 \(n = 2^{k} \ (k \geq 0)\) 。(任意多项式可以通过高位补 \(0\) 化为这个形式。)
有:

\[\begin{aligned} F(x) &= a_0 + a_1 x^{1} + a_2 x^{2} + \cdots + a_{n - 1} x^{n - 1} \\ &= a_0 + a_2 x^{2} + a_4 x^{4} + \cdots a_{n - 2} x^{n - 2} + a_1 x + a_3 x^{3} + a_5 x^{5} + \cdots a_{n - 1} x^{n - 1} \\ &= a_0 + a_2 x^{2} + a_4 x^{4} + \cdots a_{n - 2} x^{n - 2} + x (a_1 + a_3 x^{2} + a_5 x^{4} + \cdots a_{n - 1} x^{n - 2}) \\ &= G(x^{2}) + x H(x^{2}) \\ \end{aligned} \]

最终可以分治为 \(n\) 个关于 \(x^{k}\) 的多项式。

单位根
(高斯注意到)如果我们将 \(x = e^{i \frac{l}{n} 2 \pi}\) 代入多项式 \(F(x)\) ,则能以 \(O(n \log n)\) 的时间通过分治得到一组关于 \(e^{i \frac{l}{n} 2 \pi}\) 的点值表示。

由欧拉公式 \(e^{i \theta} = \cos \theta + i \sin \theta\) (正确性是你会发现左右两边的泰勒展开一样),将 \(\cos \frac{l}{n} 2 \pi + i \sin \frac{l}{n} 2 \pi\) 的形式写作 \(w_{n}^{l}\) ,即 \(n\) 次单位复根的 \(l\) 次幂。

蝴蝶变换
(高斯)注意到:

\[\begin{aligned} F(w_n^{l}) = G(w_n^{2l}) + w_n^{l} H(w_n^{2l}) &= G(w_{n/2}^{l}) + w_n^{l} H(w_{n/2}^{2l}) \\ F(w_n^{l + n/2}) = G(w_n^{2(l + n/2)}) + w_n^{l + n/2} H(w_n^{2(l + n/2)}) = G(w_{n/2}^{l + n/2}) + w_n^{l + n/2} H(w_{n/2}^{l + n/2}) &= G(w_{n/2}^{l}) - w_n^{l} H(w_{n/2}^{l}) \\ \end{aligned} \]

这展示了 \(F\) 如何由 \(G, H\) 计算。

由于分治最多有 \(\log n\) 层,每层的计算都需要 \(O(n)\) 复杂度。\(DFT\) 的时间是 \(O(n \log n)\) 。

typedef long doube db;
const double PI = acos(-1.);
void DFT(std::vector<CP> *a) {
	int n = a.size();
	for (int m = 2; m < n; m *= 2) {
		db Arg = 1.L / m * 2.L * PI;
		std::complex<db> wm(cos(Arg), sin(Arg));
		for (int l = 0; l + m - 1 < n; l += m) {
			std::complex<db> wmi = wm;
			for (int i = l, j = l + m / 2; i < l + m / 2; i++, j++) {
				std::complex<db> u = a[i], v = a[j];
				a[i] = u + wmi * v;
				a[j] = u - wmi * v;
				wmi *= wm;
			}
		}
	}
}

IDFT(离散傅里叶逆变换)

注意多项式点值表示的矩阵

\[\begin{aligned} \left [ \begin{matrix} 1 & x_{0}^{1} & x_{0}^{2} & \cdots & x_{0}^{n - 1} \\ 1 & x_{1}^{1} & x_{1}^{2} & \cdots & x_{1}^{n - 1} \\ 1 & x_{2}^{2} & x_{2}^{2} & \cdots & x_{2}^{n - 1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n - 1}^{n - 1} & x_{n - 1}^{n - 1} & \cdots & x_{n - 1}^{n - 1} \\ \end{matrix} \right ] \cdot \left [ \begin{matrix} a_{0} \\ a_{1} \\ a_{2} \\ \vdots\\ a_{n - 1} \\ \end{matrix} \right ] = \left [ \begin{matrix} F(x_{0}) \\ F(x_{1}) \\ F(x_{2}) \\ \vdots \\ F(x_{n - 1}) \\ \end{matrix} \right ] \end{aligned} \]

\(Vandermonde\) 矩阵可逆,于是:

\[\left [ \begin{matrix} a_{0} \\ a_{1} \\ a_{2} \\ \vdots\\ a_{n - 1} \\ \end{matrix} \right ] = \begin{aligned} \left [ \begin{matrix} 1 & x_{0}^{1} & x_{0}^{2} & \cdots & x_{0}^{n - 1} \\ 1 & x_{1}^{1} & x_{1}^{2} & \cdots & x_{1}^{n - 1} \\ 1 & x_{2}^{2} & x_{2}^{2} & \cdots & x_{2}^{n - 1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n - 1}^{n - 1} & x_{n - 1}^{n - 1} & \cdots & x_{n - 1}^{n - 1} \\ \end{matrix} \right ]^{-1} \cdot \left [ \begin{matrix} F(x_{0}) \\ F(x_{1}) \\ F(x_{2}) \\ \vdots \\ F(x_{n - 1}) \\ \end{matrix} \right ] \end{aligned} \]

关于范德蒙德矩阵的逆矩阵,有初等爆算和拉格朗日插值法。这里仅说明如何用构造法得到 \(IDFT\) 的逆矩阵。

注意 \(DFT\) 的点值矩阵

\[\begin{aligned} \left [ \begin{matrix} 1 & (w_{n}^{0})^{1} & (w_{n}^{0})^{2} & \cdots & (w_{n}^{0})^{n - 1} \\ 1 & (w_{n}^{1})^{1} & (w_{n}^{1})^{2} & \cdots & (w_{n}^{1})^{n - 1} \\ 1 & (w_{n}^{2})^{1} & (w_{n}^{2})^{2} & \cdots & (w_{n}^{2})^{n - 1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & (w_{n}^{n - 1})^{1} & (w_{n}^{n - 1})^{2} & \cdots & (w_{n}^{n - 1})^{n - 1} \\ \end{matrix} \right ] \end{aligned} \]

\[\begin{aligned} &\left [ \begin{matrix} (w_{n}^{-i})^{0} & (w_{n}^{-i})^{1} & (w_{n}^{-i})^{2} & \cdots & (w_{n}^{-i})^{n - 1} \end{matrix} \right ] \cdot \left [ \begin{matrix} (w_{n}^{0})^{j} \\ (w_{n}^{1})^{j} \\ (w_{n}^{2})^{j} \\ \vdots \\ (w_{n}^{n - 1})^{j} \\ \end{matrix} \right ] \\ &= \sum_{k = 0}^{n - 1} (w_{n}^{-i})^{k} \cdot (w_{n}^{k})^{j} = \sum_{k = 0}^{n - 1} (w_{n}^{-i})^{k} \cdot (w_{n}^{j})^{k} = \sum_{k = 0}^{n - 1} (w_{n}^{j - i})^{k} \\ \end{aligned} \]

当 \(j = i\) , \(\sum_{k = 0}^{n - 1} 1^{k} = n\) 显然。
当 \(i \neq j\) ,由恒等式 \(x^{n} - 1 = (x - 1)(x^{n - 1} + x^{n - 2} + \cdots + 1)\) ,令 \(x = w_{n}^{l} \ (0 < l < n)\) ,则 \((w_{n}^{l})^{n} = (w_{n}^{n})^{l} = 1\) ,显然 \(w_{n}^{l} \neq 1\) ,则 \(\sum_{k = 0}^{n - 1} (w_{n}^{l})^{k} = 0\) 。
于是

\[\begin{aligned} \frac{1}{n} \cdot \left [ \begin{matrix} 1 & (w_{n}^{-0})^{1} & (w_{n}^{-0})^{2} & \cdots & (w_{n}^{-0})^{n - 1} \\ 1 & (w_{n}^{-1})^{1} & (w_{n}^{-1})^{2} & \cdots & (w_{n}^{-1})^{n - 1} \\ 1 & (w_{n}^{-2})^{1} & (w_{n}^{-2})^{2} & \cdots & (w_{n}^{-2})^{n - 1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & (w_{n}^{-(n - 1)})^{1} & (w_{n}^{-(n - 1)})^{2} & \cdots & (w_{n}^{-(n - 1)})^{n - 1} \\ \end{matrix} \right ] \cdot \left [ \begin{matrix} 1 & (w_{n}^{0})^{1} & (w_{n}^{0})^{2} & \cdots & (w_{n}^{0})^{n - 1} \\ 1 & (w_{n}^{1})^{1} & (w_{n}^{1})^{2} & \cdots & (w_{n}^{1})^{n - 1} \\ 1 & (w_{n}^{2})^{1} & (w_{n}^{2})^{2} & \cdots & (w_{n}^{2})^{n - 1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & (w_{n}^{n - 1})^{1} & (w_{n}^{n - 1})^{2} & \cdots & (w_{n}^{n - 1})^{n - 1} \\ \end{matrix} \right ] = \left [ \begin{matrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \\ \end{matrix} \right ] \end{aligned} \]

于是得到了用于 \(IDFT\) 的逆矩阵。
注意 \(w_{n}^{-l} = \overline{w_{n}^{l}}\) 。于是 \(IDFT\) 时,只需将用共轭单位根进行 \(DFT\) 能得到答案的 \(n\) 倍。只需最后让答案除以 \(n\) 。

\(DFT + IDFT\) 即数论函数意义下的 \(FFT\) 。

FFT 精度问题
显然浮点单位根幂乘上去精度很可能爆炸,最好是使用 \(long\ double\) 的浮点数。

NTT

数论函数中,当所需要的结果很大。无法使用普通数据类型存储,需要对结果取模。

显然不能用 \(FFT\) 直接做乘法,至少单位根累乘的过程中无法在浮点数上取模。

特殊模数模意义下的 \(FFT\) 又叫 \(NTT\) ,全程快速数论变换。

考虑一个质数 \(P\) 同时也为 \(2\) 的幂次。令 \(\mod p\) 的一个原根为 \(g\) ,钦定 \(g^{\frac{P - 1}{n} \times l} = g_{n}^{l}\) 。考虑 \(n\) 为 \(2\) 的幂次。

考虑如下一些性质:

  1. \(g_{n}^{2l} \equiv g_{n / 2}^{l}\) 。
    • 证明: \(g_{n}^{2l} = g^{\frac{P - 1}{n} 2l} = g^{\frac{P - 1}{n / 2} l}\) 。由 \(n\) 为 \(2\) 的幂次显然成立。\(\square\)
  2. \(g_{n}^{0} \equiv g_{n}^{n} \equiv 1\) 。
    • 显然成立。
  3. \(g_{n}^{l + n} \equiv g_{n}^{l}\) 。
    • 显然成立。
  4. \(g_{n}^{n/2} \equiv -1\) 。
    • 证明: 显然 \(g\) 不是 \(\bmod P\) 的二次剩余,则 \(g\) 是 \(\bmod P\) 的非二次剩余。\(\square\)
  5. \((g_{n}^{l})^{k}) \equiv (g_{n}^{k})^{l} \equiv g_{n}^{lk}\) 。
    • 显然成立。
  6. \(\sum_{k = 0}^{n - 1} g_{n}^{l} \equiv 0 \ s.t. 0 < l < n\) 。
    • 证明: 考虑 \(x^{n} - 1 = (x - 1)(x^{n - 1} + x^{n - 2} + \cdots + 1)\) 。\((g_{n}^{l})^{n} = 1\) 显然,\(g_{n}^{l} \neq 1\) 显然,于是 \(\sum_{k = 0}^{n - 1} (g_{n}^{l})^{k} = 0\) 。\(\square\)

于是 \(g^{\frac{P - 1}{n} l} = g_{n}^{l}\) 具备了 \(w_{n}^{l}\) 在 \(FFT\) 中的所有性质。

void NDFT(std::vector<i64> *a) {
	int n = a.size();
	for (int m = 2; m < n; m *= 2) {
		i64 gm = g; // g^{-1} = ksm(g, P - 2);
		for (int l = 0; l + m - 1 < n; l += m) {
			i64 gmi = g;
			for (int i = l, j = l + m / 2; i < l + m / 2; i++, j++) {
				i64 u = a[i], v = a[j];
				a[i] = (u + gmi * v) % P;
				a[j] = ((u - gmi * v) % P + P) % P;
				gmi *= g;
			}
		}
	}
}

标签:begin,end,matrix,DFT,NTT,cdots,浅谈,aligned,vdots
From: https://www.cnblogs.com/zsxuan/p/18439615

相关文章

  • Svnlook使用浅谈(配置svn上传必须添加备注和删除权限)
      在配置svn上传必须添加备注和删除权限前,我先隆重介绍下今天用到的svn命令svnlook。svnlook是检验Subversion版本库不同方面的命令行工具,不会对版本库有任何修改,只是查看,包括作者信息、文件内容、更改历史、文件大小、属性等。当然它有自己的独特语法(1)语法格式:svnlookREPOS......
  • 浅谈笛卡尔树
    [介绍(百度百科)](笛卡尔树_百度百科(baidu.com))笛卡尔树是一种特定的二叉树数据结构,可由数列构造,在范围最值查询、范围\(top_k\)查询(rangetopkqueries)等问题上有广泛应用。它具有堆的有序性,中序遍历可以输出原数列。笛卡尔树结构由Vuillmin(1980)在解决范围搜索的几何数据结......
  • 浅谈数据代理
    <!DOCTYPEhtml><htmllang="en"><head><metacharset="UTF-8"><metaname="viewport"content="width=device-width,initial-scale=1.0"><title>Document</title><......
  • Manacher 算法浅谈
    \(Zero.\)\(~~\)前言杂谈认识我的人都喜欢叫我马拉车,如今,马拉车来浅谈Manacher了(不就是某天打板子的时候打错了吗,不就是啪啪打脸了吗)。首先大家需要知道,Manacher不是很常考,但是也是一项必备的算法。当遇到回文串之类的问题时,别人辛辛苦苦打一堆哈希,你用Manacher算法两个并......
  • 轻松编排工作流,浅谈DolphinScheduler如何使用Python调用API接口?
    最近,在做某大型零售企业项目时,有客户用到DolphinScheduler,并咨询是否可以用Python脚本编排工作流?该如何实现?相信有很多人会有这样的疑问,那么,本文将为我们简单分享DolphinScheduler的优势和实际使用。为什么企业数据开发要使用海豚调度?当企业在做数据开发时,任务调度平台会扮演自......
  • 轻松编排工作流,浅谈DolphinScheduler如何使用Python调用API接口?
    最近,在做某大型零售企业项目时,有客户用到DolphinScheduler,并咨询是否可以用Python脚本编排工作流?该如何实现?相信有很多人会有这样的疑问,那么,本文将为我们简单分享DolphinScheduler的优势和实际使用。为什么企业数据开发要使用海豚调度?当企业在做数据开发时,任务调度平台会扮演自动......
  • 浅谈分时电价下含电动汽车的微电网群双层多目标优化调度
    摘要:为解决大规模电动汽车无序充电导致电网出现“峰上加峰”现象,依据电动汽车充电地点的不同将配电网划分为居民区、办公区、商业区微电网,提出基于峰谷差、分时电价、用户充电满意度多目标下的电动汽车充电模式,建立了微电网内运营商峰谷差—用户充电费用少和充电满意度的双盈多目标......
  • 浅谈医院配电系统谐波分析与治理技术方案
    摘要:文章从谐波治理的危害、治理意义、谐波源组成、谐波治理等方面进行了论述。目的在于通过综合整治电网的谐波,有效地改善医院配电系统的安全、可靠、节能。关键词:医院;配电系统;谐波治理0引言配电系统中存在的谐波成分,会导致电力系统的电能品质下降,从而导致电力供应的可靠性下降。......
  • 浅谈如何处理大语言模型训练数据之三开源数据集介绍
    随着最近这些年来基于统计机器学习的自然语言处理的算法的发展,以及信息检索研究的需求,特别是近年来深度学习和预训练语言模型的研究以及国内国外许多大模型的开源,研究人员们构建了多种大规模开源数据集,涵盖了网页、图片、论文、百科等多个领域。在构建大语言模型时,数据的质量和多......
  • 浅谈软件工程
    基本概念软件工程是指导计算机软件开发和维护的一门工程学科,将合理的管理技术和前沿的技术方法结合起来,经济地开发出高质量的软件并有效地维护。软件工程是:①把系统的、规范的、可度量的途径应用于软件开发、运行和维护过程,也就是把工程应用于软件;②研究①中提到的途径。——1......