NonTrivialMIPS - 十级双发射顺序MIPS32处理器

前一段时间我和几个同学参加了名为“龙芯杯”的比赛,这个比赛是自己设计一个CPU,在其上设计SoC,运行操作系统等… 是一个系统类的比赛。我主要负责写CPU除Cache外的部分,我们最终设计了一个有十级流水的双发射顺序执行的MIPS32处理器。在龙芯的实验板上达到了123MHz的主频,同时还具有较高的IPC。代码现在已经在 github 上开源:https://github.com/trivialmips/nontrivial-mips

我来说说这个流水线的设计,除去各种数据旁路,流水线的架构大约如下图NonTrivialMIPS CPU

这十级流水基本上是从标准的五级流水扩展而来,其中有3级的取指,2级的发射,3级的访存。其实我们最开始计划的是双发射乱序处理器,并且之前在组成原理课上有一部分顺序双发射的经验,当时有一个很麻烦的地方就是取指阶段,这里需要根据到底发射一条指令还是两条指令控制PC的增加,但是在决定下一个PC值时所需要的指令还没有取到,这就需要一部分预测和旁路,导致取指阶段很复杂并且容易写错。具体可以参考我当时的计原报告

  Report of TrivialMIPS (588.0 KiB, 835 hits)

这次我基本上重写了整个CPU,在取指阶段增加了一个FIFO,把取到的指令放入FIFO。之后在发射的时候从FIFO读取指令。这样一定程度上将取指和之后的执行等阶段进行了解耦。取指有3级流水基本上是由于指令Cache的需要。由于分支要在较后才能解析,CPU必然需要一个分支预测器,我实现的是2bit的分支预测。假设分支预测完全正确的情况下,可以保证放入FIFO的指令流就是执行所需要的指令流(当然这里也处理了延迟槽)。同时,为了支持双发射,我们的指令Cache一次可以取出8字节的数据,如果其后的8字节仍然在同一个Cache行内,也会一并取出。这是为了在预测到分支指令的之后能够尽早将延迟槽放入FIFO而不需要进行一次多余的取指。也就是,我们的FIFO在一个周期最多可以放入3条指令。

在之后发射阶段,为了频率的提高,我们将读操作数单独作为一个流水段。虽然这样会造成分支预测失败的损失增大进而减小IPC,但是和频率的提升相比这部分的损失还是比较小的。双发射基本上是对称的,只要不存在数据冲突和结构冲突(例如两条除法指令或者两条访存指令)都可以同时发射。当然,为了设计简便,我额外增加了一个限制,也就是分支指令和延迟槽必须一起发射。发射之后的读操作数和执行就没有太多技术含量,和通常的五级流水线相同。我们的异常处理和CP0的读写都是在执行之后的一个阶段,也就是访存的第一个阶段进行的。

在这里有一个很严重的问题就是,由于数据Cache的流水太常,导致访存相关(即访存后立即使用数据)的暂停很高,这极大影响了性能。为了缓解这个问题,我们允许部分指令的读操作数和执行延迟到访存阶段进行。这样就可以很大地减小访存暂停。当然不是所有指令都能做到这件事,首先这需要满足一些条件。其一就是这类指令不能是访存指令,其二是这类指令不能在执行阶段出现异常(取指的TLB miss还是可以的,这在它执行前就可以判断)。当然,这两个条件其实不是很严格,很大一部分程序访存后基本上是跟着一些ALU指令或者分支指令。在增加这个优化之后,整个处理器的IPC大约有了15%的提升,同时频率仍然保持。其实我们或许还可以进一步优化,对于非分支指令还可以将其延迟到访存的第二个阶段读操作数,这样就完全不会有访存相关。至于分支指令,我们必须要能够在分支预测失败时撤销其后的指令,流水线上能够影响处理器状态的第一个阶段就是访存的第一个阶段,因此分支指令最迟必须要在访存第二个阶段执行完成。不过我并没有实现这个优化。另外还有一个较小的优化,对于存储指令的数据,可以在执行阶段再读取。

除了执行操作系统所需要的指令以外,我们还做了一定的扩展,例如实现了一个32位的FPU以及一个硬件AES加速单元同时移植到OpenSSL上。MIPS32在早期版本对单精度浮点的规范非常奇怪,有两条指令叫做ldc1和sdc1,它允许同时两个FPU寄存器和内存交换数据,也就是一次访存的数据宽度是64位。我实现的方法是在发射阶段将其拆开成为两条32位的访存指令,同时不允许中断打断这两条指令。至于AES加速,实际上是在CP2上实现的。只用到了mfc2和mtc2两条指令,AES单元具有自己的寄存器,通过这两条指令和CPU的寄存器进行数据的交互以及控制。

至于操作系统和SoC的内容,这就不是我负责的范围了,具体可以看看我们最后的报告 https://github.com/trivialmips/nontrivial-mips/tree/master/report.

其实实现的过程中还有一个乱序的版本,不过最后因为频率没法提升最终没有使用。

球坐标下Laplace算子的公式

前几天给人找了几道微积分题目,里面有一个是计算球坐标下面的Laplace算子的表达式。这当然可以暴力用复合函数的求导来算,但是可能比较复杂。事实上可以有一种积分的办法来算出这个结果。

我们知道,在标准坐标下,Laplace算子的形式为

 \Delta u = \frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2} + \frac{\partial^2u}{\partial z^2}

现在计算在球坐标

\begin{cases}
x = r\sin\theta\cos\varphi \\
y = r\sin\theta\sin\varphi \\
z = r\cos\theta
\end{cases}

下Laplace算子的表达式。首先计算其Jacobi矩阵

 \frac{\partial(x, y, z)}{\partial(r, \theta, \varphi)} = \begin{pmatrix}
\sin\theta\cos\varphi & \sin\theta\sin\varphi & \cos\theta \\
r\cos\theta\cos\varphi & r\cos\theta\sin\varphi & -r\sin\theta \\
-r\sin\theta\sin\varphi & r\sin\theta\cos\varphi & 0
\end{pmatrix}^T

对应逆矩阵为

 \frac{\partial(r, \theta, \varphi)}{\partial(x, y, z)} = \begin{pmatrix}
\sin \theta\cos\varphi & \sin\theta\sin\varphi & \cos\theta \\
\frac{1}{r}\cos\theta\cos\varphi & \frac{1}{r}\cos\theta\sin\varphi & -\frac{1}{r}\sin\theta \\
-\frac{\sin\varphi}{r\sin\theta} & \frac{\cos\varphi}{r\sin\theta} & 0
\end{pmatrix}

注意到这是一个正交矩阵。

我们利用Laplace算子的积分表达式,考虑

 \begin{aligned}\frac{ \partial{u(r, \theta, \varphi)}}{\partial(x, y, z)} \cdot \frac{\partial{v(r, \theta, \varphi)}}{\partial(x, y, z)} &= \frac{\partial{u(r, \theta, \varphi)}}{\partial(r, \theta, \varphi)}\frac{\partial(r, \theta, \varphi)}{\partial(x, y, z)}^T\frac{\partial(r, \theta, \varphi)}{\partial(x, y, z)}\frac{\partial{v(r, \theta, \varphi)}}{\partial(r, \theta, \varphi)} \\ &= u_rv_r + \frac{1}{r^2}u_\theta v_\theta + \frac{1}{r^2\sin^2\theta}u_\varphi v_\varphi \end{aligned}

对于任意的$v \in C_c^\infty$,利用分部积分

 \begin{aligned}
\int \Delta u v dxdydz &= -\int \nabla u \cdot \nabla v dxdydz \\
&= -\int \left( u_r v_r + \frac{1}{r^2} u_\theta v_\theta + \frac{1}{r^2\sin^2\theta}u_\varphi v_\varphi \right) r^2\sin\theta drd\theta d\varphi \\
&= \int \left( \left (u_r r^2\sin\theta \right )_r + ( u_\theta \sin\theta )_\theta + u_{\varphi\varphi} \right) v drd\theta d\varphi \\
\end{aligned}

再根据

 \int \Delta u v dxdydz = \int \Delta u v r^2\sin\theta drd\theta d\varphi

这样就可以得到

 \begin{aligned}
\Delta u &= \frac{1}{r^2\sin\theta}\left( \left (u_r r^2\sin\theta \right )_r + ( u_\theta \sin\theta )_\theta + u_{\varphi\varphi} \right) \\
&= \frac{1}{r} \frac{\partial^2}{\partial r^2}(ru) + \frac{1}{r^2\sin \theta} \frac{\partial}{\partial \theta}\left( \sin \theta \frac{\partial u}{\partial \theta} \right) + \frac{1}{r^2\sin^2 \theta} \frac{\partial^2 u}{\partial \varphi^2}
\end{aligned}

TrivialDB - 一个简单的SQL引擎

TrivialDB是一个简单的数据库管理系统,我们实现了大部分常见的SQL语句和类型。同时支持多表连接、复杂表达式运算、多主键约束、外键约束、CHECK约束、UNIQUE和DEFAULT约束、聚集查询、利用B+树索引的查询优化,同时,我们支持任意长度的VARCHAR类型。

项目地址是 https://github.com/miskcoo/TrivialDB,这是这学期一门课程的大作业。

编译及运行

你需要有支持C++11特性的编译器,以及Bison和Flex两个库。本项目通过CMake来构建,在根目录运行

进行编译选项的设置,之后运行

进行项目的编译,编译后的可执行程序在build/trivial_db目录下。

编译后可以选择在testcase目录下运行python3 run_test.py运行测试程序。

系统功能

数据类型

数据库支持的基本类型有:

  • 整型(INT)
  • 浮点型(FLOAT)
  • 字符串型(VARCHAR)
  • 日期型(DATE),日期格式YYYY-MM-dd

日期类型的字面值和字符串相同,在实现中如果必要可以转换为字符串。

SQL语句

我们支持的SQL语句一共有如下几种

  • 插入语句:INSERT INTO ... VALUES ...
  • 删除语句:DELETE FROM ... WHERE ...
  • 查询语句:SELECT ... FROM ... WHERE ...
  • 更新语句:UPDATE ... SET ... WHERE ...
  • 创建数据库:CREATE DATABASE ...
  • 删除数据库:DROP DATABASE ...
  • 切换数据库:USE ...
  • 显示数据库信息:SHOW DATABASE ...
  • 创建表:CREATE TABLE ...
  • 删除表:DROP TABLE ...
  • 显示表信息:SHOW TABLE ...
  • 创建索引:CREATE INDEX ...
  • 删除索引:DROP INDEX ...

复杂表达式处理

表达式大致可以分为两种:算术表达式和条件表达式。由于采用Bison进行解析,可以支持任意深度嵌套的复杂表达式。我们所支持的基本运算主要如下

  • 四则运算,针对整数和浮点数进行。
  • 比较运算符,即<=, <, =, >, >=, <>。
  • 模糊匹配运算符,即LIKE,其实现采用C++11的正则表达式库。
  • 范围匹配运算符,即IN,可以在表的CHECK约束中以及WHERE子句中使用。
  • 空值判定运算符,即IS NULL和IS NOT NULL两种。
  • 逻辑运算,包含NOT、AND和OR三种。

以下是一些复杂表达式运算的例子

聚集查询

我们实现了五种聚集查询函数COUNT、SUM、AVG、MIN和MAX。其中COUNT不支持DISTINCT关键字。例如

属性完整性约束

我们支持多种属性完整性约束,分别是

  • 主键约束。一个表可以有多个列联合起来作为主键,只有在所有主键都相同时才认为两条记录有冲突,即这种情况下主键是一个元组。
  • 外键约束,每个域都可以有外键约束,引用另外一个表的主键。
  • UNIQUE约束,该约束限制某一列的值不能重复。
  • NOT NULL约束,该约束限制某一列不能有空值。
  • DEFAULT约束,该约束可以在INSERT语句不指定值是给某列赋予一个默认值。
  • CHECK约束,该约束可以对表中元素的值添加条件表达式的检查。

下面是一个简单的例子,注意如果在多个列都指定了PRIMARY KEY,那么就认为主键是一个元组,而不是有多个主键。例如Infos表的主键为(PersonID, InfoID)。

多表连接查询

在SELECT语句中,我们支持任意多表的连接操作,例如

并且,对于多个表的连接中形如A.Col1 = B.Col2的条件,那么如果这两个列的某一个拥有索引,会利用索引进行查询优化。例如如下查询就可以优化

具体的优化方法以及何种查询可以优化见文档中"查询优化"部分。

表别名

我们在多表连接查询时支持通过别名(alias)的方式对一个表进行连接,例如

FPGA Console——硬件实现的VT220兼容终端

这是和@HarryChen在这学期数字逻辑课上做的项目。主要是利用FPGA模拟一个终端。项目放在https://github.com/miskcoo/fpga-virtual-console上。

现在大部分的操作系统都是内置有虚拟终端,就像Linux下的xterm、gnome-terminal。但是,在最初控制台和计算机是相互独立的。控制台一般是由键盘和屏幕组成,键盘将用户输入传入计算机,计算机将输出发送显示给用户。这整个通信的过程是有一套标准来规定的,例如VT100就是一个早期的标准,更进一步有VT220、xterm-256color等。由于计算机回显的数据并不是单纯地整个屏幕的像素数据,而是更高层次的类似控制命令的内容,例如表示光标移动到某个位置,在某个位置插入某个字符,删除某一行的字符等。这些标准主要是定义了计算机发送回的这些控制命令的行为,和键盘发送给计算机的控制命令的行为。

我们做的这个项目就是一个物理上的终端!通过FPGA接受键盘的输入,将输入转化为控制命令通过串口输出给计算机。同时也通过串口接受计算机传回的控制命令名且解析、执行,修改对应位置的字符,再将字符进行渲染通过VGA输出到屏幕。

主要支持的标准是VT220,以及xterm-256color下大部分影响最终显示效果的指令,对于大部分linux下的程序的输出都可以正确地解析。同时串口速度提升到了3M可以支持以大约170fps的帧率播放黑白字符视屏,和约30fps的帧率播放彩色字符视屏。简单来说,就是正常终端可以做的事情这个物理终端都可以做到。

比如下面就是一个调用笔记本上摄像头的例子(我真的不是直接把VGA连到电脑上了),这时候串口速度还没那么快所以看起来有点卡……

fpga-camera

具体的细节我就不打算说了,文档里都写得非常明白。

https://github.com/miskcoo/fpga-virtual-console/raw/v1.0/doc/report.pdf

Pólya’s Urn

Pólya’s Urn是这样一个模型:假设现在有一个盒子,里面有r个红球和g个绿球,现在不断地随机从盒子取出一个球,之后再放回去,外加放入额外的c个和取出的球颜色相同的球。现在问,盒子中绿色球的比例在最后是否会趋向与某个分布?

实际上最后这个比例会是一个参数是g/cr/c的Beta分布。

首先我们设S_n = r + g + nc是第n次操作后盒子中所有球的数目。再设G_n, R_n分别是第n次操作后盒子中绿色球、红色球的数目。以及令X_n = G_n / S_n, \mathcal F_n = \sigma(G_k, S_k, 0 \leq k \leq n)。计算条件期望

 \begin{aligned}
\mathbb E[G_{n+1} | \mathcal F_n] &= \frac{G_n}{S_n}(G_n + c) + \frac{R_n}{S_n}G_n \\
&= (c + S_n) X_n = S_{n+1} X_n
\end{aligned}

可以看出\mathbb E[X_{n+1}|\mathcal F_n] = X_n\{X_n\}是一个鞅。再由于0\leq X_n \leq 1,我们知道X_n几乎处处收敛到一个随机变量X。现在考虑其分布。

首先可以计算出在第n步及以前,有m \leq n次取出绿色球的概率,并且写为紧凑的形式如下

 \begin{aligned}
&\mathbb P\left[ X_n = \frac{mc + g}{nc + g + r} \right] \\&= \frac{n!}{m!(n-m)!}\cdot \frac{\displaystyle\prod_{k=0}^{m-1}(k + g/c) \prod_{k=0}^{n - m - 1}(k + r/c)}{\displaystyle\prod_{k=0}^{n - 1}(\frac{g+r}{c} + k)}\\&= \frac{\Gamma(n + 1)}{\Gamma(m + 1)\Gamma(n - m + 1)}\cdot \frac{\Gamma(m + g/c)\Gamma(n - m + r/c)}{\Gamma(n + (g + r)/c)\text{B}(r/c, g/c)}
\end{aligned}

其中\Gamma\text{B}分别是Gamma函数和Beta函数。

接下来利用Stirling公式\Gamma(x) = \sqrt{2\pi/x} (x/e)^x (1 + O(1/x)),可以得到

 \frac{\Gamma(x + a)}{\Gamma(x + b)} = x^{a - b} \left( 1 + O\left( \frac{1}{x} \right) \right)\text{ as } x \to \infty

那么对于固定的\alpha, \beta, m满足0 < \alpha < \beta < 1以及\alpha n < m < \beta n,可以有

 \begin{aligned}
&\mathbb P\left[ X_n = \frac{mc + g}{nc + g + r} \right] \\&= \frac{1}{\text{B}(r/c, g/c)} n^{1 - (g + r)/c} (n - m)^{r/c - 1} m^{g/c - 1} \left( 1 + O\left( \frac{1}{n} \right) \right)\\
&= \frac{1}{n \text{B}(r/c, g/c)} \left (1 - \frac{m}{n} \right )^{r/c - 1} \left( \frac{m}{n} \right)^{g/c - 1} \left( 1 + O\left( \frac{1}{n} \right) \right)
\end{aligned}

对于所有满足\alpha n < m < \beta nm求和

 \lim_{n\to \infty} \mathbb P[\alpha \leq X_n \leq \beta] = \frac{1}{\text{B}(r/c, g/c)}\int_{\alpha}^\beta (1 - x)^{r/c - 1}x^{g/c - 1} dx

最后可以得到X的分布是\text{Beta}(g/c, r/c)

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

\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的模板是很久以前写的就别看了,可以看看下面处理的部分。

Read More

程序设计训练 - 数独、国际跳棋、人物信息检索系统

这次程序设计训练小学期一共有三个大作业,分别是数独、在线国际跳棋和人物信息检索系统。这里是简单的一些介绍,代码放在 https://github.com/miskcoo/programming-training

Week 1: Sudoku Game

简介

Sudoku 是一款利用 Qt 实现的数独游戏,提供了多达 10 个难度的关卡选择,同时还有丰富的功能来帮助玩家更加高效地求解数独问题,例如候选数、高亮相同数字、高亮选中的行列、撤销当前操作以及提示等功能。玩家还可以手动输入数独题目利用 Sudoku 帮助求解。

除了传统 9x9 的数独游戏以外,还提供了更高难度的 16x16 的数独游戏。

基本功能

Sudoku提供了多个方便的按钮:

  • 新游戏:玩家可以开始一局新的游戏。
  • 重玩:玩家可以重新开始本局游戏。
  • 暂停:玩家可以暂停该局游戏(即暂停计时)。
  • 提示:如果当前已经确定的数都是正确的,玩家将会得到一个未填空格的正确数字;如果当前已经确定的数和答案矛盾,导致整个数独无解,那么所有与答案矛盾的数字将会被粗体标出。
  • 清除:清除当前选中格子的所有数字。
  • 撤销:撤销前一步的操作,以及取消撤销(最多可支持 50 步撤销)。

同时可以通过菜单来实现多达 10 种难度的游戏选择,可以求解任意用户输入的数独问题。

玩家在空格中填入的数字分为确定数(采用正常大小的字体表示)以及候选数(采用小号字体表示)。确定数只能存在一个,候选数可以存在多个,并且确定数和候选数不能同时存在。

玩家可以通过多种方式进行格子的选择:

  • 鼠标直接点击格子进行选中。
  • 键盘方向键来进行当前选中格子的切换。
  • Tab键快速切换到下一个格子。

玩家可以用鼠标右键点击格子来对其进行标记。

玩家在选中一个格子之后,如果该格子的数字已经确定,那么所有与该数相同的数将会被加粗显示,以方便确认是否满足数独的条件。同时,所有与该格在相同行或列的格子都会被高亮显示。此外,无论该格子数字确定与否,右侧数字列表中该格子的所有数都会被加粗。

玩家可以通过多种方式在空格中填入数字:

  • 键盘直接输入数字将填入对应的确定数,如果 Ctrl 键被按下,那么填入的将会是候选数。
  • 鼠标点击右侧数字列表填入。使用鼠标左键点击将会填入确定数,右键点击则会填入候选数。

玩家也可以通过多种方式在空格中填入数字:

  • 键盘删除键删除格子中的一个数字。
  • 鼠标点击右侧数字列表已经选中的数进行删除。
  • “清除”功能键来清除该格子所有的数字。

Read More

Vim里一些好用的插件

我主要是想要把一些用过的好用的插件记下来,万一哪天换了系统还可以在这里找回来。

另外这篇文章会不断更新,如果我发现新的有意思的插件也会加进来。

YouCompleteMe

绝对好用的自动补全插件,支持C/C++等各类语言的自动补全

cpp-demo-of-youcompleteme

具体请戳:http://blog.miskcoo.com/2015/12/vim-plugin-youcompleteme

Jedi-Vim

一个很好用的Python的自动补全插件

具体请戳:https://github.com/davidhalter/jedi-vim

Vim-Latex

这是一款帮助编辑Latex的插件,里面集成了很多环境的快捷输入方式,例如输入EFI就会自动变成如下

另外还能够用 \ll 编译并且 \lv 方便地查看错误信息等。

我比较喜欢的是类似 \alpha 的命令可以用 `a 来快速输入,还有就是 \frac{}{} 这样的命令可以用 `/ 或者输入 frac + <F7> 来补全在输入数学公式的时候特别方便。

具体请戳:http://vim-latex.sourceforge.net

Markdown-Preview

一个国人写的插件,支持Markdown的同步预览,感觉挺好用的。

具体请戳:https://github.com/iamcco/markdown-preview.vim

较短步数复原魔方的算法

这其实是我们上学期徐 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 算法的启发函数可以现场计算。

三维空间中的旋转:旋转矩阵、欧拉角

考虑这样一个问题:如何计算三维空间中一个点绕着某一条向量旋转一个特定角度之后的坐标?旋转矩阵、欧拉角和四元数都是用来解决这个问题的方法。关于四元数的内容可以看平面等距变换、三维空间的旋转和四元数

接下来我们来讨论一下旋转矩阵和欧拉角这两个方法,并且我们选取右手坐标系作为我们的坐标系。

旋转矩阵

首先,对于一个三维空间的点 P(x, y, z),要将其绕 z 轴旋转 \theta 角度是可以很简单地用旋转矩阵来表示的

 \begin{equation}
R_z(\theta) =
\begin{bmatrix}
\cos \theta & -\sin \theta & 0 \\
\sin \theta & \cos \theta & 0 \\
0 & 0 & 1 \\
\end{bmatrix}
\end{equation}

类似地,绕另外两个坐标轴旋转的矩阵可以表示如下(注意 R_y 有些不同)

 \begin{equation}
R_x(\theta) =
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos \theta & -\sin \theta \\
0 & \sin \theta & \cos \theta \\
\end{bmatrix}
\end{equation}

\begin{equation}
R_y(\theta) =
\begin{bmatrix}
\cos \theta & 0 & \sin \theta \\
0 & 1 & 0 \\
-\sin \theta & 0 & \cos \theta \\
\end{bmatrix}
\end{equation}

对于这三个特殊的旋转轴我们已经有了解决方案了,那么对于任意的轴 p 要怎么办呢?我们可以将这个旋转分解:

  • 将整个坐标轴旋转,使得旋转轴 pz 轴重合
  • 再将点 Pz 轴旋转 \theta
  • 再将整个坐标轴旋转回原位

3d-angle-rotation-matrix

如图,我们可以用两个角 \phi\psi 表示一个旋转轴的位置(这里认为旋转轴是个单位向量

 \begin{equation}
\left \{ \begin{aligned}
x &= \sin \phi \cos \psi \\
y &= \sin \phi \sin \psi \\
z &= \cos \phi
\end{aligned} \right.
\end{equation}

这样整个旋转就可以表示如下(最先进行的旋转它对应的旋转矩阵在最右侧)

R_z(\psi)R_y(\phi)R_z(\theta)R_y(-\phi)R_z(-\psi)

或者利用 R(-\alpha) = R^{-1}(\alpha) = R^T(\alpha),可以改写为

R_z(\psi)R_y(\phi)R_z(\theta)R_y^T(\phi)R_z^T(\psi)

经过化简,就可以得到最终的旋转矩阵

\begin{equation}
\begin{bmatrix}
\cos \theta + x^2(1 - \cos \theta) & xy(1 - \cos \theta) - z\sin \theta & xz(1 - \cos \theta) + y\sin \theta \\
yx(1 - \cos \theta) + z\sin \theta & \cos \theta + y^2(1 - \cos \theta) & yz(1 - \cos \theta) - x\sin \theta \\
zx(1 - \cos \theta) - y\sin \theta & zy(1 - \cos \theta) + x\sin \theta & \cos \theta + z^2(1 - \cos \theta)
\end{bmatrix}
\end{equation}

我们将旋转矩阵左乘需要旋转的向量就可以得到旋转后的结果了!

旋转矩阵一个很方便的地方是它可以沿着任意轴任意角度的旋转,但是,旋转矩阵缺点是它需要有 9 个元素来表示一个旋转,而且矩阵的乘法也比较慢。

欧拉角

euler-angle

欧拉角是用三个旋转角度 \alpha, \beta, \gamma 来标示旋转的。如图,图中蓝色坐标系是起始的坐标系,红色的坐标系是最后旋转完成的坐标系。整个旋转分为三个步骤:

  • 将坐标系绕 z 轴旋转 \alpha
  • 将旋转后坐标系绕自己本身的 x 轴(也就是图中的 N 轴)旋转 \beta
  • 将旋转后坐标系绕自己本身的 z 轴旋转 \gamma

由于绕不同的轴旋转所得到的欧拉角是不同的,所以欧拉角在使用的时候必须要先指明旋转的顺序,这里使用的是“zxz”的顺序。

欧拉角表示的旋转转换成旋转矩阵就是

 \begin{equation} R_z(\alpha)R_x(\beta)R_z(\gamma) \end{equation}

需要注意,这里的后面两次旋转并不是在原本固定的坐标系下的旋转。在旋转矩阵中,每一次旋转的叠加都是在左边乘上对应的旋转矩阵,然而,在这里乘法的顺序是相反的。

这可以这样来理解,假设最原始的固定的坐标系是 C_0

  • 假设有一个和 C_0 重合坐标系 C_3,先将 C_3C_0z 轴旋转 \gamma
  • 假设有一个和 C_0 重合坐标系 C_2,将 C_2 和前一步旋转后的 C_3 一起绕 C_0x 轴旋转 \beta
  • 假设有一个和 C_0 重合坐标系 C_1,将 C_1 和前一步旋转后的 C_2, C_3 一起绕 C_0z 轴旋转 \alpha

然后我们仅看这三个坐标系的关系:C_0 绕自己的 z 轴旋转 \alpha 角就可以和 C_1 重合;C_1 绕自己的 x 轴旋转 \beta 角可以和 C_2 重合,这是因为 C_1C_2 在最后一步是一起旋转的,它们的相对位置不会改变;同样可以知道 C_2 绕自己的 z 轴旋转 \gamma 角就可以和 C_3 重合。

在这里 C_1, C_2 相当于是前面欧拉角旋转的前两步得到的坐标系。

这样的话从 C_0C_3 的变换就相当于之前欧拉角的旋转变换,因此按照这个过程,旋转矩阵就是按照上面的顺序相乘了。