这篇文章的内容主要是上课的时候写的，复习的时候我又修改了一些地方。

计算机组成原理复习：第三章 指令系统

指令对应机器语言，是计算机中最基本的命令。

本章先整体介绍了指令的结构，包括了操作码、地址码两个部分，并简要介绍两者所占长度的问题。之后，本章对地址码和操作分别进行了进一步介绍。先是介绍了地址码的具体知识，因为地址码可能不表示一个具体的地址，而是通过一些间接的方式来获取地址。特别地，有一部分存储区构成了栈结构，有一些操作专门是针对堆栈的。然后是一些具体的指令操作。最后，本着从普适性到具体性的思想，本章介绍了指令集中应该具有的操作，以及一些具体的指令集的发展过程。

指令格式

指令就是一段二进制代码。指令可以分成两部分，前一部分是操作码，后一部分是地址码。我觉得可以把操作码看成函数名，地址码看成函数参数。准确来讲，地址码是表示数据的代码，用地址码获取数据的方式再寻址技术部分会提到。

地址码结构

函数的参数可以有多有少。我们默认讨论双操作数系统。在指令这个“函数”中，最多支持4个参数。可惜的是，这个函数只支持默认参数来实现重载。也就是说，所有的指令实际上都会需要4个参数，但有些参数可能存在其它位置（可以看成一个全局变量），或者某些操作不需要所有4个参数。所以指令实际传入的操作数数量可能是0~4个。

指令传入的4个参数是：第一操作数地址，第二操作数地址，运算结果地址，下一条指令的地址。和C++的默认参数一样，每少输入一个参数，就代表最后的那个参数被设为了默认值。

3参数指令默认将下一条指令的地址设为当前指令地址+1。2参数指令把运算结果到第一操作数的地址。1参数指令把另一操作数存在某一个寄存器里。0参数指令一般是堆栈操作，把栈顶元素当成参数来进行操作。

操作码结构

操作码的长度可以固定，也可以不固定。后者可以节约空间，因为多地址的操作中地址码的长度较长，可以分配比较短的操作码。

显然，和哈夫曼编码一样，不能有一个操作码是另一个操作码的前缀。

寻址技术

寻址指获取操作数或指令的地址。实际上，这一节不仅讲了如何寻址，还先讲了如何编址。

编址

编址就是给寄存器、主存或IO设备中的存储单元编号，这样通过一个编号就能唯一地找到存在某处的一个数了。

编址最重要的是编址的最小单位。如果编址最小单位是字，那么在一个字长16位的系统里，你只能通过地址区分每16位数，而不能唯一确定数据某个字中的哪个地方；如果编址最小单位是字节，那么你就可以区分每8位的数字；如果编址最小单位是位，那么你可以精确地访问每一位了。

显然，编址最小单位是位的话，地址就会变得非常长了。还按字节编址比较好一些。在C语言中，我们可以得到每一个char的地址，这大概是字节编址的一个体现。

编址长度和主存的字长不是一个概念。

寻址

主要讲怎么寻找数据的地址。

最简单的是把数据直接放在指令里，这样不需要取通过地址来访问数据了。也就是往函数直接传了个值。但这样传进去的值的位数可能不够多。

还可以传数据的地址，也就是传指针。有时为了方便，还可以寻两次地址，也就是传指针的指针。数据可以在寄存器里找，也可以在主存里找，共$2\times 2 =4$种方法。这里还要多提一句，寄存器间接寻址指的是第一次在寄存器里找地址，第二次还是在主存里找。

此外，为了方便地获取一段连续的数据，还可以用基础地址加偏移量的方式获取数据。用C语言来写的话，就是

*(p + i)。p表示基础地址，i是偏移量。如果i在寄存器中，就叫做变址寻址；如果p在寄存器中，就叫做基址寻址。理论上这两种寻址方式表示的意思都是一样的，硬件实现方式也可能是相同的，但二者使用上有一些细微的差别。由于用户希望改变对于某个基址的偏移量i，因此变址寻址是一个面向用户的操作。

此外，在基址寻址的基础上，基址p可以来自程序计数器。这种寻址方式叫相对寻址。

指令类型

数据传输

数据是在寄存器和主存之间传输的，理论上存在4种不同起点和重点的传输指令。把寄存器看成0，主存看成1，有序对<0,0>,<0,1>,<1,0>,<1,1>就可以表示4种数据传输的指令了。

此外，还有对堆栈的PUSH和POP指令。但是，这两种指令本质上也是对于主存的操作，并不能称为最基本的操作。

除了单向传输，还可以交换两个位置的数据。

运算

包含C语言中的一些基本运算。比如加减乘除模、按位逻辑运算、移位运算。但书中没有具体讲这种指令的结构（比如有几个操作数）。

流程控制

　转移

转移有两种：有条件转移和无条件转移。有条件转移用于实现条件语句，而无条件转移让程序结构更加灵活，不必受到顺序存储的束缚。

转移都要定义转移到的地址这一参数。对于有条件转移，还要把比较条件作为参数。

子程序调用

一个程序通过一些指令，可以调用另一个程序。只是根据调用关系，把被调用的程序叫做子程序。

程序本身就是被调用的，这种调用其他程序的关系中体现了递归的思想。

子程序调用也可以让程序不按顺序执行，而是跑到另一个地方去执行。在这一点上，它和转移十分类似。为了搞清这两个概念，要知道这两个概念间的区别。

子程序调用与转移的区别

1. 子程序有“返回”这一操作。在调用了子程序后，一定要进行返回，返回到调用子程序的下一条指令。
2. 机器中应该有“程序”这个概念，也就是说不同程序之间是有明确的界限的。转移只能在当前程序中转移，而调用子程序会跳转到子程序的开头。

既然调用子程序的时候要返回，那么机器必须存储返回的地址。事实上，每调用一次子程序就会产生一条返回地址信息。有多种存储返回地址的方法。

返回地址的存储

1. 存到子程序的开头。我觉得这种方法很蠢，递归调用完全不可行了。但如果没有递归的话，这样设计节约了存储空间。
2. 先存到寄存器中，再让子程序把寄存器的内容存到内存中。这样不错。
3. 存到堆栈中。这种方法和2的本质是一样的。函数调用本身就类似栈的行为，用栈模拟的话更加自然，更加易于操作。但是，我觉得这种分类还是不够清楚。方法2里，你把东西存在内存里，在内存里也可以构成一个堆栈。方法3理论上是方法2的一个改进，而不是完全没有交集的两个方法。

返回指令

返回指令明明是和子程序调用配套出现的，不懂书上为什么要单独列出来。

输入输出指令

输入输出往往要用到除主存之外的存储区域。因此，IO指令可以分成两类。一种是IO地址单独编址，从0开始命名；另一种是把IO的地址接到主存后面。前者较独立编址，后者较统一编址。

我感觉这一小节涉及的内容还是比较少的，考试应该不会考。

指令集的发展

讲了一些很虚的东西，没什么价值。

计算机组成原理复习：第二章 数据表示

实际上，这章讲的是数据是怎么编码的。大部分时候，都是在讲二进制编码。本章编码不特别说明都是在指二进制下的编码。

数值

数值是运算中最常见的成分。

非得分类的话，数字可以分成整数和小数。但小数的表示方法中也是用到了整数的。因此本章先介绍整数怎么表示。

原码 补码 反码

无符号的数直接用二进制表示。

有符号的数因为有正负，所以编码的第一个任务就是把正数和负数区分开来。如果把最高位当成符号位，0表示正1表示负，后面当成无符号数表示，那么这种表述法就是原码表示。反码和补码都是对原码的负数表示进行更改。反码是把负数的每一位取反，补码是把原码的每一位取反加一。补码和反码的原理都是取模，把一个模意义下很大的正数视为一个负数。

其中，补码的特点是0只有一种表示方法，从而负数可以多表示一个数；原码和补码都有一个+0一个-0.

定点表示

定点可以表示小数和正数。

定点小数的第一位是符号位，小数点固定在第一位后面。也就是说，所有数是零点几。定点整数第一位也是符号位，小数点固定在最后一个数后面，和之前讲的整数表示完全相同。

定点小数根本不实用。但是它在后面会被利用到。

浮点表示

计算机中，小数都用浮点表示法，即小数点的位置可以移动。在这种表示法下，一个数的数位被拆成了两部分：一部分是指数，称为阶码；一部分是底数，称为尾数。

我写完了这篇文章，写完了忘了保存，直接上传到博客。我退出文章的时候提示是否要保存，由于文章上传到博客草稿文件消失，一定会提示是否保存，我就选择了丢弃。结果我辛辛苦苦写的内容全部消失了。复习时间太少，这篇文章今天晚上就不补了。以后有缘再补完这篇文章或者删掉这篇文章。气死我了！

期末考试又到了，又有一些课要复习。

我在上课的时候曾经做过一些“学习记录”，其形式是以自己的话来重述课本的内容，并且加上很强的主观评论。我复习的时候也采用这种形式。为了节约时间，有些吐槽会少一些，但这并不代表这门课没有喷点。

计算机组成原理复习：第一章 概论

这一章做为全书的总结，都是一些很泛的东西，没有太多有价值的知识。准确来说，都是一些具体知识而没有任何理论知识在内。我照着考试复习PPT把重点列一下。

计算机的组成部分

五个部分：

• 存储器

• 运算器

• 控制器

• 输入设备

• 输出设备

这些东西大一就记过了。

总线

特点：共享、分时

类别：地址、数据、控制总线

计算机系统组成部分

软件+硬件

计算机性能指标（部分）

时钟频率和时钟周期

可以类比游戏里的FPS和每帧所需时间。

吞吐量、响应时间

这个在操作系统里提过。吞吐量指单位时间完成的请求的数量，似乎没有一个严格的单位。响应时间指一个任务等待和最终执行所耗时间之和。

CPI

Cycles per Instructions，每条指令所花时间周期数

IPC就是倒数

MIPS MFLOPS

MIPS(Million Instructions per Second)是每秒执行多少百万条指令。等于主频乘IPC

MFLOPS不是指指令而是浮点运算。

数字图像处理小作业2

代码仓库

https://github.com/SingleZombie/WienerFilteringCpp

需求分析

老师给的需求如下：

• 读入一幅灰度图像
• 生成模糊图像：高斯模糊滤波+噪声
• 利用维纳滤波去模糊

总结来说，本次作业处理的图像对象是灰度图。灰度图要经过扭曲和修复两个过程。扭曲是高斯模糊和噪声（我选择了高斯噪声），修复用的是维纳滤波的方法。

相比上次的作业，本次需求清晰了不少。

技术学习

高斯模糊

模糊一个图像，可以通过卷积核来实现。具体来说，就是把一个像素的值加上它旁边几个像素的值。这样从直观上来看，像素之间的差异被减小了，图片就变模糊了。

问题是，一个像素值该如何加上旁边几个像素的值？每个像素相加的权重是什么？有人把二维高斯分布函数当成概率密度函数，权重就是概率密度函数的值。这样的利用高斯函数来计算卷积核权重的模糊方法就是高斯模糊。

高斯噪声

噪声就是把某个位置的像素值改变多少看成一个随机变量。如果随机变量符合高斯分布，那么这个噪声就是高斯噪声。随机变量和位置无关，只和随机变量自己的性质（均值、方差）有关。

维纳滤波

图像修复就是假设一个扭曲模型（输入为原图像的一个函数），之后设法求出这个模型中的参数，再利用这些参数复原图像。

维纳滤波就是假设原图像经过了一个卷积和一个噪声添加。该方法要的思想是让复原出的图像和原图像的”差异“尽可能小，这个”差异“被定义为均方误差。由于我没有学过信号处理这门课，而本课程又更关注算法的使用而非原理，这里不加证明地给出最终的式子：

其中$\hat{F}(u,v)$是复原出来的图像的频域值，$H(u,v)$是退化函数在频域下的表示，$G(u,v)$是被扭曲过的图像。而$k$是一个可调的参数。

本次任务的难点应该就是调试$k$这个参数。

有一篇博客对维纳滤波的原理以及matlab下的实现做了讲解。

OpenCV用法

写着写着发现又碰到了很多不会用的函数。

cv::normalize

cv::normalize(src, dest, upper_bound = 1, lower_bound = 0, norm_type = cv::NORM_L2);

normalize用于归一化，避免一个向量的模变大。在图像中，卷积核做归一化可以保证图像的亮度不变。

我只会用这个函数的前5个参数，有几个参数的名字是我取的。

src：原图像

dest：输出图像

upper_bound：归一化后的上界

lower_bound：归一化后的下界

norm_type：归一化的范式

归一化的基本操作就是要让模为1，因此输出的下界是0，上界是1。归一化的操作是把向量的每个分量都除以向量的模。

范式就是模是如何计算的。第一范式把模计算成所有分量之和，第二范式把模计算成所有分量平方和开方…………无穷范式把模计算成最大分量。

cv::filter2D

1

cv::filter2D(srcMat, destMat, ddepth, kernel, anchor = (-1, -1))

这个函数用于做图像的卷积，省去了手写循环的无聊过程。

同样，我只列出了我会的参数，参数名字是随意写的。

srcMat：原图像

destMat：输出图像

ddepth：输出图像深度，-1为和原图像相同

kernel：卷积核

anchor：锚点，也就是被卷积的点在卷积核的哪个位置。(-1, -1)是卷积核正中心

cv::dft

cv::dft(srcMat, destMat, flag)

这个函数就是做离散傅里叶变换，内部实现方法肯定是用FFT做的。

srcMat：原图像

destMat：输出图像

flag：操作的一些参数。我知道的参数有：

​ cv::DFT_CCOMPLEX_OUTPUT: 输出结果以复数形式存储，通道数增加。我输入单通道的图像出来是双通道的。

​ cv::DFT_REAL_OUTPUT：仅保留结果复数的实数部分，通道数不变。

​ cv::DFT_SCALE：运算结束后除以一个数，用于IFFT中

​ cv::DFT_INVERSE：做的是逆变换。

在我的程序里，FFT的写法是cv::dft(src, dest, cv::DFT_COMPLEX_OUTPUT);，IFFT的写法是cv::dft(src, dest, cv::DFT_SCALE | cv::DFT_INVERSE | cv::DFT_REAL_OUTPUT);

C++标准库用法

本次程序中我竟然碰到了完全没见过的c++标准库函数。

std::normal_distribution

1
2
3
4
5
6
7

#include <random>

std::random_device divice;
std::default_random_engine rng(divice());
std::normal_distribution<> normDis(mu, sigma);

double randomNumber = normDis(rng);

正态分布函数在c++的随机数库中。使用c++的随机数库有一些麻烦。语句前两行貌似是一个分布函数的初始化步骤，我不求甚解地直接使用了。语句第三行用平均值和方差初始化一个正态分布。初始化了之后每次就能获取一个符合正态分布的值了。

用OpenCV完成的组合操作

本节讲的是完成某一功能的，使用到OpenCV的自定义函数。

我在完成这个作业的时候，参考了很多别人的博客。但是有些人用的是matlab实现的，里面出现了C++OpenCV里没有的函数。我只好去搜了一些这些函数的原理，参考了别人的C++实现并自己实现了一下。

circShift

circShift(mat, x, y)

该函数的作用是，把一个矩阵mat向x轴正方向移动x个单位，向y轴正方向移动y个单位。移动中溢出的数字会填补到空出来的地方，也就是整个矩阵是循环移动的。

该函数实现起来很简单，先把x,y变成模意义下的值，再用临时矩阵做比较繁琐的copyTo就行了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

cv::Mat circShift(const cv::Mat& mat, int x, int y)
{
	int mWidth = mat.cols;
	int mHeight = mat.rows;
	x = (x % mWidth + mWidth) % mWidth;
	y = (y % mHeight + mHeight) % mHeight;

	cv::Mat res(mat.size(), mat.type());
	cv::Mat m1(mat, cv::Rect(0, 0, mWidth - x, mHeight));
	cv::Mat m2(mat, cv::Rect(mWidth - x, 0, x, mHeight));
	m1.copyTo(res(cv::Rect(x, 0, mWidth - x, mHeight)));
	m2.copyTo(res(cv::Rect(0, 0, x, mHeight)));
	
	cv::Mat res2(mat.size(), mat.type());
	m1 = cv::Mat(res, cv::Rect(0, 0, mWidth, mHeight - y));
	m2 = cv::Mat(res, cv::Rect(0, mHeight - y, mWidth, y));
	m1.copyTo(res2(cv::Rect(0, y, mWidth, mHeight - y)));
	m2.copyTo(res2(cv::Rect(0, 0, mWidth, y)));

	return res2;
}

由于我对OpenCV的用法不算很熟，没有写性能更好的方法，我感觉我的写法copyTo用得有点多了。我写的时候比较匆忙，只顾着实现，没有考虑到性能的问题。

psf2otf

psf2otf(mat, outSize = Size(-1, -1))

该函数的作用是，把一个二维点扩散函数转换成指定大小频域形式。

点扩散函数也许在数学上有其他定义，但根据我对该函数效果的理解，这里的点扩散函数指的就是一个用矩阵表示的二维函数，矩阵某位置的值表示二维函数在此处的函数值，该矩阵中心点是原点$(0, 0)$。

该函数的直接作用非常不好理解，但知道该函数的应用后就容易理解了。一个图像用一个卷积核卷积在时域上的计算可能比较慢，但是把图像和卷积核转换到到频域后，只要相乘，再转换回时域就可以得到卷积的结果了。由于卷积核的原点在中心点而不在左上角，且其大小可能比图像小，因此在卷积操作前需要把卷积核调整位置，并且扩大到和原图像一样的大小。经过处理后的卷积核的频域表示就能和原图像的频域表示进行正确的计算了。

该函数作用的宏观解释在明白函数的应用后，就比较容易理解了。但要理解其具体实现方法，需要了解FFT本身原理。我是通过理解一维FFT的原理后推广到二维的。

这里不加解释地直接给出该函数的具体操作了。该函数先把卷积核扩展至指定大小，再把卷积核的中心移动到左上角（准确来说是原点，计算机里左上角是原点）。注意这一步移动是循环移动，卷积核的部分内容被循环移动到了 右上角、左下角，右下角。扩展和移动结束后，再做FFT就行了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

cv::Mat psf2otf(const cv::Mat& mat, cv::Size outSize)
{
	if (outSize.height == -1 || outSize.width == -1)
	{
		outSize = mat.size();
	}

	cv::Mat otf;
	cv::Size paddleSize = outSize - mat.size();
	cv::copyMakeBorder(mat, otf,
		paddleSize.height - paddleSize.height / 2, paddleSize.height / 2,
		paddleSize.width - paddleSize.width / 2, paddleSize.width / 2,
		cv::BORDER_CONSTANT,
		cv::Scalar::all(0));

	otf = circShift(otf, -otf.size().width / 2, -otf.size().height / 2);
	otf = grayImageFFT(otf);
	
	return otf;
}

我在写这个函数的时候参考了这篇博客。这篇博客的实现 在卷积结束以后，还判断了虚数部分是否比较小，可以忽略。但我没有多加这一步判断。

图像与卷积核在频域的计算

这部分没什么新的内容，就是把代码组合了一下：

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

// Get frequent field of origin image
cv::Mat fOrigin = grayImageFFT(origin);

// Get frequent field of kernel and rearrange it
cv::Mat h = getKernal(..........);
cv::Mat funcH = psf2otf(h, mat.size());

// Execute complex number multiply
for (int i = 0; i < fOrigin.rows; i++)
{
	for (int j = 0; j < fOrigin.cols; j++)
	{
		std::complex<double> g(fOrigin.at<cv::Vec2f>(i, j)[0], fOrigin.at<cv::Vec2f>(i, j)[1]);
		std::complex<double> h(funcH.at<cv::Vec2f>(i, j)[0], funcH.at<cv::Vec2f>(i, j)[1]);

		g *= h;

		fOrigin.at<cv::Vec2f>(i, j)[0] = g.real();
		fOrigin.at<cv::Vec2f>(i, j)[1] = g.imag();
	}
}

// Do IFFT
cv::Mat matRes = grayImageIFFT(fOrigin);
cv::Mat result = cv::Mat(matRes, cv::Rect(0, 0, origin.cols, origin.rows));

以上代码算半个伪代码，完成了普通的频域相乘操作。如果要做其它的操作（比如维纳滤波），改循环里的公式就行了。

结构设计

同上次任务。

其中图像处理分成两个模块：退化图像生成与退化图像修复。

退化图像生成又分两个模块：模糊和加噪声。

最终实现的时候又加了一个交互模块，可以动态调参。这个模块要调用退化图像生成模块和退化图像修复模块。模块化设计的好处在这里就体现出来了。

虽然我感觉我模块分得特别好，最终代码清晰简明，符合数据流的内在逻辑结构，但我还是不太想在这里炫耀了。

程序设计心路历程

我看完高斯模糊和高斯噪声的原理后，很快就把这两部分实现了。

在了解原理，知道了C++有正态分布生成器后，这部分应该非常容易实现。

先提一句，最开始高斯模糊卷积我是直接拿filter2D做的。

问题就出在维纳滤波上。

一开始，我把原图像做了FFT，卷积核直接做了FFT再扩展到和原图像同样大小，再逐像素用公式计算，最后把计算结果IFFT。最终，我得到了一幅乱七八糟的图像。

频域的操作难以可视化，难以调试。我甚至都不知道我哪一步错了。但是，我凭借着我冷静的大脑，准备把这个过程逐一调试。

首先，我打算验证FFT和IFFT是否正确。我对一个图像进行FFT再IFFT回来，发现图像没有复原。果然，FFT和IFFT这步我就做错了。为了解决这个问题，我查了好多东西，还瞎照网上写了个FFT的优化函数。最后我理解了一下cv::dft的参数的含义，才把FFT和IFFT做对。

接着，我开始怀疑我的卷积核是不是不太对。FFT的原点在左上角，而卷积核的原点在中间，这好像不太对劲。但我一时半会儿解决不了这个问题。于是我干脆把卷积核的原点改到了左上角，即每次卷积的时候只对右下角的部分卷积，忽略掉这个问题，先解决其他问题。

我现在想验证除了卷积核外，我整个维纳滤波公式的运算是否有问题。我知道，如果没有噪声的话，逆滤波可以直接把图像复原，此时维纳滤波k值取0。因此，我暂时去掉了噪声，把k值设成0，看算法效果如何。结果图像确实复原出来了，但四周有一些不和谐的条状噪声。这个结果搞得我一头雾水。

我进一步思考，逆滤波能够复原图像，是因为噪声卷积操作和逆变换互为逆操作。但我卷积是拿filter2D做的，逆变换是在频域做的，这二者会不会工作原理不一样？于是，我暂时把卷积操作也变成频域操作。

终于，经过我这次的瞎折腾后，露娜sama美丽的面庞终于清晰地出现在了屏幕上。我又对退化图像加上了噪声。这次，k值取0的时候画面一片模糊，但k值只要稍稍增大，依旧能生成较退化图像更清晰的复原图像。至此，调试第一阶段结束，我高兴了好一会儿。

顺带一提，这个阶段为了调试方便，我把调k值放进了GUI里面，这样可以动态看到调参的结果。由于Debug模式下性能捉鸡，我又把OpenCV的Release版本编译了一遍。

问题是，现在我忽略了两个问题，一个是卷积核中心位置的问题，一个是用Filter2D的问题。我准备先解决第一个问题。

我把卷积核的原点又放回了中心。这次修改后，复原图像倒是能够生成了，但是图像的位置完全乱了：图像的左上角在右下角，右下角在左上角……。整个图像的四周都颠倒了。我去写了个OpenCV的fft2shift，试图把图像移正。但移正之后发现图像复原是失败的，最终图像和退化图像长得差不多。

在反复阅读了网上的代码后，我发现我代码唯一不太对的地方就是少了个psf2otf函数。我去搜了一下这个函数，竟然惊喜地搜到了这个函数的C++实现。我赶紧学了一下，看懂了函数的工作原理。我把这个函数加入了代码中，但怎么调都调不对。无奈之下，我只好双手离开键盘，选择冷静思考一波。

吃饭的时候，我仔细想了一想FFT的原理。我通过我以前对一维FFT的认识，以及昨天搜索到的零碎知识，渐渐搞清了psf2otf的原理以及目的。原来，图像FFT溢出的部分会“往回填”，图像靠右的部分表示的是负频率。为了保证卷积操作的正确进行，卷积核也要在扩充至和图像一样大小后，循环移动到左上角。

理清楚了每一条逻辑后，我修改了代码，在卷积核原点在中心时，图像终于能够正常地复原了。

最后是Filter2D的问题。如果卷积操作是拿这个函数做，得到的复原图像四周就会有条纹。我开始分析这种问题出现的原因。从“只有四周有问题”这个事实上，我想到很可能是Filter2D对边界处理的方式和用FFT处理不同。FFT是把整幅图像当成一个循环的图像来处理，因此会出现最右边像素受到最左边的像素的影响这种问题。而Filter2D对边界的处理有多种方式。我查了一下函数说明，果然Filter2D默认是把边缘反射处理，和在频域运算的方法不相同。但我尝试把边缘循环处理的时候，函数却报错了。源代码Assert掉了循环处理，不允许使用。无奈之下，我只好不使用Filter2D。反正我已经完全理解了这些操作的细节了。

最终程序的详细内容我就不讲了。感觉这份程序结构化做得非常棒，功能实现得非常完全。由于是作业，我还不得不加上了很多注释。相信大部分人都能看懂我那清晰易懂的代码。具体代码可以去代码仓库查看。

处理结果

原图

噪声图

修复结果

附原图的彩图版本：

感想

这个作业实在是太烦人了。原来我以为一个下午就能做完的东西，硬是做了两天。

我遇到的问题实在是有些令人恼火。我既不是程序代码不会写，也不是算法原理没搞懂，而是API怎么调，某个操作的实现细节这些边界知识没有掌握。本次作业没增加多少我对复原算法的认识，倒是进一步强化了我搜索信息和理解信息的能力。

之前，包括课本在内，很多介绍模糊和复原的文章都在用一个女性照片做为原图像。我一开始还没明白为什么。后来，在写程序碰到瓶颈之余，我好好分析了一下这个问题，总算得到了一个非常漂亮的结论：科研工作者大多为单身男性，每日面对的只有无情的代码和冷冰冰的文字。在这种恶劣的条件下，偶尔出现的女性照片，能给你黯淡的内心点上一丝救赎的光芒。所以，人们更多地选择那张女性照片。

可惜我的内心对于这种照片不为所动。我想，凭什么只用三次元的照片？于是，我随便在电脑里找了一张露娜sama的图片。果然，看着这幅CG被一步一步复原，我被各种BUG蹂躏得受伤的心灵得到了安慰。我为自己的机智和随机应变能力而赞叹不已。

数字图像处理小作业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函数非常简明，getInput和outputImage隐藏了输入输出细节，方便修改。事实上，图像读入时像素是用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的时候碰到了这个问题，写这个程序的时候就跳过了这个坑。

处理结果

原图

RGB风格转换结果

CMY风格转换结果

HSI风格转换结果

感想

写博客真浪费时间。

数独求解软件

Github链接

Sudoku Solver

PSP耗时表

具体任务划分

• 开发 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个数独。如果要更多的终局的话，只需要在求解的时候多求几组解就行了。

数独求解

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

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

待求解数独生成

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

设计实现过程

面向对象分析与设计

用例图

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

基本类图

基本的类图建立如上。

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

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

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

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

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

（以上两个任务共花费27分钟）

活动图

数独求解

（花费40分钟）

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

终局生成

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

（花费24分钟）

精化类图

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

数独矩阵

数独状态

求解算法

（以上三个任务共花费38分钟）

生成算法

(花费5分钟)

残局生成算法

(花费5分钟)

终局生成器

求解器

残局生成器

(共花费14分钟)

用户交互

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

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

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

（花费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)的。

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

残局生成器

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

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

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

求解器

（测试求解算法正确性花费30分钟）

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

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

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

(之后共花了100分钟)

集成测试

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

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

• 命令个数错误

• 命令错误（无‘-’）

• 命令错误（有’-‘，第二个及后面的字符不对）
• 命令参数错误

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

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

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

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

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

（花费56分钟）

代码复审

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

一开始，有很多警告：

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

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

（花费25分钟）

性能分析及程序优化

数独生成

性能分析结果

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

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

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

优化方法及结果

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

输出一次性用fwrite输出。

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

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

数独求解

初次结果

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

优化方法及结果

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

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

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

现在程序能90秒跑完了。

单元测试

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

最终程序代码说明

求解算法

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的操作。不管是怎样的区间操作、查询，所有modify和query都可以调同一个格式的函数，只要修改node中函数的实现细节就可以了。

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

markMult和preMove似乎是相关的两个变量，但多用一个变量也不会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”，心里就很不是滋味。我想的事情并不是写出好的文章，而是不希望文章中有漏洞，有缺点。我把写文章当成了一种自我展示、自我实现的工具。文章中有纰漏，文章写得不好，就会显得我这个人有缺点。为了让我自己显得更完美，我把这种对于完美的追求代入了文章中。我希望用一篇质量高的文章，来反映一个能力强的自己。我产生这些想法的原因和其他很多表现都是一样的。有些人在乎穿着打扮，认为外貌能反映自己；有些人在乎言行，以自谦来衬托自己的高尚。我虽然不在意别人的看法，但我也在用我写的文章，或者我的一些其他成果，来作为自己好与差的评判标准。对于事物的缺点的掩饰，不过是出于自卑，出于对自己的不足的掩饰罢了。

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

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

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

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

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

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

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

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

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