Chengyuan Ma's Blog

幽居默处而观万物之变,尽其自然之理,而断之于中

0%

投影

高中的时候就学过,向量\boldsymbol b\boldsymbol a上的投影为

\boldsymbol p = \frac{\boldsymbol a^{\mathrm T}\boldsymbol b}{\boldsymbol a^{\mathrm T}\boldsymbol a} \cdot \boldsymbol a
当时推导主要利用了点积的几何意义$a b = |a||b|$。在代数意义上上式是怎么推导的呢?

其实也不难:注意到\boldsymbol b和其投影\boldsymbol p的差一定是垂直于基向量\boldsymbol a的。因为\boldsymbol p\boldsymbol a平行,设\boldsymbol p = x\boldsymbol a

\begin{aligned}
    \boldsymbol a^{\mathrm T}\left(\boldsymbol b - x\boldsymbol a\right) &= 0 \\
    \Rightarrow \boldsymbol a^{\mathrm T}x\boldsymbol a &= \boldsymbol a^{\mathrm T}\boldsymbol b \\
    \Rightarrow x\boldsymbol a^{\mathrm T}\boldsymbol a &= \boldsymbol a^{\mathrm T}\boldsymbol b \\
    \Rightarrow x &= \frac{\boldsymbol a^{\mathrm T}\boldsymbol b}{\boldsymbol a^{\mathrm T}\boldsymbol a}
\end{aligned}
我们就得到了一开始的结论。如果把\boldsymbol a写到左边(因为是向量的实数积所以可以交换)
p = \boldsymbol a\cdot \frac{\boldsymbol a^{\mathrm T}\boldsymbol b}{\boldsymbol a^{\mathrm T}\boldsymbol a} = \frac{\boldsymbol a\boldsymbol a^{\mathrm T}}{\boldsymbol a^{\mathrm T}\boldsymbol a}\boldsymbol b
我们这下看清楚了:将\boldsymbol b转变为其在\boldsymbol a上的投影的是一个矩阵!这个矩阵P = \frac{\boldsymbol a\boldsymbol a^{\mathrm T}}{\boldsymbol a^{\mathrm T}\boldsymbol a}被称为投影矩阵

代数推导和投影矩阵有啥用?它们可以让我们把投影从投影到一个向量推广到投影到向量空间。

考虑如何把一个向量\boldsymbol b投影到矩阵A = \begin{bmatrix} \boldsymbol {a_1} &\boldsymbol {a_2} &\cdots &\boldsymbol {a_n}\end{bmatrix}的列空间上(假设A的列向量线性无关)。我们一样考虑\boldsymbol b\boldsymbol p的差,其一定垂直于C(A),自然也垂直于\boldsymbol {a_1} \cdots \boldsymbol {a_n},于是就有:

\boldsymbol {a_k} ^{\mathrm T}\left(\boldsymbol b - \boldsymbol p\right) = 0, \quad k = 1,2,\cdots,n
或者写在一起:
\begin{bmatrix}
    \boldsymbol {a_1}^{\mathrm T} \\
    \boldsymbol {a_2}^{\mathrm T} \\
    \vdots \\
    \boldsymbol {a_n}^{\mathrm T} \\
\end{bmatrix}
\left(\boldsymbol b  - \boldsymbol p\right) = \boldsymbol 0 \Rightarrow 
A^{\mathrm T} \left(\boldsymbol b  - \boldsymbol p\right) = \boldsymbol 0
因为\boldsymbol p \in C(A),所以\boldsymbol pA列向量的线性组合,可以写作\boldsymbol p = A\boldsymbol x。于是我们就得到了方程
A^{\mathrm T} \left(\boldsymbol b  - A\boldsymbol x\right) = \boldsymbol 0 \Rightarrow A^{\mathrm T}A\boldsymbol x = A^{\mathrm T}\boldsymbol b
因为A的列是线性无关的,所以A^{\mathrm T}A是可逆的(这个我们待会证明),于是我们可以解得
\boldsymbol x = \left(A^{\mathrm T}A\right)^{-1}A^{\mathrm T}\boldsymbol b
\boldsymbol p = A\left(A^{\mathrm T}A\right)^{-1}A^{\mathrm T}\boldsymbol b
此时我们很清楚地看到,投影矩阵是P = A\left(A^{\mathrm T}A\right)^{-1}A^{\mathrm T}


我们接下来证明之前我们用到的一个结论:A^{\mathrm T}A可逆的充分必要条件是A的列向量线性无关

我们发现这个结论等价于:A^{\mathrm T}AA有相同的零空间。

怎么证明?一个方向是简单的:

A\boldsymbol x = \boldsymbol 0 \Rightarrow A^{\mathrm T}\left(A\boldsymbol x\right) = \boldsymbol 0 \Rightarrow A^{\mathrm T} A\boldsymbol x = \boldsymbol 0
反方向的证明需要一些技巧:
A^{\mathrm T} A\boldsymbol x = \boldsymbol 0 \Rightarrow \boldsymbol x^{\mathrm T} A^{\mathrm T} A\boldsymbol x = 0 \Rightarrow \left(A\boldsymbol x\right)^{\mathrm T}A\boldsymbol x = 0 \Rightarrow A\boldsymbol x = \boldsymbol 0
倒数第二步写得有点花里胡哨其实就是向量自己和自己的点积,如果是零的话显然原向量就是零。

证毕。

最小二乘拟合

投影和最小二乘拟合是怎么联系起来的呢?

线性拟合的本质是用给定的若干列向量的线性组合表示另一个向量。

比如我有一组数据\boldsymbol x,\boldsymbol y,我的模型是y = a + bx,拟合的过程就是用向量\boldsymbol 1\boldsymbol x的线性组合表示\boldsymbol y的过程:

\begin{bmatrix}
    1 &x_1 \\
    1 &x_2  \\
    \vdots &\vdots\\
    1 &x_n 
\end{bmatrix}
\begin{bmatrix}
    a \\
    b
\end{bmatrix} \stackrel{?}{=} 
\begin{bmatrix}
    y_1 \\
    y_2 \\
    \vdots \\
    y_n 
\end{bmatrix}
“线性拟合”的线性指的是这种组合的线性,和模型本身的线性是没有关系的。模型里面有二次项不要紧,只要二次项和其他项是线性组合的,那也不过是在左边的矩阵里多了一个列向量而已:
\begin{bmatrix}
    1 &x_1 &x_1^2 \\
    1 &x_2 &x_2^2 \\
    \vdots &\vdots &\vdots \\
    1 &x_n &x_n^2
\end{bmatrix}
\begin{bmatrix}
    a \\
    b \\
    c
\end{bmatrix} \stackrel{?}{=} 
\begin{bmatrix}
    y_1 \\
    y_2 \\
    \vdots \\
    y_n 
\end{bmatrix}
那么我们发现线性拟合的本质其实就是解一个形如A\boldsymbol x = \boldsymbol b的方程组(这里的\boldsymbol x表示模型中各项的系数,和上面表示自变量的\boldsymbol x没有关系),这个我们老熟悉了。

但是事情不太妙的就是A又高又瘦,方程组比未知数多,所以几乎肯定不存在准确解。

什么时候存在准确解?\boldsymbol b要在A的列空间C(A)里。

所以如果只能求近似解,我们可以考虑把\boldsymbol b“近似到”C(A)中的一个最近的向量,然后就能求解了。

这个“近似”的过程其实就是投影。在C(A)\boldsymbol b的投影肯定最接近\boldsymbol b,因为只有二者的差是和C(A)垂直的,这很符合我们的代数直觉。

所以联系之前的投影的推导,假设这个投影向量是A\boldsymbol{\hat{x}},我们直接写出方程和解:

A^{\mathrm T} \left(\boldsymbol b  - A\boldsymbol{\hat{x}} \right) = \boldsymbol 0 \Rightarrow \boldsymbol{\hat{x}} = \left(A^{\mathrm T}A\right)^{-1}A^{\mathrm T}\boldsymbol b
这就叫“最小二乘拟合”。


如果非要较真,从微积分的角度理解一波也是可以的:定义拟合的误差函数为

L= \left|\boldsymbol b - A\boldsymbol{\hat x}\right| = \sum_{i = 1}^n \left(\boldsymbol b_i - \left(A\boldsymbol{\hat x}\right)_i\right)^2
如果要最小化L,那么L关于\boldsymbol{\hat x}的所有分量的导数都为0
\frac{\mathrm d L}{\mathrm d\boldsymbol{\hat x}_k} = 2\sum_{i = 1}^n A_{i, k} \left(\left(A\boldsymbol{\hat x}\right)_i - \boldsymbol b_i\right) =  0, \quad \forall k = 1, 2, \cdots, n
Ai列的列向量为\boldsymbol a_i,那么上面也可以写成
\frac{\mathrm d L}{\mathrm d\boldsymbol{\hat x}_k} = 2\boldsymbol a_k^{\mathrm T} \left(A\boldsymbol{\hat x} - \boldsymbol b\right) =  0, \quad \forall k = 1, 2, \cdots, n
并在一起得
2A^{\mathrm T}\left(A\boldsymbol{\hat x} - \boldsymbol b\right) =  0
这就是上面我们通过投影推出的方程。

L中的平方项也是最小二乘中“二”的由来。

矩阵的Gram-Schmidt正交化

我们发现无论是投影还是最小二乘拟合,计算A^{\mathrm T}A和它的逆都是绕不过去的一步。A^{\mathrm T}A本身就有一些特殊的性质,比如说其一定是正方阵,一定是对称的(通过转置等于自身易证)。但是这些性质在计算的时候帮助不大。如果我们能够优化A的形态从而简化A^{\mathrm T}A和它的逆的计算,那么各种算法的效率都会高上不少。

注意到A^{\mathrm T}A中的每一个位置都是A两个列向量的点积。如果A当中的列向量是两两正交的,那么A^{\mathrm T}A就变成了一个对角阵。对角阵的逆不要太好算!更进一步,如果A的每一个列向量都是单位向量,那么A^{\mathrm T}A就是单位阵,连逆都不用求了

所以说啊,正交是个好东西。那对于一个矩阵A,通过什么样的方式将其列向量完成正交化呢?

比较简单的一个算法称为Gram-Schmidt正交化。设A的列向量为\boldsymbol {a_1},\cdots, \boldsymbol{a_n},则向量正交化之后的结果\boldsymbol{q_i}可以这么计算:

\boldsymbol{q_i} = \boldsymbol{a_i} - \sum_{j = 1}^{i - 1} \frac{\boldsymbol{q_j}^{\mathrm T}\boldsymbol{a_i}}{\boldsymbol{q_j}^{\mathrm T}\boldsymbol{q_j}} \boldsymbol{q_j}
即对于第i个列向量通过减去其在所有之前的(完成正交化的)向量上的投影来确保其和前面的向量都正交。

最后,把所有的\boldsymbol{q_i}除以其模完成归一化,我们就可以得到一个标准正交列向量组成的矩阵Q,满足Q^{\mathrm T}Q = I。(线性代数中常使用Q表示列向量标准正交的矩阵)

A变到Q的过程没有改变列空间,所以不会改变投影矩阵,但计算投影矩阵的开销大大降低了:

P = Q\left(Q^{\mathrm T}Q\right)^{-1}Q^{\mathrm T} = QQ^{\mathrm T}

QR分解

Gram-Schmidt正交化的过程可以看作是在原矩阵上进行列变换,而列变换是可以用矩阵表示的。类比从高斯消元导出LU分解的思路,我们可以从Gram-Schimidt正交化导出一个矩阵的QR分解:

A = QR
我们有A,可以通过Gram-Schmidt正交化算出Q,但R怎么计算呢?因为Q^{\mathrm T}Q=I,所以只要在等式两边同时乘以Q^{\mathrm T},就得到R=Q^{\mathrm T}A。或者用QA的列向量表示:
R = \begin{bmatrix}
    \boldsymbol{q_1}^{\mathrm T} \\
    \boldsymbol{q_2}^{\mathrm T} \\
    \vdots \\
    \boldsymbol{q_n}^{\mathrm T}
\end{bmatrix}
\begin{bmatrix}
    \boldsymbol{a_1} &\boldsymbol{a_2} &\cdots &\boldsymbol{a_n}
\end{bmatrix}
= \begin{bmatrix}
    \boldsymbol{q_1}^{\mathrm T}\boldsymbol{a_1} &\boldsymbol{q_1}^{\mathrm T}\boldsymbol{a_2} &\cdots &\boldsymbol{q_1}^{\mathrm T}\boldsymbol{a_n} \\
    \boldsymbol{q_2}^{\mathrm T}\boldsymbol{a_1} &\boldsymbol{q_2}^{\mathrm T}\boldsymbol{a_2} &\cdots &\boldsymbol{q_2}^{\mathrm T}\boldsymbol{a_n} \\
    \vdots &\vdots &\ddots &\vdots \\
    \boldsymbol{q_n}^{\mathrm T}\boldsymbol{a_1} &\boldsymbol{q_n}^{\mathrm T}\boldsymbol{a_2} &\cdots &\boldsymbol{q_n}^{\mathrm T}\boldsymbol{a_n}
\end{bmatrix}
注意到,因为\boldsymbol{q_i}正交于所有\boldsymbol {q_j}(i < j),而\boldsymbol{a_j}又是由\boldsymbol{q_1},\cdots,\boldsymbol{q_j}表示的,所以i > j时,\boldsymbol{q_i}^{\mathrm T}\boldsymbol{a_j} = 0。所以R其实是一个上三角阵
R = \begin{bmatrix}
    \boldsymbol{q_1}^{\mathrm T}\boldsymbol{a_1} &\boldsymbol{q_1}^{\mathrm T}\boldsymbol{a_2} &\cdots &\boldsymbol{q_1}^{\mathrm T}\boldsymbol{a_n} \\
    0 &\boldsymbol{q_2}^{\mathrm T}\boldsymbol{a_2} &\cdots &\boldsymbol{q_2}^{\mathrm T}\boldsymbol{a_n} \\
    \vdots &\vdots &\ddots &\vdots \\
    0 &0 &\cdots &\boldsymbol{q_n}^{\mathrm T}\boldsymbol{a_n}
\end{bmatrix}
再观察R的各个非零位置的值以及Gram-Schimidt的公式,其实R是可以在计算Gram-Schmidt正交化的时候顺便算出来的。

QR分解有什么用呢?之前我们说正交化不改变投影矩阵,原因是“正交化不改变列空间”,QR分解能让我们在代数上证明这一点:

\begin{aligned}
    P &= A\left(A^{\mathrm T}A\right)^{-1}A^{\mathrm T} \\
    &= QR\left(R^{\mathrm T}Q^{\mathrm T}QR\right)^{-1}R^{\mathrm T}Q^{\mathrm T} \\
    &= QR\left(R^{\mathrm T}R\right)^{-1}R^{\mathrm T}Q^{\mathrm T} \\
    &= QRR^{-1}\left(R^{\mathrm T}\right)^{-1}R^{\mathrm T}Q^{\mathrm T} \\
    &= QQ^{\mathrm T}
\end{aligned}
注:因为R是一个方阵,所以我们这里能够把\left(R^{\mathrm T}R\right)^{-1}改写成R^{-1}(R^{\mathrm T})^{-1}\left(A^{\mathrm T}A\right)^{-1}是不能这么转换的,因为长方阵不存在常规意义上的逆。

类似地,QR分解还简化了最小二乘拟合的计算:

\boldsymbol{\hat{x}} = \left(A^{\mathrm T}A\right)^{-1}A^{\mathrm T}\boldsymbol b = R^{-1}Q^{\mathrm T}\boldsymbol b

又一次在网上看到中文编程的日经贴后,我觉得这似乎会是一个不错的博文题目。

这篇文章我计划讲的内容在标题里写得很明白——“迷思”。

迷思,经常被用作英文myth的音译,衍生出“神话”、“传言”的意项。在这个定义下,“中文编程”或许可以称得上是一个迷思了,初学编程的萌新和温饱思扯淡的键盘程序语言设计家(我或许可以忝列后者)一天到晚在讲,Github上偶尔有个小项目也会有不少的曝光率,但是也就止步于笑谈中了,现实的应用几乎没有(唯一勉强成气候的易语言常年处于鄙视链的最低端)。这样一种现状,似乎和动漫当中的“都市传说”颇有共同点。可是动漫的“都市传说”最后都是真的,“中文编程”则未必然。

若纯粹地拆字组词,“迷思”又有一种“胡思乱想”的意项在,这是我对于我想法的总结——胡思乱想。这篇文章就是胡思乱想,想着想着觉得有暴论出现可以写一下,仅此而已。我会尽力组织自己的逻辑,但是要这篇文章变成如同高考作文一般逻辑严密,既非我愿,也非我力所能及。

对中文编程我一直保有比较浓厚的兴趣,记得小学刚毕业的时候简短地用过所谓易语言。后来自己以前闲着无聊的时候写过自己脚本语言的解释器所以对于语言设计实现大约还算有一点点心得。恰好最近手痒又想写parser,我就想:为什么不试试看写一个中文的脚本语言呢?

然而思前想后,哪怕只是纸上谈兵的幻想,结果都不算乐观。

做个暴论吧:目前的所谓中文编程语言不过花瓶,真正实用的中文编程(如果真的有必要的话)任重而道远。


我一直认为,要大致把握一个编程语言的特征,写一段快速排序的代码是一个非常好的选择。快排是一个不简单且不平凡的算法,在传统命令式语言的理论框架内,顺序,分支和循环结构都用到了;如果更特殊的特性,快排的复杂度也勉强允许一些语言施展一下其特殊的语言构造。例如Haskell经典的四行快排我觉得就是对函数式语言哲学的一种不错的体现,C++或许也可以炫一波模板元编程和新出的concept。与此同时快速排序也不是特别复杂,让人能够多花心思在语言的特征而不是辨认算法本身上。最后的最后,它还能让我复习一下快排咋写,何乐而不为?

先放一段快排的C代码,代码是我从网上随手拉来的,显然还有很大的优化空间,但这并不重要:

void quick_sort(int *arr, int left, int right) {
    if (left >= right) return; 
    int p = arr[left];
    int i = left, j = right;
    while (i < j) {
        while (arr[j] >= p && i < j) j--;         
        while (arr[i] <= p && i < j) i++;
        int temp = arr[j];
        arr[j] = arr[i];
        arr[i] = temp;
    }
    arr[left] = arr[i]; 
    arr[i] = p;
    quick_sort(arr, left, j - 1);
    quick_sort(arr, j + 1, right);
}

那么问题来了:人们心目中与上面这一段C对应的“中文代码”应该是什么个样子呢?

是这样吗?

定义 函数 快速排序(输入:整数 数组,左边界:整数,右边界:整数):
    如果(左边界 >= 右边界):返回
    令 基准 为 输入【左边界】
    令 甲 为 左边界
    令 乙 为 左边界
    当(甲 < 乙)时循环:
        当(输入【乙】 >= 基准 且 甲 < 乙)时循环:
            乙 自减
        当(输入【甲】 <= 基准 且 甲 < 乙)时循环:
            甲 自增
        交换(输入【甲】,输入【乙】)
    令 输入【左边界】 为 输入【甲】
    令 输入【甲】 为 基准
    快速排序(输入,左边界,乙 - 1)
    快速排序(输入,乙 + 1,右边界)

中文编程是为了接地气,不是接地府

以上文法的所谓“中文编程语言”在现实开发中是极为不便的。

但很遗憾,这就是现在网络上通行的“中文编程”的现状。

有错吗?乍一看似乎是没有的,在某种意义上这还是“要素完全”的:

  1. 上面的代码一个拉丁字母都没有用到——“完全汉化”
  2. 尽量使用全角字符——“贴合国人使用习惯”
  3. 可读性勉强令人满意。

然而上面这坨东西我写的时候血压直线飙升。是哪里出了问题?

不妨再看看下一种可能的情况:

定义快排(整数组输入,整数左界,整数右界):
    若左界>=右界:返回
    基准赋输入于左界
    甲赋左界;乙赋右界
    当甲<乙:
        当输入于乙>=基准且甲<乙:
            乙自减
        当输入于甲<=基准且甲<乙:
            甲自增
        交换(输入于甲,输入于乙)
    输入于左界赋输入于甲
    输入于甲赋基准
    快排(输入,左界,乙-1)
    快排(输入,乙+1,右界)

看起来似乎比最上面一个例子要好一点是不是?

如果在这个基础上进行一定的代码高亮,那感觉就更上一个台阶:

虽然一开始的那一版“中文编程”也可以通过适当的代码高亮在效果上进行优化,但各位试一试就知道,效果应当还是不如上图的。那么,两种文法,都是“中文编程”,是什么导致了舒适度上的巨大差异?


回顾编程语言的历史,我们发现,编程语言正逐渐以语言本身的定位为基础,探寻如下三者的最优平衡:

  1. 对于机器来说好读:文法无二义性,且运行效率高
  2. 对于人类来说好读:文法大致贴合阅读习惯,或者说,在训练之后可以快速阅读
  3. 对于人类来说好写:文法“废话少”,开发效率高。

在文法解析技术、编译器优化技术、以及硬件技术如此发达的当下,第一点的重要性正在逐渐被弱化。那么剩下的,就只有对人好读好写的要求了。既然是“对人”,就自然需要考虑编程语言所依附的自然语言的语言特性。所有人都知道英语和中文无论是在直觉上还是在语言学意义上都有着巨大的差异,因此基于中文的编程语言,又怎么可以生搬硬套基于英语的编程语言的结构呢?以上示例中的两门语言之所以在观感上具有差距,就是因为前者是可以通过直接的单词替换预处理还原为另一门语言的“套皮”,而后者则做了一定程度上的“本地化”。

我们首先注意到,中文当中是没有空格分词的。在输入大段中文的时候,空格键往往只作为输入法选词的快捷键而本身不进入文本。这是很多中文编程语言所忽视的。例如

每 个 中文 词 之间 都有 空格, 写 起来 很 别扭, 读 起来 也 不 习惯。

每个中文词之间没有空格读写起来会更为自然。

这是很多中文编程语言所忽视的一点。然而也必须承认,空格的存在给词法分析提供了极大的便利。联系到当下对于大段中文的分词往往采用,也只能采用启发式的算法的现状,可以说在一门中文编程语言当中删去空格,虽然可以让读写更为顺畅,但对文法的设计与解析提出了更大的挑战。(同时,空格的去除让代码开起来更为紧凑,在汉字本身就结构相对复杂的情况下就更会对可读性产生一定的影响,因此在这种情况下,代码高亮是必要的

我们还发现,当下的编程语言在英语的基础上对于关键词进行大量的简化。例如,define在Python中被简化为def;function在Go中被简化为func,在Kotlin中被简化为fun,在Rust中更是被简化到fn;integer在很多语言当中都被简化到int,structure被简化到struct……对于英语等拼音文字,这种简化往往是通过去除辅音,删去末音节等完成的,较为简单。但是对于中文象形文字,这种简化就很难了,因为一个词往往就双音节两个字,去了哪个都不行,而仿照去除辅音搞去除笔画之类的毫无意义:结构变成吉勾?还是两个汉字不说,词本身的辨识性已经被破坏殆尽。

然而,我们关注的是目的而不是手段,抛开如此简单粗暴的类比,事情就豁然开朗了。编程语言简化关键词拼写的目的何在?一方面是让代码更加简洁,阅读效率更高。另一方面则是提升输入效率。中文或许在阅读效率方面已经无可挑剔,但是在输入效率上一直差强人意——我不是说中文本身的输入效率,而是写代码本身的效率。二者有何区别?

如今不用代码补全写代码的人已经不多了。代码补全的本质,是在用户敲击键盘进行原始输入的同时,结合上下文,给出符合编程语言语法以及(理想情况下)代码语义的代码建议。仔细想想,这和中文的输入法何其相似!然而现在并没有专为写代码而生的输入法。要写中文代码,就必须先敲键盘,再输入法选词,等到中文字实实在在地落到了编辑器里,代码补全才会出现,然后用户就要再次选词,如图

显然,在上例中,如果能把输入法整合进自动补全(或者反过来),使得在键入ks时就自动跳出快速排序的提示项,对于开发效率和流畅性会有可观的提升(在没有空格分词的中文编程语言中,这样的融合自动补全可以综合考虑到关键词以及文法结构,这种流畅性提升会更为显著)。否则,对于针对自然语言输入而设计的输入法而言,编程的上下文不免会让提词非常别扭;对于代码补全而言,变量名一般就三四个字个字,在输入法打完两个字后再补全一般似乎也如同鸡肋。两者没有沟通各自为战,折磨的就是开发者。理想情况下,中文编程中的输入法和代码补全不说是合二为一也应当是互相合作的。相比于独立地弹出一个下拉框,代码补全应当使用一个输入法提供的接口访问用户在键盘上的原始输入并给出建议,这样用户在输入法的选词界面就可以直接看到代码提示,效率何止提升了两倍。然而,现实中的各输入法都热衷于闭门造车,即使是以开源闻名的Rime都不提供这种动态提词接口。因此,如果要实现前文的创想,倒是IDE的开发者自己从零写一个兼容自动补全的输入法更现实一点。而后者可就造轮子造大了去了。

简单的总结一下,如果要实现我心目中理想的中文编程语言,应当有这两个特点:

  1. 本土化的语言文法:在不引入二义性的前提下,降低空格分词的必要性,或更多地从中文本身借鉴一部分文法结构,使得汉字的输入和阅读更符合自然语言的习惯。不满足这一点,则不免对于阅读的流畅性产生影响,产生生硬之观感,被好事者讥为“套皮”。
  2. 高效率的开发环境:在开发环境上完成输入法和自动补全的整合,提升写代码的流畅性,同时提供靠谱的代码高亮方案以辅助阅读。不满足这一点,则读起来再赏心悦目的编程语言写起来都有不便,难免花瓶一个,华而不实。

这两点的难度不是一般的大,因此我觉得实用的中文编程任重而道远。当然,考虑到中文编程除了看着一乐以外未必有其他特殊的吸引力,实际上这种编程语言大概是遥遥无期了。

迷思迷思,思到这大概也就结束了。整篇文章想到哪里写到哪里,估计很乱吧(笑),如果有人读着这个觉得摸不着头脑,实在是抱歉。

这篇文章是在经历\LaTeX改格式改得死去活来后对于以后可能经常会用到的格式调整小技巧的一个汇总。许多片段的来源是TeX StackExchange, 源链接就不一一表明了。这些技巧在XeTeX上编译通过,对于其他的引擎尚未测试。

该列表预计会随着我被\LaTeX虐的次数增加而不断更新。

页边距设置

使用geometry宏包:

\usepackage{geometry}
\geometry{a4paper, top = 34.4mm, bottom = 45mm, left = 23mm, right = 23mm}

更高级的用途可以参阅其官方文档。

字体设置

Times New Roman

为了全局设置Times New Roman字体,我们需要fontspec宏包:

\usepackage{fontspec}
\setmainfont{Times New Roman}

如果希望数学公式的字体也使用Times New Roman,可以使用newtx宏包:

\usepackage{newtxmath}

注:newtx宏包下还有newtxtext可以代替fontspecsetmainfont的功能,但是不知为何在Windows的XeTeX上编译错误。归根结底,XeTeX以及LuaTeX的一个设计理念就是使用户可以直接调用系统字体而不用拘泥于宏包,因此我个人还是偏向于setmainfont的解决方案。

宋体粗体

CTeX在Windows系统下宋体采用的是中易的字库,而不幸的是中易的宋体并不包括粗体。Word之所以可以给中易宋体加粗,是因为它有在标准字重基础之上自动加粗的算法。而高贵的\LaTeX显然是不屑于用这种歪门邪道加粗字体的,所以你惊喜地发现在Windows下你似乎无论怎么\textbf/\bfseries都得不到加粗的宋体(于是CTeX就机智地把黑体作为粗体时的字体)。如果因为种种原因迫不得已一定要使用粗体,CTeX给出了两个解决方案:

  1. 伪粗体:通过对于一个汉字堆叠若干个微距平移标准字号的宋体实现视觉上的粗体效果。这种伪粗体效果肯定不如专门的宋体效果好,而且在PDF内选词的时候偶尔会遇到问题,因此我是没有使用过。
  2. 使用带有粗体的宋体:可以使用\setCJKmainfont命令转而使用例如思源黑体等设计时考虑粗体的字体。这在Windows下其实也不简单。最傻瓜的做法是在\documentclass{ctexart}前加上[fontset = fandol]选项,使得CTeX采用Fandol的字库。这家的字库可以在GPL协议下自由使用,但是缺点是缺的字不少。
  3. 不是解决办法的办法:装一台Ubuntu的虚拟机,Ubuntu下的默认宋体是带有粗字号的,问题直接解决。

标题格式的修改

一般来说,LaTeX内修改标题格式最合适的宏包是titlesectitlesec的文档相当详细,SE上也有非常多的教程,所以在这类就不多赘述了。我想记一笔的是CTeX内部自带的格式设置功能:这个功能在CTeX的官方手册上确有记载,但是在网上却不是很好找。个人认为用法比titlesec简洁明了,对于一般的中文文档已是相当够用了。这里举一个例子:

\ctexset{
    section/format = \songti\sffamily\zihao{-4}\bfseries, % 中文宋体,英文无衬线体,字号小四,加粗
    section/afterskip = 0pt, % 标题下方不留空,直接接段落,更紧凑
    subsection/format = \songti\sffamily\zihao{5}\bfseries, % 中文宋体,英文无衬线体,字号五号,加粗
    subsection/beforeskip = 0pt, % 标题上方不留空,直接在上一段下方,更紧凑
    subsection/afterskip = 0pt, % 标题上方不留空
    subsubsection/format = \songti\sffamily\zihao{5}\bfseries, % 中文宋体,英文无衬线体,字号五号,加粗
    subsubsection/beforeskip = 0pt,
    subsubsection/afterskip = 0pt   
}

“参考文献”标题格式的修改

和上节相同,只需要在\thebibliography之前修改section的格式就行了(如果还有附录的话不要忘记改回来)。

如果要修改“参考文献”这四个字本身,CTeX提供了refname选项。

自动引号配对

使用csquotes宏包:

\usepackage[autostyle=false, style=english]{csquotes}
\MakeOuterQuote{"}

中文的定理环境

CTeX会自动将amsthm包当中的证明环境从“Proof”改成“证明”,但在ctexart类型下不会对定理,引理,定义等环境进行汉化,因此需要手动加入如下代码:

\usepackage{amsthm}
\newtheorem{theorem}{定理}[section]
\newtheorem{definition}{定义}[section]
\newtheorem{corollary}{推论}[section]
\newtheorem*{remark}{注}
\newtheorem{lemma}{引理}[section]

列表环境下的间距

一般列表去除间距

使用enumitem宏包的nosep选项让enumerateitemize中相邻项之间的垂直间距为0:

\usepackage{enumitem}

\begin{enumerate}[nosep]

\end{enumerate}

文献列表去除间距

通过以下代码去除文献列表中相邻项的垂直间距:

% 去除参考文献一节文献项之间过大的间隙
\let\OLDthebibliography\thebibliography
\renewcommand\thebibliography[1]{
    \OLDthebibliography{ #1 }
    \setlength{\parskip}{0pt}
    \setlength{\itemsep}{0pt plus 0.3ex}
}

图注格式修改

最为通用,最为灵活的方法:

% 宋体五号粗体,#1#2#3分别是“图xxx”,“:”,以及图注文字本身
\DeclareCaptionFormat{mycaptionformat}{\songti\zihao{-5}\bfseries#1#2#3\par}
\captionsetup{format = mycaptionformat}

页眉页脚相关

首页特殊

\usepackage{fancyhdr}
\fancypagestyle{firststyle}{
    \fancyhf{}
    \renewcommand{\headrule}{\hrule height 0.5pt \vspace{0.8pt}\hrule height 0.5pt}
    \renewcommand{\footrule}{\hrule height 0.5pt width 0.3\textwidth}
    ...
}

% 在\begin{document}之后
\thispagestyle{firststyle}

装饰线

% 页眉双线
\renewcommand{\headrule}{\hrule height 0.5pt \vspace{0.8pt}\hrule height 0.5pt}
% 页脚部分装饰线
\renewcommand{\footrule}{\hrule height 0.5pt width 0.3\textwidth}

参数可以视情况自行调整;可以写在\fancypagestyle内部。

在一行内同时包括左中右对齐的文本

```latex \newcommand{}[3]{
\par

\makebox[][s]{

3月15日清晨6:30分

某人迷迷糊糊地刚刚从被窝里爬出来

迷迷糊糊地唤醒一夜没关的电脑

迷迷糊糊地点进MIT的portal

啊,果然status更新了呢,虽然不该抱什么期待的

自己这个水平申请纯属做分母哦,大概一定是thank you for applying to MIT这种套路话罢(确信)

为了彻底让自己死心点了进去并开始阅读

On behalf of the Admissions Committee, it is my pleasure to offer you admission to the MIT Class of 2025! You stood out ...

嗯?(要素察觉)

嗯?!?!?

千言万语在早晨上学的急忙中汇成一句话

草,居然录了?

然后简单地发了个朋友圈,吃完饭,开开心心上学去了。

学校申请的竞争是残酷的,所谓成王败寇莫不如是。抢救性写感言,各路学弟开始问经验,然后我其实一时半会并答不出来。另一方面,我在学姐的帮助下倒是在一个小时内就找到了组织见到了前几届MIT的学长,然后发现群里面个个都是集训队级别的大佬,区区国铜瑟瑟发抖。

一切都引向了一个我必须思考的问题:

我这么菜,为啥会被MIT看上呢?

论标化,自己的标化成绩在国内传统认知内对于MIT的申请只能称得上是堪堪够用。

论文书,因为自己原本对MIT的申请不抱特别大的希望,文书写的是中规中矩,没有特别多的灵光一现和奇妙构思。

论活动,自己没有上过许多国内奉为圭臬的夏校。没有辩论,没有模联,没有特别多社会活动。

论竞赛,自己国铜的最高水准大概只到了最低线?

论申请策略,自己直接说喜欢EECS的直球做法在许多人看来也是头铁。

那是为什么呢?

转了一圈,思绪回到1月6号申请最后一天在电脑前的奋笔疾书。

是了,原来如此,答案或许见诸于MIT申请系统的首页:

The key is to be yourself.

我在最后一天才意识到MIT的申请居然允许递交一份maker portfolio,那时的兴奋是无以言表的。

自己就把自己平时乱搞的小项目连简介带Github链接填了上去,填了九个,简介写了一千多个词,自己从未如此畅快。

想必让招生官高看我一眼的,大概就是这么多乱搞项目当中体现出来的对摸鱼计算机的一种热情?(暴论)

记得自己曾经做第一个小项目的时候,一方面是觉得自己开心,另一方面为未尝没有“如果项目牛逼了会不会被人赏识”的小功利心和小幻想。做了那么多小项目,快乐是总有的,但牛逼的项目是没有的,所以在高中的几年觉得大学申请愈发迫近也曾焦虑迷茫过,尝试过竞赛,但也难以完全投入精力,最终屈服于自己强大的摸鱼欲之下重操旧业。也不止一次地心中有“自己摸鱼做项目的时间是不是浪费了呢?”的疑问。

但是现实证明,所有的努力都不会被辜负。

自己的气质或许也和MIT很般配呢(大雾)。

当然,这也可以理解为是某人摸鱼主义哲学的大胜利(大雾)。

上图:当你摸鱼的contribution都状似摸鱼,你就达成了二阶摸鱼的至高境界,物极必反,这是MIT的隐藏通道

人啊,还是要有理想的,有理想地摸鱼,这就达成了一组矛盾的对立统一,矛盾摸鱼,乃摸鱼之大境界。

最后还是要说一句,Cornell你拒的好啊!

录了技校,人很开心,加上之前录了剑桥,所以现在就是一个去Cambridge还是去Cambridge的问题了,我选择去Cambridge。

今天数学课讲了一点组合数的性质,因为太简单了所以稍微划了一会水。

结果发现自己整出来了一个算\pi的式子?

化简以后发现是传说中的Wallis公式?

虽然原理一样的证明Wiki上有提到但是还是当场惊了。

简单记一笔吧。

一切的开端是今年寒假做夏校申请的时候证明过的一个极限

\lim_{n\to\infty} 2^{-2n}\binom{2n}{n} = 0
当时是用斯特林近似暴力代换进行证明(现在看来并不严谨)
\begin{aligned}
    \lim_{n\to\infty} \binom{2n}{n} \frac{\left[\sqrt{2\pi n} \left(n \over e\right)^n \right]^2} {\sqrt{4\pi n} \left(2n \over e\right)^{2n}} &= \lim_{n\to\infty} \left[\frac{\sqrt{2\pi n} \left(n \over e\right)^n}{n!}\right]^2 \frac{(2n)!}{\sqrt{4\pi n} \left(2n \over e\right)^{2n}} \\
    &= 1 \\     
    \Rightarrow 
    \lim_{n\to\infty} 2^{-2n}\binom{2n}{n} &= \lim_{n\to\infty} 2^{-2n}\frac {\sqrt{4\pi n} \left(2n \over e\right)^{2n}}{\left[\sqrt{2\pi n} \left(n \over e\right)^n\right]^2} \cdot \lim_{n\to\infty} \binom{2n}{n} \frac{\left[\sqrt{2\pi n} \left(n \over e\right)^n\right]^2} {\sqrt{4\pi n} \left(2n \over e\right)^{2n}} \\
    &=  \lim_{n\to\infty} 2^{-2n}\frac {\sqrt{4\pi} 2^{2n}n^{2n+\frac{1}{2}}e^{-2n}}{2\pi n^{2n+1}e^{-2n}} \\
    &= \lim_{n\to\infty} \frac{1}{\sqrt{\pi n}}\\
    &= 0
\end{aligned}
今天划水的时候发现由以上过程,这个极限可以加强为
\lim_{n\to\infty} \frac{\sqrt{\pi n}\binom{2n}{n}}{2^{2n}} = 1
稍作整理即得
\lim_{n\to\infty} \frac{2^{4n}}{n\binom{2n}{n}^2} = \pi
拿卡西欧摁了一下,发现式子没有假,虽然收敛得有亿点点慢但确实是收敛到了\pi

关键是这个式子我似乎没见过啊?内心直接膨胀,可把我牛逼坏了。

因为形式看起来较复杂,接着我就想可不可以通过邻项作比的方式变换一下形态

\begin{aligned}
\frac{2^{4n}}{n\binom{2n}{n}^2} \Bigg / \frac{2^{4(n-1)}}{(n-1)\binom{2n-2}{n-1}^2} &= \frac{16(n-1)}{n} \left[\binom{2n-2}{n-1} \bigg/ \binom{2n}{n}\right]^2 \\
&= \frac{16(n-1)}{n} \left[\frac{n^2}{2n(2n-1)}\right]^2 \\
&= \frac{16n(n-1)}{(4n-2)^2} \\
&= \frac{2n}{2n-1}\cdot \frac{2n-2}{2n-1}
\end{aligned}
这个形式似曾相识,结合上式把式子化为连乘积的形式,要素察觉!
\begin{aligned}
\frac{2^{4n}}{n\binom{2n}{n}^2} &= \frac{2^{4}}{1\binom{2}{1}^2} \cdot \frac{2}{3} \cdot \frac{4}{3} \cdot \frac{4}{5} \cdot \frac{6}{5} \cdots \frac{2n-2}{2n-1}\cdot \frac{2n}{2n-1} \\
&= 4\cdot \frac{2}{3} \cdot \frac{4}{3} \cdot \frac{4}{5} \cdot \frac{6}{5} \cdots \frac{2n-2}{2n-1}\cdot \frac{2n}{2n-1}
\end{aligned}
这不就是Wallis乘积公式嘛?翻出Wikipedia一看:
\begin{aligned}
\frac{\pi}{2} &=\prod_{n=1}^{\infty} \frac{4 n^{2}}{4 n^{2}-1}=\prod_{n=1}^{\infty}\left(\frac{2 n}{2 n-1} \cdot \frac{2 n}{2 n+1}\right) \\
&=\left(\frac{2}{1} \cdot \frac{2}{3}\right) \cdot\left(\frac{4}{3} \cdot \frac{4}{5}\right) \cdot\left(\frac{6}{5} \cdot \frac{6}{7}\right) \cdot\left(\frac{8}{7} \cdot \frac{8}{9}\right) \cdots
\end{aligned}
完全一致,直接得证。

哇,也就是说我划着水就把Wallis公式不严谨地整了一遍?

果然还是要膨胀.jpg

以前觉得这个式子很高端的,现在有种莫名的幻灭感。

但是写到这的时候多看了一眼,发现斯特林逼近的一个推导里用到了Wallis公式?

突然有点不确定这算不算是循环论证了。

It turns out that writing a STL container from scratch is mostly a tedious physical labor.

STL is created to save C++ programmers the time of reinventing wheels. Unfortunately, many STL data structures, most notably self-balancing BSTs (aka. std::(multi)set), are not extendable and are by itself too limited to be used in the context of competitive programming, forcing us to write our own BSTs again and again in competitions.

There comes the fact I find really interesting: A quick and dirty self-balancing BST implementation written during a competition is only about 50-60 lines long, while the STL implementation of std::set and std::multiset is usually more than a thousand-lines long in total.

It does make me wonder: what makes this huge difference? And, will the code of our BST bloat as well if we write it the STL way -- with generics, iterators, and all the necessary bits and pieces as specified in the reference?

To answer my question, I have here tried creating my implementation of multiset using treap as the underlying data structure. The result, treap_multiset, is almost a drop-in replacement to std::multiset. The few places where it does not conform to the C++ standard are:

  1. It is currently not allocator-aware, so all allocator-related features are not implemented.
  2. emplace and emplace_hint are not implemented.
  3. All operations that have logarithmic time complexity in std::multiset still have logarithmic time complexity here, but only in the average sense (because treap is a randomized data structure), and could have linear worst-case time complexity (though very, very, very unlikely).
  4. void erase(iterator) takes amortized logarithmic time instead of constant time.
  5. A few uncommonly-used member types are missing.

treap_multiset also supports two new operations that are not supported in the original std::multiset:

  1.    size_type rank(iterator it) const;
       size_type rank(const_iterator it) const;

    Both take average logarithmic time and return the rank / position of the iterator.

  2.    iterator at(size_type index);
       const iterator at(size_type index) const;

    Both take average logarithmic time and return the iterator at the specified index/position.

The code is shown below:

#include <stdexcept>
#include <random>
#include <algorithm>

static std::random_device random_device;
static std::mt19937_64 random_engine(random_device());

template <typename T>
struct treap_node {
    using rand_weight_type = decltype(random_engine)::result_type;
    using size_type = std::size_t;

    treap_node *left, *right, *parent;
    rand_weight_type weight;
    size_type size;
    T value;

    treap_node(const T &value): left(nullptr), right(nullptr), parent(nullptr),
        weight(random_engine()), size(1), value(value) {}

    treap_node(treap_node *left, treap_node *right, treap_node *parent, 
        rand_weight_type weight, size_type size, const T &value): 
        left(left), right(right), parent(parent), weight(weight),
        size(size), value(value) {}

    void update_size() {
        size = 1 + (left ? left->size : 0) + (right ? right->size : 0);
    }
};

#define IMPL_ITERATOR_MOVE_NEXT do { \
    if (!node) break; \
    if (node->right) { \
        node = node->right; \
        while (node && node->left) node = node->left; \
    } else { \
        bool from_right = true; \
        while (from_right) { \
            if (!node->parent) { node = nullptr; break; } \
            from_right = node->parent->right == node; \
            node = node->parent; \
        } \
    } \
} while (0)

#define IMPL_ITERATOR_MOVE_PREV do { \
    if (!node) break; \
    if (node->left) { \
        node = node->left; \
        while (node && node->right) node = node->right; \
    } else { \
        bool from_left = true; \
        while (from_left) { \
            if (!node->parent) { node = nullptr; break; } \
            from_left = node->parent->left == node; \
            node = node->parent; \
        } \
    } \
} while (0)

#define TREAP_ITERATOR_DECL(name, qualifier, inc, dec) \
template <typename T> struct name { \
    qualifier treap_node<T> *node; bool past_the_end; \
    name(): node(nullptr), past_the_end(true) {} \
    name(qualifier treap_node<T> *node, bool past_the_end) \
        : node(node), past_the_end(past_the_end) {} \
    name(qualifier treap_node<T> *node): node(node), past_the_end(false) {} \
    qualifier T &operator *() qualifier { \
        if (!node || past_the_end) \
            throw std::runtime_error("dereferencing null/past-end iterator"); \
        return node->value; \
    } \
    bool operator ==(const name<T> &b) const { \
        return node == b.node && past_the_end == b.past_the_end; \
    } \
    bool operator !=(const name<T> &b) const { \
        return node != b.node || past_the_end != b.past_the_end; \
    } \
    name<T> &operator ++() { \
        qualifier treap_node<T> *backup = node; \
        IMPL_ITERATOR_MOVE_##inc; \
        if (!node) past_the_end = true, node = backup; \
        return *this; \
    } \
    name<T> operator ++(int) { \
        name<T> ret(*this); \
        return ++(*this), ret; \
    } \
    name<T> &operator --() { \
        if (past_the_end) past_the_end = false; \
        else IMPL_ITERATOR_MOVE_##dec; \
        if (!node) \
            throw std::runtime_error("can't decrement at the beginning"); \
        return *this; \
    } \
    name<T> operator --(int) { \
        name<T> ret(*this); \
        return --(*this), ret; \
    } \
}

TREAP_ITERATOR_DECL(treap_iterator, /* NO QUALIFIER */, NEXT, PREV);
TREAP_ITERATOR_DECL(reverse_treap_iterator, /* NO QUALIFIER */, PREV, NEXT);
TREAP_ITERATOR_DECL(const_treap_iterator, const, NEXT, PREV);
TREAP_ITERATOR_DECL(const_reverse_treap_iterator, const, PREV, NEXT);

template <typename T, typename Compare = std::less<T>>
class treap_multiset {
public:
    using key_type = T;
    using value_type = T;
    using size_type = typename treap_node<T>::size_type;
    using key_compare = Compare;
    using value_compare = Compare;
    using node_type = treap_node<T>*;
    using iterator = treap_iterator<T>;
    using reverse_iterator = reverse_treap_iterator<T>;
    using const_iterator = const_treap_iterator<T>;
    using const_reverse_iterator = const_reverse_treap_iterator<T>;

    treap_multiset(): root(nullptr) {}

    treap_multiset(const treap_multiset &b): root(deep_copy(b.root)), comp(b.comp) {}

    treap_multiset(treap_multiset &&b): root(b.root), comp(b.comp) {}

    ~treap_multiset() { if (root) recursive_free(root); }

    bool empty() const { return root == nullptr; }

    size_type size() const { return root ? root->size : 0; }

    size_type max_size() const { return 0x7FFFFFFF; }

    key_compare key_comp() const { return comp; }

    value_compare value_comp() const { return comp; }

    iterator begin() 
        { return iterator(leftmost(root)); }

    const_iterator begin() const 
        { return const_iterator(leftmost(root)); }

    const_iterator cbegin() const 
        { return const_iterator(leftmost(root)); }

    iterator end() 
        { return iterator(rightmost(root), true); }

    const_iterator end() const 
        { return const_iterator(rightmost(root), true); }

    const_iterator cend() const 
        { return const_iterator(rightmost(root), true); }

    reverse_iterator rbegin() 
        { return reverse_iterator(rightmost(root)); }

    const_reverse_iterator rbegin() const 
        { return const_reverse_iterator(rightmost(root)); }

    const_reverse_iterator crbegin() const 
        { return const_reverse_iterator(rightmost(root)); }

    reverse_iterator rend() 
        { return reverse_iterator(leftmost(root), true); }

    const_reverse_iterator rend() const 
        { return const_reverse_iterator(leftmost(root), true); }

    const_reverse_iterator crend() const 
        { return const_reverse_iterator(leftmost(root), true); }

    iterator insert(const value_type &value) {
        node_type left, right;
        split_le(root, value, left, nullptr, right, nullptr);
        node_type temp = new treap_node<T>(value);
        root = join(join(left, temp), right);
        return iterator(temp);
    }

    iterator insert(iterator position, const value_type &value) 
        { return insert(value); }

    template <typename II>
    void insert(II first, II last) { 
        for (; first != last; first++) 
            insert(*first); 
    }

    size_type rank(iterator it) const { return rank(it.node); }
    
    size_type rank(const_iterator it) const { return rank(it.node); }

    iterator at(size_type index) {
        if (index < 0 || index > size()) return end();
        return iterator(at_internal(index));
    }

    const_iterator at(size_type index) const {
        if (index < 0 || index > size()) return end();
        return const_iterator(at_internal(index));
    }

    void erase(iterator pos) {
        assert_valid(pos);
        node_type a, b, c;
        size_type rank = this->rank(pos.node);
        split_size(root, rank, a, nullptr, c, nullptr);
        split_size(a, rank - 1, a, nullptr, b, nullptr);
        root = join(a, c);
        // assert(b == pos.node);
        delete b;
    }

    size_type erase(const value_type &key) {
        node_type a, b, c;
        split_le(root, key, a, nullptr, c, nullptr);
        split_re(a, key, a, nullptr, b, nullptr);
        root = join(a, c);
        if (b) {
            size_type ret = b->size;
            recursive_free(b);
            return ret;
        }
        return 0;
    }

    void erase(iterator first, iterator last) {
        size_type rank_first = rank(first);
        size_type rank_last = rank(last);
        node_type a, b, c;
        split_size(root, rank_last, a, nullptr, c, nullptr);
        split_size(a, rank_last - 1, a, nullptr, b, nullptr);
        root = join(a, c);
        if (b) recursive_free(b);
    }

    void clear() { if (root) recursive_free(root); }

    void swap(treap_multiset &b) { swap(root, b.root); }

    iterator find(const value_type &key) {
        node_type ret = find_internal(key);
        return ret ? iterator(ret) : end();
    }

    const_iterator find(const value_type &key) const {
        node_type ret = const_cast<treap_multiset<T>*>(this)->find_internal(key);
        return ret ? const_iterator(ret) : end();
    }

    size_type count(const value_type &key) const {
        node_type a, b, c;
        treap_multiset<T> *thiz = const_cast<treap_multiset<T>*>(this);
        thiz->split_le(root, key, a, nullptr, c, nullptr);
        thiz->split_re(root, key, a, nullptr, b, nullptr);
        size_type ret = b ? b->size : 0;
        thiz->root = thiz->join(thiz->join(a, b), c);
        return ret;
    }

    iterator lower_bound(const value_type &key) {
        node_type ret = lower_bound_internal(key);
        return ret ? iterator(ret) : end();
    }

    const_iterator lower_bound(const value_type &key) const {
        node_type ret = const_cast<treap_multiset<T>*>(this)->lower_bound_internal(key);
        return ret ? const_iterator(ret) : end();
    }

    iterator upper_bound(const value_type &key) {
        node_type ret = upper_bound_internal(key);
        return ret ? iterator(ret) : end();
    }

    const_iterator upper_bound(const value_type &key) const {
        node_type ret = const_cast<treap_multiset<T>*>(this)->upper_bound_internal(key);
        return ret ? const_iterator(ret) : end();
    }

    std::pair<iterator, iterator> equal_range(const value_type &key) {
        return std::make_pair(lower_bound(key), upper_bound(key));
    }

    std::pair<const_iterator, const_iterator> equal_range(const value_type &key) const {
        return std::make_pair(lower_bound(key), upper_bound(key));
    }

    size_type depth(node_type node) {
        if (!node) return 0;
        return 1 + std::max(depth(node->left), depth(node->right));
    }
    
private:
    Compare comp;
    node_type root;

    void recursive_free(node_type root) {
        if (root->left) 
            recursive_free(root->left);       
        if (root->right) 
            recursive_free(root->right);
        delete root;
    }

    void split_le(node_type root, const value_type &key, 
        node_type &left, node_type left_parent,
        node_type &right, node_type right_parent) {
        if (!root) { left = right = nullptr; return; }
        if (!comp(key, root->value)) {
            left = root; root->parent = left_parent;
            split_le(root->right, key, root->right, root, right, right_parent);
        } else {
            right = root; root->parent = right_parent;
            split_le(root->left, key, left, left_parent, root->left, root);
        }
        root->update_size();
    }

    void split_re(node_type root, const value_type &key, 
        node_type &left, node_type left_parent,
        node_type &right, node_type right_parent) {
        if (!root) { left = right = nullptr; return; }
        if (comp(root->value, key)) {
            left = root; root->parent = left_parent;
            split_re(root->right, key, root->right, root, right, right_parent);
        } else {
            right = root; root->parent = right_parent;
            split_re(root->left, key, left, left_parent, root->left, root);
        }
        root->update_size();
    }

    void split_size(node_type root, size_type size, 
        node_type &left, node_type left_parent,
        node_type &right, node_type right_parent) {
        if (!root) { left = right = nullptr; return; }
        size_type left_size = 1 + (root->left ? root->left->size : 0);
        if (left_size <= size) {
            left = root; root->parent = left_parent;
            split_size(root->right, 
                size - left_size, root->right, root, right, right_parent);
        } else {
            right = root; root->parent = right_parent;
            split_size(root->left, size, left, left_parent, root->left, root);
        }
        root->update_size();
    }

    node_type join(node_type left, node_type right) {
        if (!left) return right;
        if (!right) return left;
        if (left->weight <= right->weight) {
            node_type temp = join(left->right, right);
            if (temp) temp->parent = left;
            left->right = temp;
            left->update_size();
            return left;
        } else {
            node_type temp = join(left, right->left);
            if (temp) temp->parent = right;
            right->left = temp;
            right->update_size();
            return right;
        }
    }

    node_type leftmost(node_type x) const {
        if (!x) return nullptr;
        node_type ret = x; while (ret->left) ret = ret->left;
        return ret;
    }

    node_type rightmost(node_type x) const {
        if (!x) return nullptr;
        node_type ret = x; while (ret->right) ret = ret->right;
        return ret;
    }

    node_type find_internal(const value_type &key) {
        node_type a, b, c;
        split_le(root, key, a, nullptr, c, nullptr);
        split_re(root, key, a, nullptr, b, nullptr);
        root = join(join(a, b), c);
        return b;
    }

    node_type lower_bound_internal(const value_type &key) {
        node_type left, right;
        split_re(root, key, left, nullptr, right, nullptr);
        node_type ret = leftmost(right);
        root = join(left, right);
        return ret;
    }

    node_type upper_bound_internal(const value_type &key) {
        node_type left, right;
        split_le(root, key, left, nullptr, right, nullptr);
        node_type ret = leftmost(right);
        root = join(left, right);
        return ret;
    }

    node_type at_internal(size_type index) const {
        node_type temp = root;
        while (true) {
            size_type left_size = 1 + (temp->left ? temp->left->size : 0);
            if (index == left_size) return temp;
            else if (index < left_size) temp = temp->left;
            else temp = temp->right, index -= left_size;
        }
        return nullptr; // UNREACHABLE
    }

    size_type rank(node_type node) {
        bool from_right = true;
        size_type ret = 0;
        while (node) {
            if (from_right)
                ret += 1 + (node->left ? node->left->size : 0);
            if (node->parent)
                from_right = node == node->parent->right;
            node = node->parent;
        }
        return ret;
    }

    node_type deep_copy(node_type node) {
        if (!node) return nullptr;
        node_type left = deep_copy(node->left);
        node_type right = deep_copy(node->right);
        node_type ret = new treap_node<T>(
            left, right, nullptr,
            node->weight, node->size, node->value
        );
        if (left) left->parent = ret;
        if (right) right->parent = ret;
        return ret;
    }

    void assert_valid(iterator it) {
        if (!it.node || it.past_the_end)
            throw std::runtime_error("invalid iterator");
        node_type temp = it.node;
        while (temp->parent) temp = temp->parent;
        if (temp != root)
            throw std::runtime_error("invalid iterator");
    }
};

(The code above hasn't been thoroughly tested yet and could still contain bugs).

引入

给定一串n个珠子组成的环链,每个珠子可以被染成黑白两色,求本质不同的染色方案的个数?

群,置换和置换群

群的不严谨定义

如果我们对于一个集合G定义一种运算\cdot ,该运算满足以下性质:

  1. 封闭性:\forall a, b \in G, a \cdot b \in G
  2. 结合律:\forall a, b, c \in G, (a \cdot b) \cdot c = a \cdot (b \cdot c)
  3. 单位元:\exists e \in G, \forall a \in G, e \cdot a = a
  4. 逆元:\forall a \in G, \exists a^{-1} \in G, a \cdot a^{-1} = e

则称(G, \cdot)一起构成了一个

置换的概念

我们抽象并拓展“变换”这个概念,定义集合G \to G的一个双射为一个置换,举个例子:

\left(\begin{matrix}
1 &2 &3 &4 \\
2 &3 &4 &1
\end{matrix}\right)
这个置换表示将将集合中的1变成22变成33变成44变成1

置换群

稍微想一想可知,置换是可以成群的,称为置换群,对于这个置换群我们定义的运算就是置换的复合,对于这个运算我们考察其是否满足群的性质:

  1. 封闭性:大部分情况下应该是成立的
  2. 结合律:显然
  3. 单位元:\left(\begin{matrix} 1 &2 &3 &\cdots \\ 1 &2 &3 & \cdots \end{matrix}\right)
  4. 逆元:置换是集合到集合的双射,显然

置换群的引入允许我们使用群论来处理一开始的问题,说抽象了我们刚才的问题就是:求一个序列的集合X,在一个置换群G作用下本质不同的元素个数

Burnside引理

序列集合X在置换群G的作用下不同的序列数等于不动点的平均数

|X \setminus G| = \frac{1}{|G|} \sum_{g \in G} |X^g|
其中X^g = \{ x \in X | g(x) = x\}表示在置换g作用下不动点的集合,何为不动点呢?

例如对于置换\left(\begin{matrix} 1 &2 &3 &4 \\ 3 &4 &1 & 2 \end{matrix}\right)而言,序列(1, 2, 1, 2)就是它的一个不动点,因为它在置换后还是原来的样子。

我们接下来考察对于给定置换如何数出不动点的数目,从例子入手,比如说我们要求的是置换,\left(\begin{matrix} 1 &2 &3 &4 \\ 4 &3 &2 & 1 \end{matrix}\right)的不动点数目,我们挨个考虑这样的序列拥有什么样的性质:

  1. 第一位与第四位相同
  2. 第二位与第三位相同
  3. 第三位与第二位相同
  4. 第四位与第一位相同

显然我们最后得到了一个集合\{ \{1, 4\}, \{2, 3\} \},其中的每一个集合中的元素必须相等,考虑染色问题,对于每一组“相等”,我们可以把它整个染成两种颜色,因此这个置换共有2^2 = 4个不动点。

由于等号的传递性,如果我们逐位查看一个置换,把上下的两位并入同一个“相等集合”,那么最后我们一定会得到若干个这样的相等集合,每个“相等集合”都可以被独立地染成不同的颜色,由乘法定理,我们就得到了——

Polya定理

|X \backslash G| = \frac{1}{|G|} \sum_{g \in G} m^{c(g)}

其中m是颜色个数,c(g)是置换g形成的“循环数”,一个“循环”是指某一位经过不断地置换之后最后回到自身。可以看到Polya定理其实就是我们上面的推论,不需要更多解释。

循环同构序列计数

终于进入了我们的正题,循环同构序列计数,到这一步我们终于可以解答一开始的问题了:

Polya定理的精髓就在于求出一个置换的循环数,而朴素的枚举显然是不够的,我们考察所谓“循环同构”背后依赖的哪些置换,比如说对于一个长度为4的序列,以下置换后它与自身是循环同构的:

\left(\begin{matrix} 
    1 &2 &3 &4 \\ 
    2 &3 &4 &1 
\end{matrix}\right),
\left(\begin{matrix} 
    1 &2 &3 &4 \\ 
    3 &4 &1 &2 
\end{matrix}\right),
\left(\begin{matrix} 
    1 &2 &3 &4 \\ 
    4 &1 &2 &3 
\end{matrix}\right),
\left(\begin{matrix} 
    1 &2 &3 &4 \\ 
    1 &2 &3 &4 
\end{matrix}\right)
我们发现一个循环同构的置换就是把它的后面全部移到前面来,前面几位补到后面去

为了接下来的推导方便我们让位数起始于0

那么,对于一个偏移了k的置换,第a位被置换之后就到了(a + k) \bmod n的位置,其中n为序列长度,而对于一个长度为l的完整循环来说,它必须满足:

\begin{aligned}
a + kl &\equiv a &\pmod{n} \\
kl &\equiv 0 &\pmod{n}
\end{aligned}
显然当kl = \operatorname{lcm}(k, n)的时候l有最小值\frac{\operatorname{lcm}(k, n)}{k}。由于每一位都偏移了k位,因此置换中每个循环的长度都是相等的,而因为每个循环都cover了l个元素,因此循环的个数为:
\begin{aligned}
\frac{n}{l} &= \frac{n}{\frac{\operatorname{lcm}(k, n)}{k}} \\
&= \frac{kn}{\operatorname{lcm}(k, n)} \\
&= \gcd(k, n)
\end{aligned}
这个发现让我们免于枚举每个循环置换,相反,我们枚举\gcd(k, n)的每一个取值d,即n的因数,统计出\gcd(k, n) = d的所有k的个数,依Polya定理统一累加到答案上。

那怎么统计呢?枚举?

我们发现\gcd(k, n) = d等价于\gcd({k \over d}, {n \over d}) = 1,即小于n \over d的与n \over d互质的数个数,这显然等于\varphi({n \over d})

因此最后的答案就是:

\sum_{d|n} \varphi\left(\frac{n}{d}\right)m^d

完整的代码在此

背景

最近我在用Rust写一个ICMP隧道,因为ICMP包本身是不可靠的,于是需要在ICMP之上写一个可靠协议。一个完整的TCP协议栈显然过于臃肿了(何况也并没有现成的无IO的轮子),所以我就看上了skywind3000大佬的KCP协议——轻量、简洁、代码连我这种网路萌新都看得懂,实在是再好不过了。

一开始是我直接使用原版的C实现+FFI封装,在不开流控的情况下效果不错,但可惜原版的拥塞控制用的是最朴素的TCP Tahoe(作者也表明出于简洁性考虑不准备在标准实现当中使用复杂的流控算法),其在国际网络环境下的表现实在差强人意。想改进,但自己对于自己C的编程水平实在是不抱信心。既然Rust也是性能一流的系统编程语言,我最终还是决定用Rust再写了一个实现。本实现具有以下特点:

  • 相较于C实现进行了架构上的些许调整。
  • 在C实现的基础之上,使用链表+滚动数组优化大窗口下的发送性能。
  • 在C实现的基础之上,使用小根堆优化RTO计时器的效率,提升重传性能。
  • 将著名的BBR拥塞控制算法进行一定修改后试验性地运用到KCP中。

依赖的包

为了使编写更加简便,我们的实现依赖以下Rust crates:

  • bytes——简便的字节处理(代替原来C实现当中的encode_xxx/decode_xxx)。
  • num_enum——简化Rust枚举与字节的互相转换。
  • derivative——简化一些trait的实现。
  • thiserror——简化错误类型的定义。
  • rand——用于BBR随机相位初始化。

架构上的调整

常量与配置

相对于C实现,本实现大幅减少了常量的数量。最后仅剩的常量有五:

/// KCP包头大小
const OVERHEAD: u32 = 24;
/// 最大分段
const MAX_FRAGMENTS: u16 = 128;
/// KCP段类型
#[derive(Debug, Clone, Copy, TryFromPrimitive, IntoPrimitive)]
#[repr(u8)]
enum Command {
    Push = 81,
    Ack = 82,
    AskWnd = 83,
    TellWnd = 84,
}
/// BBR各阶段的增益
const BBR_GAIN_CYCLE: [usize; 8] = [5, 3, 4, 4, 4, 4, 4, 4];
/// BDP增益的分母,见后文
const BDP_GAIN_DEN: usize = 1024;

常量少了,变量自然就多了,原来C实现的常量在本实现中成为可配置项:

/// 大部分配置的意思如字面
#[derive(Clone, Debug, Deserialize, Derivative)]
#[derivative(Default)]
pub struct Config {
    #[derivative(Default(value = "536"))]
    pub mtu: u32,
    #[derivative(Default(value = "200"))]
    pub rto_default: u32,
    #[derivative(Default(value = "100"))]
    pub rto_min: u32,
    #[derivative(Default(value = "6000"))]
    pub rto_max: u32,
    #[derivative(Default(value = "7000"))]
    pub probe_min: u32,
    #[derivative(Default(value = "120000"))]
    pub probe_max: u32,
    #[derivative(Default(value = "1024"))]
    pub send_wnd: u16,
    #[derivative(Default(value = "1024"))]
    pub recv_wnd: u16,
    #[derivative(Default(value = "40"))]
    pub interval: u32,
    /// 若一个包重传dead_link_thres次后依然失败,则视作底层链路失效。
    #[derivative(Default(value = "20"))]
    pub dead_link_thres: u32,
    /// nodelay模式下, rto_min = 0且rto在重传失败后不指数增长。
    #[derivative(Default(value = "false"))]
    pub nodelay: bool,
    /// stream模式下, 多个数据包可以被合并在同一段内从而减少开销。
    #[derivative(Default(value = "false"))]
    pub stream: bool,
    /// 如果指定,则一个包在fast_resend_thres个在其之后的包ACK之后会直接重传
    #[derivative(Default(value = "None"))]
    pub fast_resend_thres: Option<u32>,
    /// 快速重传的次数上限
    #[derivative(Default(value = "None"))]
    pub fast_resend_limit: Option<u32>,
    /// 是否启用BBR控制算法
    #[derivative(Default(value = "false"))]
    pub bbr: bool,
    /// BBR中RTprop(往返时间)滑动窗口的时间长度(单位:毫秒)
    #[derivative(Default(value = "10000"))]
    pub rt_prop_wnd: u32,
    /// BBR中BtlBw(瓶颈带宽)滑动串口的长度(单位:RTT)
    #[derivative(Default(value = "10"))]
    pub btl_bw_wnd: u32,
    /// BBR中一次RTT/RTprop探测的时间(单位:RTT),减少该值可以减轻RTT探测对于流量的影响。
    #[derivative(Default(value = "200"))]
    pub probe_rtt_time: u32,
    /// BDP增益,见后文
    #[derivative(Default(value = "1024"))]
    pub bdp_gain: usize,
}

impl Config {
    pub fn mss(&self) -> usize {
        (self.mtu - OVERHEAD) as usize
    }
}

impl ControlBlock {
    pub fn new(conv: u32, config: Config) -> ControlBlock {
        ...
    }
}

异常类型

KCP原本的C实现仅使用负数表达异常值,虽简介但其含义并不明晰,在本实现中我们对于异常进行了清晰定义:

#[derive(Debug, Error)]
pub enum Error {
    #[error("packet to be sent too large to be fragmented")]
    OversizePacket,
    #[error("incomplete KCP packet")]
    IncompletePacket,
    #[error("invalid KCP command: {0}")]
    InvalidCommand(u8),
    #[error("empty queue (try again later)")]
    NotAvailable,
    #[error("wrong conv. (expected {expected}, found {found})")]
    WrongConv { expected: u32, found: u32 },
}

pub type Result<T> = std::result::Result<T, Error>;

以上异常类型还有更精确的空间,但是目前应该已经堪堪够用了。

发包方式

在原来的C实现在发包时直接调用callback,其优点是简洁,但其缺点在于callback的运行时间不定以及异常处理不明对运行产生的影响。何况在Rust当中安全存储callback需要和borrow checker拼命。在本实现中,我们将flush出去的包暂存在一个队列中,然后通过外部不断poll的方式拉出去。一方面,主动poll的方式和底层收到包时的push呼应;另一方面,这有助于分离底层发包和KCP逻辑本身,是“无IO/Sans IO”理念的一种体现。缺点是缓存队列可能会膨胀得厉害。当中tradeoff见仁见智。

impl ControlBlock {
    ...
    /// 底层收包push
    pub fn input(&mut self, mut data: &[u8]) -> Result<usize> { ... }
    /// 底层发包poll
    pub fn output(&mut self) -> Option<Vec<u8>> { ... }
}

去除checkupdate

这是一个比较大胆的改动,未必适合所有情形。去除的原因是在数据结构的优化下计算重传、更新发送窗口的开销大幅度减小,已经可以在每一次调用inputsend的时候进行一次,没有必要去不断checkupdate。上层只需要按照固定的时间间隔调用flush就行了。我进行这样的设计是为了简化上层的代码,而且我的应用情形恰好是高流量的反正都要一直flush,也无所谓。

flush的代码也其实很简单:

pub fn flush(&mut self) {
    self.sync_now(); // 更新now
    self.flush_probe(); // 更新窗口探测
    self.flush_push(); // 计算重传以及更新发送窗口
    self.flush_ack(); // 发ACK
    if !self.buffer.is_empty() {
        let mut new_buf = Vec::with_capacity(self.config.mtu as usize);
        std::mem::swap(&mut self.buffer, &mut new_buf);
        self.output.push_back(new_buf);
    }
}

既然sync_nowflush_pushinputsend当中都可以廉价地调用,那为什么还需要不断checkupdate呢?直接调用flush了事。

如果要参考有checkupdate的实现,可以参照早些时候的commit

窗口数据结构的改进

KCP原版的实现中发送/接收的队列/窗口全部使用队列作为数据结构,这固然使得代码变简单了,但也一定程度上降低了性能:在队列中查找KCP段最差需要线性时间,这在某些情形下未必是最优的。在本实现中,我们优化数据结构,以最优的复杂度实现发送/接受窗口需要的若干操作:

  • \mathcal{O}(1)以分段的序号为键查找分段。
  • \mathcal{O}(1)插入分段。
  • \mathcal{O}(1)以分段的序号为键查找分段。
  • \mathcal{O}(1)以分段的序号为键删除分段。
  • \mathcal{O}(1)查询/弹出最早插入的分段(在发送窗口中,最早插入的分段自然是序号最小的分段)。
  • \mathcal{O}(k)遍历以插入顺序为序,某分段的所有k个前驱(在发送窗口中,分段插入顺序即序号顺序,因此该操作可直接用于快速重传的计算)。
  • \mathcal{O}(1)查询大小。

考虑到任何时刻窗口内分段序号之差不会大于窗口大小这一常数,符合上述要求的数据结构就可以用链表+滚动数组高效实现。代码不长,百行左右:

struct Element<T> {
    /// 前驱下标
    prev: usize,
    /// 后继下标
    next: usize,
    data: T,
}

pub struct Window<T> {
    size: usize,
    entry: Vec<Option<Element<T>>>,
    end: Option<usize>,
    len: usize,
}

impl<T> Window<T> {
    pub fn with_size(size: usize) -> Self {
        Self {
            size,
            entry: (0..size).map(|_| None).collect(),
            end: None,
            len: 0,
        }
    }

    pub fn is_empty(&self) -> bool {
        self.end.is_none()
    }

    pub fn get_mut(&mut self, index: usize) -> Option<&mut T> {
        match self.entry[index % self.size].as_mut() {
            Some(elem) => Some(&mut elem.data),
            None => None,
        }
    }

    pub fn push(&mut self, index: usize, data: T) {
        let index = index % self.size;
        if self.entry[index].is_some() {
            return;
        }
        self.entry[index] = Some(match self.end {
            Some(prev) => {
                let prev_elem = self.entry[prev].as_mut().unwrap();
                let next = prev_elem.next;
                prev_elem.next = index;
                self.entry[next].as_mut().unwrap().prev = index;
                Element { prev, next, data }
            }
            None => Element { prev: index, next: index, data },
        });
        self.end = Some(index);
        self.len += 1;
    }

    pub fn remove(&mut self, index: usize) -> Option<T> {
        let index = index % self.size;
        let elem = self.entry[index].take()?;
        let (prev, next) = (elem.prev, elem.next);
        self.entry[index] = None;
        self.len -= 1;
        if index == self.end.unwrap() {
            if prev == index {
                self.end = None;
                return Some(elem.data);
            } else {
                self.end = Some(prev);
            }
        }
        self.entry[prev].as_mut().unwrap().next = next;
        self.entry[next].as_mut().unwrap().prev = prev;
        Some(elem.data)
    }

    pub fn front(&self) -> Option<&T> {
        self.end.map(|end| {
            let head = self.entry[end].as_ref().unwrap().next;
            &self.entry[head].as_ref().unwrap().data
        })
    }

    pub fn pop_unchecked(&mut self) -> T {
        let end = self.end.unwrap();
        let head = self.entry[end].as_ref().unwrap().next;
        self.remove(head).unwrap()
    }

    pub fn len(&self) -> usize {
        self.len
    }

    pub fn for_preceding(&mut self, index: usize, mut action: impl FnMut(&mut T)) {
        let mut index = index % self.size;
        index = match self.entry[index].as_ref() {
            Some(elem) => elem.prev,
            None => return,
        };
        while index != self.end.unwrap() {
            let elem = self.entry[index].as_mut().unwrap();
            action(&mut elem.data);
            index = elem.prev;
        }
    }
}

因为滚动数组是连续空间,在内存布局上相较于链表对于缓存更加友好,所以速度应该还可以再快一点。唯一的不足是unwrap有点多看着心惊肉跳,并且用指针可能会比用下标快一丁点,但是用Rust写数据结构大约就是这个尿性。

有这个打底,窗口大小开到8192实测是一点问题都没有的,更大的没试过。

但KCP的设计本质上是不适合大流量的,因为快速重传无论如何优化数据结构最坏的线性复杂度就在那里无法消除,除非可以限制快速重传向后看的范围,但后者又削弱了快速重传的意义与效用。

重传计时器的改进

原版的KCP实现在check的时候需要遍历发送窗口来确定最近的重传时间,在flush的时候又要遍历才能重传,这在窗口较大的时候显然是比较吃性能的。原作者记得在issues里的讨论中提过可以用时间轮进行优化。诚然,时间轮是最好的方案,但是实现起来较为复杂。因此,本实现使用借助Rust的标准库实现起来相对简单的小根堆进行优化:

use std::cmp::Reverse;
use std::collections::BinaryHeap;

pub struct Timer(BinaryHeap<Reverse<u64>>);

impl Timer {
    pub fn with_capacity(capacity: usize) -> Self {
        Self(BinaryHeap::with_capacity(capacity))
    }

    pub fn schedule(&mut self, ts: u32, sn: u32) {
        self.0.push(Reverse(((ts as u64) << 32) | sn as u64));
    }

    /// 获取截止到now发生的 一个 事件,应该重复调用
    pub fn event(&mut self, now: u32) -> Option<(u32, u32)> {
        let key = (now as u64 + 1) << 32;
        match self.0.peek() {
            Some(&Reverse(val)) if val < key => {
                let sn = val & (u32::max_value() as u64);
                let ts = val >> 32;
                self.0.pop();
                Some((ts as u32, sn as u32))
            }
            _ => None,
        }
    }
}

计时器只需要存时间和分段序号即可。调用的代码如下:

fn flush_push(&mut self) {
    // ... 省去流控以及把队列里的分段加入发送窗口的部分
    let mut send_buf = std::mem::take(&mut self.send_buf);
    while let Some((ts, sn)) = self.timer.event(self.now) {
        if sn < self.send_una || sn >= self.send_nxt {
            continue; // 分段被ACK于是不在发送窗口里了,自然跳过
        }
        if let Some(seg) = send_buf.get_mut(sn as usize) {
            if ts == seg.ts {
                seg.ts = self.prepare_send(seg); // 更新RTO并计算下一次重传的时间
                seg.ts_last_send = ts;
                self.dead_link |= seg.sends >= self.config.dead_link_thres;
                self.flush_segment(Command::Push, seg.frg, seg.sn, ts, seg.payload.len());
                self.buffer.extend_from_slice(&seg.payload);
                self.timer.schedule(seg.ts, seg.sn); // 安排下一次重传
            }
        }
    }
    self.send_buf = send_buf;
}

由于查看小根堆堆顶是\mathcal{O}(1)的,因此在没有重传的时候flush_push的开销确实很小。足以在inputsend时都调用一次。真的要重传时,更新小根堆的时间复杂度也是对数级别的,这就给去除checkupdate提供了基础。

BBR

最后的改进是用BBR取代了KCP原版实现中朴素的基于丢包的流控算法。

我试图参照原论文实现BBR,但因为计时精度的问题packet pacing是做不到了。实现的部分有

  • 基于单调队列的滑动窗口BtlBw max-filter。
  • 基于单调队列的滑动窗口RTprop min-filter。
  • BBR状态机。
  • 基于以上三者计算inflight limit进行流控。

和BBR有分歧的一点在于在ProbeRTT状态采用BDP的一半作为拥塞窗口而不是原文的4个包。

此外,本实现只对只传输一次的分段计算BBR的各项参数,如RTT,带宽,更新各个filter等。原因是实际上大部分的分段都在看到包头的UNA之后就被ACK掉了而不是被单独的ACK包ACK的。ACK包带有分段的序号与时间戳,所以可以清楚知道ACK的是哪一次传输,但被UNA ACK掉的就不清楚,唯一的例外是分段只被传输了一次。如果对于多次传输的包仍然直接计算BBR,那么万一ACK恰好在重传之后到达,那么误算出的RTT就非常小,导致RTprop非常小,进而BDP非常小,整个BBR就堵住了。诚然,可以把每一次传输的时间戳都存起来,然后在被UNA时选择最接近now - srtt的传输,但这就增加了代码的复杂度。考虑到丢包的毕竟是少数,如此未必会有特别大的优势(实现这个的代码在这个commit,也可作为参考)。

实测效果还行。一个比较重要的问题是在间歇性flush的情况下对于带宽的计算并不准确甚至有低估的倾向,往往导致BDP过于保守,即使ProbeBW状态有一个1.25x的激进phase也解决不了问题。我想到的一个解决方案是将拥塞窗口不是简单的设置为BDP,而是乘上一个增益,也就是配置里面的bdp_gain。为了避免浮点数运算,bdp_gain使用1024为基数。一般来说设置成1280就差不多了——其实就是主动创造轻微的拥塞来确保占有带宽。这对于其他的TCP连接固然有些不公平的,是否采用见仁见智。

真·背景

其实一开始,只是为了编译的时候能够不带着unsafe和C编译器而选择把原版实现移植到Rust,当时代码大部分几乎一模一样。

后来引入BBR魔改了一回。

为了代码更Rust魔改了一回。

再后来做配置分离魔改了一回。

再后来数据结构优化又魔改了一回。

来来回回地改,到最后除了架子还和原版实现相似,内部的代码已经大变样了。

但是在ICMP隧道上试验下来仍然不是最令人满意,CPU占用仍然不少,带宽仍然不能跑满,goodput仍然不高。

我不知道是我应用层以及底层的代码写的有问题,还是KCP本身就不是很适合高流量大窗口的应用场景。

或许二者兼有之?

前两天突然想起了QUIC,找了一下,Cloudflare有一个优秀的QUIC实现,是Rust的,而且是Sans IO的。大概QUIC才是最适合我的应用情形的吧。我准备这几天试验一下。

或许之后就转QUIC了呢?(笑)

那我魔改的KCP就放在这吃灰?

于是乎,就有了这篇文章。

前几天借着Github的学生优惠在name.com上嫖了一个域名,这个域名是带SSL的,而恰巧家里的主机因为443端口被封的缘故无法通过Let's Encrypt获得证书,所以把这个域名给家里的主机升级HTTPS再合适不过了。

唯一的问题是name.com只提供静态DNS服务,而服务器放在家里自然是动态IP的。虽然IP不是经常换但是如果换了IP没有及时更新记录就会出问题,何况手动更新记录也有点烦。有两个方案:

  1. 我目前使用的是花生壳的DDNS,这个的记录是动态更新的。我可以在新域名下面新建一条CNAME记录指向动态域名。这样的好处是省事,坏处是可能会增加DNS解析的时间,我目前还不清楚HTTPS要不要求CNAME指向的域名也有证书,如果要求的话这个方法就更不行了。
  2. name.com作为一家比较大的域名商有自己的API以及完备的文档,可以自己写一个定时更新脚本来实现类似DDNS的功能。

经过考虑之后我选择后者。自己写的脚本如下:

#!/bin/sh

domain='<DOMAIN>'
credential='<ACCESS TOKEN>'
ttl=300
interval=10

echo 'Querying type A record ID...'
rec_id=$(curl -su $credential "https://api.name.com/v4/domains/${domain}/records" | jq '.records|map(select(.type=="A"))|.[0].id')
echo 'Found type A record ID:' $rec_id

while true; do
    res_ip=$(host $domain | grep -ohP '\b(?:\d{1,3}\.){3}\d{1,3}\b')
    real_ip=$(curl -s myip.ipip.net | grep -ohP '\b(?:\d{1,3}\.){3}\d{1,3}\b')
    echo 'Resolution IP:' $res_ip 'Real IP:' $real_ip
    if [[ $real_ip != $res_ip ]]; then
        echo "Resolution mismatch! Updating record..."  
        req_data="{\"type\":\"A\",\"fqdn\":\"${domain}.\",\"answer\":\"${real_ip}\",\"ttl\":${ttl}}"
        echo "Request data:" $(echo $req_data | jq '.')
        req_res=$(curl -m 30 -su $credential "https://api.name.com/v4/domains/${domain}/records/${rec_id}" \
            -X PUT -H 'Content-Type: application/json' --data $req_data | jq '.')
        echo "Request result:" $req_res 
        sleep $ttl
    fi
    sleep $interval
done

其实很简单,主要分为一下几个部分:

  1. 记录id的获取。name.com对每一个记录都分配了一个id以便于API操作,这个id在网站管理面板上是不可见的,因此需要在运行时查询,命令为curl -su $credential "https://api.name.com/v4/domains/${domain}/records"。查询之后需要解析JSON,这里我使用的是jq这个第三方JSON parser。
  2. 当前DNS解析的IP。这个使用host结合grep提取IP字符串即可。
  3. 获取本机真实IP。这个接口就多了,我这里用的是myip.ipip.net的接口。据我所知还有接口是直接返回IP字符串的,还可以后处理的功夫。
  4. 如果解析IP和真实IP不符,那就调用API更新记录。注意在更新完之后最好等待TTL的时间以避免更新生效前多次更新。

写完这个脚本再写一个配套的systemd service file,然后systemctl挂在后台运行就好了。运行到现在效果非常好。如果有name.com的域名同时也有类似需求的或许可以参考一下。

Intro

Redstone has been a core element in the game Minecraft for quite some years. It is presumably the most untrivial one as well: while anyone could master nearly all Minecraft mechanics through experiences, it takes not only experience but also ingenuity to design a good redstone circuit. Few of us are bold enough to claim "I master redstone", even after playing Minecraft for a decade.

So here comes the question: Can the design of redstone circuits, the core of Minecraft automations, be automated? and If so, how? (appreciate how meta this is :)

Theoretically, the answer is "yes...but". Minecraft has a finite world, and each position has a finite number of possible blockstates. We can write a program to enumerate all possible placements of blocks until we find some placement corresponding to the desired redstone circuit. However, this needs exponential time and we may need to wait for a century before it could give us, for instance, a decent piston door. Moreover, if a circuit involves manipulation of entities (which could have infinite many states), then we are easily screwed.

Well, perhaps it is difficult to let a program design any redstone circuit. But there is indeed a subset of redstone circuits whose design can very likely be automated -- computational redstone circuits, aka. logic gates, calculators, CPUs etc. Why? Because software that design their real world electronic counterparts are readily available, i.e. EDA software.

As a high school student I, of course, know little about the inner workings of real-world EDA applications (and there doesn't seem to be a lot of resources out there). I am convinced that this problem is NPC (further articulated below), so designing an efficient polytime combinatorial algorithm doesn't seem plausible. That said, what about reducing this problem to some other NPC problems which we can solve relatively quickly with optimized algorithms / heuristics -- say, ILP? This is what I am trying to do here.

Formulating the problem

"Designing computational redstone circuit automatically" is a vague idea, so it is necessary that we know what this truly means.

What's the input?

The input should describe the intended functionality of a circuit. Recall how we usually describe a circuit: we draw a circuit diagram. I here characterize a redstone circuit diagram by the assumptions and constraints below:

  1. A circuit contains two parts: wires and components.
  2. Components are the primitives of a circuit. E.g. A torch or a wire junction.
  3. A component has interfaces, either in or out, as where the component receives signals from and sends signals to.
  4. A wire connects an out-interface from a component ("source") to an in-interface of another component ("target").
  5. Wires are directed.
  6. Components are independent. i.e. they do not interfere with other components in any way other than being connected by wires from interfaces.

A circuit diagram can be represented in a directed graph, with components as vertices and wires as edges. Source/target interfaces as extra information stored on edges.

What's the output?

We want our program to tell us how the circuit we described in the input can be built in the Minecraft world. Therefore, we could define the the output to be a set of position - blockstate pairs, (which, in implementation, can be stored in a schematic file).

However, we don't want to jump straight from a circuit diagram to a detailed Minecraft schematic because that means taking interference between components, quasi connectivity, update order -- basically everything that makes redstone engineering complex -- into consideration in the first place.

Instead, we could first build our circuit in an ideal world, in which we forget about all those factors above, and then convert the ideal placement into an actual Minecraft schematic.

What's an idea world?

  1. A circuit consists of multiple ideal blocks.
  2. A component fully occupies a set of ideal blocks, some of which are its interfaces. How many and which blocks a certain type of component occupies depend on its size in Minecraft and how we plan to convert the ideal placement to a real schematic.
  3. A wire is a chain of ideal blocks, where any adjacent two share a face. The first block is always the source interface and the last is always the target interface.
  4. Exclusiveness: All components and all wires (ignoring their first and last block) mustn't overlap.
  5. Mutual Independence: Unless both blocks are occupied by the same component / wire, anything in two adjacent cells do not interfere with each other. To meet this constraint, we could say (rather conservatively) that an ideal block should map to 3\times 2 \times 3 Minecraft blocks.
  6. Wire junctions are special components and are exceptions to rule 2 and 3. A wire junction always have three interfaces (1 in & 2 outs, or 2 ins & 1 out). Multiple junctions can overlap and they can overlap with an interface of some component.
  7. There are times when we want to fix the location of some components in the input. These components are usually just placeholders that mark the position of IO. (We don't want to produce a circuit with an unreachable input/output in the center of everything else, right?)

The Objective

  1. The circuit represented by the output must have the same functionality as described by the input circuit diagram.
  2. The delay of the circuit should be minimized.

Example: the AND gate

Let's see how we design a simple AND gate.

Suppose the only primitive components we have are NOT gate (torch), wire junction, and IO placeholder. The circuit diagram of AND gate is:


Reduction to Integer Programming

Once the problem statement is clear, we can proceed to convert the problem into an integer programming. Let's start by defining some notations. First, assume there are n components and m wires in a circuit.

The circuit we are designing obviously must meet some space constraints. Hence, let F \subset \mathbb{Z}^3 be the set of all feasible ideal block coordinates.

Then we have different types of components, which can be placed with different orientations (at most 6). Some component may not be placed in a certain direction due to limitations in Minecraft mechanisms (such contraptions are often known as "directional"). Let D_i be the set of feasible orientations of the ith component.

For each component i and a feasible orientation d\in D_i, let C_{i,d} \subset \mathbb{Z}^3 be the set of ideal block coordinates this component occupies when placed in this direction. For convenience, let (0,0,0) be a tight lower bound of the three components in this set. We assume that the shape of a component, after fixing its type and orientation, is translation-invariant.

It's finally time to introduce some variables. Let binary variables x_{i,d,\mathbf{u}} (d\in D_i, \mathbf{u} \in F) denote whether the ith component is placed with orientation d at coordinate \mathbf{u}.

Every component should be placed only once. This gives the constraint

\sum_{d\in D_i}\sum_{\mathbf{p} \in F} x_{i,d,\mathbf{p}} = 1, \quad \forall i = 1,\cdots,n.
And components should not overlap! Consequently,
\sum_{i=1}^n\sum_{d\in D_i}\sum_{\substack{\mathbf{c}\in C_{i,d} \\
\mathbf{p}-\mathbf{c} \in F}} x_{i,d,\mathbf{p}-\mathbf{c}} \le 1,\quad \forall \mathbf{p} \in F.
This constraint applies for each location \mathbf{p} in the feasible space. We enumerate all components and their possible placements that would occupy \mathbf{p}, making sure that their corresponding variables sum to at most 1.

Now we consider wires. For convenience, let N(\mathbf{u}) denote \mathbf{u}'s neighbors in F, N(\mathbf{u})\le 6.

Let binary variables y_{i,\mathbf{u}, \mathbf{v}}(\mathbf{v} \in N(\mathbf{u})) denote whether the ith wire has a segment from \mathbf{u} to \mathbf{v}, and let a_i and b_i be the indices of the source and destination component of wire i.

Furthermore, let \mathbf{a}_{i,d}(d\in D_{a_i}) be the coordinate of the source interface relative to the source component's (a_i) location when the source component of the ith wire is placed with orientation d. Define \mathbf{b}_{i,d}(d\in D_{b_i}) in a similar manner, denoting the coordinate of the destination interface relative to the destination component's location.

With these notations, the constraint that enforces the connectivity of wires can be written as

\sum_{\mathbf{v} \in N(\mathbf{u})} \left(y_{i,\mathbf{u}, \mathbf{v}} - y_{i,\mathbf{v}, \mathbf{u}}\right) = \sum_{\substack{d\in D_{a_i}\\ \mathbf{u} - \mathbf{a}_{i,d} \in F}} x_{a_i,d,\mathbf{u} - \mathbf{a}_{i,d}} - \sum_{\substack{d\in D_{b_i}\\ \mathbf{u} - \mathbf{b}_{i,d} \in F}} x_{b_i,d,\mathbf{u} - \mathbf{b}_{i,d}}, \quad \forall \mathbf{u}\in F, i=1,2,\cdots,m.
Check: If \mathbf{u} is the starting location of wire i. both LHS and RHS are 1; If \mathbf{u} is the ending location of wire i. both LHS and RHS are -1; Otherwise both sides equal 0.

Meanwhile, given that a location is either occupied by a component or a wire but never both, and that any two wires shall not overlap, we have the following constraint

2\sum_{i=1}^n\sum_{d\in D_i}\sum_{\substack{\mathbf{c}\in C_{i,d} \\
\mathbf{p}-\mathbf{c} \in F}} x_{i,d,\mathbf{p}-\mathbf{c}} + \sum_{i=1}^m\sum_{\mathbf{v} \in N(\mathbf{p})} \left(y_{i,\mathbf{p}, \mathbf{v}} + y_{i,\mathbf{v}, \mathbf{p}}\right) \le 2,\quad \forall \mathbf{p} \in F.
The intuition: Usually if a wire (say, i) passes through \mathbf{p}, \sum_{\mathbf{v} \in N(\mathbf{p})} y_{i,\mathbf{p}, \mathbf{v}} + y_{i,\mathbf{v}, \mathbf{p}} = 2, so the RHS should be 2. At the same time, we expect that when some component has occupied \mathbf{p}, no wires should occupy \mathbf{p} again. We then need to use indicator
\sum_{i=1}^n\sum_{d\in D_i}\sum_{\substack{\mathbf{c}\in C_{i,d} \\
\mathbf{p}-\mathbf{c} \in F}} x_{i,d,\mathbf{p}-\mathbf{c}},
which is at most 1 by our second constraint, so we give it a coefficient of 2 to balance it with the wire term.

The above constraint, however, isn't quite right: Wires and components do over lap at the interface block. To fix this, we amend the constraint to

2\sum_{i=1}^n\sum_{d\in D_i}\sum_{\substack{\mathbf{c}\in C_{i,d} \\
\mathbf{p}-\mathbf{c} \in F}} x_{i,d,\mathbf{p}-\mathbf{c}} + \sum_{i=1}^m\sum_{\mathbf{v} \in N(\mathbf{p})} \left(y_{i,\mathbf{p}, \mathbf{v}} + y_{i,\mathbf{v}, \mathbf{p}}\right) \le 2 + \sum_{i=1}^m\left(\sum_{\substack{d\in D_{a_i}\\ \mathbf{p} - \mathbf{a}_{i,d} \in F}} x_{a_i,d,\mathbf{p} - \mathbf{a}_{i,d}} +\sum_{\substack{d\in D_{b_i}\\ \mathbf{p} - \mathbf{b}_{i,d} \in F}} x_{b_i,d,\mathbf{p} - \mathbf{b}_{i,d}} \right),\quad \forall \mathbf{p} \in F.
The complicated summation term on the right will be greater than 0 if \mathbf{p} is an interface block for some component and wire.

These are all the constraints that I had thought of. Together, the integer programming instance is

\begin{aligned}
    \text{minimize}\ &\sum_{i=1}^m\sum_{\mathbf{u}\in F}\sum_{\mathbf{v}\in N(\mathbf{u})} y_{i, \mathbf{u}, \mathbf{v}}, \\
    \text{subject to}\ &\sum_{d\in D_i}\sum_{\mathbf{p} \in F} x_{i,d,\mathbf{p}} = 1, &\forall i = 1,\cdots,n;\\
    &\sum_{i=1}^n\sum_{d\in D_i}\sum_{\substack{\mathbf{c}\in C_{i,d} \\
\mathbf{p}-\mathbf{c} \in F}} x_{i,d,\mathbf{p}-\mathbf{c}} \le 1,&\forall \mathbf{p} \in F;\\
&\sum_{\mathclap{\mathbf{v} \in N(\mathbf{u})}} \left(y_{i,\mathbf{u}, \mathbf{v}} - y_{i,\mathbf{v}, \mathbf{u}}\right) = \sum_{\mathclap{\substack{d\in D_{a_i}\\ \mathbf{u} - \mathbf{a}_{i,d} \in F}}} x_{a_i,d,\mathbf{u} - \mathbf{a}_{i,d}} - \sum_{\mathclap{\substack{d\in D_{b_i}\\ \mathbf{u} - \mathbf{b}_{i,d} \in F}}} x_{b_i,d,\mathbf{u} - \mathbf{b}_{i,d}}, &\quad \forall \mathbf{u}\in F, i=1,\cdots,m;\\
&2\sum_{i=1}^n\sum_{d\in D_i}\sum_{\substack{\mathbf{c}\in C_{i,d} \\
\mathbf{p}-\mathbf{c} \in F}} x_{i,d,\mathbf{p}-\mathbf{c}} + \sum_{i=1}^m\sum_{\mathbf{v} \in N(\mathbf{p})} \left(y_{i,\mathbf{p}, \mathbf{v}} + y_{i,\mathbf{v}, \mathbf{p}}\right) \\
&\quad \le 2 + \sum_{i=1}^m\,\sum_{\mathclap{\substack{d\in D_{a_i}\\ \mathbf{p} - \mathbf{a}_{i,d} \in F}}} x_{a_i,d,\mathbf{p} - \mathbf{a}_{i,d}} +\sum_{i=1}^m\,\sum_{\mathclap{\substack{d\in D_{b_i}\\ \mathbf{p} - \mathbf{b}_{i,d} \in F}}} x_{b_i,d,\mathbf{p} - \mathbf{b}_{i,d}},& \forall \mathbf{p} \in F.
\end{aligned}
The objective is simply the sum of all y variables because it denotes the total length of all wires, which we aim to minimize if we want the circuit to be compact.

Keep in mind, however, that these constraints are merely the necessary conditions of a valid circuit. Given the complexity of the constraints, it is difficult to reason if all circuit configurations conforming to these constraint are valid circuits. I suppose the best way to check is to write some code and see if it gives what we what.

Preliminary Experiments

I did some very rudimentary test of my ILP formulation above with the Gurobi solver and ~1kLoC of Kotlin. The goal was to route the simple circuit described below:

val in1 = ComponentNode(FixedLocation(Vec3(0, 0, 0)))
val in2 = ComponentNode(FixedLocation(Vec3(1, 0, 0)))
val out = ComponentNode(FixedLocation(Vec3(0, 0, 3)))
val a = ComponentNode(NotGate)
val b = ComponentNode(NotGate)
val c = ComponentNode(NotGate)
val d = ComponentNode(NotGate)
val e = ComponentNode(NotGate)

val circuit = MutableCircuitGraph()
// Connect the IO interface of in1 to the IN interface of a
circuit.addWire(in1, FixedLocation.IO, a, NotGate.IN)
circuit.addWire(in2, FixedLocation.IO, a, NotGate.IN)
circuit.addWire(in1, FixedLocation.IO, b, NotGate.IN)
circuit.addWire(in2, FixedLocation.IO, c, NotGate.IN)
circuit.addWire(a, NotGate.OUT, b, NotGate.IN)
circuit.addWire(a, NotGate.OUT, c, NotGate.IN)
circuit.addWire(b, NotGate.OUT, d, NotGate.IN)
circuit.addWire(c, NotGate.OUT, d, NotGate.IN)
circuit.addWire(d, NotGate.OUT, e, NotGate.IN)
circuit.addWire(e, NotGate.OUT, out, FixedLocation.IO)

This describes an XOR gate with NOT gate as the only primitive. It has eight components and 10 wires -- decently complex. A FixedLocation component takes up 1 block of space and a NOT gate takes up 2 blocks, with the input interface in one block and the output in another. NOT gates can only be placed horizontally (i.e. 4 possible directions.)

Time to solve it:

val solver = MIPLayoutSolver()
val result = solver.solve(circuit, Vec3.inBox(Vec3(3, 2, 10)))

I gave the solver a lot of slack: the feasibility space is 3\times 2 \times 10, considerably larger than what's actually needed. Anyway, the code ran very fast (thanks to Gurobi). To visualize the result I wrote some extra code to generate an HTML which internally calls THREE.js:

Complicated, isn't it? The rightmost magenta block is the output component out. The leftmost magenta block on the bottom level is the input component in1. The cyan block next to in1 is in2. And then blocks of the same color represents a NOT gate, and arrows represent wires. The visualization is a mess but I couldn't do better without writing significantly more code. If you try hard to follow the the arrows then you will find that this routing is actually correct and does the exact thing the Kotlin description above describes -- and indeed, the circuit is amazingly compact.

Yay!

(Note that here a component's interface is connected to multiple wires. This might not be a good thing because by our constraints, if the fan in/out of any component is greater than 6, the constraint can not be satisfied. In an attempt to work around this I wrote a preprocessor that automatically introduces OR gates in circuit graphs when necessary. However, it turns out that with this modification things no longer work and the solver starts giving bogus results for some reason. I didn't have time to investigate into this issue, but this could be an error in my ILP formulation somewhere!)

Afterwords -- Two Years Later

Aha, I am finally here!

I did this project in 2020 but I somehow left this article unfinished for two years. What a shame! I spent an entire evening doing some archaeological work today and dug out my code and notes for this project (and confirmed they did work). The article is complete -- finally -- but please tolerate some inconsistencies in the writing. Everything after the AND gate example was written in 2022. Hopefully this article is helpful.

In retrospect, this small project is by no means complete. The test was too simple and I didn't even make a Minecraft circuit out of the result! Probably too much handwaving. I could have of course went much much deeper if I hadn't been busy with my college applications back then. I wrote the first part of this article on Nov. 25, 2020. A couple of days later I would have found myself rejected by Cornell. That was very frustrating then and I had to write more and more essays -- I guess that was when this unfinished article got buried somewhere.

Anyway. it had been great fun doing this project and I learned a ton of interesting stuff despite the shallow progress on the surface. I still vividly recollect the moment I found in my inbox an email from the Gurobi China sales team because I was the first high school student in China to apply for an academic license and they really wanted to confirm I was not joking. I said that I was serious and several days later I got my license -- that reply email literally made my day back then. Things have changed. When I re-install Gurobi to re-run my code last night, the Gurobi server recognized my MIT IP and gave me a license instantly. That was fast, but felt... different, certainly not as rewarding as before. I no longer play Minecraft as much as I did now, or do these fun projects as much as I had done. The decline in motivation on these things doesn't seem to be a good sign to me. I can't blame everything on the course load or UROP or anything else. It is a problem I am actively seeking a solution to.