首页 > 其他分享 >《不一般的 DFT》阅读随笔

《不一般的 DFT》阅读随笔

时间:2022-11-27 06:33:06浏览次数:79  
标签:frac DFT sum 点值 阅读 多项式 随笔 omega

感觉上前置知识是毛啸 16 年的论文?
我手头也有,到时候发现有 at 到的地方就插一嘴说一句
srds 先这篇是因为有纸质版的这篇

感觉上大篇幅在讲复杂度模数大小相关的做法。

1 引言

我这写个啥?
概括一下。FFT nb,但是光会板子不够,给你整点有用的新玩意。我先扯点别人写过的,再给你讲讲任意长度 FFT。反正不要钱多少看一点~

2 定义及说明

2.1 约定

  1. \(|f(x)| :=\) 多项式 \(f\) 的次数/系数。
  2. \(f(p) :=\) 多项式 \(f(x)\) 在 \(p\) 处的点值。
  3. \(\omega_n :=\) \(n\) 次单位根,即 \(\cos\left(\frac {2\pi} n\right) + \sin\left(\frac {2\pi} n\right)i\)。
  4. \(DFT(f):=\) 多项式 \(f(x)\) 在进行 DFT 后得到的系数序列。
  5. \(f_i :=\) 多项式 \(f(x)\) 的第 \(i\) 项系数。一般地,\(0\le i < |f(x)|\)。

2.2 DFT

离散傅里叶变换(DFT)是将 \(\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}\) 依次带入多项式 \(f(x)\) 中得到的长度为 \(n\) 的点值序列 \(a\)。称 \(n = 2^k\) 为模长。一般有 \(n > |f(x)|\)。

在已知长度为 \(n\) 的多项式 \(f(x)\) 的点值序列 \(a\) 的情况下,我们构造多项式 \(g(x) = \sum_{i=0}^{n-1} a_i x^i\),并求出依次将 \(\omega_n^0, \omega_n^{-1},\cdots,\omega_{n}^{-(n-1)}\) 带入得到的点值序列 \(b\),则 \(f(x) = \frac 1n\sum_{i=0}^{n-1} b_i x^i\)。

证明:

\[\begin{aligned} b_k &= \ \sum_{i=0}^{n-1} a_i\left(\omega_n^{-k}\right)^i \\ & = \ \sum_{i=0}^{n-1} \left(\sum_{j=0}^{n-1} f_i\left(\omega_n^i\right)^j \right)\left(\omega_n^{-k}\right)^i \\ & = \ \sum_{i=0}^{n-1} f_i\left(\sum_{j=0}^{n-1} \left(\omega_n^{i-k}\right)^j \right) \\ & = \ \sum_{i=0}^{n-1} f_i\times n[i - k = 0] \\ & = \ n f_k \end{aligned}\]

\(3\to 4\) 的转化考虑 \(i-k=0\) 时 \(\omega_{n}^{i-k} = 1\),反之通过等比数列求和得到总和为 \(0\)。

2.3 MTT

为啥它会把一般模数 FFT 叫做 MTT 呢?

首先有一种拆系数 FFT。
假设我们需要求解 \(A(x)\times B(x)\) 的各项系数对 \(1e9 + 7\) 等模数取模后的值。如果直接干上去,运行多项式乘积时产生值的绝对值是 \(O(p^2|A(x)|)\) 的,不能接受。
我们将多项式的各项系数写作 \(xC + y(0\le y < C)\) 的形式,其中 \(C\) 是一个正整数。那写 \(A(x) = CA_1(x) + A_2(x)\),\(B(x) = CB_1(x) + B_2(x)\)。则

\[A(x)\times B(x) = C^2 A_1(x)B_1(x) + CA_1(x) B_2(x) + CA_2(x) B_1(x) + A_2(x)B_2(x) \]

我们取 \(C = \left\lfloor \sqrt p\right\rfloor\),则运行多项式乘积时产生值的绝对值是 \(O(p|A(x)|)\) 的,可以接受。

然后说一下三次变两次的优化。srds这部分论文里没有
就是说啊,我们构造多项式 \(P(x) = A(x) + B(x)i,\ Q(x) = A(x) - B(x)i\),则根据定义由 \(DFT(A)_ i = \frac{DFT(P)_ i + DFT(Q)_ i}2,\ DFT(B)_ i = -i \frac{DFT(P)_ i - DFT(Q)_ i}2\)。并且,我们有 \(DFT(P)_ i = DFT(Q)_ {n-i}\)。证明考虑在写出 \(DFT(Q)_ i\) 的表式时展开两个复数的相乘部分,并取共轭。
容易发现这样只需要两次 \(DFT\),分别是取得 \(DFT(P),\ DFT(Q)\) 和 \(A(x)\times B(x)\)。

然后通过这样的优化,我们能得到一种四次 DFT 的 MTT 写法。这里有一份常数不大的实现。

其次有一种做两/三次 NNT 模数下的 DFT 后使用 CRT 合并的写法。常数太大了,没人写。

3 任意模长 DFT

3.1 引入

例题 1:单位根反演练习题。

\[\sum_{i=0}^{\infty}\binom{k}{in+m} = \sum_{i=0}^{k}\binom{k}{i}[n | i - m] = \sum_{i=0}^k\binom ki \frac 1n\sum_{j=0}^{n-1}\omega_{n}^{(i-m)j} \]

例题 2:【模板】Chirp Z-Transform

没讲做法。因为下面:

3.2 快速傅里叶变换

讲了一下 FFT 的基础分治部分。为下面做准备用的。
具体内容左转oiwiki

3.3 基于分治的快速傅里叶变换

类比 FFT,我们推广。一般地,我们需要计算 \(V(j) = \sum_{i=0}^{n-1} a_i \omega_n^{ij}\)。
仍然是拆 \(n\)。设 \(n\) 的 \(lpf\) 为 \(d\),我们有这样的做法的复杂度为 \(T(n) = dT(\frac nd) + O(dn)\)。容易看出来 \(T(n) = O(n^2)\)。
换种做法。

3.4 基于卷积的快速傅里叶变换

嗯,正主来了。Chirp Z-Transform,简称 czt。又称 Z 算法、bluestein 算法。
套路地,\(ij = \binom{i + j}2 - \binom i2 - \binom j2\)。

对于 \(DFT\),我们有式子

\[V(j) \omega_n^{\binom j2} = \sum_{i = 0}^{n-1} a_i\omega_n^{-\binom{i}{2}}\omega_n^{\binom{i + j}{2}} \]

等号右侧可以看作构造多项式 \(F(x) = \sum_{i=0}^{n-1} \omega_n^{-\binom{i}{2}}x^i\) 对 \(G(x) = \sum_{i=0}^{n-1}\omega_n^{a_i\binom{i}{2}}x^i\) 进行了减卷积后的结果。然后右侧的部分不需要模长的限制,可以直接做不限制长度的 DFT。

对于 \(IDFT\),我们有式子

\[V(j) \omega_n^{-\binom j2} = \sum_{i = 0}^{n-1} a_i\omega_n^{\binom{i}{2}}\omega_n^{-\binom{i + j}{2}} \]

直接减卷积即可。

总时间复杂度 \(O(n\log n)\)。

[这里放一个 czt 的代码,待更]

优化技巧:减卷积,模长需要开到 \(3n\) 级别,但是我们只需要中间 \(n\) 个元素。因此开 \(2n\) 长度就够了。

事实上,对于任意的非零数字 \(a\),我们都可以通过此类思路来计算 \(V(j) = \sum_{i=0}^{n-1} x_i a^{ij}\)(此处的上下标问题存疑)的值。我们只需要保证存在 \(a^{ij}\) 及其逆元即可。也就是说,对于任意复数该算法仍然适用。这点在下方有实际应用。

3.5 例题

略(
U498 新年的追逐战
CF901E Cyclic Cipher

4 多点求值

4.1 思路

大家都知道多项式多点求值是啥了,我就不说了。
由于 DFT 是将一系列点值 \(\omega_n^i\) 依次带入 \(f(x)\) 得到的长度为 \(n\) 的点值序列,因此在解决一类点值有特殊性质的多点求值问题时,我们可以将点值重排列后使用 czt 加速多点求值的进程。

例题:点值性质相关多点求值
给一个多项式 \(f(x)\)。\(q\) 次询问,第 \(i\) 次询问 \(f(q_i)\) 对 \(p = 998244353\) 取模后的结果。对于 \(q_i\),我们按 \(q_i = aq_{i-1} + b\) 生成。

做法:
假设我们已知了 \(f(x)\)。则我们可以在 \(O(|f(x)|)\) 的复杂度内找到唯一一个次数不超过 \(|f(x)|\) 的多项式 \(g(x)\) 满足 \(g(x) = f(x \times k)\),也可以在 \(O(|f(x)|\log |f(x)|)\) 的复杂度内找到唯一一个次数不超过 \(|f(x)|\) 的多项式 \(g(x)\) 满足 \(g(x) = f(x + k)\)。
归纳得到 \(q_i = q_0a^i + \sum_{j=0}^{i-1} ba^j\)。乘 \((a-1)\) 再加 \(b\),再除 \(q_0(a - 1) + b\)。推一下式子你会发现这个过程就是在将等号右侧消成单纯的 \(a^i\)。
因此构造 \(g(x) = f(\frac{x(q_0(a - 1) + b) - b}{a - 1})\)。我们只需要求得 \(\forall\ 1\le i \le q,\ g(a^i)\) 的点值,即求解 \(V(j) = \sum_{i=0}^{n-1} g_ia^{ij}\)。使用 czt 求解即可。
总时间复杂度 \(O((n+q)\log (n+q))\)。

例题:模数值域相关多点求值
给一个多项式 \(f(x)\)。\(q\) 次询问,每次询问给定 \(y\),求 \(f(y)\) 对 \(p = 786433(3\times 2^{18} + 1)\) 取模的值。\(0\le y < p\)。

做法:
考虑 DFT 是个啥。其实就是把 \(\omega_n^i\) 分别带入 \(f(x)\) 后的点值。
因为 \(p\) 是质数,因此其存在一个原根 \(g\)。同时不难发现 \(g^0,g^1,\cdots,g^{n-2}\) 在模 \(p\) 意义下互不相同,且组成了一个 \(1\sim p-1\) 的排列。
因此对 \(f\) 进行长度为 \(p-1\) 的定长度 DFT,取各位点值重排后即可得到所有可能的询问的答案。\(0\) 处的点值显然是 \(DFT(f)_ 0\)。
总时间复杂度 \(O(p\log p)\)。

4.2 推广

模数值域相关多点求值的模数其实可以被推广到 \(p = q^c\),其中 \(q\) 是奇质数。由于 \(p\) 存在原根 \(g\),因此任意 \((y,p) = 1\) 的询问 \(y\) 都可以经一次定长度 DFT 求得。现在考虑 \((y,p) \neq 1\) 的情况。显然 \(y^c \equiv 0 \pmod p\)。因此我们对于这类值只需要考虑小于 \(c\) 的项的贡献。因此直接暴力计算前 \(c\) 项即可。由于 \(c = O(\log p)\),因此这部分的总时间复杂度为 \(O(p \log p)\)。
因此总时间复杂度 \(O(p\log p)\)。

特殊的,对于 \(p = 2^c\) 的特殊模数,虽然没原根了,但是论文里说能证明 \(\pm 3^0, \pm 3^1, \cdots, \pm3^{2^{c-2}}\) 各项在模意义下彼此不同且与 \(2\) 互质。然后总时间复杂度仍然是 \(O(p\log p)\) 的。一堆奇怪的 corner cases,但是我不是很认为有人会闲的出自己论文里的这种分讨题。所以没有应用,就这样吧。

这同样表示我们可以通过 CRT 合并答案,因此能做到 \(O(p\log p)\) (复杂度存疑)进行任意模数多项式多点求值。

多项式多点插值

由于 DFT 能在多点求值问题上获得很好的应用,我们可以考虑将 IDFT 应用到多点求值的逆运算多点插值上。

5.1 连续点值的多点插值

对。给 \(0\le i < n\) 的 \(f(i)\),求 \(f(x)\) 的系数模 \(p\) 意义下的值。\(p\le 5\times 10^5\),\(p\) 是质数。

首先点值形式转下降幂形式。然后计算 \(0\sim p-1\) 处的点值。随后我们可以通过一次定长度 IDFT 取得一个多项式 \(g(x)\)。构造 \(h(x) \equiv [x = 0]\pmod p\),容易(?)发现 \(f(x) = g(x) + h(x)(a_0 - g(0))\)。
这里的多项式 \(h(x)\) 可以构造为 \(\sum_{i=0}^{p-2}x^i\)。总时间复杂度 \(O(p\log p)\)。

5.1 非连续点值的多点插值

对。给 \(0\le i < n\) 的 \(f(a_i)\),求 \(f(x)\) 的系数模 \(p\) 意义下的值。\(p\le 3\times 10^5\),\(p\) 是质数。

上来先插一波。

\[f(x) = \sum_{i=1}^n f(a_i) \prod_{j\neq i} \frac {x - a_j}{a_i - a_j} \]

设 \(g(x) = \prod{i=1}^n (x - a_i)\)。则 \(f(x) = g(x)\sum_{i=1}^n \frac{f(a_i)}{(x - a_i) g'(a_i)}\)。
首先是 \(g(x)\) 怎么求。你当然可以直接分治求,总时间复杂度是 \(O(n\log^2 n)\) 的。
但是两边除以 \(\prod -a_i\),再取 \(\ln\) 并求导后得到

\[\ln' \left(\frac{g(x)}{\prod -a_i}\right) = \sum_i \sum_j x^j \frac 1{a_i ^{j+1}} \]

复杂度 \(O(p \log p)\)。
然后是形如 \(F(x) = \sum_{i}\frac{v_i}{x - a_i}\) 的式子怎么求。这法子挺高级的。

\[\begin{aligned} F(x) & = \ \sum_{i}\frac{v_i}{x - a_i} \\ & = \ \sum_{i} v_i \sum_{j} x^j \left(\frac 1 {a_i}\right)^{j+1} \\ & = \ \sum_{j} x^j \sum_{i} v_i \left(\frac 1 {a_i}\right)^{j+1} \\ & = \ \sum_{j} x^j \sum_{i} \sum_{k} v_i \omega_{p-1}^{-k(j+1)}[\omega_{p-1}^k ={a_i}] \\ & = \ \sum_{j} x^j \sum_{k} \sum_{i} v_i[\omega_{p-1}^k ={a_i}] \omega_{p-1}^{-k(j+1)} \end{aligned}\]

设 \(G(k) = \sum_{i} v_i[\omega_{p-1}^k ={a_i}]\)。

\[\begin{aligned} F(x) & = \ \sum_{j} x^j \sum_{k} \sum_{i} v_i[\omega_{p-1}^k ={a_i}] \omega_{p-1}^{-k(j+1)} \\ & = \ \sum_{j} x^j \sum_{k} G(k) \omega_{p-1}^{-k(j+1)} \\ & = \ x^{-1} \sum_{j} x^j \sum_{k} G(k) \omega_{p-1}^{-kj} \end{aligned}\]

最后的式子的右半边两个求和号很符合 czt。算出 \(G\) 的点值后可以直接 \(O(p\log p)\)。

最后只需要将 \(g\) 求导后多点求值。

容易做到 \(O(p\log p)\)。

当然,上面的推导用到了 \(\frac 1x\) 和求导,因此在 \(a_i = 0\) 的情况下会产生错误。可以构造 \(g(x) = f(x + k)\) 使得点值非 \(0\),算完平移回来就好了。另一种方法是忽略 \(a_i = 0\) 处的点值,最后加入 \(k\prod_{a_i \neq 0}x - a_i\) 得到 \(f(0)\) 处的点值。第二种方法也被称作部分拉格朗日插值。

上面求 \(F(x)\) 的方法等价于给定一个列向量,左乘一个 Vandermonde 矩阵。因此我们也可以将上述算法用于解 Vandermonde 矩阵。
设我们需要解 \(F(i) = \sum_{j}a_jx_j^i\)。仍然设 \(G(x) = \sum_{i} [\omega_{p-1}^k = x_i]v_i\)。方程组转化成解 \(F(i) = \sum_{i} H(i) \omega_{p-1}^{ij}\)。不难发现这是 czt 形式。解决该问题需要 \(O(p \log p)\) 的时间复杂度。

03 年提出了一个 iczt。

感谢大家的阅读!没了!

标签:frac,DFT,sum,点值,阅读,多项式,随笔,omega
From: https://www.cnblogs.com/joke3579/p/paperessay221127.html

相关文章

  • [论文阅读] 颜色迁移-N维pdf迁移
    [论文阅读]颜色迁移-N维pdf迁移文章:N-DimensionalProbabilityDensityFunctionTransferanditsApplicationtoColourTransfer,[paper][code]1-算法原理简单......
  • C语言随笔5
    分支与循环(一)语句:由一个分号隔开的就是一个语句分支语句:if    switch/breakif 语法结构: 1.单分支语句   if(表达式)       语句; 2.......
  • 人工智能文献阅读1
     一、2021年人工智能领域科技发展综述                        总结:从各个的研发力度来阐述人工智能的发......
  • 阅读 OSTEP 第一、二章笔记
    虚拟化不等于虚拟机不等于Hypervisor,定义如下TheprimarywaytheOSdoesthisisthroughageneraltechniquethatwecallvirtualization.Thatis,theOStakesa......
  • 阅读开源代码时:idea中如何使用todo标记、活动模板 (史上最全)
    接下来,尼恩要带大家完成一个超级牛逼的大厂offer收割机项目——100Wqps三级组件实操,实操中,用到caffeine并且,尼恩要带大家穿透式、起底式的学习caffeine的......
  • 数据挖掘理论与算法,随笔1
    资源:b站本系列课程主要是启发为主,不会介绍很多很多的算法,适合初学者。一、学习资源书籍:国际会议:InternationalConferenceonDataMiningInternationalConference......
  • [Object-C语言随笔之一]Mac os 下搭建iOS开发环境
    ​​ 李华明Himi ​​​原创,转载务必在明显处注明 从这一章开始,Himi将一步一步的带大家走进Iphone4的开发,当然开发语言则不再是Java,而是Objective-C,简单来说是C的变种......
  • [Object-C语言随笔之四]创建视图并绘制简单图形
    ​​ 李华明Himi ​​​原创,转载务必在明显处注明这段时间N忙,没办法,创业公司,当然抽时间也仍然再自学ios~OK,基础的语言基础,我就不多说了,从今天开始直接写游戏开发部分......
  • 西门子PLC指令学习新想法随笔
    1,上升下降沿的用法:(用来记录运动中的物品数据,变化中的信息数据)(1),与set和rst配合,记录报警信息,放置用线圈记录时,报警信息同步消失;(2),与光电开关配合用于输送线体上的物品计数,一......
  • 【Cocos2d-X(1.x 2.x) 】iOS6与iphone5适相关设置随笔(解决第三方类库无法通过armv7s编
    本站文章均为​​ 李华明Himi ​​原创,转载务必在明显处注明:​一、 很多项目使用>=4.5version的Xcode无法,发现很多第三方库,比如SWavesSDK、AdmobSDK、91SDKMobage......