实序列离散傅里叶变换的快速算法

\newcommand{\dft}[1]{\text{DFT}\left[#1\right]}从多项式乘法到快速傅里叶变换中我们介绍了快速傅立叶变换的算法。由于在计算多项式乘法或者大整数乘法时,我们所使用的都是实数,而快速傅立叶变换并不仅适用于实序列,对于复序列也是可以计算的,这让我们考虑是否可以利用虚部来保存额外的信息而减少计算的次数。

事实上,对于实序列x(n),其DFT是共轭对称的,也就是X^*(n) = X(N - n),如果将两个实数序列x(n)y(n)合并称为新的z(n) = x(n) + i\cdot y(n),对其进行FFT,我们可以从中复原出分别对两个序列进行FFT的结果。也就是,可以通过两次FFT就进行多项式乘法的计算(原先是需要三次的)。

沿用其中的记号,先回忆一下离散傅里叶变换的公式。

考虑长度为N的序列x(n)及其N点离散傅里叶变换X(n) := \dft{x(n)}

\begin{equation} X(n) = \sum_{k=0}^{N - 1} x(k)\cdot\omega_N^{kn},~ n = 0, 1, \cdots, N - 1 \end{equation}

其中\omega_N = e^{-2\pi i/N},即N次单位根。

现在来考虑X(n)X(N - n)之间的关系,我们用X^*(n)\overline{X(n)}来表示X(n)的共轭

\begin{equation} X^*(n) = \sum_{k=0}^{N - 1} x^*(k)\cdot\omega_N^{-kn} = \sum_{k=0}^{N - 1} x^*(k)\cdot\omega_N^{k(N-n)} \end{equation}

这是因为 \overline{\omega_N^{kn}} = e^{2\pi ikn/N} = e^{-2\pi ik(N - n)/N} = \omega_N^{k(N - n)}

那么,对于实序列,x^*(n) = x(n),我们可以得到 X^*(n) = X(N - n)。也就是说,实序列的DFT是共轭对称的。

现在来考虑两个实序列x(n)y(n),其对应的DFT为X(n)Y(n),定义新的序列z(n) = x(n) + i\cdot y(n),其对应DFT为Z(n),那么

\begin{equation} Z(n) = \dft{x(n) + i\cdot y(n)} = X(n) + i\cdot Y(n) \end{equation}

注意这里X(n)Y(n)本身并不是一个实数,因此它们也不是Z(n)的实部和虚部,我们设Z(n)的实部和虚部分别为A(n)B(n)

\begin{equation} \begin{aligned} \overline{X(N - n) + i\cdot Y(N - n)} &= \overline{A(N - n) + i\cdot B(N - n)} \\ \Rightarrow X(n) - i\cdot Y(n) &= A(N - n) - i\cdot B(N - n) \end{aligned} \end{equation}

这样就可以得到

\begin{equation} \begin{aligned} 2X(n) &= [A(n) + A(N - n)] + i[B(n) - B(N - n)] \\ 2Y(n) &= [B(n) + B(N - n)] - i[A(n) - A(N - n)] \end{aligned} \end{equation}

这样,对于两个实序列的DFT可以只利用一次FFT以及一点额外的运算进行计算而不必分别计算两次FFT。

下面来考虑一个长度为 2N 的实序列 x(n)。我们将其按照下标的奇偶分解为两个长度为N的序列

\begin{equation} \begin{aligned} x_1(n) &= x(2n)\\ x_2(n) &= x(2n + 1) \end{aligned} \end{equation}

现在可以利用刚刚的方法通过一次N点的FFT来计算它们分别的DFT序列X_1(n)X_2(n),通过类似于FFT的办法将其合并

\begin{equation} \begin{aligned} X_1(n) + \omega_{2N}^n X_2(n) &= \sum_{k=0}^{N-1}x_1(k)\omega_N^{kn}+\omega_{2N}^n\sum_{k=0}^{N-1}x_2(k)\omega_N^{kn} \\ &= \sum_{k=0}^{N-1}x(2k)\omega_{2N}^{2kn}+\sum_{k=0}^{N-1}x(2k + 1)\omega_{2N}^{(2k+1)n} \\ &= X(n) \end{aligned} \end{equation}

这里只有前N项,对于其余项可以利用共轭对称这个性质直接计算。 对于大整数乘法,我们也可以利用两个实序列DFT的快速算法进行计算,这样只需要两次FFT即可,下面我提供一个多项式乘法的代码(UOJ#34),FFT的模板是很久以前写的就别看了,可以看看下面处理的部分。

(more…)

Read More

较短步数复原魔方的算法

这其实是我们上学期徐 BOSS 程设课的一个大作业,我本来很早就想把它记录下来了,可是实在太懒了加上拖延症严重于是就拖到了这个学期。本来作业是要求复原魔方就好了,但是后来发现没有图形界面的话似乎调试很困难,就顺便写了一个 UI,刚好还学了一下 OpenGL。最开始我是写了能够找到最优解的 Krof 算法,但是估计是我写得比较丑再加上状态空间实在太大了(大约在搜索到第 15 步还是 16 步的时候状态数就到百亿级别了,而我 8 个线程一秒只能搜几千万个状态好像),根本没有办法在可以接受的时间内得到解。于是退而求其次,改为去搜索一个比较优的解,最后基本上一秒以内就可以得到解了。

整个工程的代码我放在了我的 Github 上:https://github.com/miskcoo/rubik-cube

魔方的表示

如何表示一个魔方呢?一个很直接的方法就是用六个面一共 54 个色块来表示它,我在最开始写的时候就是用这种方法来表示一个魔方的,但是这样存在着很多的冗余。事实上,我们可以按照魔方的各个棱块和角块的位置和方向来存储。对于魔方的一个角块来说,我们只需要存储这个块在 8 个角的哪一个地方以及它的方向如何(关于方向接下来会进行规定)。对于棱块来说也是相同的,只需要位置和方向就可以决定了。

角块方向的表示

首先对于一个角块来说它可能会有 3 个不同的方向(如下图,暂不考虑这个魔方是否合法,仅关注角块)

corner-0corner-2

想象一下,从左到右相当于单独把面对我们的这个角块(按照从魔方中心指向该角块顶点的向量)旋转了两次。我们可以按照这样来规定角块的方向:首先将对应角块像上图中一样面对我们,并且以绿色或蓝色作为顶面。现在以绿色或蓝色作为参照色(因为一个角块有且仅有这两种颜色之一)。现在来规定角块的方向:

  • 如果参照色位置在如同第一张图一样的位置,就规定方向为 0
  • 如果参照色位置在如同第二张图一样的位置,就规定方向为 1
  • 否则就规定方向为 2

棱块方向的表示

对于棱块方向的规定就稍微有些复杂了,因为没有办法像角块这样找到对立的两个面使得每个棱都包含这两种颜色之一。那么对于棱块的参照色我们规定:如果一个棱存在蓝色或绿色,则该颜色为参照色,否则白色或者黄色为参照色。规定了参照色之后就可以规定棱块的方向了:

  • 如果棱块位于顶层或者底层(也就是绿色或者蓝色面所在的层),那么当参照色在这两个面中的一个时方向为 0,否则为 1
  • 如果棱块位于中间层,那么当参照色在白色或者黄色面时方向为 0,否则为 1

edge-0edge-1

在第一张图中,两个棱块的方向都是 0。这是因为对于红白色棱块,白色为参照色并且其位于顶层绿色面;对于绿白色棱块,绿色为参照色并且位于中层白色面。

在第二张图中,两个棱块的方向都是 1。原因与之前的刚好相反。

Kociemba 算法

首先我们规定一种魔方的中间状态:仅仅通过 U、D、L2、R2、F2、B2 这几种旋转就可以将魔方复原的状态。也就是不允许 90^\circ 地旋转魔方的前后左右这四个面。这种状态有什么特点呢?

首先很明显可以观察到按照我们刚刚规定的方向,旋转顶面和底面都是不会改变棱块和角块的方向的,并且仅按照 180^\circ 来旋转周围四个面也是不会改变棱块和角块的方向。那也就是说,在这种中间状态下,魔方的棱块和角块的方向应该都是 0,这是因为目标状态的魔方如此。除此之外,魔方中间层的四个棱块必然也在中间层。

可以简单地计算出这种中间状态一共只有 8! \times 8! \times 4! = 3.9 \times 10^{10} 种可能,事实证明,在不超过 18 步以内就可以从这种状态复原魔方。我们用 IDA* 来进行搜索,当我们设计启发函数的时候分别针对角块(一共 8! 种可能)、顶层和底层的棱块(一共 8! 种可能)、中间层的棱块(一共 4! 种可能)计算复原它们所需要的最少步数(这可以直接进行搜索),之后的启发函数就是这三者的最大值。

剩下的问题就是如何将魔方旋转到中间状态,我们同样采用 IDA* 来进行搜索。对于启发函数,分别针对角块的方向(一共 3^8 种可能)、棱块的方向(一共 2^8 种可能)、中间层棱块的位置及方向(一共 A^4_{12}\times 2^4 种可能)计算复原它们所需要的最少步数。经过计算,不超过 12 步就可以将一个魔方旋转到中间状态。这样总共在不超过 30 步的情况下就可以复原一个魔方了!

刚刚的状态表示方式在这个算法里显示出了很大的优势,在进行 Hash 的时候可以直接使用这个状态表示,而不必从其它表示方法(例如存储 54 个面的方法)计算出棱和角的方向。事实证明,在通常情况下,这个算法很快就能够复原一个魔方。

Krof 算法

这个算法就是我前面说的搜索最优解的算法,思路和 Kociemba 算法几乎一样,不过 Krof 算法没有划分阶段,启发函数是分别对于角块的方向和位置(一共 8!\times3^7 种可能)、两组每组 6 个梭块的位置和方向(每组都是 A^6_{12}\times2^6 种可能)计算复原它们所需要的最少步数。由于这个启发函数的计算也是十分耗时的,所以可以先计算一次之后记录下来之后直接读取就可以了,而前面 Kociemba 算法的启发函数可以现场计算。

Read More

CodeChef上一道有趣的网络流题目

这是一道2016年二月CodeChef上的一道题目,叫做Call Center Schedule

题目大意是这样的:

一个电话公司有 P 个员工,一个星期有 D 天,每天需要工作 H 个小时。每个员工在每个小时只能做参加会议和与客户通话两件事情中的一件,或者不做。

现在要为安排一个和客户通话的时间表,要求满足

  • 每个人每天最多花费 N 个小时在参加会议和通话上
  • 第 i 个人每周最多花费 Li 个小时和客户通话
  • 每个人每天在 LTbegin 到 LTend 的时间内至少要有 1 小时没有被安排任务
  • 对于第 i 天第 j 个小时,恰好有 Ri, j 个人可以和客户通话

另外,开会的时间表已经安排完成,如果 Fk, i, j 是 0 则表示第 k 个人在第 i 天的第 j 个小时开会,如果是 1 则表示他可以被安排去和客户通话。

时限是 2s,最多 5 组数据,并且 1 \leq N \leq H \leq 70, 1 \leq D, P \leq 70, 1 \leq L_i \leq ND, 0 \leq R_{i, j} \leq 15, F_{k, i, j} \in \{0, 1\}, 1 \leq LT_{begin} \leq LT_{end} \leq N

(more…)

Read More

Tarjan算法寻找有向图的强连通分量

在图论中,一个有向图被成为是强连通的(strongly connected)当且仅当每一对不相同结点 u 和 v 间既存在从 u 到 v 的路径也存在从 v 到 u 的路径。有向图的极大强连通子图(这里指点数极大)被称为强连通分量(strongly connected component)。

tarjan-graph

比如说这个有向图中,点 1, 2, 4, 5, 6, 7, 8 和相应边组成的子图就是一个强连通分量,另外点 3, 9 单独构成强连通分量。

Tarjan 算法是由 Robert Tarjan 提出的用于寻找有向图的强连通分量的算法。它可以在\mathcal O(|V|+|E|)的时间内得出结果。

Tarjan 算法主要是利用 DFS 来寻找强连通分量的。在介绍该算法之前,我们先来介绍一下搜索树。先前那个有向图的搜索树是这样的:

search-tree

有向图的搜索树主要有 4 种边(这张图只有 3 种),其中用实线画出来的是树边(tree edge),每次搜索找到一个还没有访问过的结点的时候就形成了一条树边。用长虚线画出来的是反祖边(back edge),也被叫做回边。用短虚线画出来的是横叉边(cross edge),它主要是在搜索的时候遇到了一个已经访问过的结点,但是这个结点并不是当前节点的祖先时形成的。除此之外,像从结点1到结点6这样的边叫做前向边(forward edge),它是在搜索的时候遇到子树中的结点的时候形成的。

(more…)

Read More

CodeChef February 2016. Weird Sum 题解

这是 CodChef February Challenge 2016 上的一题,题目链接在这里:WRDSUM

给定一个整数 n(n \geq 2),考虑其质因子分解 n = p_1^{k_1}p_2^{k_2}\cdots p_r^{k_r}

g = \text{gcd}(k_1, k_2, \cdots, k_r)m_i = \frac{k_i}{g}

定义函数 F(n) = p_1^{m_1}p_2^{m_2}\cdots p_r^{m_r}

现在给定一个整数 N(N \leq 10^{2016}),你需要计算

 S = F(2) + F(3) + \cdots + F(N)

答案对 998244353 取模

(more…)

Read More

特殊多项式在整点上的线性插值方法

首先我来解释一下题目是什么东西,假如给你一个 k 次多项式 P(x),并且已知了 P(0), P(1), \cdots, P(k),要求计算出 P(n), n \in \mathbb N 的值,由于已经知道了 k + 1 个点值,那么要插值出系数是可以在 \mathcal O(k^2) 完成,如果使用 FFT 的话是可以在 \mathcal O(k \log^3 k) 的时间内完成的[1],但是对于这个特别的问题,我们有一个 \mathcal O(k) 的算法!这好像是杜教先提出来的,然后还有这篇文章

(more…)

Read More

Codeforces 上的一些组合计数问题

Codeforces 15E. Triangles

Problem:给出一个有规律的金字塔(如图),黑色的边是可以走路径,要求求出包含 H 点的简单回路(不经过一个顶点两次)的条数,并且这个回路围住的部分不可以包含灰色的三角形,答案对 10^9+9 取模,其中 n \leq 10^6,并且一定是一个偶数

5945210be972aa5fe947f3b8a3a0378a4cade844

Example:当 n = 2 时,一共有 10 种方案,其中 5 种是这样,另外的是顺时针顺序地走

5E-triangles

Solution:首先第 1 层是最特殊的,因为这里没有灰色的三角形,然后注意到灰色的部分把整个金字塔分成了左部和右部

5E-triangles-level1

路径现在可以分成三个部分

  • 在第一层移动,这一共 10 种(见样例)
  • 进入左部或右部后直接回到 H
  • 进入左部(右部)后经过 E 点进入右部(左部)

对于第二种情况,在第一层的部分走法可以是(只有逆时针部分,实际上有 8 种)

  • A \rightarrow B \rightarrow D \rightarrow \text{LeftPart} \rightarrow E \rightarrow C \rightarrow A
  • A \rightarrow B \rightarrow D \rightarrow \text{LeftPart} \rightarrow E \rightarrow F \rightarrow C \rightarrow A
  • A \rightarrow B \rightarrow E \rightarrow \text{RightPart} \rightarrow F \rightarrow C \rightarrow A
  • A \rightarrow B \rightarrow D \rightarrow E \rightarrow \text{RightPart} \rightarrow F \rightarrow C \rightarrow A

对于第三种情况,在第一层(原来图中的两层现在算一层,也就是下面的 n 是题目中 n\frac{1}{2})的部分走法有 2 种,左部到右部和右部到左部

现在要统计的就是在左部和右部内有多少种方案,因为左部右部对称,所以答案是一样的,我们来考虑左部

S_n 表示在有 n 层的金字塔中,走左部的方案数,因为走进左部后每一层都会遇到凹进去的部分,这部分只要枚举一下走了多远可以统计出第 n 层一共有 2\sum_{i=1}^{n-2} 2^i + 1 = 2^n - 3 这么多的情况

然后要走到第 n 层后就是 \prod_{i=2}^{n} (2^i - 3) 种情况,然后走回去的路径有 4 种选择,因此

 S_n = \sum_{i=2}^{n} 4\prod_{j=2}^i (2^j - 3)

答案是 2S_n^2 + 8S_n + 10 (more…)

Read More

牛顿迭代法在多项式运算的应用

总算是快把 FFT 和生成函数的各种东西补了好多,膜拜策爷的论文和 Picks 的博客 QAQ

这篇文章大概就是说如何用牛顿迭代法来解满足 G(F(z)) \equiv 0 \pmod {z^n}F(z)

然后这东西可以比较方便地计算 \sqrt{A(z)}e^{A(z)},也就是多项式开根、求指数对数之类鬼畜的东西,在生成函数计数中十分有用

顺便一提,这里说的“多项式”实际上你可以直接理解为生成函数或者形式幂级数

(more…)

Read More

平面图转对偶图及点定位

平面图(planar graph)就是一个可以在平面上画出来的图,并且所有的边只在顶点处相交。平面图的对偶图(dual graph)是将这个平面的每个区域看成点,原图每一条边所属的两个相邻的区域对应在对偶图中的点有连边

比如下面这张图(来源于 Wikipedia)

300px-Duals_graphs.svg

蓝色的部分就是一个平面图,红色的就是它的对偶图

平面图的点定位(point location)就是找出一个点属于这个图的哪个区域

这篇文章介绍如何将一个平面图转换为其对偶图,和基于扫描线的离线的点定位算法(对着代码看了一天终于会了QAQ)

(more…)

Read More