0%

数字图像处理小作业1

最近软件工程课要我们用博客记录项目。我感觉这样做挺好的,一来可以让项目的开发更规范化,二来可以给博客和github添加一些内容,不让大学生活看起来很空虚。以后我所有项目作业尽量都建立仓库和博客。

代码仓库

https://github.com/SingleZombie/RGB-CMY-HSI-transformer-homework

需求分析

  • 原图像一张RGB(彩色图像)
  • 分别在RGB,CMY,和HSI三个空间通过调整颜色实现一种风格化转换(风格自选)
  • 需要提交源码以及不超过半页的代码说明(说明使用的库和参数)
  • 压缩包内文件命名方式为:
    -图像命名:
    0(此原图像)
    1RGB,2CMY,3HSI(此为处理结果)
    -文档命名:学号_姓名_1
    -代码包命名:学号_姓名_1

以上是老师给的作业要求。据说老师上课的时候还提到,三种通道的风格化结果要看起来相等。

我对风格化的理解是,改变某种彩色表示下某个分量后得到的结果。

经总结,需求有以下内容:

  • 实现RGB到CMY和HSI的转换
  • 在某种表示下改变图片风格,并调整另两种表示下的结果使得图片看起来相同
  • 写说明文档

技术学习

OpenCV入门

OpenCV安装与配置

在本作业中,我打算用OpenCV来读、写、显示图像文件。在一些任务开始前,我要先学OpenCV的用法。

我是在Win10下用VS2017编程,生成的是x86程序。

  • https://opencv.org/releases/ 下载Windows版本的库文件。下载后能得到一个exe文件,把该文件“解压”到一个路径即可。解压的文件夹里有一个opencv子文件夹。记此子文件夹为$(OPENCV)
  • x86的OpenCV库编译(可选)

    1. https://opencv.org/releases/ 下载Sources,这样编译出来的库可以满足当前编译器的配置。(直接安装的话没有x86的文件)
    2. 用cmake配置源代码。注意源代码的路径最好不要包含中文。(我第一次包含了中文,结果vs编译到一半出错了)
    3. 用vs打开cmake生成的目录下的OpenCV.sln,编译源代码。
    4. 编译结束后在目录的lib文件夹中找到编译的结果文件夹(Debug或Release),里面所有.lib文件就是编译的结果,把它们放到上一步的build文件夹里。我按照格式放到了$(OPENCV)\build\x86\vc15\lib文件夹里。
    5. 编译出的bin文件夹需要加入系统环境变量,因为里面的dll会被用到。改完环境变量重启电脑。
  • 在VS中的属性管理器新建配置,添加用户宏OPENCV,宏值就是开始那个解压出的opencv子文件夹。$(OPENCV)现在就表示库目录了。
  • 在属性的包含目录里加上$(OPENCV)\build\include,库目录改成对应版本的库的目录,比如我的就是$(OPENCV)\build\x86\vc15\lib
  • VS链接器-输入-附加依赖项中加入需要的.lib文件。我的第一个程序加入了opencv_highgui420d.lib;opencv_core420d.lib;opencv_imgcodecs420d.lib这三个库
  • 我随便到网上找了下函数,依葫芦画瓢地写了以下代码。如果一切都配置好了,程序会显示main.cpp目录下wallpaper.jpg这张图片。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>

const std::string windowName = "window1";

int main()
{
cv::Mat mat = cv::imread("wallpaper.jpg");

cv::namedWindow(windowName);
cv::imshow(windowName, mat);

cv::waitKey();

return 0;
}

OpenCV 像素获取与修改

OpenCV的操作很多都基于Mat的对象,也就是图像矩阵。

1
2
3
4
unsigned char R, G, B;
B = mat.at<cv::Vec3b>(i, j)[0];
G = mat.at<cv::Vec3b>(i, j)[1];
R = mat.at<cv::Vec3b>(i, j)[2];

对矩阵的at操作可以get或set颜色分量的值。

注意RGB是反的。

OpenCV输入输出图像

1
2
cv::Mat mat = cv::imread("in.png");
cv::imwrite("out.png", mat);

输入输出都是基于Mat的对象。

RGB转CMY、HSI

RGB转CMY

看了一下,RGB符合人对颜色的认知。但现实中,为了方便打印出黑色,用CMY表示颜料的颜色更加方便。

CMY就是RGB取反,即:

1
2
3
C = 255 - R;
M = 255 - G;
Y = 255 - B;

RGB转HSI

我对HSI的大致理解是:H是颜色的属性,比如是红色、绿色、蓝色还是其他颜色;S是颜色被稀释的属性,也就是颜色的深浅;I是光强,是颜色向量的模的大小。

https://blog.csdn.net/yangleo1987/article/details/53171623 里有好几个RGB转HSI方法的介绍。其中第一种是上课提到的方法。

风格转换

我问了一下别人,又看了一下书。感觉风格转换就是对某个颜色通道做一个函数映射,使得图片整体风格改变。

结构设计

由于程序过于简单,该程序用结构化的设计方法。

输入经过输入预处理模块,进入处理模块,最后进入输出模块输出。整个结构非常简明而无趣。

1
2
3
//         ******************    ***********************    *******************
// data -> ** input module ** -> ** processing module ** -> ** output module ** -> output
// ****************** *********************** *******************

由于要提交源代码文件,我把代码全部放到一个main.cpp里面了。但在我心中,还是为每个模块各建了一个头文件和cpp文件的。

程序设计

程序概览

代码放在Github上:代码仓库

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
int main()
{
cv::Mat mat = getInput("0.jpg");

auto matArr = processImage(mat);

outputImage(matArr[0], "1RGB.jpg");
outputImage(matArr[1], "2CMY.jpg");
outputImage(matArr[2], "3HSI.jpg");

cv::imshow("w1", matArr[0]);
cv::imshow("w2", matArr[1]);
cv::imshow("w3", matArr[2]);

cv::waitKey();

return 0;
}

main函数非常简明,getInputoutputImage隐藏了输入输出细节,方便修改。事实上,图像读入时像素是用unsigned char存储的,我在输入的时候把它转换成了float,输出的时候又转了回去。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
for (int i = 0; i < mat.rows; i++)
{
for (int j = 0; j < mat.cols; j++)
{
// tranformation
matCMY.at<cv::Vec3f>(i, j) = bgrToYmc(matCMY.at<cv::Vec3f>(i, j));
matHSI.at<cv::Vec3f>(i, j) = bgrToHsi(matHSI.at<cv::Vec3f>(i, j));

// tonal function
// ........


// inverse transformation
matCMY.at<cv::Vec3f>(i, j) = ymcToBgr(matCMY.at<cv::Vec3f>(i, j));
matHSI.at<cv::Vec3f>(i, j) = hsiToBgr(matHSI.at<cv::Vec3f>(i, j));
}
}

processImage的主体是遍历像素,先把原图像做格式转换,再用自己的风格转换函数瞎搞,最后转换回正常的格式。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
inline cv::Vec3f bgrToYmc(const cv::Vec3f& vec)
{
return cv::Vec3f(1.0f - vec[0], 1.0f - vec[1], 1.0f - vec[2]);
}
inline cv::Vec3f ymcToBgr(const cv::Vec3f& vec)
{
return cv::Vec3f(1.0f - vec[0], 1.0f - vec[1], 1.0f - vec[2]);
}
inline cv::Vec3f bgrToHsi(const cv::Vec3f& vec)
{
float R = vec[2], G = vec[1], B = vec[0];
float RG = (R - G), RB = (R - B);
float tmp = sqrt(RG * RG + RB * (G - B)); // 分母不能为0!!!
float theta = tmp != 0 ? acos(0.5f * (RG + RB) / tmp) : 0;
theta = B > G + EPS ? CV_2PI - theta : theta;
return cv::Vec3f(
theta / CV_2PI,
1.0f - 3.0f * std::min(std::min(R, G), B) / (R + G + B),
(R + G + B) / 3.0f
);
}
inline cv::Vec3f hsiToBgr(const cv::Vec3f& vec)
{
float H = vec[0] * CV_2PI, S = vec[1], I = vec[2];
float H2 = fmod(H, CV_2PI / 3);
float v1 = I * (1.0f - S);
float v2 = I * (1.0f + S * cos(H2) / cos(CV_PI / 3 - H2));
float v3 = 3 * I - v1 - v2;
if (H + EPS < CV_2PI / 3.0f)
{
return cv::Vec3f(v1, v3, v2);
}
else if (H + EPS < 2.0f * CV_2PI / 3.0f)
{
return cv::Vec3f(v3, v2, v1);
}
else
{
return cv::Vec3f(v2, v1, v3);
}
}

格式转换函数完全是在照搬公式而已,没有什么价值。

写程序过程中碰到的问题

  1. 为了公式处理方便,我把像素格式转换成了float。我一开始以为直接就能通过mat.at<Vec3f>来访问浮点格式像素,但是发现cv::Mat的机制好像是只能允许像素按一种格式存。要把图像矩阵用.convertTo(mat, CV_32FC3, 1 / 255.0)来转换成浮点表示,最后输出的时候还要用.convertTo(tmpMat, CV_8UC3, 255.0f)转换回去。
  2. 我还好提前看了一下cv::Mat的工作原理。该类只存了图像数据的指针,要复制图像数据的话要调用src.copyTo(dest)
  3. 图片转HSI再转回来后,我发现图像中出现一些绿色方块。经调试发现,H分量在运算的时候变成了nan。计算H的时候一定要判断分母是否为0!
  4. 图像的RGB是倒着存的。我之前测试OpenCV的时候碰到了这个问题,写这个程序的时候就跳过了这个坑。

处理结果

0.jpg

原图

1.jpg

RGB风格转换结果

2.jpg

CMY风格转换结果

3.jpg

HSI风格转换结果

感想

写博客真浪费时间。

数独求解软件

Github链接

Sudoku Solver

PSP耗时表

PSP2.1 Personal Software Process Stages Expected time (min) Actual time (min)
Planning 计划 20 50
· Estimate · 估计这个任务需要多少时间 20 50
Development 开发 980 1119
· Analysis · 需求分析(包括学习新技术) 50 120
· Design Spec · 生成设计文档 150 27
· Design Review · 设计复审(和同事审核设计文档) - -
· Coding Standard · 代码规范(为目前的开发指定合适的规范) 30 5
· Design · 具体设计 210 161
· Coding · 具体编码 300 575
· Code Review · 代码复审 60 25
· Test · 测试(自我测试,修改代码,提交修改) 180 206
Reporting 报告 170 100
· Test Report · 测试报告 80 60
· Size Measurement · 计算工作量 30 10
· Postmortem & Process Improvement Plan · 事后总结,并提出过程改进计划 60 30
合计 1170 1269

具体任务划分

  • 开发 980
    • 需求分析 50
      • 列下任务的具体项 20
      • 评估自己是否缺少相关知识,如果缺少则去收集资料 30
    • 生成设计文档 150
      • 用例模型、基本类图 30
      • 动态模型设计 120
    • 设计复审
    • 代码规范 30
      • 查现有代码规范 20
      • 定义代码规范 10
    • 具体设计 210
      • 精化类图 60
      • 精化动态模型 120
      • 构件图 30
    • 具体编码 300
    • 代码复审 60
    • 测试 180
      • 单元测试 120
      • 集成测试 60
  • 报告 170
    • 测试报告 80
      • 用例设计 20
      • 测试 20
      • 报告撰写 40
    • 计算工作量 30
    • 总结 60

解题思路

需求分析

阅读作业要求文档后,我对需求的总结如下:

任务主要分两个部分:数独生成和数独求解。

数独生成要求通过可处理异常输入的命令行参数传入生成终局的数量,按照格式要求生成指定数量的终局。格式要求是数字用空格隔开,行末无空格,文件末无换行,终局间空一行。终局第一个数字为学号后两位之和模9加1。

数独求解要求从命令行读取待求解数独文件名,从文件中读取所有终局并按和数独生成同一格式输出求解结果。

由于数据求解需要数据,我还需要加入一个生成待求解的数独。参考附加任务的需求,待求解数独满足挖空数在30~60之间,每个宫不少于2个。

解法思考

数独终局生成

由于数独的每一行都是不重复的,因此每一行可以看成一个9个数的排列。一行一行地生成数独或许是一个比较好的方法。

数独的第一行就可以区别出一类数独。如果能得到一个终局的话,只需要修改第一行的映射,就能快速得到大量终局。同时,由于数独的限制比较强,搜索算法不会浪费很多时间在不合法终局上。如果知道了前几行,可以对剩下的部分进行dfs暴力搜索。这一步和数独求解完全相同,将在下一节分析。现在的问题就是如何生成前几行。

假设现在我们得到了第一行的一个合法排列。对于给定的一行,考虑有多少新行,使得新行和给定的行的三部分都能放在同一个宫里。这是一个排列组合的问题,我用分类讨论的方法对它进行了求解。

对于给定的行,行内的每个数只需要考虑它们来自哪个宫。不妨把给定的行表示为[1 1 1 2 2 2 3 3 3],新行表示为[0 0 0 0 0 0 0 0 0]。0表示未填,1、2、3表示这个数来自旧行的哪一个宫。我们先考虑位置信息,最后再区分三个1、2、3。

分类讨论1 1 1被怎么放到了新行。根据这三个数有几个放到了第二个宫、第三个宫,可分成四种情况:

1
2
3
4
[0 0 0 0 0 0 1 1 1]
[0 0 0 1 0 0 1 1 0]
[0 0 0 1 1 0 1 0 0]
[0 0 0 1 1 1 0 0 0]

类似地,可以分类讨论,把2、3填入第二行,使得第二行合法。通过排列组合的计算,总共能得出$C_3^3C_3^3+C_3^2C_3^2C_3^1+C_3^1C_3^1C_3^1+C_3^0C_3^0=56$种不同的结果。

这样分类的好处是,对于一个给定的第一行,每种生成的序列的不重复的。改变第一行的映射就能改变终局。而第一行由于第一个数固定,有$8!=40320$种可能。因此,用这种方法至少可以生成$56\times40320=2257920‬$种终局,这个数满足$10^6$的需求。

因此终局生成的算法总结如下:

  1. 固定第一行第一个数,随机生成一个排列

  2. 按上述方法枚举出第二行的所有可能,把对应的1、2、3转换成任意一个合法的序列。

  3. 使用求解方法得到一个终局。

  4. 交换数字映射来得到额外的终局。

如果仅仅要满足需求的话,最多只需要求解56个数独。如果要更多的终局的话,只需要在求解的时候多求几组解就行了。

数独求解

我在玩数独的时候发现了两个性质:某个宫的某个数字只能填入一个位置;某一个位置只能被填入一个数。

因此搜索时,可以加入下列策略:对于每个宫的每个数建立一个表,表示它的可选位置;对于每个位置建立一个候选数字表。如果数字的位置唯一或位置的候选数唯一,就填上这个数,更新所有的表。如果没有唯一的数,就选一个候选数最少的位置,尝试填入一个数字。如果某个数的可选位置为空或者某个位置的候选数字为空,就回溯,直到数独填完。

待求解数独生成

待求解数独要从终局中生成。根据需求,只需先随机总空数,再每个宫随机挖掉两个空,再随机挖空直至满足空格数即可。

设计实现过程

面向对象分析与设计

用例图

用例图

根据需求分析能够快速得到用例图。

基本类图

classdiagram

基本的类图建立如上。

用户交互类负责处理输入输出,并且根据输入参数来调用执行该任务的类。

为了使算法能够复用,每项任务的算法单独建立一个类。

数独矩阵只存储了数独81个数的信息,该类被所有的任务类使用。在终局生成算法和求解算法中,会用到数独的其他信息,这些附加信息用一个新的数独状态类来存储。

classdiagram

后来修改了一下类图的设计,三个任务器不应该做为用户交互类的子对象,应该是一种关联关系。

classdiagram

后来写终局生成算法的时候发现终局生成器和终局生成算法的关联设计得不太对。

(以上两个任务共花费27分钟)

活动图

数独求解

(花费40分钟)

根据之前分析问题时的算法,得到以下活动图:

solve

终局生成

根据之前分析问题时的算法,得到以下活动图:

generate

(花费24分钟)

精化类图

感觉visio建模很不方便,没有提供数组的定义。以下类图表示我在建类时的大致想法,与实际的类定义应该存在不小的出入。

数独矩阵

用例图

数独状态

用例图

求解算法

solve

(以上三个任务共花费38分钟)

生成算法

Generate

(花费5分钟)

残局生成算法

GetEndGame

(花费5分钟)

终局生成器

Generator

求解器

Solver

残局生成器

EndGameGenerator

(共花费14分钟)

用户交互

用户交互部分需要对数独矩阵进行输入输出。这一部分应该再用一个类来完成。

用户的命令行命令可以分成两个部分:命令和参数。应该有检查命令和执行对应命令的函数也应该分成两大部分。

执行命令的具体细节都交给对应的命令器完成。

CmdInteraction

(花费35分钟)

程序设计

(已花费1小时)

(又花费了2小时)

(又花费了1小时40分钟)

(又花费了40分钟)

(又花费了30分钟)

(又花费了45分钟)

(又花费了60分钟)

(优化花了120分钟)

求解算法

根据之前的面向对象分析与设计,我实现了4个类:SudoMatrix, SudoChoice, SudoState, SudoSolveAlgorithm。和设计的时候相比,每个类有一些变化。

SudoMatrix多写了很多方便调用的静态函数,比如位置坐标和第几个宫第几个数坐标的互相转换。

SudoChoice所有的操作都是封装好的,比如加入一个选项,去掉一个选项。

由于大致的结构都已经预先设计好了,代码实现得比较顺利,Debug时间很少。

最终求解算法定义如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
class SudoSolveAlgorithm
{
public:
// return true if the sudoku has solution
static bool solve(SudoMatrix& mat);
private:
SudoSolveAlgorithm();
~SudoSolveAlgorithm();

static bool dfs();
static bool fillState(const SudoMatrix& mat);

static SudoState _sudoState;
};

算法只维护了一个变量:当前的数独状态。

算法主要函数solve实现如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
bool SudoSolveAlgorithm::solve(SudoMatrix& mat)
{
_sudoState = SudoState();
if (!fillState(mat))
{
return false;
}
if (!dfs())
{
return false;
}
mat = _sudoState.getMat();
return true;
}

和之前的活动图设计的一样,算法先把已有内容填充进数独矩阵中,之后进行dfs。如果有任意一个步骤中发现数独无解,就立刻返回false。如果数独成功求解出来,就把结果数独进行拷贝。

dfs只对数独当前状态这个变量进行操作。不断找可能最小的下一状态,如果碰到了无解的情况就让状态回溯。

终局生成算法

该类定义如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
class SudoGenerateAlgorithm
{
public:
// To get the sudoku, you should generate the template first.
static void calTemplate(int count = MAX_TEMPLATE_COUNT);
// Search valid diffrent initial sudoku passed to solve algorithm
static void dfs(SudoMatrix& crtMatrix, int pos, int cnt1, int cnt2, int cnt3);

// The permutation indicates the first row of the generated sudoku.
static SudoMatrix getSudokuFromPermutation(std::vector<int>& permutation, int templateId);


const static int MAX_TEMPLATE_COUNT;
private:
SudoGenerateAlgorithm();
~SudoGenerateAlgorithm();

static int _templateCount;
static std::vector<SudoMatrix> _templateSudokus;
};

和设计的时候相比,实现时又加了一个dfs函数,以搜索前两行的合法状态。

算法仅提供生成模板的函数和从模板和排列生成最终终局的函数。具体生成多个矩阵的工作在终局生成器中完成。

残局生成算法

该类定义如下:

1
2
3
4
5
6
7
8
9
10
11
12
class SudoEndGameAlgorithm
{
public:
static SudoMatrix getEndGame(
const SudoMatrix& sudokuFinalState,
int minEmptyEntryPerPalace,
int totalEmptyEntryCount,
int seed = 123456);
private:
SudoEndGameAlgorithm();
virtual ~SudoEndGameAlgorithm();
};

和设计的时候相比,算法的主要函数中多了一个参数,表示随机种子。

交互类与命令器

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
class CmdInteraction
{
public:
CmdInteraction();
virtual ~CmdInteraction();

bool inputCmd(int argc, char* argv[]);
private:
Solver _solver;
EndGameGenerator _endGameGenerator;
FinalStateGenerator _finalStateGenerator;

bool processCmd(char* cmdText, char* argPtr[]);
bool generateSudoku(char* arg[]);
bool solveSudoku(char* arg[]);
}

交互类在最终实现的时候还是把每种命令器当成了子对象。

1
2
3
4
5
6
7
8
9
10
11
12
class FinalStateGenerator
{
public:
FinalStateGenerator();
virtual ~FinalStateGenerator();

// Generate matrices of specific count and stored in pre-allocated array.
// The max count is 2e6 if the first number is fixed.
// Otherwise the max count is 1.8e7.
// Return the actual count of generated matrices.
int generateFinalState(int count, Sudo::SudoMatrix* matrixArr, int fixFirstNumber = 0) const;
};

终局生成器和设计的时候差不多,只有一个调用函数。做为输出的数组需要提前分配内存。根据需求,可以设置固定在第一行第一个的数字。

1
2
3
4
5
6
7
class Solver
{
public:
Solver();
virtual ~Solver();
bool solve(SudoMatrix& matrix) const;
};

求解器仅仅在求解算法上面加了一层调用,和设计时几乎一样。

1
2
3
4
5
6
7
8
9
10
11
12
13
class EndGameGenerator
{
public:
enum Difficulty
{
EASY,
MIDDLE,
HARD
};
EndGameGenerator();
virtual ~EndGameGenerator();
SudoMatrix generateEndGame(const SudoMatrix& finalState, Difficulty difficulty = MIDDLE, int seed = 123456) const;
};

残局生成器和设计时也差不多,只不过难度是用枚举表示的。对应生成算法,生成器也需要输入一个随机种子。

难度之间只有总空数和每个宫空数的不同。简单难度总空数20,每个宫至少2个空;中等难度总空数40,每个宫至少2个空;困难难度总空数60,每个宫至少3个空。

单元测试

由于VS2017 Community查看不了代码覆盖率,单元测试仅能测试代码的大致正确性。

终局生成器

由于生成器只有一个函数generateFinalState,我仅对这一个函数做了测试。

我测试了三个指标:生成的数独是否正确(满足行列宫不重)、第一个数是否正确、是否有重复的数独。

其中是否有重复的数独直接比较是O(n^2)的,测试起来十分慢。我把每个数独的每一行转换成一个数字,一整个数独就可以看成是一个9个数字的数组。我把每个数独转换成数独,扔到一个set里。如果在set发现了重复的数组,就说明有数独矩阵重复。这样测试重复的复杂度是O(nlogn)的。

FinalStateTest

如图,我总共运行了两次测试。两次测试只有生成数独的数量不同,一次是1000,一次是100000.

残局生成器

我用终局生成器生成了终局,之后把终局放入残局生成器了生成残局。

由于终局生成器已经验证过正确性,现在只对残局生成器验证两个指标:每个宫的最小空格数、所有空格数之和。

我对3种难度的残局生成器各设置了3个种子进行测试。生成数独的数量都是100个。

EndGameTest

求解器

(测试求解算法正确性花费30分钟)

Result1-1

我在main函数里添加了一些临时代码,以验证求解算法的正确性和效率。

我在http://www.sudoku.name/index-cn.php 这个网站上随便找了一个难度为困难的数独,在Release模式下运行代码。从结果中可以看出,算法高效且正确的完成了任务。

最后我把残局生成器单元测试中得到的9组数据放入求解器进行测试,仅考虑求解出来的数独是否完全合法。

SolverTest

(之后共花了100分钟)

集成测试

我直接运行最终程序,做集成测试。

交互模块的测试用例可分成以下等价类:

  • 命令个数错误

  • 命令错误(无‘-’)

  • 命令错误(有’-‘,第二个及后面的字符不对)
  • 命令参数错误

命令参数错误根据具体的命令可以细分:

  • -c 之后是字母
  • -c 之后数字不合法
  • -s 之后不是路径
  • -s 之后文件不存在

我先对以上七种情况设计了10个测试用例,运行结果如下:

ITest1

最后我输入了两个正确的命令。(sudokuInput.txt里有预先生成好的20个数独题目)

ITest2

在集成测试的时候,我发现我输出的最后一行有空行。经过修改后这个BUG被排除了。

(花费56分钟)

代码复审

我在VS中启用了代码分析,规则集为Mircrosoft本机建议规则。

一开始,有很多警告:

CodeReview1

经过修改,警告被改掉了。

比较特殊的是一个和迭代器有关的警告。迭代器的偏移量我提前计算好不会报错,如果直接把偏移量的计算式和迭代器放到一起运算就会报错。这个可能和64位平台有关的警告我没有太理解。

(花费25分钟)

性能分析及程序优化

数独生成

性能分析结果

使用vs自带的性能分析,生成1000000个数独,得到的结果如下:

Performance1

可以发现,主要的时间都花在输入上。

Performance2

其中,最花时间的是fputc函数。

优化方法及结果

我把输出的字符先存到一个缓冲区里,每存100个矩阵的信息就输出一次。

输出一次性用fwrite输出。

另外,由于生成算法需要的模板较少,我直接把最终的模板放进了文件里面,而不是每次都计算模板。

Performance3

从结果可以看出,运行效率有所提升。

数独求解

初次结果

最开始的时候,求解1000000个数独的速度非常慢。我没有等程序跑完就关掉了程序。

优化方法及结果

Performance4

我用性能分析工具检查,发现很多时间都浪费在了数据结构操作上。

我把存状态的set改成了bitset,把部分操作改成了位询问、修改。

同时,我稍微优化了状态初始化部分的算法。在初始化一个数独的时候,不会再写入修改日志了。

Performance5

现在程序能90秒跑完了。

单元测试

UTest1

经过一些修改后,单元测试又能通过了。

最终程序代码说明

求解算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
bool SudoSolveAlgorithm::solve(SudoMatrix& mat)
{
// Fill the matrix with known message
if (!fillState(mat))
{
return false;
}
// Use dfs to find the solution
if (!dfs())
{
return false;
}
// Get result
mat = _sudoState.getMat();
return true;
}

求解算法是整个软件的核心。在求解时,算法先更新当前已知的信息,之后对其他部分做DFS。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
bool SudoState::fill(const SudoMatrix& mat)
{
// It's the first(0) step.
_step = 0;

// Create a new matrix
_mat = SudoMatrix();

// For each entry and number, initialize their state
for (int i = 0; i < SudoMatrix::SUDO_SIDELENGTH; i++)
{
for (int j = 0; j < SudoMatrix::SUDO_SIDELENGTH; j++)
{
_entryChoises[i][j].init(SudoChoice::Type::ENTRY, i, j);
_numberChoises[i][j].init(SudoChoice::Type::NUMBER, i, j);
if (mat(i, j) == 0)
{
_entryChoises[i][j].setAllOne();
}
_numberChoises[i][j].setAllOne();
}
}

// For each known number, update the message through setNumber()
for (int i = 0; i < SudoMatrix::SUDO_SIDELENGTH; i++)
{
for (int j = 0; j < SudoMatrix::SUDO_SIDELENGTH; j++)
{
if (mat(i, j) != 0)
{
if (!setNumber(i, j, mat(i, j), false))
{
return false;
}
}
}
}
return true;
}

FillState会调用SudoState::fill函数。该函数先初始化整体的信息:当前已经运行到第几步,当前的数独矩阵。并且把每个位置的候选数和每个数字在每个宫的侯选位置初始化。最后,对于每个已经存在的数,用setNumber()填入并更新状态。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
bool SudoSolveAlgorithm::dfs()
{
if (_sudoState.getStep() >= SudoMatrix::SUDO_ELEMENTS_CNT)
{
return true;
}

const auto& op = _sudoState.getMinForkOperation();

if (op.getType() == SudoChoice::ENTRY)
{
// ENRTY
const auto& set = op.getChoices();
for (auto num : set)
{
bool yes = false;
if (_sudoState.setNumber(op.getPosIOrNumber(), op.getPosJOrPalace(), num + 1))
{
if (dfs())
yes = true;
}
if (!yes)
_sudoState.recall();
else
return true;
}
}
else
{
// NUMBER
.......// do the same thing
}

return false;
}

dfs函数是用于搜索解的函数。该函数每次从数独当前状态中获得候选情况最少的位置或数字,之后对这个地方进行搜索。搜索完全部解就会返回成功;setNumber返回false说明搜索失败,数独当前状态会进行回溯。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
bool SudoState::setNumber(int i, int j, int num, bool recordLog)
{
assert(_mat(i, j) == 0); // the entry must be empty
assert(num >= 1 && num <= 9); // number must in 1~9

int palaceId = SudoMatrix::getPalaceId(i, j);
int idInPalace = SudoMatrix::getIdInPalace(i, j);

bool valid = true;
_mat(i, j) = num;

// record currentLog
.....


// Ban all the choices in (i, j)
......

// Update the choices of entries with same column
......

// Update the choices of entries with same row
......

// Update the choices of entries with same palace
......

_step++;
return valid;
}

数独状态的setNumber函数是算法中的十分重要的函数。该函数在填入一个数字后,更新相关行、列、宫的候选数字信息及每个数的候选位置信息。同时,这一步的操作会被记录下来,以进行回溯。

生成算法

在模板是本地储存而不是实时计算后,生成算法只剩了从排列构造这一项任务。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
SudoMatrix SudoGenerateAlgorithm::getSudokuFromPermutation(std::vector<int>& permutation, int templateId)
{
assert(permutation.size() == 9);
assert(templateId >= 0 && templateId < _templateCount);

SudoMatrix result = _templateSudokus[templateId];
for (int i = 0; i < 9; i++)
{
for (int j = 0; j < 9; j++)
{
int newIndex = result(i, j) - 1;
result(i, j) = permutation[newIndex];
}
}
return result;
}

生成算法类提供了排列构造的函数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
int FinalStateGenerator::generateFinalState(int count, Sudo::SudoMatrix* matrixArr, int fixFirstNumber) const
{
SudoGenerateAlgorithm::calTemplate();

int generateCount = 0;
std::vector<int> permutation;
if (fixFirstNumber == 0)
{
permutation = { 1, 2, 3, 4, 5, 6, 7, 8, 9 };
}
else
{
permutation.push_back(fixFirstNumber);
for (int i = 1; i <= 8; i++)
{
if (i >= fixFirstNumber)
permutation.push_back(i + 1);
else
permutation.push_back(i);
}
}
do
{
if (fixFirstNumber != 0 && fixFirstNumber != permutation[0])
break;

for (int i = 0; i < SudoGenerateAlgorithm::MAX_TEMPLATE_COUNT && generateCount < count; i++, generateCount++)
{
matrixArr[generateCount] = SudoGenerateAlgorithm::getSudokuFromPermutation(permutation, i);
}
} while (std::next_permutation(permutation.begin(), permutation.end()) && generateCount < count);

return generateCount;
}

按照需求生成数独的工作交给生成器来完成。生成器生成许多固定第一个数字的排列,利用排列生成数独,直到生成的数独满足数量要求为止。

心得体会

从我对整个项目的评价和我学到的东西两个方面来谈我的体会。

我对项目的评价可以分成我的开发时间管理评价和项目结果评价。我这次时间管理得非常不好。我一开始觉得这个项目难度较低,花费不了什么时间。我很早就想好整个项目的框架和具体算法,设计和编码的过程也还算顺利。但我错估了测试和项目优化的时间,项目后期的时间非常紧,我也没来得及完成附加题。

项目完成得也不够好。我在设计数独状态类的时候没有设计好,返回给数独求解类的信息和数独状态类自己保存的信息应该分成两个类。我在后来代码优化的时候由于这两个类没有分开,碰到了一些麻烦。项目的性能还存在优化的空间,求解算法可以进一步压缩时间。

在这个项目中,我学习到了一些工具的使用方法、软件工程知识、软件工程实践经历。我学会了VS单元测试工具、性能分析工具、代码分析工具的使用。我在画用例图、类图等UML图的过程中对于软件工程的理解进一步加深。在本次项目中,我不局限于编码,而是完成了整个软件的开发流程,体会到了软件开发的不易。

我本身就有比较多的大型程序开发经验。在学过了软件工程课并经历了这次的开发过程后,我能把知识和我的经验结合起来,以后写出模块划分得更好、测试做得更全面、文档更完整的软件。

更新记录

  • 19.12.24 11:36

    1. 添加博客地址。
    2. 完成PSP耗时表。
  • 19.12.27 12:03

    1. 分析题意,想出算法。
  • 19.12.28 23:03

    1. 设计用例图。
    2. 设计基本类图。
  • 20.1.1 21:41

    • 项目进展

      1. 画完求解活动图。
      2. 精化求解模块的类图。
      3. 求解模块的程序设计中。
    • 项目更新

      1. 修改基本类图。
    • 博客进展

      1. 博客跟进项目进展。
      2. 添加未来任务的大标题。
  • 20.1.2 16:30

    • 项目进展
      1. 求解模块的程序设计中。
    • 博客进展
      1. 记录了花费在程序上的时间。
  • 20.1.10 15:48

    • 项目进展
      1. 完成求解算法。
      2. 求解算法基本测试完毕。
    • 博客进展
      1. 更新了求解算法的代码简要说明。
      2. 更新了程序设计花费时间。
      3. 更新了求解算法的正确性测试结果。
  • 20.1.15 23:30

    • 博客进展
      1. 添加了终局生成的活动图。
  • 2020.1.16 20.39

    • 项目进展
      1. 设计并完成了终局生成算法模块。
      2. 给SudoMatrix类增加比较函数与复制函数。
    • 博客进展
      1. 再次修改了基本类图。
    1. 添加了终局生成算法类图。
  • 2020.1.17 10:35

    • 项目进展
      1. 设计并完成残局生成算法模块。
    • 博客进展
      1. 添加了残局生成算法类图。
  • 2020.1.17 15:10

    • 项目进展
      1. 设计并完成终局生成器、残局生成器、求解器、用户交互、数独矩阵IO类,但没有测试。
    • 博客进展
      1. 添加项目进展中的设计类图。
  • 2020.1.18 21:26

    • 项目进展
      1. 添加了单元测试的所有代码。
    • 博客进展
      1. 统计了设计文档、具体设计、编码的时间。
      2. 更新了单元测试有关内容。
  • 2020.1.18 22:33

    • 项目进展
      1. 完成集成测试的BUG修改。
    • 博客进展
      1. 更新集成测试有关内容。
  • 2020.1.19 10:52

    • 项目进展
      1. 完成代码复审
      2. 进行优化的分析
      3. 完成优化及优化后的测试
      4. 发布最终版本
    • 博客进展
      1. 更新代码复审内容
      2. 更新优化内容
      3. 更新最终代码说明
  • 2020.1.19 21:14

    • 博客进展
      1. 完成了心得体会。

今天——准确来说应该是发这篇文章的好几天前(发博客的时候竟然已经过去半个月了)——我过得十分愉快,而且我做的事情和平常也大不相同,十分有趣。我打算把今天发生的事情稍微记录一下。

昨天晚上,我躺在床上辗转反侧,过了两点才勉强睡着。于是,早上我理所应当地在九点半醒来,并且赖床赖到了十点多。如果是平常的我,肯定会想:“哎呀,半天都快过去了,下午肯定干不了什么事了。今天就过得随便一点吧。”但是,这周末我有三门作业要写,无数项目要开坑,理智与责任阻止了我怠惰的想法。我充满自信地对自己说道:“不慌。虽然一个上午的时间没了,但是问题不大。今天我来展示一下大局观选手是怎么样学习的。”

首先,我进行了自身能量管理。考虑到没吃早饭,我塞了几块威化饼干,稍微补充了一下糖分,并且立刻点了外卖,准备11:20的时候吃午饭。之后,我开始规划今天的任务。我计划中午写完一门作业,下午在电脑上装个东西,去学校完成项目,晚上再写一门作业。一切准备就绪,我正式进入了工作模式。

在我正在写第一门学科的作业时,外卖到了。考虑到如果我立刻进食,肯定会血糖升高,头脑无法转起来,我毅然决然地不顾自己健康,选择边吃饭边看书。我每写一题,就吃一两口。饮食的诱惑令我写起题目来动力十足。最终,下午一点的时候,我圆满地写完了第一门学科的作业,并且完全理解了整个章节,顺便吃完了饭。

接下来,我准备休息一下,让电脑继续编译之前编译了一半的linux内核。很不幸,经过了长时间的编译后,我发现这个linux版本太高,虚拟机加载不出来。我又去打了一把炉石酒馆乱斗,倒数第二轮没有升级,导致伤害少了一点,给别人留了1点血,最后一轮被1血翻盘,以第二的成绩含恨结束游戏。我带着怨气和一肚子的葡萄糖趴在桌上睡了一觉。

下午四点,我准备去学校编程。为什么我非得去学校用电脑呢?这里要解释一下事情的前因后果了。上周末,我不小心再次把水弄进笔记本电脑里了。虽然我已经多次见到这种事了,但我还是抱有侥幸心理地按了一次电脑电源。电脑拿去维修后,其他地方勉强能用,“只是”充不进电。电脑最后被拿去返厂维修了。这一周,我只能拿着陪伴我多年,配置却已经被时代抛弃的老电脑进行工作了。

老电脑慢就慢点吧,我那颗波澜不惊的大心脏对电脑运行速度也没有太多的要求。可是,这个笔记本的显卡竟然不支持OpenGL 3.3。OpenGL 3.3开始支持核模式,我的OpenGL程序都是在核模式下写的。我现在有个OpenGL的项目急着要做,但现在我的代码在我自己电脑上能编译,却无法运行,编程工作开展不了。我想了好多办法,最终决定每次在这个电脑上编译,然后把可执行文件放到网盘上,用学校机房的电脑来运行程序。这样麻烦是麻烦了一点,但起码我还可以正常地工作。周五我已经证明了这种方法的可行性。

但是,今天下午,当我跑到机房后,却发现机房的门全关了。机房门口的“自习室”三个字,完全融入了图书馆的背景。这下我无处可去了。

我再次动用我智商300的大脑思考问题的解决方案。我排除了找同学借电脑等65535个方案,决定使用第65536个方案——去网吧编程。这周我实在被我老电脑的性能气坏了,不用一用好电脑就不解气。我很快开始计划起了晚上的行动。

现在四点多,离吃饭还早。我决定先回去尝试写掉一门作业,之后自己煮个饺子吃,吃完再去网吧。可惜,这门学科当前章节的内容过于复杂,我还没有看完这章内容饺子就煮好了。吃自己做的东西的体验还是挺棒的,只是不知道从新奇到厌烦的过程究竟有多快。不管怎样,吃完晚饭后的我精神饱满,斗志昂扬,接下来的冒险深深吸引着我。

只看了一眼地图后,我凭借多年走迷宫经验,顺利来到了一家比较好的网吧,噢,现在的网吧要叫作网咖。我本身都只玩一些配置不高的游戏,对硬件条件不怎么需求,平时根本不会想去网吧。但我也偶尔会去,都是和同学一起开黑。而今天来到网吧,我凭借着以前的经验,一切流程也算是轻车熟路。不过现在的网吧都很有套路,最便宜的上网区都要每小时10元,而办了会员是5~6元(节假日6元)。很明显,这种折扣就是逼着你要去办会员。我向来是个意志坚定的人,不会因贪图便宜而在大局上失亏本。但听店员说办会员只要50块钱,而且网吧还是全国连锁时(以后还能用),我就立刻真香地办了个会员。本身今天晚上就要待个4小时左右,不办会员显然亏了。不贪白不贪啊。

开了机子,身躺舒适的沙发,面对宽大的屏幕,手敲噼里啪啦的机械键盘,舒适感充斥了我的全身。作为非富人(我比较喜欢程序员描述事情的方法),我们一定要培养花钱买服务的思想。人生苦短,及时行乐。在网吧花钱,不仅是为了上网本身,更多的是享受这种舒适的上网环境。以前网吧还全是抽烟的人,网吧的空气极其恶劣,未成年人可以提前享受成年人的糟糕气味。现在的网吧这点倒好多了,大家抽烟只会躲到厕所抽,游戏里潜行规避敌人伤害,现实里躲藏逃离法律制裁(北京室内是禁烟的)。

花了约半小时下好VS、从老电脑上迁移文件后,我正式开始工作了。非常幸运,网吧的特殊性没有影响到我:本来VS装好是要重启电脑来修改电脑配置的,不然有些库文件可能找不到,而在网吧重启电脑会让下载的东西清零。我手动给项目多添加了一个C盘上的包含目录,解决了这一问题。不过VS用起来还是有一些问题,程序在调试的时候看似停在了断点上,实际上程序只是停下来了,我只好用祖传的std::cout输出调试。经历千辛万苦的调试,当看到程序顺利运行的那一刻,我感动得快要流出热泪:我发现我的6块钱没了。

网吧说吵也吵,但这种环境下我反而能集中精神。所谓大隐隐于市,就是指我这样心如止水的高人。听着大家一句一句的粗口,听着一些我很熟悉的游戏名词,我的心理十分踏实。大概是人的优越感在作怪吧。大家都在玩游戏,就我在“学习”;我也不是只会学习,你们玩的游戏我也会玩,我也不死板。在图书馆,别人看书,你也看书。你恨不得眼睛不离开书本,好证明自己比周围的人更认真,最后反而忘了自己在干嘛;而在网吧,你不管写代码认不认真,哪怕是空着屏幕不玩游戏,你就已经与众不同,超过周围所有人了,心情反而放松。而且,图书馆里过于安静,反而放不开手脚;在网吧,程序出了BUG,调试了10分钟还调不出来,我可以一边自言自语地念着代码,一边入乡随俗地骂着“fuck”。烦起来敲一敲桌子,代码成功运行起来后大笑几声,这些举动在网吧都是那么地正常。对我来说,网吧绝对是一个比安静的图书馆,或者只有我一个人的自习室更好的地方。

三小时过后,我顺利完成了这一阶段的工作。这段时间里,我基本没有摸鱼,注意力高度集中,这在平时是不太可能的。我再一次见识到了自己工作能力的上限。我十分希望我能每天都保证这样的状态。但每次的情况都在表明,人的自我干预能力,远小于外界条件的干预能力。有时间的话我想写一篇有关这个理论的文章,但我也可能不写,因为我没有足够的例子,要说的话并不是很多。

工作过后,自然是娱乐。连我自己都被我劳逸结合,堪称完美的工作安排折服了。经理论分析,用脑过后不宜玩智力游戏,因此我选择去玩《CupHead》。打了几盘,我突然觉得来网吧不玩平时玩不到的游戏就亏了,于是想去玩《城市:天际线》。结果这破游戏教程都没有,严重违反周弈帆游戏第八定理之游戏必须设法让玩家上手。网上说,要学会玩这游戏,去B站上看个教学视频就好了。看视频哪里都能看,我来网吧看什么视频呢?

此时,十点半到了,我在把网吧待了4个小时。看着我花掉的网费,我深深体会到了时间就是金钱这一道理,于是速速结账下机。这个网吧的消费制度还真是不错,送新用户3个10块红包,每次上网只能用一个。我一边体会着白嫖十块钱的愉快感,一边试图拉出已经落入消费陷阱的自己。

网吧离我住的地方挺远,网吧门口也没有自行车了。正好,一天的辛苦劳动后,我需要进行一些真正的放松活动了。我开始了我的巡游演奏会。在晚上街上没什么人,口哨声异常清楚的时候,我一路高歌的爱好总算能得到满足了。在欢快的乐曲声中,我这美妙的一天也步入了尾声。

今天过得确实很有意思。有懒有勤,有怒有笑,有作有玩。一成不变的人生被一些有趣的事情所扰乱,总会给人带来新奇的体会。正如伟大思想家周弈帆在他本来应该写却没写的游记中的话:旅行不是为了去一个地方游览名胜、品味民俗,只是为了打破现有的生活规律,享受新的环境而已。人生创造快乐的方式很多,我得好好记住这种打破习惯以产生快乐的方式。


update:这篇文章估计是我快睡着的时候写的,语句不通顺、逻辑不清楚、缺乏上下文,还没有注意保护个人隐私信息。天知道我写的时候在想什么。新版本我修正了文章中的逻辑错误和连贯性问题,文章起码读起来通顺了。

Preface

本题是2018CCPC吉林站的H题。吉林站是我参加的第一场正式比赛,当时我们写的最后一题就是这道题。我在赛场上想出了正确的大致解法,却因为细节没有想全,没有在赛场中A掉这道题。也是从这次比赛开始,我所在的队伍在拿银的时候,几乎每次离金牌就差一题,而且这一题的算法已经被我们想到,只是没有写完。每场比赛打完,我的感受大多是“要是再写出一题就好了”。

时光飞逝,如今我都快退役了。我还有两场正式比赛要打,而其中一场比赛将在这个周末举办。我想在比赛前临阵磨枪,找一些代码量比较大的题来提升写代码的能力。我想起了这道当时认为比较繁琐的题目。

没想到WA了一次后,这道题就被我十分轻松的A掉了。思路清晰,大脑能量充足的情况下,这道题对我来说竟然算不上一道难题,甚至谈不上一道有代码量的题。我都不知道该怎样评价自己的水平了。该写出题的时候写不出题,训练的时候又看得过去。

不管怎样,我要给这道对我来说意义重大的题写一篇题解。

Description

数轴上有$n$个字符串,每个字符串初始为空。现在要进行$m$次以下操作:

  • 操作1:把一个$[l, r]$区间里所有字符串前面加上一位十进制数字,后面加上一位十进制数字。比如”44”这个字符串“加上”5后,就变成了”5445”.
  • 操作2:求$[l, r]$区间里所有字符串代表的数字之和。

样例组数$1\leq t\leq5$,长度$1 \leq n \leq 10^5$,操作次数$1 \leq m \leq 10^5$

Sample Input & Output

1
2
3
4
5
6
7
8
9
2
3 2
wrap 1 3 1
query 1 2
4 4
wrap 1 3 0
wrap 2 4 3
query 1 4
query 2 3
1
2
3
4
5
Case 1:
22
Case 2:
6039
6006

Solution

看到关于区间的操作、求和,第一反应肯定是线段树。

问题的关键在于,在一个数字前面和最后加了一位,这个操作的实质是什么。

假设当前某个数字是$x$,有$d$位,在它的前后放了一个$y$,那么这步操作的结果是$x \leftarrow y * 10^{d+2} + x \times10 + y$。可以发现,后面的乘法和加法都十分好维护。用两个lazy变量来标记区间该乘和该加的值。每次乘的时候还要改变加的标记。问题的关键是,前面加的那一堆东西和当前数无关,而是和当前数的位数有关。

既然操作和数位数有关,那么维护每个数的位数一定是一个可行的做法。单由于操作是对区间进行操作,一个区间里不同的数有不同的位数,这个东西哪怕可以用vector存下,合并的时间复杂度都会爆炸。有没有一种办法能够合并区间所有数的位数这个信息呢?

观察操作的式子我们可以发现,我们做的所有运算都是加法和乘法。我们不关心区间里每一个数有几位,而是要让这些数的最高位乘上一个操作数$y$再乘100,最后再加到一起。这样想的话,我们可以把每个数的最高位加起来,以获取这个区间所有数的最高位信息。比如一个区间的数是$11,2222,4444$,那么最高数位之和就是$10 + 1000 + 1000 = 2010$。那么对这个区间信息操作时,不再是$x \leftarrow y 10^{d+2} + x \times10 + y$,而是$x \leftarrow y 2010^{d+2} + x \times10 + y$。如果我们能维护一个最高数位和信息,我们就能计算出当前区间的数字和,也就是操作2询问的内容了。

现在问题就变成了如何维护区间的最高数位和及相关lazy变量。为了计算出新的$x$,我们要用$y$来乘最高数位和。由于区间操作要在线段树上下放,而父区间和子区间$y$乘最高数位和的结果是不同的。因此,我们只能每次下放$y$这个操作数。我们需要一个新的lazy变量来储存我们下放的$y$。

光存了$y$不够,因为子区间可能有没有计算过的操作数。比如子区间进行了操作数3的操作,子区间的数字应当是$3x3$,但子区间没有被访问,标记还没有计算过。现在又对整个区间进行了操作数4的操作,子区间的数字现在应当是$43x34$。此时子区间依然没有被访问。也就是说,之前lazy变量是3,经过操作后lazy变量变成了43。

每次进行操作,lazy变量前面都会加一个数。往前面加数这个操作可不好做。但是,我们可以记录这个lazy变量的最大位数。比如lazy变量现在是43,那么最大位数就是2。下次进行操作数5的操作时,$5 \times 10^2 + 43 = 543$,就能正确计算出新的lazy变量了。这个lazy变量的最大位数也是一个新的lazy变量。

所以,所有的变量我们都可以用线段树来维护,这道题就写完了。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
struct node
{
node *l, *r;
int preSum;
int preMult;
int preMove;
int sum;
int markSum;
int markMult;
int len;
node()
{
l = r = nullptr;
preSum = 1;
len = 0;
preMult = 0;
markMult = 1;
preMove = markSum = sum = 0;
}
void update()
{
sum = ((preMult * (ll)preSum) % MOD * p10[preMove] % MOD + (ll)sum * markMult % MOD + markSum * (ll)len % MOD) % MOD;
preSum = preSum * p10[preMove * 2] % MOD;
}
void pushDownTo(node* l)
{
l->preMult = (((ll)preMult * p10[l->preMove] % MOD) + l->preMult) % MOD;
l->preMove += preMove;
l->markSum = ((ll)l->markSum * markMult % MOD + markSum) % MOD;
l->markMult = (ll)l->markMult * markMult % MOD;
}
void pushDown()
{
update();
if (l)
{
pushDownTo(l);
pushDownTo(r);
}
clearMark();
}
void pushUp()
{
if (l)
{
l->pushDown();
r->pushDown();
preSum = (l->preSum + r->preSum) % MOD;
sum = (l->sum + r->sum) % MOD;
}
}
void setMark(int x)
{
preMult = x;
preMove = 1;
markMult = 10;
markSum = x;
}
void clearMark()
{
preMult = preMove = markSum = 0;
markMult = 1;
}
};

/*………………a lot of code…………*/

void modify(int l, int r, int v, node* nd = root, int cl = 1, int cr = n)
{
nd->pushDown();
if (l == cl && r == cr)
{
nd->setMark(v);
nd->pushDown();
}
else
{
int mid = (cl + cr) / 2;
if (r <= mid)
{
modify(l, r, v, nd->l, cl, mid);
}
else if (l > mid)
{
modify(l, r, v, nd->r, mid + 1, cr);
}
else
{
modify(l, mid, v, nd->l, cl, mid);
modify(mid + 1, r, v, nd->r, mid + 1, cr);
}
nd->pushUp();
}
}

ll query(int l, int r, node* nd = root, int cl = 1, int cr = n)
{
nd->pushDown();
if (l == cl && r == cr)
{
return nd->sum;
}
else
{
int mid = (cl + cr) / 2;
if (r <= mid)
{
return query(l, r, nd->l, cl, mid);
}
else if (l > mid)
{
return query(l, r, nd->r, mid + 1, cr);
}
else
{
return (query(l, mid, nd->l, cl, mid) +
query(mid + 1, r, nd->r, mid + 1, cr)) % MOD;
}
}
}

我一般是不放代码的,但我觉得这段代码写的非常有条理,所以就打算展示一下。

线段树操作的精髓在于node的操作。不管是怎样的区间操作、查询,所有modifyquery都可以调同一个格式的函数,只要修改node中函数的实现细节就可以了。

len变量不是一个动态的变量,我是在建树的时候把len算了一下。

markMultpreMove似乎是相关的两个变量,但多用一个变量也不会MLE会TLE,用两个变量能够让写代码时思路更加清晰。

Review

写多了线段树的区间操作后,你就能对线段树有更深的认识。线段树本质是维护一些满足幺半群的运算,也就是一些不可逆,互相叠加的运算,比如加法、乘法、求最值。

在对区间进行操作的时候,要想清楚两件事。一个是lazy标记怎么对固有变量产生影响的,一个是操作怎么对lazy标记产生影响的。如果只是简单的加法和乘法,这两个步骤实在是过于显然,以至于问题一变复杂,我们就搞不清应该如何维护lazy变量了。只要认识了区间操作的本质,任何单纯线段树区间操作都是很简单的题。

前几天在写HDU6619的时候,问题最终转化成了一个形如下式的DP:

这个式子很简单,也就是当前状态的值由前面某个状态的值加上两个状态之间的数组值的和得到,且当前状态的值要尽可能大。按照式子直接枚举当前状态的上一状态的话,复杂度是$O(n^2)$。但本题的$n$是$10^5$,直接写肯定会T掉。我一时半会没有想出优化的方法,看了题解才想起来,这个DP可以用斜率优化来加速转移。于是,我准备快速复习一下斜率优化的用法。

斜率优化DP用法示例

还是回到刚刚那个式子:

右边一堆求和实在是太碍眼了。我们可以先对数组求前缀和,这样求和直接被前缀和做差代替:

如果没有那个$(n - i)$就好了。这样上一状态的贡献永远是$dp[j] + sum[j]$,只要对每个状态求这个值再取最大就可以直接$O(1)$获取之前答案的最大贡献了。问题是,前面乘了一个$(n - i)$后,之前状态的贡献就和当前的$i$有关系了。之前所有状态的部分贡献$dp[j],sum[j]$已经是完全确定的东西了,加上当前的$i$就能算出之前状态的全部贡献。那么,有没有一个仅通过当前的$i$来衡量之前哪个方案更优的方法呢?

如果之前存在两个状态$a,b$,不妨设$a < b$。那么,当从$a$比从$b$转移更优的时候是怎么样的呢?我们不妨列出此时两个状态的贡献值不等式,并且把含$i$的部分移到不等号的一边去。

由于前缀和数组递增,$sum[a] > sum[b]$,不等式不用变号。

也就是说,通过计算两个状态的$dp$数组差除以$sum$数组差,并将这个结果与$n - i$比较,我们就能知道之前两个状态孰优孰劣。这个计算出的值是反映状态好坏的本质属性,因为这个值的计算与当前的$i$无关。由于这个计算值和数学上斜率的计算十分相似,因此才有斜率优化这个名字。

如果把之前所有状态都看成一个$(sum[j], dp[j])$的二维点,现在我们知道了$n - i$的值,那么之前哪个状态最优呢?

我们从左到右从两两相邻的点连线。如果两两相邻的点的斜率是递减的,那就好办了:这就说明,前面有很多点对的斜率大于$n - i$,后面有很多点对的斜率小于$n - i$。那么中间的某个点就是最优的,那个点和它前一个点的斜率大于$n - i$,和它下一个点的斜率小于$n - i$。它比左边的点优,也比右边的点优。

问题是,如果出现了某三个从左到右相邻点$a, b, c$,且$ab$斜率小于$bc$呢?也就是相邻点对的斜率不是增减的。这种情况下,若斜率$ab$和$bc$都比$n - i$小,则最左边的$a$点最优;反之,若两个斜率都比$n - i$大,则最右边的点$c$最优。剩下的情况只有斜率$ab$小于$n - i$,斜率$bc$大于$n - i$这种情况。但是,此时$a$比$b$优,$c$也比$b$优。可以看到,无论什么情况下,$b$都不是最优的。因此,$b$点永远不会是上一个最优状态,可以把它“抛走”,再也不考虑它了。

也就是说,可能作为上一个最优状态的点,它们的$(sum[j], dp[j])$构成的二维点的斜率永远是单调递减的。这个单调性,是斜率优化DP的根本原因。

有了上述的分析,选上一个最优状态,不再需要枚举每一个之前的状态,而是维护一个可能作为上一状态的集合就可以了。这个集合从左到右存了一系列的可能的状态。每次要找最优状态的话,就像之前分析的那样,找到第一个斜率小于的$n - i$的左边的点,那个点对应状态就是最优的上一状态。由于单调性的存在,这个过程可以用二分完成。但由于我们每次的$n - i$是一个递减的量,也就是说我们每次需要找的斜率只会越来越小。那么那些斜率很大的点都可以扔掉去。因此,在本题中,可以不去二分查找,而去删掉前面那些已经斜率过大的点。最后剩下来的第一个点,它和下一个点的斜率小于$n - i$,它就是本次要寻找的最优上一状态。这个集合可以用双端队列Deque来维护。

同时,这个集合还需要进行维护。每次新得到了一个dp值后,我们就要试图把新的状态加入这个集合中。为了保证集合中点斜率的单调性,我们要取出集合最后两个点,拿它们和当前的点算出两个斜率,看看这个斜率是否满足单调递减的性质。如果不满足,就把中间的那个点去掉,也就是把原来集合中的最后一个点去掉。

Deque实现斜率优化DP方法简述

其实上面那道题并不适合用来学习斜率DP,因为它实际上牵扯的细节很多,dp式子并没有上面写得那么简单。HDU3507应该是最适合写的第一道斜率优化DP题。

这道题状态转移的式子比上面那个例子要复杂一点,而且是取最小值,需要维护的斜率应该满足单调递增而不是单调递减,但本质还是对前缀和什么的做一些运算,再进行斜率的比较,最终优化的方法还是完全一样的。

这道题中,不等式右边的与当前状态$i$有关的值是$sum[i]$,它是一个单调递增的量,正好与我们维护的是一个单调递增的斜率集合“相符合”(这里“符合”的意义比较微妙,也就是说,如果这里$sum[i]$不是单调递减,是递增或是无规律的,那么就不能用Deque来做了)。因此,我们在找最优状态的时候,不用去二分的找第一个斜率大于某个值的点,而是把前面斜率小的点都删掉去,把剩下来的第一个点作为上一状态。这个东西可以用Deque维护。

因此,在写代码时,应该按照下列步骤。

  1. 读入数据,初始化dp等数组,该求前缀和求前缀和。
  2. 声明一个deque\< int >,作为上一状态的集合。如果0位置的状态一开始就能被转移,就往deque里push一个0。
  3. 开始for循环,遍历每一个当前状态$i$。
  4. for循环体中:如果deque的大小大于等于2,说明deque中可能存在斜率不合法的点,循环寻找是否存在这样的点。每次循环得到front(),再pop_front(),再得到front(),也就是取出队列中前面两个点。判断这这两点的斜率,若斜率小于(大于)关键值(本例是2 * sum[i]),则continue掉,也就是扔掉第一个点,继续循环判断;否则,把取出的第一个点push_front()放回去,退出循环。这样,第一个斜率大于(小于)关键值的点就变成了que.front()。
  5. 不妨设j = que.front(),这个j就是最优的上一状态。用推好的dp式子计算dp[i]。
  6. 把状态i加入队列并维护队列的单调性。和开始的操作类似,循环判断deque大小是否大于等于2,若是,则取出队列队尾的两个点,计算它们和新的点$i$的两个斜率是否满足应该具有的单调性。若不满足单调性,则扔掉队尾的点,循环继续。
  7. 把状态i用push_back()扔进队列的队尾。
  8. for循环结束
  9. 输出dp数组中的答案。

二分查找实现斜率优化DP方法简述

CodeForces 674C也是一道斜率优化DP的题目。这道题递推式子相较之前两题更复杂了一些,而且它进行比较的和$i$有关的关键值并不是单调的。因此,这题不能再用Deque来维护最优上一答案,而是要用一个vector来存可能的状态,并进行二分查找。

写代码的步骤也类似,只是寻找上一最优答案的过程有所变更。

  1. 读入数据,初始化dp等数组,求好所需变量。
  2. 声明一些vector\< int >,作为上一状态的集合。如果dp某一维的0位置的状态一开始就能被转移,就往对应vector里push一个0。
  3. 开始for循环,遍历每一个当前状态$i$。如果dp数组和本题一样是多维的,那么再加一层for循环,按照顺序遍历dp数组的每一维。
  4. for循环体中:如果对应vector数组是空的,那么说明该状态不存在上一状态,直接跳过;否则,声明一个存储最优上一个状态的变量ii。ii的初始值是vector数组的第一项,因为我们将要二分查找的是最后一个斜率满足不等式的右端点,若二分查找什么都没找到,那么说明最优的点就是数组中的第一项。
  5. 二分查找最后一个斜率满足不等式的右端点。实现时把二分查找的mid对应状态套入不等式判断即可,若满足不等式则更新临时变量ii为vector中mid处对应的状态。
  6. 二分查找结束,ii就是最优的上一状态。用推好的dp式子计算dp[i]。
  7. 把状态i加入数组并维护数组的单调性。循环判断deque大小是否大于等于2,若是,则取出数组尾的两个点,计算它们和新的点$i$的两个斜率是否满足应该具有的单调性。若不满足单调性,则扔掉数组尾的点,循环继续。
  8. 把状态i用push_back()扔进数组尾。
  9. for循环结束
  10. 输出dp数组中的答案。

写代码的注意事项

写代码的时候会发现,不等式左边分子的那个量会被反复地用到。用一个函数来计算这个值会更加方便可靠。

在比较相邻两个斜率,而分子分母都是整数的时候,得用交叉相乘来比较斜率,规避除法。

总结/感想

算法原理

我算是大概理解了斜率优化DP的原理,但对什么时候可以用这个算法还有一点模糊。斜率优化DP本质上是在寻找最优上一状态时,找到衡量每一个状态的本质量,也就是那个二维点。对于每一对二维点,都有同一个判断它们两个谁更优的标准,那就是比较二者斜率是大于还是小于某个值(等于的情况可以忽略)。由于判断条件相同,只有斜率单调的点才可能成为最优上一状态。又由于有了单调的性质,才能在状态中进行二分查找。如果进行判断的斜率值还相对应是单调的,那么二分查找都可以省略,直接就可以得到最优的上一个状态。

能用斜率DP优化问题的核心原因就在于,对每一对二维点,判断它们谁更优的标准是一样的。但是,如果对于一个dp递推式子推导出的不等式中,作为分母的不是$(sum[a] - sum[b])$,而是另一组值,这一组值并不满足单调。这样的话,哪怕对于$a > b$的点,算出来的分母可能是个负数,除过去的话不等式要变号。对于相邻的两个点,有的时候斜率大于某个值更优,有的时候是斜率小于某个值更优。好像在这种情况下就不能用斜率优化DP了。

算法适用情况

能够用斜率优化DP的题都有类似的地方。这些题的dp在计算新值时,是在旧dp值的基础上加上一个和当前状态$i$有关的贡献。这个贡献往往是和上一位置到这一位置的区间和有关,以保证不等式分母是永远是正的。当看到一个类似这样的DP,且数据范围不允许$O(n^2)$时,就要考虑用斜率优化DP了。

和大部分DP题目一样,这种题目推式子是最重要的,代码实现上只要写个2、3题就能学会套路,很难出错。

我之所以开通博客,是因为我感觉自己有很多想法,想把它们整理、发表出来。但在我开通了博客之后,我写的文章却非常少。我仔细分析了一下自己的心理,发现我并不是没有东西可写;相反,我经常会想到一些想写的话题。但是,当我准备去构思文章该怎么写的时候,我总不自觉地会开始想一些问题:“我要写的东西会不会太短内容不够?”“我写的东西会不会内容涉及得不够全面?”“我不去查阅资料就来写会不会失去严谨性甚至是正确性?”想的问题多了起来后,我就会发现自己还没有做好写这篇文章的准备。于是乎,一篇文章就这么搁置下来了。

再仔细一想,我这种有想法却没写出来的现象是有原因的。原因分客观原因和我自己的主观原因。

先来控诉一下这个世界的不公。文章的质量的越高,所需要的投入也越高。正如我想的那些问题那样,一篇高质量的文章,需要先对所谈的话题进行信息的收集,保证全面性、严谨性;再把自己之前的思考同收集的信息整合起来,整理总结自己的思想,确定表述的内容;最后确定好文章的结构,按照合理的顺序,用文字把想法表达出来。文章的质量越高,花费的精力肯定就越多。但是,写一篇好文章,对我来说没有任何利益上的回报。不管我再怎么欣赏自己的才华,再怎么愿意给未来的自己多留一份回忆,也掩盖不了我的内心对于做一件没有收益的事情的抵制。投入和回报不成比例,我写高质量的文章注定会缺乏动力。要不是实在是时间有多,精力有剩,我是无法具备足够的动力的。因此,从客观上来看,对我来说,一直写质量较高的文章是不可能的。

现在可以谈我自己的问题了。我为什么非得要写那么好的文章呢?想到什么就说什么不行吗?非得把博客当成高考作文来写?我开始分析起了自己的心理。我也有想过快速完成一些短一点的文章,但我一想到这些文章可能会写得不够好,像程序一样充满“BUG”,心里就很不是滋味。我想的事情并不是写出好的文章,而是不希望文章中有漏洞,有缺点。我把写文章当成了一种自我展示、自我实现的工具。文章中有纰漏,文章写得不好,就会显得我这个人有缺点。为了让我自己显得更完美,我把这种对于完美的追求代入了文章中。我希望用一篇质量高的文章,来反映一个能力强的自己。我产生这些想法的原因和其他很多表现都是一样的。有些人在乎穿着打扮,认为外貌能反映自己;有些人在乎言行,以自谦来衬托自己的高尚。我虽然不在意别人的看法,但我也在用我写的文章,或者我的一些其他成果,来作为自己好与差的评判标准。对于事物的缺点的掩饰,不过是出于自卑,出于对自己的不足的掩饰罢了。

显然,客观的原因不可改变。除非有人说:“你写的文章太好了,让我五体投地,感慨万分。从今天起,你也不用干其他事了。你写一篇好的文章,我就给你十万块钱。”不然的话,缺乏客观上的利益,我没有足够的动力去写东西的。

但是,我还是有很多东西要表达出来,我得实现我开通博客的初衷。那么,我只能设法改变自己的想法了。我的问题就在于,过分追求完美。

说来也有点讽刺。我对于过分追求完美这个问题,也有过很多思考,也有把这些想法整理出来的打算,但也因为上述种种顾虑迟迟动不了笔。我今天只能说对完美主义进行一个简单的分析,并把如何改变自己而不是探讨理论作为文章的重点。

完美主义的具体定义我是不知道。但很明显,我上述的表现,都是完美主义的常见表现:想把事情做到最好,却最后什么也没有做。

从本质来看,完美主义是不成立的,因为人本身就不是完美的。人的精力有限,人的能力有限,人注定不可能把每件事做到极致。而完美主义所追求的”完美“,不是事物质量的完美,而是能否通过事情的完美展示人自己的“完美”。由于人的限制,这种完美主义注定会失败。

完美主义,错就错在动机。人生来自卑,人想尽可能让自己变得更好,或者说仅仅是看上去很好,这种希望本身没有错。但是,把事情做得很好,做得“最好”,并不等价于你这个人就是“最好”。这两者没有任何关系。甚至反过来,不去急功近利地追求成果,而是潜心雕琢,反而能达到结果的完美。一个的人价值,不能被任何外在的事物所体现,所代表。认识到了完美主义出错的那个关键点,纠正事物完美等于自我完美这个想法,就能克服这种完美主义了。

现在谈回我自己。既然我对自己的想法、自己的问题都深入思考到这个程度了,我以后该做出哪些改变呢?写长的文章没有动力,我就写一些短的文章啊!文章短又如何?质量差一点又如何?只有让写一篇的文章的成本低到一定程度,直到低于我写文章的回报,我自然就有动力去写了。

那又回到我第一篇博客谈的问题了。我为什么要写博客?我能得到什么回报呢?我也不放眼未来,把这个问题展开来谈。我就谈一下近期我写博客能获得的回报。如果是写一篇介绍算法的博客,我能利用写博客这个机会,去把一个算法彻底学会。写的文章好不好,有没有人看,都和我的收益毫不相关。如果是写一篇像这篇文章一样,吐诉内心想法的文章,那就想到什么写什么,管他写的是什么东西,写得自己爽就好了。如果我想把我的想法一些整理出来,给别人看,认同我提出的观点,那我就得提升文章质量,增强说服力,让他人被我的文笔和人格魅力所折服。如果我写的是题解,那么我只能稍微提升一下对于题目的理解,并且略微提升一下我本来就很强大的语言组织能力,再期待这篇题解能被别人看到,增加自己的知名度,收益就不高了。如果我复习时没有动力,我就像上学期那样,通过写博客来敦促自己系统地复习知识。如果我吃饱了没事,想吹吹牛展示一下自己,就可以针对某项成果单独开一篇博客来吹嘘一番。

这样一看,不同类型的博客,收益的种类也是不同的。我之前只是把博客当成一种证明自己完美的手段,而忽略了写博客还有如此多其他方面的收益。只要我能把写博客的成本降低,我自然而然地就能获得写博客的收益,就像这篇文章一样。我需要去掉的成本,就是那些名字好听,被称为“完美”的成本。

题目按个人感受的难度升序排列。

J Fraction Comparision

Description

给$\frac{x}{a},\frac{y}{b},0\leq a,b\leq 10^9,0\leq x,y\leq 10^{18}$,问这两个分数比较大小的结果(大于、小于或等于)。

Solution

由于浮点数会有误差,不能除了以后再判断。

用有大整数的语言的话,可以直接通分比较。

用C++的话,可以把两个分数划为带分数。比较带分数时先比较整数部分,再通分比较真分数部分。由于分母在1e9范围内,不会爆long long。

F Random Point in Triangle

Description

给平面上一个三角形,顶点为ABC。现在随机再三角形内部取一点S,S把三角形划分为三个小三角形。问这三个小三角形面积的最大值的期望乘36(保证这样算出来的结果是整数)。

Solution

我在知乎上看到过类似的题目,原理没记住,结论倒记住了:这种三角形随机取点求面积期望的题,与三角形形状无关,只与面积有关。原理好像是坐标变换什么的,不管三角形形状是怎么样的,都可以通过坐标变换变成同一个三角形。这个坐标变换的过程面积也是按照比例变化的。大概是这个意思,我也说不清楚具体的原理。

有了这个结论,只需要知道面积为1的三角形的答案就行了,这个直接在本地模拟就行。我把三点设为(0, 0), (0, 1), (1, 0),取两个随机数后把随机数都映射到[0, 1]中,再判断这个点是否在三角形内。如果这个点在三角形内,就用叉积算三个小三角形面积。不断地大量随机取点后再除一下有效点数量就是答案。

判断点是否在三角形内可以考虑判断三个小三角形的面积是否“正负号相同”(即叉积结果符号是否相同)。

最后发现面积乘22就是答案,也就是叉积的结果乘11。

B Integration

Description

已知$\int_{0}^{\infty}\frac{1}{1+x^2}dx=\frac{\pi}{2}$,
给定$a_1, a_2…a_n(1\leq n\leq 10^3,1\leq a_i\leq 10^9,a_i互不相同)$,求$\frac{1}{\pi}\frac{1}{\Sigma_{i = 1}^{n}(a_i^2+x^2)}dx$模$10^9+7$意义下的值。

Solution

这道题只需要微积分的基础知识,即积分的线性性质:积分里面如果是相加的话可以拆开来积分再相加。但题目里面的积分是相乘,只能考虑把积分内的乘积化为简单的分式相加,再套入题目给的公式了。

分式的乘拆分成分式的和则是初等数学的内容了。不妨先考虑$n=2$的情况,考虑$\frac{1}{(a^2+x^2)(b^2+x^2)}$的拆分结果。

设$\frac{1}{(a^2+x^2)(b^2+x^2)} = \frac{n}{a^2 + x^2} + \frac{m}{b^2+x^2}$,则:

由于上式恒成立,可得:

解得:

用同样的方法可以求得n为任意值时的拆分结果。对于分母是$a_i^2+x^2$的项,它的分子是$\frac{1}{\Sigma_{j\neq i}(a_j^2-a_i^2)}$。

再由微积分的知识,$\int_{0}^{\infty}\frac{C}{a^2+x^2}dx=\frac{C}{a^2}\int_{0}^{\infty}\frac{1}{1+(x/a)^2}dx =\frac{C}{a}\int_{0}^{\infty}\frac{1}{1+(x/a)^2}d(x/a) =\frac{C\pi}{2a}$,

因此,对于分母是$a_i^2+x^2$的项,积分的结果是$\frac{\pi}{2a_i\Sigma_{j\neq i}(a_j^2-a_i^2)}$。算出了公式,转换成程序就行了。

C Euclidean Distance

Description

给一个n维空间的点$A(a_1/m, a_2/m…a_n/m)$,求它到点$P(p_1, p_2, …p_n)$的最短距离的平方,即$\Sigma_{i = 1}^{n}(a_i/m - p_i)^2$。其中点$P$要满足下列条件:
$p_i\in \Reals$
$p_i \geq 0$
$\Sigma p_i = 1$

答案用分数的形式表示。

数据范围:
$1 \leq n \leq 10^4,1\leq m \leq 10^3, -m\leq a_i \leq m$

Solution

为了方便说明,先把点$P$的每个坐标分量乘一个m,我们要求的东西变成了$\Sigma_{i = 1}^{n}(a_i - p_im)^2$。只要把这个结果除以一个$m^2$就是题目的答案。

由于点$P$的分量都大于等于0,且和为$m$,我们可以考虑把点$P$的坐标看成是实数$m$“分配”到每一个坐标分量$p_i$上。即我们可以让每一个$p_i$变大一点。

我们如何让距离尽可能短呢?如果在某一维上,比如第$i$维,若$a_i < 0$,那么此时令$p_i$增长的话,那么距离只会越来越长;如果$a_i > 0$,那么$p_i$越靠近$a_i$,那么最终的距离就会越短。当$p_i = a_i$时,我们就不用让$p_i$增长了。

由不等式$a^2+b^2\geq2ab(仅当a=b时等号成立)$我们发现,数越大,其平方也增长得越快。如果此时我们去缩小一个本身就很大的数,就能让平方和减少得更多。具体来说,对于两个点的第$i, j$两个维的坐标分量,若$|a_i-p_i|>|a_j-p_j|$,那么让$a_i,p_i$更接近一点,比让$a_j,p_j$更接近一点,更能让总距离减小。那么在分配$m$到$P$的每一维时,应该先分配给$a_i$最大的地方,缩小此处$a_i和p_i$的差距。等这一维的距离和$a_i$中第二大的数时,就要同时缩小这两个维上的距离了。

比如$A(2,3),m = 3$,我们发现$a_2 > a_1$,于是先从$m$中分配$1$给$p_2$,此时$P(0,1)$,$|a_2 - p_2| = |a_1 - p_1|$。我们再让这两个距离同步减小,从$m$再分配1分别给$p_1,p_2$,此时$m=0$,分配结束,而$P(1,2)$,$|a_2 - p_2| = |a_1 - p_1|$依然成立。这种情况下得到的$A,P$两点的距离是最短的。

有了上述思想,算法也就呼之欲出了。把所有$a_i$降序排序,按照上述方法分配$m$的值,尽可能让最大的几个距离相等,再同步减少这几个最大的距离。这种分配方式对于$a_i$是负数的情况也是成立的,不用多加考虑。注意一下有理数分子分母的处理方式即可。

Review

有时对于一道看上去复杂的数学题目,用一些比较简单的、定性的方法就能解出来。

A Equivalent Prefixes

Description

在两个数组$a, b$的前$m$项上,定义一个等价关系:若两个数组前$m$项任意一个区间$l,r$的最小数的位置相同,则称两个数组前$m$项等价。现在给两个数组$a,b$,问它们最多前几项是等价的?

数组长度为$n$

数据范围:
$1 \leq n \leq 10^5$
$1 \leq a_i,b_i \leq n, 所有a_i,所有b_i互异$

Solution

这个等价关系的意思好像就是说,每个数作为最小值的影响范围相同,即每个数碰到的左边、右边第一个比它小的数的位置都相同。这个东西好像就是笛卡尔树啊!

二分答案,对前$mid$项构建笛卡尔树,若笛卡尔树完全相同就check成功。

Review

从知识和写代码的复杂性来看,这道题肯定不是难度第二低的题目。但是,从全场的通过情况来看,这是大部分队第二道通过的题,也是通过人数第二多的题。

可以得出结论:很多题目,不是难到不会写,而是难到不敢写;有本来比较简单的题目没来得及写,不是你没有AC,而是别人没有AC。

话说笛卡尔树的定义和用途其实我现在都没有弄清,但我知道这个东西用单调栈就可以构造,而且所涉及的所有知识都建立在单调栈方面的知识上。

E ABBA

Description

给定$n,m$,问有多少个长为$2(n+m)$的字符串,’A’和’B’各$(n+m)$个,使得该字符串能够被拆成$(n+m)$个长度为2的子序列,其中有$n$个子序列是”AB”,$m$个子序列是”BA”。

数据范围:
$0 \leq n,m \leq 10^3$

Solution

本题最大的障碍在于,题目对于n个AB子序列和m个BA子序列的限制条件,等价于说:对于任意的$i$,字符串前$i$中不能出现A比B多$n$个以上的情况,也不能出现A比B少$m$个以上的情况。证明如下:

充分性:
即证逆否命题:如果字符串中出现了A比B多$n$个以上的情况,则不能满足题目的条件。假设在某个位置,A比B多$n+1$个,则前面这些多出来的$n+1$个A永远不能放在B的后面,也就是最多只有剩下的$(m+n)-(n+1)=m-1$个A能够构成BA,但是$m-1$个BA是不能满足题目的限制条件的。

不能有A比B少$m$个以上的情况的证明类似的。

必要性:
即证满足AB的数量条件即能满足题目的限制条件。考虑构造出题目中要求的$(n+m)$个子序列。从左到右读一个满足AB数量条件的栈,用一个栈来构造这些序列。首先,构造$n$个AB。像括号匹配一样,A当成左括号,B当成右括号。一旦栈中出现了AB,就弹出AB并且令AB子序列数+1,直到构造完$n$个AB。由于在任何时刻不会出现A比B少$m$个的情况,栈中最多有$m$个B,剩下的B全部可以构造成AB。所以这个$n$个AB一定可以构造出来。

栈中可能会剩下很多B,B后面又跟了很多A。而现在我们要构造出$m$个BA。先把栈中已经配对成功的BA弹出,直到栈中全部是A或者全部是B为止。继续从左向右读字符串,如果当前栈顶是B,又读入一个A,则弹出这个BA;若当前栈顶是A,又读入了一个B,则从前面我们构造好的$n$个AB中随便挑一个出来,则之前的AB和现在栈顶的两个AB构成的子序列是ABAB,注意到这个子序列可以拆成A(BA)B,因此我们依然可以构造出一个BA来。后面这种替换是有限度的,如果连续的A碰到了连续的B,则需要之前有足够的AB才行。比如栈里面现在是AAAA,下面4个字符是BBBB,则需要连续替换4次。为了保证替换成功,前面一定要有4个AB才行。但我们现在有前提条件,不存在A比B多$n$个以上的情况,即栈中不可能有$n$个以上的A,也就是说连续替换的次数一定小于等于$n$。而我们已经构造出了$n$个AB,这种替换的一定可行的。

综上,按照上述方法,我们可以构造出题目要求的$(n+m)$个子序列。

有了上面的证明,题目的要求就变成我们自己的等价描述。那么只需要一个dp数组来统计方案就行了。$dp[i][j]$表示长度为$i$的字符串,A的个数减B的个数等于j的方案数。这个dp十分巧妙,dp的范围自动删除掉了A与B数量不符合要求的字符串。那么转移方程为:

$dp[n][0]$就是题目的答案。数组下标不能是负数,稍微处理一下就好了。

Review

很久没打比赛,我的思维现在非常死板。我想这道题的时候误入歧途,一直想用排列组合算这道题。明明我在暴力枚举验证我猜的结论的时候,已经用到了栈来判定字符串是否合法了,但我还是没有想到只要保证A与B之间的数量关系就行了。

当我看到这个结论的时候,觉得这个结论十分简单。但证明起来却发现,充分性很容易想到,必要性的构造却是要花一番功夫才能证明出来的。比赛中我们往往想到了某个结论的充分性,就会匆匆忙忙地猜结论。这样有时候或许侥幸能过题,但一旦过不了就会很伤士气。我平时就要养成严格证明的习惯,并且同时加快证明的速度,对于一个想到的结论要简单地从两个方面证明等价,不要想到一个结论就立刻跃跃欲试,脑袋都空掉了。

H XOR

Description

给$n$个数$a_1, a_2…,a_n$构成的集合,从该集合找出子集,这些子集的异或和为0,求这些子集的集合大小之和。

答案模$1e9+7$

数据范围:
$1 \leq n \leq 10^5$
$0 \leq a_i \leq 10^18$

Solution

首先,求出所有的集合并统计它们的大小是十分困难的。但所有集合的大小和,等于分别讨论每个元素的贡献,再求出每个元素在多少个集合中出现,最后对每个元素出现集合数求和。

集合的异或信息可以用线性基来得到。得到了集合中的一组线性基后,那些不是基的元素无论怎么异或,都可以用这组基得到,也就是和某些基在一起是线性相关的,也就是异或和为0。假设基的大小是$r$,那么还剩$n - r$个不是基的元素。对于每一个剩下的元素分别讨论贡献。$n - r$个元素中,某个元素一定存在的集合数有$2^{n - r - 1}$2^{n - r - 1}$$个。对于每个这样的集合,我都可以选取一些基,使得最后的集合异或和为0。因此,这些元素的贡献都是$2^{n - r - 1}$。

现在要单独讨论剩下的$r$个作为基的元素了。算这$r$个元素的贡献的方法是类似的,如果除掉其中某个元素外,剩余的$n - 1$个元素又构成了一个空间,又选出了一些基,则这个单独的元素又可以被新的基所表示,这个元素的贡献同样是
$2^{n - r - 1}$。

但是,如果运气不好,$r$个基中的某个元素是必不可少的一个基。少了这个基,整个空间就少了一维。换句话说,谁和它的异或和都不能为0。在这样的情况下,这个元素就毫无贡献了。

那么对于作为基的$r$个元素,我们的任务就是判断它们是否可以被任意一个其它的基所表示。对于某个基,有一种很简单的方法能构造出其它$n - 1$个元素的基。先构造出$n - r$个元素的基,它们代表了剩余$n - r$个数的空间。再把$r$个基中除了当前这个基外的$r - 1$个基也加入这个新的空间里面。现在,这个空间就能代表剩下$n - 1$个数了。判断选出的这个基是否在新的空间里即可。

Review

在写这题之前我是不会线性基的,但稍微学了一下就发现,只要学过了线性代数的话很快就能理解线性基。只要理解了线性基,那么这题就非常简单了。

不过这题的一个关键是把问题从求集合大小转换成求每个元素出现集合数。这个转换十分关键。我由于是看着题解写的,没有在这个问题上碰到障碍。但若是要在比赛中快速把难求的问题转换成容易求的问题,就需要极高的解题灵活度了。这种能力只有通过大量写题目才能获得。这也是大量刷题培养手感的意义所在吧。

I Points Division

Description

给$n$个点,每个点有$a,b$两种价值。现在要把所有点分成两个部分$A,B$,要保证$\neg \exists p_i \in A, p_j \in B,x_i \geq x_j \bigwedge y_i \leq y_j$。在$A$里面的点算$a$价值,在$B$里面的点算$B$价值。求最大价值和。

数据范围:
$1 \leq n \leq 10^5$
$1 \leq x_i, y_i, a_i, b_i \leq 10^9$
$所有点互异$

Solution

显然,我们要做的第一件事就是把划分$A,B$的条件弄得更直观一些。题目的要求是:没有两个点,使得左上角的点属于$B$,右下角的点属于$A$。也就是说,如果一个点被分配进了$B$,那么它右下角的所有点都得进入$B$;如果一个点被分配进了$A$,那么它左上角所有的点都得进入$A$。

也就是说,空间被分成了左上角和右下角两部分。那么是谁把空间划分的呢?仔细一想,划分空间的应该是一条折线。折线左上角的点属于$A$,折线右下角的点属于$B$。不妨再规定一下,折线上面的点都属于$B$(同一种划分下,如果折线上的点属于$A$的话,会得到另一条不同的折线,但效果是一样的)。这条折线是从左下向右上的,也就是随着x的增长y是不减的。

分$A,B$算价值太麻烦了。不妨假设所有点最开始都属于$A$,然后我们试图让一些点加入$B$,使得总价值能够有所增长。这样的话,每一个点的价值就唯一了,新价值是$b - a$。

之后要统计结果并保证不重复,不妨尝试对所有点从左到右进行统计。

如果选择了一个点,那么$x$坐标大于等于它且在它下方的点一定也要被加入,也要统计答案。而只要唯一确定了一条折线,我们就确定了一种方案。因此,我们只要考虑折线上的关键点有哪些,让那些右下角的点在计算关键点的时候被计算到。我们从左到右考虑每个点作为关键点的结果。但如果在把当前点加入的同时就想统计右边所有的点,并立刻算出答案,那么再考虑下一个点时,它就不好知道右边有哪些点已经被左边之前的点计算过了。所以,在考虑加入一个关键点的时候,不能立刻算出它的最终贡献,只能考虑在所有已经访问过的点中,它的临时答案。而经过后面那些点时,后面的点会对之前所有$y$坐标大于等于它的临时答案都加上自己的贡献,因为这些点被强制加入了。在讨论当前进入$B$的点时,不直接计算所有右下角的点,而是反过来,先计算一个临时答案,等碰到了右下角的点再加上它们的贡献。这种统计方式是本题的关键。

对于一个关键点,我们已经能统计它右下角的点了,现在来考虑如何计算当前点的临时答案。如果把一个点放到折线上,那么之后下一个点无论怎么选点,都与折线上左边的点无关,而只与现在选的这一个点有关。这表明了选择折线时的子结构性质。既然如此,我们把当前点放到折线上时,也不用管以前的点,只要选一个最优的上一个点就行了。这里的最优,指的是临时答案的最优。对每个点选择上一个点的过程,和动态规划的过程是一样的。发现最优子结构性质,是本题的第二个关键。

现在我们已经可以得到一个正确的算法了:从左到右,从上到下(这个很关键)地枚举每一个点,讨论它作为折线上的点时的情况。对每个点,考虑所有之前走过的,且$y$坐标小于等于它的点(要保证能和这个点构成不降折线),从中选一个临时答案最大的点作为这个点的上一个点。此时,新的这个点的临时答案和上一个点是一样的,因为我们还没有考虑新的这个点自己还有它右边和下面的点。我们把这个点的临时答案也写入记录里面。最后,把已经访问过的,且$y$坐标大于等于它的所有点的临时答案都加上新的这个点的价值,表示这个点对之前所有$y$坐标大于等于它的点的贡献。只要对每个点都进行上述的操作,等处理完最后一个点后,每个点的临时答案就变成了最终答案。从所有最终答案选一个最优的就是题目的答案了。

但是,这样的话时间复杂度似乎不能接受。遍历每个点肯定要右$O(n)$,讨论之前所有点并修改临时答案是$O(n)$,总的时间复杂度是$O(n^2)$,肯定会超时。我们必须优化一下从之前点找最优值以及修改值的过程。

我们可以发现,我们在找之前所有点和修改答案时,都是在访问$y$坐标连续的点。找最大临时答案是找所有$y$小于等于当前点的点,修改是修改$y$大于等于当前点的点。那么,我们可以把所有点的临时答案按照$y$坐标的顺序存储。我们存的不再是某个点答案,而是某个$y$坐标的答案,这样我们依然不会漏掉某些情况:如果后面出现了一个和前面$y$坐标相同的点,如果选了之前的点的话,后面这个点一定要选到。因此,不存在不选新的这个点而只选前面那个点的情况,我们可以放心大胆地覆盖掉之前那个临时答案。现在所有答案都存到了$y$轴上。对于连续区间的修改、查询,我们就可以用线段树来维护了。线段树完成区间加一个值、找到区间最大值这两个操作就行了。

$y$坐标的范围是$[1,10^9]$,读完数据后离散化一下就可以了。

最终时间复杂度是$O(nlogn)$

Review

这道题考察的知识不难,但对问题的思考、分析是本题的难点,在修改哪些点、访问哪些点、如何选择枚举方向等一些细节上也需要多加注意。我感觉这道题出得非常好。

我自己写的时候犯的错误只有搞错了从左到右,从上到下的枚举方向,其它部分都很顺利。由于我是看了题解后再写的,少了很多对于问题的思考、分析,我不能保证下次碰到同样难度的题目也能写出。只能说,很多题目用的算法、数据结构大家都会,但是要把问题转换到能用熟悉的方法来解决的那一步,不需要太多的经验积累,而需要对问题做出正确的分析。碰到难题要敢去想。

在过去的半年里,由于自己和外部的种种原因,我没有认真地进行ACM训练。用愧疚来形容我的感受的话,或许有一些夸张。更多的感受或许是无奈与遗憾。

但是,暑假ACM集训即将开始了。这或许是我人生中最后一段能够把全部精力投入算法竞赛的时光。每天不需要考虑太多,只是享受思考与解题的快感就好。

为了能够让自己更有动力,更有目标性,我决定开一篇文档来记录自己的刷题成果。我在这个暑假计划写完200道题。每场比赛都能补7~8题的话,就已经能够完成150题了,可见这个目标的达成并不困难。但是,我希望我写每一道题都能认真地去思考,有质量地完成,并顺便补习题目中出现的不会的知识。

9.7:虽然暑假过了,但是过题清单统计到所有网络赛打完吧

9.17:虽然我没有完成我的目标,但我这个暑假还算是投入了很多精力在ACM上,学到了不少了东西的。开学后做的题目会另开一个文章来存。

下面是每道题完成的情况及题解链接(如果有题解的话):

下划线表示这道题很有价值/我很喜欢

进度(120/200)

7.18

共计4题
2019牛客1A(笛卡尔树)
2019牛客1B(数学)
2019牛客1F(找规律)
2019牛客1J(水)

7.19

共计4题
2019牛客1C(贪心、数学)
2019牛客1E(贪心、dp)
2019牛客1H(线性基)
2019牛客1I(二维点,线段树)

7.20

共计2题
2019牛客2F(蛮力搜索)
2019牛客2H(单调栈)

7.21

共计8题
2019CCPC江西省赛1
2019CCPC江西省赛4
2019CCPC江西省赛6
2019CCPC江西省赛7
2019CCPC江西省赛8
2019CCPC江西省赛9
2019CCPC江西省赛10
2019CCPC江西省赛11
(今天题比较水)

7.22

共计3题
2019HDU多校1-2(区间异或和最大)
2019HDU多校1-4(思维)
2019HDU多校1-9(序列自动机、贪心)

7.24

共计5题
2019HDU多校1-1(dp)
2019HDU多校1-12(组合数学、卷积)
2019HDU多校2-I(回文树)
2019HDU多校2-J(签到)
2019HDU多校2-K(可持久化线段树)

7.25

共计1题
2019牛客3-F(st表、双向队列维护区间最值)

7.26

共计4题
2019牛客3-G(笛卡尔树上区间二分)
2019HDU多校1-6(后缀自动机)
2019HDU多校2-H(最小割)
2019HDU多校2-L(线段树维护信息)

7.27

共计3题
2019牛客4-A(树形DP)
2019牛客4-I(后缀自动机、回文树)
2019牛客4-J(分层图(令边权为0))

7.28

共计2题
2019HDU多校2-F(FWT)
2019牛客4-C(笛卡尔树、区间最值)

7.29

共计2题
2019HDU多校3-4(二分、线段树/平衡树)
2019HDU多校3-7(二分、线段树/平衡树)

7.30

共计2题
2019HDU多校1-11(莫比乌斯反演)
2019HDU多校2-E(期望、简单dp)

7.31

共计3题
2019HDU多校4-3(思维)
2019HDU多校4-7(结论题)
2019HDU多校4-8(二分、可持久化线段树)

8.1

共计3题
2019牛客5-G(中等dp)
2019牛客5-F(思考、补图、二分图最小覆盖方案)
2019牛客5-H(字符串前后关系转化成对图求环)

8.3

共计4题
2019牛客6-D(暴力模拟)
2019牛客6-E(构造自补图)
2019牛客6-J(简单dp)
2019HDU多校4-6(思维、斜率优化)

8.5

共计1题
2019HDU多校5-6(exkmp签到)

8.7

共计5题
2019HDU多校6-4(两颗权值线段树)
2019HDU多校6-6(曼哈顿距离模意义下性质)
2019HDU多校6-8(暴力)
CodeForeces 674C(斜率优化DP)
2019HDU多校5-1(分母最小的给定范围下的分数)

8.8

共计4题
2019牛客6-A(最小表示法暴力DP)
2019牛客6-B(三次以上多项式必可因式分解结论)
2019牛客6-D(签到)
2019牛客6-H(繁琐的数位DP)

8.9

共计5题
Comet8-A(签到)
Comet8-B(签到)
2019HDU多校6-2(随机数据、LIS、暴力重构)
2019HDU多校6-5(线段树维护最大子段和)
2019HDU多校6-12(签到)

8.10

共计1题
2019牛客8-C(找规律)

8.11

共计7题
Comet8-C(dp或者线段树维护最大子段和)
Comet8-D(思考、树状数组)
2019牛客8-A(单调栈、最大矩形去重)
2019牛客8-B(签到、单独考虑贡献)
2019牛客8-D(暴力BFS重构最短曼哈顿距离)
2019牛客8-G(签到)
2019牛客8-J(组合数学、容斥原理)

8.12

共计1题
2019HDU多校7-8(思考)

8.13

共计5题
2019HDU多校7-1(高精度模拟)
2019HDU多校7-6(思维)
2019HDU多校7-7(思考、dp、二分转移)
2019HDU多校7-10(思维、map排序)
2019HDU多校7-11(期望dp)

8.14

共计3题
2019HDU多校8-4(构造)
2019HDU多校8-9(离散化后DFS)
2019HDU多校8-11(思维)

8.15

共计2题
2019牛客9-B(模意义方程、二次剩余)
2019牛客9-E(并查集)

8.17

共计6题
2019牛客10-B(递归)
2019牛客10-D(Java写的中国剩余定理)
2019牛客10-F(线段树维护最值)
2019百度之星初赛-1(签到)
2019百度之星初赛-2(签到)
2019百度之星初赛-5(找规律)

8.18

共计3题
cf-A (水题)
cf-B (无向图最短环)
2019HDU多校8-6(树形dp)

8.19

共计1题
2019HDU多校9-3(硬核折半搜索+基数排序)

8.21

共计2题
2019HDU多校10-3(简单概率)
2019HDU多校10-9(签到DFS)(全场一血)

8.23

共计1题
2019CCPC网络赛1(签到)

8.31

共计6题
2019百度之星复赛-A(树形dp)
2019百度之星复赛-B(思维)
2019百度之星复赛-C(贪心)
2019银川网络赛C(签到)
2019银川网络赛D(猜结论概率)
2019银川网络赛H(贪心排序)

9.1

共计3题
2019南京网络赛B(欧拉降幂)
2019南京网络赛F(区间最值)
2019南京网络赛G(组合数学)

9.7

共计4题
2019徐州网络赛G(回文树)
2019徐州网络赛I(区间询问,固定右端点)
2019徐州网络赛L(打表,TSP问题)
2019徐州网络赛M(序列自动机贪心)

9.8

共计2题
2019银川网络赛重赛E(水)
2019南昌网络赛E(约瑟夫问题)

9.9

共计3题
min_25筛入门题X3

9.14

共计2题
2019沈阳网络赛F(二分查找)
2019沈阳网络赛G(物理)

9.15

共计3题
2019上海网络赛B(离散化、区间异或)
2019上海网络赛C(FFT、简单容斥、不同数据量不同算法)
2019上海网络赛E(指数生成函数)

递推关系与生成函数

递推关系

斐波那契数列通项

俗话说,斐波那契数列的通项是:

$\frac{1}{\sqrt{5}}(\frac{1 + \sqrt{5}}{2})^n - \frac{1}{\sqrt{5}}(\frac{1 - \sqrt{5}}{2})^n (n \geq 0)$

看上去如此玄妙的式子是怎么推出来的呢?

我们都知道,斐波那契数列是:$0,1,1,2,3,5,8,13,21,34,55…$。规律是,从第三项开始,后一项等于前两项之和,即$h_n = h_{n - 1} + h_{n - 2} (n >= 2), h_0 = 0, h_1 = 1$。通过观察,我们发现这个数列增长得十分快,速度似乎和指数函数差不多。因此我们产生了大胆的假设:假设$h_n = a^n$。

由通项公式:$h_n = h_{n - 1} + h_{n - 2}$,$a^n = a^{n - 1}+a^{n - 2}$。由于此时$n\geq2,a_n\neq0$,式子化简为:

$h_n$有多个解,这些解的线性组合同样满足条件。因此可以将解表示为$h_n = c_1(\frac{1 + \sqrt{5}}{2})^n + c_2(\frac{1 - \sqrt{5}}{2})^n$。注意到还有条件$h_0 = 0, h_1 = 1$没有使用。等式有两个,未知数也有两个,因此可以解出$c_1, c_2$来:

解得:

所以就得到了开头说的那个通项公式。

可以发现,解通项公式的过程和求常系数线性微分方程的方法是类似的,碰到一些的问题的处理方式也是类似的。

常系数线性齐次递推关系

通过上述的斐波那契数列的通项的求解,我们已经可以感受到递推关系的解法了。

若某序列$h_i$满足递推关系$h_n = a_1 h_{n - 1} + … + a_k h_{n - k} + b$,其中$a_i,b$可能依赖于$n$,$a_k \neq 0$,则称序列满足k阶线性递推关系。若$b = 0$,则还称它是齐次。若$a_i$都是常数,则称它还是常系数。正如标题所讲,本节讨论常系数线性齐次递推关系

由递推关系得到的方程$x^k = a_1x^{k - 1} + … a_k$被叫做特征方程。特征方程的k个根叫特征根,一个特征根对应一个原递推关系的解。如果p是根,则$p^n$是$h_n$的一个解。如果p是二重根,则$p^n,np^n$是$h_n$的解,依此类推。$h_n$解的最终表达形式是每个特征根对应的解的线性组合。如果此时数列$h_n$有初始值,那么再解个多元一次方程,就可以完全求出$h_n$的通项了。

综上,求常系数线性齐次递推关系就是先解一个一元多次的特征方程,求出数列的表达形式,才根据初始条件解一个多元一次方程,得到数列的通项。

常系数线性非齐次递推关系

由于特征方程有多个根,因此递推关系有多个解,这些解要以通解的形式写出来。而某一个单独能满足条件的递推关系叫做特解。又根据叠加原理,$h_n + a_1 h_{n - 1} + … + a_k h_{n - k} = b$和$h_n + a_1 h_{n - 1} + … + a_k h_{n - k} = c$的解加起来,就是$h_n + a_1 h_{n - 1} + … + a_k h_{n - k} = b + c$的解。因此对于一个常系数线性非齐次递推关系,不用再像之前一样直接求出通解,而是只要求出一个特解,再加上齐次情况的通解,就可以得到该递推关系的特解。我们现在要完成的任务就是求常系数线性非齐次递推关系的特解。

求非齐次的特解我只知道公式,不理解公式的由来,这里就直接给公式了。若:

则将

代入$h_n$,解出这些$c$来,就可以得到一个特解。如果$r$还是特征方程的根,且是$m$重根,则要代入

可以看出,求特解时,就是把等式右边的非齐次项按$r^n$的指数分个类。再顺便看一下$r$是否是特征方程的解,若是的话要乘上$n$的次方项。特解的未知数个数取决于非齐次项中的多项式的次数,把这个式子套回去可以解出所有的未知数来。具体看下面的例子。

例:求$h_n = 2h_{n - 1} + n^2$的特解。
解:
特征方程:

代入得:

上式恒成立,因此:

解得:
$
\left\{
\begin{aligned}
c_0 &= -6 \\
c_1 &= -4 \\
c_2 &= -1
\end{aligned}
\right.
$

得到一个特解:
$h_n = -6 - 4n - n^2$

生成函数

生成函数的用处

若有一个数列$h_0,h_1,h_2…$,则它的生成函数定义为$g(x) = h_0 + h_1x + h_2x^2…$。若数列是无穷的,则这个级数也是无穷的。

生成函数乍看上去非常奇怪,因此必须要先理解生成函数是用来做什么的。生成函数正是为了解决一些问题而定义出来的。

设$h_n$为从集合$\{2a\}$中取出n个字母的方案数,则数列为$1,1,1,0,0…$,生成函数$g_1(x) = 1 + x + x^2$。同理,从集合$\{3b\}$中取出n个字母的方案数得到的生成函数$g_2(x) = 1 + x + x^2 + x^3$。那么现在,从集合$\{2a,3b\}$中取出n个字母的方案数得到的生成函数是什么呢?

还是先求出方案数序列。$h_0 = 1,h_1 = |\{a,b\}| = 2, h_2 = |\{aa,ab,bb\}| = 3,h_3 = |\{aab,abb,bbb\}| = 3, h_4 = |\{aabb,abbb\}| = 2, h_5 = |\{aabbb\}| = 1$。则其生成函数$g(x) = 1 + 2x + 3x^2 + 3x^3 + 2x^4 + x^5 = (1 + x + x^2)(1 + x + x^2 + x^3) = g_1(x)g_2(x)$。“巧合”发生了。一个问题的解的生成函数等于子问题的解的生成函数相乘。

生成函数的本质,就是把“总和”为$n$的方案数放到$x^n$的系数上。多项式在进行乘法时,会很自然地进行排列组合,排列组合计算的结果依然保留在系数上。我目前的理解是,生成函数就是用来计算排列组合的。

生成函数既然是用来计算排列组合的,但我们现在都是在通过排列组合反推生成函数。如果是用多个生成函数的乘积来计算总的生成函数,再算出方案数,那么我们还是要手动进行多项式乘法,这个过程和我们算排列组合的过程是一样的,计算量完全没有减少。那么,生成函数到底有什么用呢?

当生成函数是无穷级数时,生成函数的形式就美观了许多,比如:

这样,再求一个复杂的排列组合时,可以先求出简单的子问题的生成函数,才把这些生成函数相乘。当生成函数是无穷级数的时候,往往能进行某些化简。求最终的生成生成函数的某个系数,就可以得到对应的方案数。

看一道经典例题:
例:确定苹果、橘子、香蕉和梨的袋装水果的袋数$h_n$的生成函数,其中各袋要有偶数个苹果,最多两个橘子,3的倍数个香蕉,最多一个梨。然后从该生成函数求出$h_n$的公式。
解:

$\therefore h_n = n + 1$

指数生成函数

在求多重集合的排列问题,比如从集合$\{2a,3b\}$中能拼出多少个长度为3的单词这样的问题时,我们要先求出各种组合的方案,除以每个数出现次数的阶乘,最后再乘上总数的阶乘。生成函数既然是解决排列组合问题的,自然也能解决此类多重集排列问题。在刚才的生成函数的基础上,加上我们平时求多重集排列的方法,就得到了指数生成函数的概念。

类似地,对于序列$h_n$,其指数生成函数定义为$g^{(e)}(x) = h_0 + h_1\frac{x}{1} + h_2\frac{x^2}{2}…$。由于在分母上除了一个阶乘,生成函数很容易变成$e^x$之类的形式,大概指数生成函数的名字就是这样得来的。$x^n$前面的系数再乘一个$n!$就是多重集n排列的方案数。

用法看例题即可快速理解:
例:如果偶数个方格被涂成红色以及奇数个方格被涂成白色,试确定用红、白、蓝和绿为1行$n$列棋盘的方格着色的方案数$h_n$。
解:

$\therefore h_0 = 0, h_n = 4^{n - 1}(n \geq 1)$
无穷级数带了一个常数项,因此通项在0的时候要单独写出来。

ACM中也有一些题要用到生成函数的思想,但解题时最后还是要用FFT。等以后碰到了相关的题目,我会把题解写在ACM分类中。

递推关系与生成函数

用生成函数解递推关系

好吧,生成函数还是有其他用途的,比如求解递推关系。不然为什么课本把这两个概念放在同一章呢?

只用一道例题来说明用生成函数求常系数线性齐次递推关系的方法:
例:
求$h_n - 5h_{n - 1} + 6h_{n - 2} = 0,h_0 = 0, h_1 = 1, h_2 = 5$的解。
设数列的生成函数为$g(x)$,可以得到:

把上面三个式子左右相加,由$h_n - 5h_{n - 1} + 6h_{n - 2} = 0$:

$\therefore h_n = 3^n - 2^n$

虽然感觉这种方法很高大上,但是比起直接解特征方程来,这种方法还是太麻烦了。既然如此,我们是不是可以不用生成函数来解递推关系的问题了?

当然不是。仔细看一下上面的解题过程,我们求解递推关系完全是在对生成函数$g(x)$做各种各样的操作,好像和递推关系是否线性、齐次什么的一点关系也没有。对于最开始特征方程法解不了的一些递推,用生成函数可能可以解出来。不妨看下面一个问题:

例:(Catalan数)对于一个有$n$个节点的二叉树,在保持左右子树的有序性,即不考虑任何变换的情况下,有多少种可能的形态?(比如$n = 2$时有两种,一个节点是根,另一个节点是根的左儿子或者右儿子)

解:
为了方便,记$n$个节点的二叉树的形态数为$h_n$。

考虑把问题转换成子问题。当现在有$n$个节点时,我们需要先选一个节点为根,再把剩下的节点放到根的左右子树中。如果根的左子树放了$a$个节点,那么右子树就要放$n - 1 - a$个节点,也就是说左右子树的节点和一定是$n - 1$。在这种情况下,问题就变成了两个子问题:$a$个节点的二叉树和$n - 1 - a$个节点的二叉树有多少个?我们把这两个数算出来,再一乘,就是当左边放$a$个节点时,$n$个节点的二叉树的形态数。

当然,我们不能只算$a$个节点,因为根左右子树的节点数可能是$0, 1,2,3…n - 1$。因此我们要讨论所有情况,并且把每种情况下的总情况数算出来,再进行累加。因此可以得到递推公式:

$h_n = \Sigma_{i = 0}^{n - 1}h_{i}h_{n-1-i}$

特别地,这里要令$h_0 = 1$。

由于这个递推关系根本不是线性的,我们无法用特征函数来解。考虑用生成函数$g(x)$来表示数列的每一项,即:

$g(x) = \Sigma_{i = 0}^{\infty}h_ix^i$

突然,我们发现:

$g^2(x) = h_0 + (h_0h_1 + h_1h_0)x + (h_0h_2 + h_1h_1 + h_2h_0)x^2 + …$

每一项中x的系数在递推公式中都可以找到。把递推公式代入式子:

$g^2(x) = h_0 + h_2x + h_3x^2 + …$

由于$h_0 = h_1 = 1$:

$
\begin{aligned}
g^2(x) &= h_1 + h_2x + h_3x^2 + … \\
&= \frac{1}{x}(h_1x + h_2x^2 + h_3x^3 + …) \\
&= \frac{1}{x}(g(x) - h_0) \\
&= \frac{1}{x}(g(x) - 1) \\
xg^2(x) &= g(x) - 1 \\
\end{aligned}
\\
\begin{aligned}
xg^2(x) - g(x) + 1 &= 0 \\
g(x) &= \frac{1 \pm \sqrt{1-4x}}{2x}
\end{aligned}
$
由于$h_0 = 0$,$g(x) = \frac{1 - \sqrt{1-4x}}{2x} = \frac{1}{x}(\frac{1}{2} - \frac{1}{2}(1 - 4x)^{1/2})$

根据牛顿二项式定理:
$
\begin{aligned}
(1 + x)^{1/2} = 1 + \Sigma_{n = 1}^{\infty}\frac{(-1)^{n - 1}}{n\times 2^{2n-1}}\tbinom{2n - 2}{n - 1}x^n
\end{aligned}
$

所以:
$
\begin{aligned}
(1 - 4x)^{1/2} &= 1 + \Sigma_{n = 1}^{\infty}\frac{(-1)^{n - 1}}{n\times 2^{2n-1}}\tbinom{2n - 2}{n - 1}x^n(-4)^n \\
&= 1 + \Sigma_{n = 1}^{\infty}\frac{(-1)^{1}}{n\times 2^{-1}}\tbinom{2n - 2}{n - 1}x^n \\
&= 1 - 2\times \Sigma_{n = 1}^{\infty}\frac{1}{n}\tbinom{2n - 2}{n - 1}x^n
\end{aligned}
$

套回开始的式子中:
$g(x) = \Sigma_{n = 1}^{\infty}\frac{1}{n}\tbinom{2n - 2}{n - 1}x^{n - 1}$

所以,$h_n = \frac{1}{n + 1}\tbinom{2n}{n}$。

如果我们把$h_0$看成第一项的话,那么第$n$项的值就是$\frac{1}{n}\tbinom{2n - 2}{n - 1}$。这个数就是著名的卡特兰(Catalan)数。

在这个过程中,我们完全是在对生成函数进行求解,而没有管递推关系的形式是怎么样的。这其中最妙的一步是发现了$g^2(x)$正好可以得到递推关系中的右半部分,只能说这种解法不是推算出来,而是突发奇想的。

总结

有时候推出来看似花里胡哨的递推关系,只要式子是线性的,就可以利用特征方程的方法求出通项来,从而把O(n)的时间复杂度变成O(1),把编程问题变成数学问题。

生成函数是用来解决组合问题的,当问固定数量的方案数时,就要想到生成函数了。在求排列时,要利用之前的方法,求出指数生成函数,再在结果前乘一个阶乘。

递推关系对应一个生成函数,可以用生成函数来解递推关系。对普通的线性方程,这种解法太麻烦了,还不如用普通的解法。对于一些比较复杂的递推关系,生成函数才能发挥出它的用途。

计算是什么呢?

数学,一门从数字与计算中出现的学科,在发展到一定阶段后,开始追问起了计算的本质。

数学家们用数学的语言,把“计算”这个早就存在的词语严格定义。用自身定义自身,这正是数学之美。

更令人开心的是,数学家在研究计算的过程中,开创了计算理论这门学科。人们在计算理论的基础上,制造了计算机。再随着计算机的应用不断变广,计算机成了一门单独的学科——计算机科学。

又扯远了。计算理论研究了什么是计算的机器,哪些问题是人们可以计算出来的,可以计算的问题怎么才能计算得更快。为了开始计算理论的学习,我们需要从最简单的计算模型,来一步一步理解计算理论研究的内容。有穷自动机及正则语言,就是我们的第一个学习对象。

第1章 正则语言

事物的状态

每时每刻,世界的万物都在发生变化。比如,昼夜不断交替着,早晨太阳升起,傍晚太阳落下;今天的我,总是比昨天的我更加帅气一点。

有些变化,它变来变去就是那么几种情况。比如,天要么是亮的,要么是暗的。有些变化则不然,是会一直进行下去的。我每天都在变帅,我的帅气程度永远在递增,不会有两个相同的值。

为了更好地掌控那些状态有限的事物,人们用一个有向图来表示这些事物。在有向图中,顶点表示状态,边表示转移条件。比如

黑夜——————-太阳升起————————> 白天
<—————太阳落下———————————-

状态,就是事物本质的情况;状态发生转移,就是外界条件对事物的本质产生了改变。

有限自动机

人们在研究一个数学问题是否可以解决时,想到数学问题涉及的内容都可以转化成字符串。问题解决,就是问题对应的字符串在一系列验证过程后,该字符串被认为是正确的。比如问一个十进制数字是否是偶数,所有问题的输入都转换成一个十进制数字字符串。”1”被认为不是符合要求的字符串,”666”被认为是一个符合要求的字符串。而问题的所有解,则是一个字符串集合,也就是一个语言

于是人们把数学问题转换成字符串是否“正确”的问题。字符串的每一个字符,都会对字符串是否正确产生影响,都会改变字符串的整体性质。联系开始我。们对于状态有向图的分析,状态转移正是事物本质发生变化的过程。我们似乎可以用开始的状态有向图,来判定一个字符串。因此,人们把字符串的每个字符做为状态转移的条件,把一些状态设为“正确状态”,意味着能通过字符一步一步走到这里的字符串是正确的。当然,还有一个初始状态,表示在什么字符也没有读取时的状态。

这样讲还是抽象了一点,还是举个例子吧。假设我们用的字符全部是0(也就是说字母表是$\{0\}$),我们现在碰到的问题是:给定一串字符串,问其长度是否为偶数。我们该怎么样用一个状态有向图来解决这个问题?或者说,如何把所有解表达出来呢?

从本质上来看这个问题,一串只有0的字符串只有长度为奇数和长度为偶数这两种可能。每多加一个0,就会改变字符串长度的奇偶性。也就是说,我们有两个状态:“偶数长度状态”、“奇数长度状态”。在碰到0的时候,两个状态间会互相转移。最开始时,字符串处于“偶数长度状态”,因为我们读入的字符串长度为0,是个偶数。如果我们读完了字符串后,发现我们停留在了“偶数长度状态”,那么这就说明该字符串的长度是偶数。状态图画出来是这样的:

Sample

左上角的箭头表示最开始进入的是q0状态,也就是“偶数长度状态”。右边箭头上的0表示转移碰到0就往箭头的方向转移。左边q0状态里面有一个小圆圈,表示这个状态是最终我们能够认可的状态。

我们构建的这个状态有向图,有一个十分大气的名字“有穷自动机”。这个名字为什么这么叫呢?大概有穷,指的是状态数是有限个。自动机,指的是我们只要把字符串到这个有向图里,按照规则在里面走迷宫,我们最终就可以知道这个字符串是否符合我们的要求。整个有向图就像一个自动工作的机器一样。虽然这个名字看上去很厉害,但正如我们分析的一样,有穷自动机就是一个有向图,它的概念十分容易理解。

每一个有穷自动机都可以由5部分确定:$(Q,\Sigma,\delta,q_0,F)$。$Q$是状态(点)集合,$\sigma$是字母表,$\delta$是转移函数(边)集合,每个转移函数接受的参数是当前状态与碰到的结果,输出的是下一状态,也就是说转移函数集合是$Q\times\Sigma\to Q$。$q_0$是初始状态,也就是我们走迷宫的起点。$F$是接受状态集合,$F\subseteq Q$,也就是迷宫的所有终点。我们需要字母表,是因为无论碰到什么字母,我们都得确切地知道我们下一步该往哪里走。我们在转移函数中,必须写清楚每个字母的转移情况。因此也可以得出,转移函数的集合有$|Q|$行$|\Sigma|$列,代表每个点在碰到每个字母转移到的下一个点。

显然这个定义是数学家给出来的。如果是程序员发明这个东西的话,一定会给每个量取一个好听的变量名,以更好地记忆和理解每个量的意思。不过和其他所有的数学定义一样,这些东西是不用背的,只要理解了它们的意思就好了。这些希腊字母和带下标的字母看起来确实让人头疼。

正则语言与正则运算

有一个很重要的名词,我们明明见过它,却只有现在才能介绍它。正则语言,就是某个有穷自动机可以识别的语言。既然有这个定义,就暗示世界上还有很多语言,很多字符串的集合,是有穷自动机表示不了的。话说回来,作为本章的标题,正则语言竟然是通过有穷自动机定义的,真是没有牌面。

但是,正则语言之所以不叫有穷自动机语言,是因为还有一些概念是“姓”正则的。我们之前讲过,所有数学问题都可以被转换成字符串,问题间的运算就转化成了字符串集合运算,也就是语言运算。比如找到两个问题任一的解,就是这个字符串满足问题一的解或者问题二的解。也就是说,问题转化成了语言求并的过程。类似的,语言还可以互相连接(符号$\circ$)、自我重复(符号$\ast$,所以也被叫做星号)。$A\circ B = \{xy|x\in A \bigwedge y\in B \}$,$A\ast= \{x_1x_2..x_k |k \geq 0 \bigwedge x_i\in A\}$。当然,求并和连接是二元运算,自我重复是一元运算。这三种运算都叫做正则运算。为什么这些也“姓”正则呢?因为正则语言在做了正则运算后,还是正则语言,也就是正则语言在正则运算下封闭。

正则语言在正则运算下封闭这个定理是可以证明的。要证明此定理,要分别证明正则语言在3种运算下封闭。不过在当前条件下,我们比较方便证明的只有求并运算。

定理:正则语言在并运算下封闭。

证明:
设原来的语言为$A_1$,$A_2$,状态集合$Q_1,Q_2$。现在我们构造一个新的自动机的状态集合$Q$,它有$|Q_1|\times|Q_2|$个元素,不妨把它们排成$|Q1|$行$|Q2|$列。其中,第$i$行第$j$列个状态表示处于原来自动机$A_1$的第$i$个状态和$A_2$的第$j$个状态。也就是说,我们在新的自动机上完全模拟出了之前两个自动机的状态。

新自动机字母表是之前字母表的并。在新自动机上,每一列状态都按$A_2$的规则向左右转移,每一行状态都按$A_1$的规则上下转移,如果是碰到某个自动机之前不存在的字母,就转移到一个失败状态——不管碰到任何字符都回到自己,且不是接受状态的状态。初始状态是既处于$A_1$初始状态也处于$A_2$初始状态对应的状态。$A_1$接受状态对应行上、$A_2$接受状态对应列上的所有状态都是新自动机的接受状态。

新自动机表示的语言就是前面两个语言的并。因此正则语言在并运算下封闭。

为了证明另外两种运算,我们还需要一个新的工具。

非确定有限自动机

在一个偌大的迷宫中,路径错综复杂,你会怎么办呢?在毫无办法的情况下,我们只能选择像无头苍蝇一样随便选择下一步了。毕竟,不停地往前走总比站在原地强。

有一类比较任性的有穷自动机,它们在碰到了一个字符后,可能选择往多个地方走,甚至拒绝接受这个字符。它们在某个状态的时候,还没碰到下一个状态,就有可能急不可耐地往其它一些状态触发。这样的有限自动机叫做非确定有限自动机(Nondeterministic Finite Automation,NFA),之前我们熟悉的有限自动机叫做确定有限自动机(Deterministic Finite Automation, DFA)。为了节约宝贵的时间,后文用简称来称呼它们。

NFA可以看成是一个会影分身的人在自动机上走迷宫。读取到下一个字符后,他可能会召唤多个自己的分身,和自己走不一样的路;可能拒绝接受这个字符,让自己就此消失。没有读取到字符,或者说读取到空字符$\epsilon$时,它也有可能召唤一个分身。只要有一个分身到达了终点,那么他就胜利了,这个字符串就算是接收。如果怎么也走不到终点,或者他的所有分身都消失,那么字符串就算是拒绝。NFA中的非确定性,就是指在碰到某个字符后,下一步的状态是非确定,可能有多个的。

要严谨地定义NFA的话,只需要稍微修改一下开始DFA的定义。DFA的状态转移集合$\delta$的类型是$Q\times\Sigma_{\epsilon}\to P(Q)$。$\Sigma_{\epsilon}$是原字母表中附加一个空字符$\epsilon$,$P(Q)$是$Q$的幂集,也就是$Q$的所有子集的集合。换言之,在每个状态,碰到了某个字符或不碰到字符后,下一步得到的状态是一个集合,而不是单个状态。

配上图的话,理解NFA会更方便。但由于图不好放,这里就不贴了。随便翻开一本计算理论的课本,或者去网上搜索非确定有限自动机,都能找到一些很直观、易于理解的NFA运算过程图。

NFA是如此强大,它用了一种很赖皮的方式来走迷宫:我尝试当前字符串表示的所有路径,有一条能走出去,就算我能走出迷宫。而在原来的DFA上,我们只能按照规则,一边读字符,一边往前走一步。但令人惊奇的是,NFA和DFA是等价的,每一个NFA都可以找到一个和它对应的DFA。我们还是使用构造性发发来证明这个定理,也就是证明每个NFA都可以转换成DFA。DFA可以转换成NFA是显然的,因为根据定义,DFA就是NFA。

我们再次回顾一下NFA的定义。NFA转移函数的结果,从单个状态,变成了状态的集合。转移函数的结果集合,不是状态集合$Q$,而是幂集$P(Q)$。仔细一想,$P(Q)$的大小也只是$2^{|Q|}$个啊!它的大小是有限的。我们可以构造一个新的DFA,它有$2^{|Q|}$个状态,每个状态对应$P(Q)$的一个元素。要得到这个对应的话,只需考虑到计算机科学常用的二进制就行了。用一个$|Q|$位二进制数来表示子集,某一位是1就代表NFA中这一个状态里有一个分身。NFA是$Q$子集与子集之间的转移,但如果我们把子集看成单个元素的话,那么子集转移就等价于状态转移了。

上面这段话其实不是用来读的,是用来启发思考的。对于一个证明,往往要通过自己的思考来理解,看别人的证明思路一般是很难看懂的。相信看到了二进制,二进制每一位对应原NFA中的一个状态,你就能灵光一闪地想出整个证明过程来了。

现在,我们知道NFA和DFA等价,它们只是名字不同的同一事物罢了。DFA的有关性质,NFA都有。比如,能被NFA识别的语言都是正则语言

嗯?这说明我们以后可能通过NFA来证明语言的正则性了。正好,我们还差两个正则运算的定理没有证明呢!这两个定理可以用一些直观但不严谨的方法证明。

还是把自动机想象成走迷宫,我们将构造出一台NFA,能在新NFA中走到终点的字符串,就是运算过后的字符串。

先证明连接运算的封闭性,即证明有这样一个DFA或NFA,它可以识别连接后的字符串。但我们手里有的,只有识别之前两种语言的DFA或NFA。如果想用DFA来判定连接后的字符串的话,我们首先要面对一个问题:这个字符串应该在哪个位置拆成两半,使得前一半属于第一种语言,后一半属于第二种语言呢?但NFA就没有这个烦恼。先构造两个语言的DFA,取名为A,B。在机器A的所有接受状态中,连一条$\epsilon$的边到B的初始状态,表示任何一个到了A终点的人都可以影分身抵达B的起点。A的初始状态为新机器的初始状态,B的所有接受状态为新机器的接受状态。NFA没有在哪个位置拆开字符串的烦恼——反正我试遍所有可能就行了。

再证明星号运算的封闭性。根据上面的思路,我们可以构造一台NFA,它的所有接受状态,都连一条$\epsilon$的边到初始状态,表示我们在任何时候从可以尝试把这个字符串“拆开”。但是,星号运算中空字符串一定会被接受,因此我们要额外建立一个初始状态,它一个被接收的状态,有且仅有一条$\epsilon$的边指向原DFA的初始状态。

总算,我们证明了正则运算在正则语言下封闭,它们是一家人,有一样的姓也是理所当然的了。

正则表达式

有了正则运算这一新武器,我们有了一个新的表示语言的工具——正则表达式。如果用一个DFA来表示语言,每次都要画一张图,实在是太麻烦了。但是,使用正则表达式的话,我们只需要用一行文字就可以表达语言了。正则表达式和我们的数字表达式一样,运算符号就是并、连接、星号,“数字”就是空字符、空集和单个字符。数学表达式里有1+1=2,正则表达式里有0*0 = {至少有一个0且全部都是0的字符串}。

等等!DFA可以表示语言,正则表达式也可以表示语言,还没有人说它们是等价的呢!但正则表达式也姓正则,暗示它表达的语言就是正则语言,正则表达式等价于DFA。我们接下来又要用构造来证明这一定理。

唉,要是世界上有这样一台机器就好了。机器只有两个状态,第一个状态是初始状态,第二个状态是接收状态。第一个状态到第二个状态的转移条件是一个正则表达式。我要验证一个正则表达式的字符串,就等价于把字符串放入这个机器中判定。可惜我们的DFA只允许在转移条件上写字符啊!

数学家们向来比较开放,喜欢扩大定义。没有这样的机器,我们就来创造这样的机器。广义非确定有限自动机(GNFA)是一个可以把正则表达式当成转移条件的自动机。为了证明正则表达式和DFA的等价,我们只需要一步一步地把DFA的转移条件变成正则表达式,最后变成我们开始说过的那个梦想中的机器——只有两个状态和一行正则表达式的机器。如果能构造出这样的转换方法,就能证明定理了。

我们这次构造的GNFA,是在原来的DFA上一步一步改造过来的。改造的第一步,是新建两个状态——新起始状态$q_s$,新结束状态$q_e$。在保证正确性的前提下,$q_s$对其它每个点连边,每个点向$q_e$。由于正则表达式包括空集,所以连边总是可行的。我们新建的这两个状态就是保留到最后的状态。我们要删掉其它所有点,并修改转移条件,使得整个GNFA依然正确。

对于任意要删的点$q_{rip}$,对于任意其它点$q_i$,$q_j$,设$q_i\to^{R1} q_{rip},q_{rip}\to^{R2} q_{rip},q_{rip}\to^{R3} q_{j}, q_i\to^{R4} q_{j}$,则$q_i$到$q_j$的正则表达式修改为$(R1)(R2)\ast(R3)\bigcup(R4)$。千言万语,胜不过下面这张图:
Sample

按照这种方法,我们总能在保证正则表达式的意思不变的情况下把中间点删掉。总有一天,我们会删到只剩$q_s,q_e$两个点,这两个点靠一条正则表达式来转移。

我们说明了DFA可以转换成正则表达式。为了证明等价,我们还得证明正则表达式可以转换成DFA。不过这一步要比开始简单得多。

借助证明NFA与DFA等价的方法,我们用如下方法构造NFA,使之与正则表达式等价。一个字符或者是空字符,就连一条边。空集就不连边。并集就是一个点可以通过空字符$\epsilon$移动到并集所表示的两个部分。星号就是结束部分连$\epsilon$连回初始状态,同时再初始状态前加入一个接收状态。同样,说了这么多,不如放一张图:(声明:该图片来自《计算理论导引》(Michael Sipser)第三版 42页)

Sample

DFA和正则表达式互相转化,那它们肯定是等价的了。

非正则语言与泵引理

开始我们提过,正则语言是可以被DFA表达的语言。换言之,还有许许多多的语言无法用DFA表达。举一个经典的例子:设语言$B = \{0^n1^n|n\geq0\}$,也是说该语言表示0和1个数相同,且先出现0再出现1的字符串。仔细一想,你哪怕使出浑身解数,也构造不出识别这种语言的DFA——为了构造一个这样的DFA,我们必须用状态来存储0的数量,但0的数量可以是无穷大,而状态数的有限的。

但有些愣头青喜欢钻牛角尖,他们偏要说道:“我就不管!你构造不出一台DFA,万一别人构造出来了呢?你凭什么说世界上不存在一台DFA识别这种语言?!”

科学家们又想出了一种证明语言不是正则语言的方法,来应对这些“对知识刨根问底”的热心青年。这种方法用了一个引理,叫做泵引理

如果一种语言$A$是正则语言,那么$\exists (int)p,\forall(string)s\bigwedge|s|\geq p\bigwedge s\in A \Rightarrow \exists xyz = s \bigwedge xy^iz \in A (i \geq 0)\bigwedge |y| > 0\bigwedge |xy| \leq p$。这一行一阶逻辑,能够让人充分复习离散数学的知识。但是,你很可能看完了也看不懂这一行话要干什么,或者干脆跳过了这行话。

泵引理是说,如果A是正则语言,我们可以随便选一个泵长度p。对于A中的每一个很长很长,长度至少是p的字符串,我们都可以把它拆成3份xyz。其中y部分一定非空,且xy加起来很短很短,长度必须小于p。如果我们在中间不断插入y部分,也就是对于任意字符串$xy^iz,i\geq0$,这个字符串还是A的语言。这个引理中有很多存在和任意,需要仔细地多看几遍。

泵引理的正确性的证明十分诡异。我们需要证明,一个正则语言对应的DFA,它满足泵引理的条件。引理说我们可以随便选一个泵长度p,那么不妨令p为$|Q|$,也就是这个DFA的状态数。对于语言中任意一个长度大于等于p的字符串s,它经过的状态大于等于p+1,因为p+1个状态需要通过p次状态转移,也就是需要读取p个字符。s经过了p+1个状态,但我们总共就只有$p = |Q|$个状态啊!这是怎么回事呢?这说明我们至少经过某个状态两次。我在走迷宫的时候,两次走到了同一个地方说明什么?说明我绕路了!我绕了一圈,又返回了原地。既然我经过了某个状态两次,就必然存在一个状态序列,通过这个状态序列可以回到同一个状态。状态的转移需要读取字符,也就是说,读取了一些字符后,我们又回到了之前的某个状态,我们绕路了。这些字符,我把它重复若干遍,我还是会回到这个状态,我永远会在这个状态绕不出去了。这个字符序列,就是我们拆成三部分xyz中的y。在第一次碰到重复字符时,也就是第一次绕路结束时,我们至多有一个状态走了两次,其它每个状态走了一次,也就是说最多走了p+1个状态,也就是最多读取p个字符,这一部分就是xy,$|xy| \leq p$。至此,泵引理得到一个描述性的证明。

可以发现,泵引理中“存在一个p”这句话一点用也没有,因为我们在证明泵引理时,直接把p钦定为状态数|Q|了。

等等,泵引理有什么用啊?泵引理说明如果A是正则语言,则可以干嘛干嘛。我都知道正则语言可以通过构造一个DFA来判断了,我要泵引理干嘛?

泵引理给出的正则语言的必要条件。也就是说,泵引理的逆否命题,得到的是判断一个语言不是正则语言的充分条件。如果一个语言满足泵引理的逆否命题,那么很遗憾,这个语言一定不是正则语言。

或者再从另一个角度上来讲,我们可以通过反证法来证明一个语言不是正则语言。我们把这个语言套进泵引理中,发现它无论如何都会导出矛盾,那么就可以得出这个语言不是正则语言了。

泵引理反过来是这样说的:对于一个语言A,如果对于任意泵长度p,A中存在一个很长很长,长度至少是p的字符串,我们无论如何都不能把它拆成3份xyz,满足下列所有条件:y部分一定非空,且xy加起来很短很短,长度必须小于p。如果我们在中间不断插入y部分,也就是对于任意字符串$xy^iz,i\geq0$,这个字符串还是A的语言,那么A不是正则语言。我们把整个引理取反,所有的存在和任意都反了过来。这就是开始说要注意存在和任意的原因。现在我们也发现,泵长度p不是没有用的。正过来时,p可以任意取值;但反过来说时,p就是任意取值了。泵长度p在证明泵引理时是完全没有用途的,因为它是一个存在的条件,可以任意取值;而在说明一个语言不是正则语言的时候,就要考虑泵长度p为任意值的情况了。

就考试而言,证明一个语言不是正则语言都是一个套路。因为要求证明的命题一定是正确的,所以我们可以无脑套入泵引理来导出矛盾反证。以经典的例子$B = \{0^n1^n|n\geq0\}$来说明一下证明语言不是正则语言的方法。

证明:
假设$B = \{0^n1^n|n\geq0\}$是正则语言。令p为泵长度。随便选择一个字符串$0^p1^p$,它属于B,且长度大于等于p。泵引理保证,它可以拆成3部分xyz,且满足之前的三个条件。由于$|xy|\leq p$,所以$xy$肯定全部由0组成。$xy^iz(i\geq0)$也一定属于B。但很明显,假设我把y重复100遍,由于y非空,xyz的前半部分一定出现了许许多多的0.这样的字符串是不属于B的,矛盾产生。因此,我们假设错误,B不是正则语言。

如果能仔细品味泵引理背后涉及的数学、逻辑原理,我相信人们都会被数学之美、思维之美所打动。这种美看不见摸不着,却深深埋藏在每个人脑中。对于那些爱思考的人来说,发现这些美,就能收获到一种无比的喜悦。能成为一名热爱自己的学科,并能在自己的学科里有所成就的科学家,真是一件令人羡慕的事情。可惜不是人人都适合科研,不是人人都有机遇把自己的一生都献给知识荒原的开拓。

总结一下。这一章讲的是正则语言,也就是可以被有限自动机识别的语言。这一章介绍了许多在后面都会介绍的概念:机器、状态、接收与拒绝、非确定性等。这些概念,是研究后面一些更强大的计算机器的基础。从全书的角度来看,这一章为后续的知识学习奠定了最低层次的基础。同时,这一章也围绕有限自动机,讲了许多实际的内容。有限自动机本身是最低级的计算机器,通过学习它我们能隐隐体会到计算的本质。

从考试的角度来看的话,这一章难点只有两个。一个是DFA向正则表达式的转换,一个是泵引理的理解与证明。而其他部分都不难。自动机名字看上去结构复杂,实际上就是一个有向图,哪怕不知道什么是图,会走迷宫,就能学会有限自动机的概念。非确定有限自动机的概念比较难理解,但一旦理解,一旦迈过那一道坎后,这个概念就会显得十分清晰。所有的证明都不会考,学习它们对于考试毫无益处。


后记

我写这篇文章的本意和之前复习概率论一样,是站在讲述者的角度来叙述这些知识,更好地把知识梳理清楚。从结果来看,复习效果是有,但是我花费了大量的时间,而且很多时间都花在了吹牛和一些我已经熟悉的知识的讲解上。而且,考虑到时间不是很够的原因,文章本身的质量也不够高,讲得也不够清楚。也就是说,我花了很多时间,就稍微复习了一些不太熟的知识,文章也没有写好来,总体上来看非常亏。我得到的结论是,复习还是不要写文章的好。如果想清楚地介绍一个知识,就在闲暇的时候慢慢地写文章;如果要复习,就用效率更高的方式来看书。这样复习的话,竹篮打水一场空,什么都没得到。