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

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

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

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

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

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

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

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

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

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

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

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

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

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

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（指数生成函数）