KUbuntu安装CIscoPacketTracer

注意：这是正版教程，需要你有Cisco账号。

第一步注册账号

先去思科官网注册账号：Cisco

可以先尝试这个链接，如果可以的话就跳过第二步，直接看第三步，如果链接失效了请继续第二步。PackeTracer

第二步下载PacketTracer

思科规定下载PacketTracer需要先免费注册任意一门课程，登陆好的页面如下：
image

image
等待跳转或者点击这个Skills for all
image

然后点击链接开始
image

image

image

进入课程后往下滚动
image

找到下载链接：
image

第三步选择需要的版本

选择对应版本并下载，安装。
image

Ubuntu||KUbuntu||Debain

如果出现缺少libgl1-mesa-glx依赖，并且apt 无法直接安装的时候，那么从官方源下载deb安装即可。libgl1-mesa-glx
选择对应版本，然后找到下面的链接，下载到本地用apt安装。记得换成你对应版本的包的名字。

1
sudo apt install ./libgl1-mesa-glx_22.3.6-1+deb12u1_amd64.deb

装好依赖后就可以正常安装CiscoPacketTracer了。

LinuxUbuntu

PVE学生自用记录

PVE记录

这篇博客主要记录自己大二阶段配置和使用PVE的过程。

什么是PVE

说到PVE，大家可能会想到Playsers Vs Environment，但是这里肯定不是指的游戏中的模式了，而是一个操作系统。

它的全称为：Proxmox VE，是一个运行虚拟机和容器的平台。基于 Debian Linux 完全开源。最大的灵活性，实施了两种虚拟化技术 （1）基于内核的虚拟机 (KVM) （2）基于容器的虚拟化 (LXC)。

我为什么会接触到这个东西？

事情还要从许的这篇论文说起，那天他在朋友圈发了论文的消息，然后我便迫不及待地把代码拉下来想要跑一下，结果就有了这个朋友圈：
image
是的，我这古老的笔记本无法支撑强大的Pytorch-3D编译所需要的内存，多次omm，于是我痛定思痛，想要换电脑。
但是我看了看我的钱包：
image
还是先将就着吧…………
然后！——————

成为垃圾佬

我在向longhao chen请教如何捡垃圾的时候，我发现了它：
image

分析一下：

  • CPU : E5 2680V4*2
    image
    • 一颗100 ，两颗200
  • 主板：Z10PA-D8
    image
    • 600左右
  • 内存条 ： 3200 * 16g *2条
    • 一条320 ，两条640
  • 电源
    • 240左右
  • 机箱
    • 50
  • 总计
    • 1730
  • 到手
    • 1030

要求不高，只要能点亮，不少配件就赚麻了。

激情下单购买配件

128G nvme 固态当缓存（但最后因为速度太慢放弃了）
image
5*500G sata机械盘组RAID5阵列

Raid5：至少需要3块硬盘
raid5优势：以上优势，raid5兼顾。任意N-1快硬盘都有完整的数据。
缺点：只允许单盘故障，一盘出现故障得尽快处理。有盘坏情况下，raid5 IO/CPU性能狂跌，此时性能烂到无以复加。
建议：盘不多，对数据安全性和性能提示都有要求，raid5是个不错选择

image

装机

配件完好，机箱像新的，用户手册非常详细。
image
外挂机械测试是否可正常组阵列。这里组阵列也是用到了PVE的软组，非常的方便。
image
Bios版本
image
CPU正确
image
操作系统安装
image
意外之喜：BMC
image

PVE

再通过zerotier进行路由，直接远程连接容器，垃圾佬也是用上自己的服务器了。
image

zerotier 路由设置

image

使用

longhao chen给我的建议是：不要在节点上拉屎，工作全部在容器中进行。
节点就是我的PVE服务器（因为PVE服务器可以存在多个一起组网），然后容器就是运行在PVE上的容器。
于是我当即立下创建了我人生中的第一个CT容器，取名为Main。开始了我垃圾佬的一生。
image

windows

PVE除了可以创建容器以外还可以创建虚拟机！而我刚好有有一些只能在windows上运行的软件，于是便可以通过远程桌面直接控制windows：
image

Ubuntu远程桌面

由于图形学的显示+计算需要，使用ubuntu的远程桌面貌似是最好的解决办法了，于是我找了一些方法，最后找到了这位大佬的帖子：
PVE下安装LXC创建桌面环境

LinuxPVE

第？课——基于矩阵快速幂的递推解法

第？课——基于矩阵快速幂的递推解法
由于中间的数论部分我自己学的很差，没有办法写出清晰的博客来，所以这里跳过了数论部分的博客，来到矩阵快速幂。

递推

递推是一个非常常用的工具。比如经典的斐波那契数列：

\[f(x)= \left\{ \begin{array}{**lr**} 1 &, 0\leq x\leq 1 \\ f(x-1)+f(x-2)&, 2 \leq x \\ \end{array} \right. \]

很明显，\(f(n)\)依赖于\(f(n-1)\)和\(f(n-2)\)，所以我们需要先计算\(f(n-1)f(n-2)\)才能计算\(f(n)\)。

假设我们现在需要求f(1e9+7)，你很快发现了，这个数字非常大。所以我们要求只需要结果对\(MOD\)取模就好了，而\(MOD=1e9+7\)。问题是：我们的迭代算法是\(O(n)\)的。那如何快速的求解这样一个递推问题呢？

矩阵与递推的联系

让我们站在巨人的肩膀上来看这个递推问题的第二项以及之后：

\[\begin{pmatrix}f(n)\\? \\\end{pmatrix}= \begin{pmatrix}1&1\\?&?\\\end{pmatrix}\dot\\ \begin{pmatrix}f(n-1)\\f(n-2)\end{pmatrix} \]

数学家们把这个问题用矩阵的形式表现了出来。但是矩阵上还有一行是空的。打个比方，我们在求\(f(2)\)的时候，矩阵可以写成：

\[\begin{pmatrix}f(2)\\? \\\end{pmatrix}= \begin{pmatrix}1&1\\?&?\\\end{pmatrix}\dot\\ \begin{pmatrix}f(1)\\f(0)\end{pmatrix} \]

那如果我们需要继续计算\(f(3)\)呢？

\[\begin{pmatrix}f(3)\\? \\\end{pmatrix}= \begin{pmatrix}1&1\\?&?\\\end{pmatrix}\dot\\ \begin{pmatrix}f(2)\\f(1)\end{pmatrix} \]

我们需要知道

\[\begin{pmatrix}f(2)\\f(1)\end{pmatrix} \]

可是\(f(2)f(1)\)哪里来呢？不妨把我们的通项改写一下，在计算\(f(n)\)的同时顺带计算\(f(n-1)\)？

\[\begin{pmatrix}f(n)\\f(n-1) \\\end{pmatrix}= \begin{pmatrix}1&1\\1&0\\\end{pmatrix}\dot\\ \begin{pmatrix}f(n-1)\\f(n-2)\end{pmatrix} \]

！！wow，我们得到了一个强大的新递推式子。这时候有人要不乐意了，这不是复杂化了么，我们本来只要两个数加一加，现在还要算一个2$\times$4的矩阵乘法？可是，我们的矩阵是常数。 再稍微改写一下这个式子：

\[\begin{pmatrix}f(n)\\f(n-1) \\\end{pmatrix}= \begin{pmatrix}1&1\\1&0\\\end{pmatrix}^{n-1} \\ \dot\\ \begin{pmatrix}f(1)\\f(0)\end{pmatrix} \]

相信你看到这里已经顿悟了，因为这个矩阵的高次幂我们可以大做文章！

快速幂与矩阵快速幂

快速幂

快速幂很简单，这里直接给出代码：

1
2
3
4
5
6
7
8
int quickpow(int a,int b,int x) {
	while(x>0) {
		if(x&1) a=a*b;
		b=b*b;
		x>>=1;
	}
	return a;
}

快速幂原理给一个友情链接：
快速幂总结

矩阵乘法

先放一个市面上常见的通用矩乘：

1
2
3
4
5
6
7
8
for(int i=1; i<=n; i++) {
	for(int j=1; j<=n; j++) {
		for(int k=1; k<=n; k++) {
			tmp.mat[i][j]+=(a.mat[i][k]%mod*b.mat[k][j]%mod)%mod;	\\拖慢速度
			tmp.mat[i][j]%=mod;
		}
	}
}

这个矩阵乘非常符合人的想法，但是有一个点让它比较慢，我们矩阵快速幂只需要进行幂次乘法，可以做一些优化：

1
2
3
4
5
6
7
8
for(int i=1; i<=n; i++) {
	for(int j=1; j<=n; j++) {
		for(int k=1; k<=n; k++) {
			tmp.mat[i][k]+=(a.mat[i][j]%mod*b.mat[j][k]%mod)%mod;
			tmp.mat[i][k]%=mod;
		}
	}
}

看出不一样的地方了吗？唯一的区别就在下标处。把k放在第一维会大大增加寻址的计算量，所以需要把k放在第二维上就能减少非常多的计算量。其实矩阵小的时候也没什么区别

矩阵快速幂

其实和普通快速幂没有多大区别，只需要重载一个乘法运算符就好了。这里给一份结构体加快速幂的代码：

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
#include <iostream>
#include <string.h>
using namespace std;
const int mod = (int)1e9+7;
struct MyMat{
    MyMat(int n_=24,int m_=24):n(n_),m(m_){
        memset(mat,0,sizeof(mat));
    }
    int mat[25][25];
    int m,n;
    friend MyMat operator*(const MyMat&a,const MyMat&b){
        MyMat tmp;
        if(a.m!=b.n)throw("Matrix size mismatch");
        for(int i=1; i<=a.n; i++) {
            for(int j=1; j<=b.m; j++) {
                for(int k=1; k<=a.m; k++) {
                    tmp.mat[i][j]+=(a.mat[i][k]%mod*b.mat[k][j]%mod)%mod;
                    tmp.mat[i][j]%=mod;
                }
            }
        }
        tmp.n = a.n,tmp.m=b.m;
        return tmp;
    }
    friend istream& operator>>(istream&is,MyMat& M){
        for(int i=1;i<=M.n;i++){
            for(int j=1;j<=M.m;j++){
                is>>M.mat[i][j];
            }
        }
        return is;
    }
    friend ostream& operator<<(ostream&os,const MyMat& M){
        for(int i=1;i<=M.n;i++){
            for(int j=1;j<=M.m;j++){
                if(j!=1)os<<" ";
                os<<M.mat[i][j];
            }
            os<<"\n";
        }
        return os;
    }
};
MyMat quickpow(MyMat a,MyMat b,int x) {
	while(x>0) {
		if(x&1) a=b*a;
		b=b*b;
		x>>=1;
	}
	return a;
}
int main(){
    MyMat a(2,3),b(3,2);
    cin>>a>>b;
    auto ans = a*b;
    cout<<"a:\n"<<a<<"*\nb:\n"<<b<<"=\n"<<ans;
    ans = quickpow(a,ans,5);
    cout<<"a*ans^5=\n"<<ans;
}

感兴趣的可以当一个参考。下面是本地运行测试输出

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
g++ -o test test.cpp
./test
1 2 1 2 1 2 
1 2 3 1 2 3
a:
1 2 1
2 1 2
*
b:
1 2
3 1
2 3
=
9 7
9 11
a*ans^5=
2480048 2480080 2480048
3188656 3188624 3188656

解题

斐波那契数列

让我们回到经典的斐波那契数列，写出这个矩阵递推式后，加上矩阵快速幂

\[\begin{pmatrix}f(n)\\f(n-1) \\\end{pmatrix}= \begin{pmatrix}1&1\\1&0\\\end{pmatrix}^{n-1} \\ \dot\\ \begin{pmatrix}f(1)\\f(0)\end{pmatrix} \]

我们就能快速地写出斐波那契数列的代码了，和上面不同的地方只有quickpow函数的x改为了long long：

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
MyMat quickpow(MyMat a,MyMat b,long long x) {
	while(x>0) {
		if(x&1) a=b*a;
		b=b*b;
		x>>=1;
	}
	return a;
}
int main(){
    long long n;
    cin>>n;
    if(n<2)cout<<"1"<<endl;
    else {
        MyMat A(2,2);
        A.mat[1][1]=1;
        A.mat[1][2]=1;
        A.mat[2][1]=1;
        MyMat x(2,1);
        x.mat[1][1]=1;
        x.mat[2][1]=1;
        auto ans = quickpow(x,A,n);
        cout<<ans.mat[1][1]<<endl;
    }
    return 0;
}

A Simple Math Problem

image
重点：
image
这个矩阵我就不写了，自己动手尝试一下吧～

Count

image
重点：
image

\[\left\{ \begin{array}{**lr**} f(n) = f(n-1)+2f(n-2)+n^3 \\n^3 = (n-1)^3-3n^3-3n-1\\n^2 = (n-1)^2+2n-1\\n = n-1 \end{array} \right. \]

\[\begin{pmatrix}f(n)\\ f(n-1) \\ n^3\\ n^2\\ n\\\ 1 \end{pmatrix}= \begin{pmatrix}1&2&1&3&3&1\\1&0&0&0&0&0\\0&0&1&3&3&1\\0&0&0&1&2&1\\0&0&0&0&1&1\\0&0&0&0&0&1 \end{pmatrix}\\ \dot\\ \begin{pmatrix}f(n-1)\\f(n-2)\\(n-1)^3\\(n-1)^2\\n-1\\1\end{pmatrix} \]

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
//AC代码
#include <iostream>
#include <string>
#include <sstream>
#include <string.h>
using namespace std;
const int mod = (int)123456789;
struct MyMat{
    MyMat(int n_=24,int m_=24):n(n_),m(m_){
        memset(mat,0,sizeof(mat));
    }
    int mat[25][25];
    int m,n;
    friend MyMat operator*(const MyMat&a,const MyMat&b){
        MyMat tmp;
        if(a.m!=b.n)throw("Matrix size mismatch");
        for(int i=1; i<=a.n; i++) {
            for(int j=1; j<=b.m; j++) {
                for(int k=1; k<=a.m; k++) {
                    tmp.mat[i][j]+=1ll*a.mat[i][k]*b.mat[k][j]%mod;		//1ll防止爆int
                    tmp.mat[i][j]%=mod;
                }
            }
        }
        tmp.n = a.n,tmp.m=b.m;
        return tmp;
    }
    friend istream& operator>>(istream&is,MyMat& M){
        for(int i=1;i<=M.n;i++){
            for(int j=1;j<=M.m;j++){
                is>>M.mat[i][j];
            }
        }
        return is;
    }
    friend ostream& operator<<(ostream&os,const MyMat& M){
        for(int i=1;i<=M.n;i++){
            for(int j=1;j<=M.m;j++){
                if(j!=1)os<<" ";
                os<<M.mat[i][j];
            }
            os<<"\n";
        }
        return os;
    }

};
MyMat quickpow(MyMat a,MyMat b,long long x) {
	while(x>0) {
		if(x&1) a=b*a;
		b=b*b;
		x>>=1;
	}
	return a;
}
int main(){
    std::string instr("1 2 1 3 3 1  1 0 0 0 0 0  0 0 1 3 3 1  0 0 0 1 2 1  0 0 0 0 1 1  0 0 0 0 0 1  2 1 8 4 2 1");
    std::stringstream is;
    is<<instr;
    MyMat A(6,6);
    MyMat x(6,1);
    is>>A>>x;
    long long n;
    cin>>n;
    for(int i=0;i<n;i++){
        long long m;
        cin>>m;
        if(m<=2)cout<<m<<endl;
        else {
            auto ans = quickpow(x,A,m-2);
            cout<<(ans.mat[1][1])%mod<<endl;
        }
    }
    return 0;
}

END

ubuntu无法进入桌面的一种情况

问题描述

系统环境：
image

  • 无法进入桌面
  • 可以进入锁屏
  • 输入密码后黑屏，并返回锁屏
  • tty能进入startx
  • startx中部分软件无法打开
    无法进入桌面最直接的错误，非常严重不可原谅。用户登陆输入密码，黑屏，然后回到用户登陆。

误导我的一个报错：

image
起初我认为是xdroid-server出毛病了，但是后来发现它无法启动，于是尝试使用startx来寻找错误。

tty如何进入

在锁屏页面使用Ctrl+Alt+F4，即可进入tty。
进入后输入用户名和用户密码。

进入startx进行debug

登陆好后使用startx进入。
后续操作使用tty执行startx展现。

大量的.desktop无法打开

执行code竟然出现Node.js的报错

image

打开code可以看到顶部的解释器选择，如果直接指定解释器，发现能够打开VSCode，说明不是可执行文件的错误。

image

指定解释器：

image

如果去掉code中的指定解释器行：

image

成功打开：

image

发现了Node.js

观察到VScode的解释器是/usr/bin/env sh

为什么不是/bin/sh或者/usr/bin/sh呢？

于是我尝试运行了/usr/bin/env发现！：

image

我天呢？这是什么道理？

问题解决

综上所述，我认为是/bin/env和/usr/bin/env被NodeJs覆写了，于是我决定删除/bin/env，然后重装coreutils问题解决：

1
2
sudo rm /bin/env
sudo apt reinstall coreutils

另外因为安装了gnome-tweaks和原来的某些东西有冲突，卸载后恢复原状：

1
2
sudo apt remove gnome-tweaks 
reboot

总结

这是一个非常奇怪的问题，因为我的/bin/env被Node.js覆写了，导致大量依赖于/bin/env寻找环境的批处理无法运行，这才导致了桌面无法打开。
debug期间我干了：

重装显卡驱动

卸载nodejs

删除所有与nodejs有关的文件

重装gnome

切换内核

修改用户权限

修改profile、~/.bashrc、/etc/environment……

LinuxUbuntu

第三课——最短路

最短路

最短路的算法是学了非常多次的一个算法了，但是这次学到的算是更加全面的一次。以前从数据结构中学到的两种最短路的算法，Dijkstra和Floyd。这两个算法在这篇文章中也会提到，最为最基础的两种最短路算法，后续的算法也都是在他们的基础上展开的。文章的最后，还提到了最短路的一个变种（故切算是？）差分约束。

Dijkstra

Dijkstra的想法其实有一种局部搜索的感觉，它做的事情是：“按照中转次数递增的顺序找最短路。”。和人的思想很像，看下面这张图：
image
假如我从\(C_{1}\)出发，那么人会先看到自己一步能走到的地方：
image
不管它最终要到哪里，都应该先走到最近的地方，也就是\(C_{3}\)，那么\(C_{1}\)到\(C_{3}\)的最短路就找到了，也就是\(edge_{1,3}\)。然后我们把\(C_{3}\)一步能到的点加入视线中：
image
这里我们通过中转可以到达\(C_{4}\)，但是发现从\(C_{3}\)中转再到\(C_{4}\)的距离大于到\(C_{2}\)的距离。所以我们把\(C_{2}\)加入候选集。同样的，\(C_{2}\)能到的点也要加入视线。
image
后续的操作我就用图片展示了：
image
image
这里有个点，提一下。可以发现当前\(edge{1,2,4}\)是比\(edge{1,3,4,6}\)短的，但是因为\(C_{4}\)已经存在最短路了，而且那条路一定比新的路更短，所以不需要考虑在候选集中的点的最短路了。接下来应该更新\(C_{6}\)：
image
到这里，从起点到所有点的最短路就都算出来了。这是人脑的求路过程，而Dijkstra也是非常相似的做法，这里给出代码描述：

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
/*数据定义*/
#define N 101  			//节点个数、
//inf 不能赋值为0x7fffffff，执行inf+inf时会溢出。
#define inf 0x3FFFFFFF
int map[N][N];
int via[N];
int dis[N];				//起点到某一个点的最短距离
/*初始化*/
for(int i=0;i<=n;i++){
	via[i] = 1;
	dis[i] = inf;
	for(int j=0;j<=n;j++)map[i][j]=inf;
}
/*算法部分*/
int start = 1, target = n;
int now = start;
dis[start]=0;
via[start]=0;
while(now!=target){
	int MIN = inf,next=now;
	for(int i=1;i<=n;i++){		//注意这里，下面要讲
		if(map[now][i]!=inf)	//没有这一句会上溢出，因为涉及了0x7FFFFFF+0x7FFFFFFF
			dis[i] = min(dis[i],map[now][i]+dis[now]);
		if(via[i]&&dis[i]<MIN){
			next = i;
			MIN = dis[i];
		}
	}
	if(MIN == inf)break;
	now = next;
	via[now]=0;
}
if(dis[target]==inf)printf("-1\n");
else printf("%d\n",dis[target]);

Dijkstra有一个非常明显的问题，通过比较与人脑的算法就能看出来，Dijkstra算法在判断下一个要进候选集的点的时候，需要遍历所有的点，这也使得复杂度来到了\(O(n^2)\)。想像一下，如果一个图比较稀疏，比如城市道路，可能有上千个十字路口，但是每个十字路口只有四条边左右。这个时候，如果也需要搜索每一个点的话，就会显得非常的多余。下面就是一个对比：
image
image
所以我可以使用链表来存储点的邻边，这样就能非常方便的访问各个点边了，但是有人会说了链表太慢了，我才不用，于是我们用数组来模拟链表——链式前向星。

链式前向星

这是一种数据结构，它用数组来存储邻接表。
邻接表我就不介绍了，这里直接介绍链式前向星：
image

图上模拟的就是内存中的数据存储。这个数据结构主要由两个表组成，第一个是head表，其中from作为下标而存在，实际上它只有一个值，就是next，表示从该下标对应的顶点出发的点存储位置。比如，我们上面提到的图中，\(C_{1}\)点有两条出边，分别是\(C_{2}，C_{3}\)，那么我们可以通过head[1].next找到edge表中对应的边。再看edge表，首先是一个下标index实际并不存储，然后是to它是边的另一个顶点的编号，再是value也就是该边的权值，最后是next它和head表中的next有着相同的属性，他们都是edge表的下标。通过edge[1].next我们可以找到\(C_{1}\)的下一条边。也就是\(edge_{1,3}\)。学过链表的同学想必已经理解了。没有学过的同学可以再思考一下。
最后是这个数据结构的一个基础实现，主要是一个插入操作。

1
2
3
4
5
6
7
8
9
10
11
12
13
#define N 100005		//顶点数
#define M 200005		//边数
long long head[N];
struct Edge{
    long long to,dis,next;
}edge[M];

void add_dege(long long from,long long to,long long value){
    edge[++cnt].to = to;
    edge[cnt].dis = value;
    edge[cnt].next = head[from];
    head[from] = cnt;
}

堆优化的Dijkstra

学会了链式前向星，我们来尝试一下思考如何用在Dijkstra上吧。
我们需要优化的点是：Dijkstra每次更新候选集的时候，需要遍历所有的点。

优先队列的设计

那么我们可以使用一个优先队列来维护我们当前访问到的点，并且排序的依据应该是该点到起始点的当前最短路长度。回到我们的人工计算最短路上来：
image
我们起初候选集里只有起点，但是我们能访问到的点有两个，并且这两个点都不在候选集里面。于是我们将他们连同他们的当前最短路加入优先队列中，下次从队列中取出最短的最短路就可以了。
我们可以设计一个放入队列中的数据包like this，其中，重载了小于运算符是为了让它在优先队列中能自动排序。因为我们需要构造最小堆，所以应该重载为a.dis<dis。因为优先队列是默认最大堆，而默认使用小于符号排序，所以我们需要重载反向的小于符号。

1
2
3
4
struct node{
    long long id,dis;
    bool operator <(const node&a)const{return a.dis<dis;}
};

点的重复

点的重复指的是一个点可能多次入队列，但是只要我们取出了某一个点，就能判断：后续出现的该点的距离一定比当前大，也就是后续的点都是多余的部分。 所以我们维护一个数组来标记顶点是否出队列，一旦出队列，就在后续的计算中不再接受该顶点。看代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
void Dijkstra(){
    priority_queue<node>q;
    q.push(node{s,0});
    for(long long i = 1;i<=n;i++){
        dis[i] = INF;
    }
    dis[s]=0;
    while(!q.empty()){
        node a = q.top();
        q.pop();
        long long now = a.id;
        if(vis[now])continue;		//一旦出队列，就在后续的计算中不再接受该顶点
        vis[now]=1;
        for(long long i = head[now];i;i=edge[i].next){
            long long j = edge[i].to;
            if(dis[now]+edge[i].dis<dis[j]){
                dis[j]=dis[now]+edge[i].dis;
                q.push(node{j,dis[j]});
            }
        }
    }
}

一个完整的例子

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
#include <iostream>
#include <queue>
#include <cstring>
using namespace std;

#define INF 0x3FFFFFFFFFFFFFFF
#define N 100005
#define M 200005

long long cnt,head[N];
long long n,m,s;
long long dis[N],vis[N];
struct Edge{
    long long to,dis,next;
}edge[M];
void add_dege(long long from,long long to,long long value){
    edge[++cnt].to = to;
    edge[cnt].dis = value;
    edge[cnt].next = head[from];
    head[from] = cnt;
}
struct node{
    long long id,dis;
    bool operator <(const node&a)const{return a.dis<dis;}
};
void Dijkstra(){
    priority_queue<node>q;
    q.push(node{s,0});
    for(long long i = 1;i<=n;i++){
        dis[i] = INF;
    }
    dis[s]=0;
    while(!q.empty()){
        node a = q.top();
        q.pop();
        long long now = a.id;
        if(vis[now])continue;
        vis[now]=1;
        for(long long i = head[now];i;i=edge[i].next){
            long long j = edge[i].to;
            if(dis[now]+edge[i].dis<dis[j]){
                dis[j]=dis[now]+edge[i].dis;
                q.push(node{j,dis[j]});
            }
        }
    }
}
int main(){
    while(scanf("%lld%lld%lld",&n,&m,&s)!=-1){
        //初始化
        cnt = 0;
        memset(head,0,sizeof(long long)*N);
        memset(dis,0,sizeof(long long)*N);
        memset(vis,0,sizeof(long long)*N);
        memset(edge,0,sizeof(Edge)*M);
        long long u,v,w;
        for(long long i=0;i<m;i++){
            scanf("%lld%lld%lld",&u,&v,&w);
            add_dege(u,v,w);
        }
        Dijkstra();
        for(long long i=1;i<=n;i++){
            i==1?1:printf(" ");
            printf("%lld",dis[i]);
        }
        printf("\n");
    }
}

Floyd

相比Dijkstra，Floyd真是再简单不过了，它使用的是一个非常清晰的动态规划的思想，使用三次循环直接算出任意两点之间的最短路。这里先给出核心代码：

1
2
3
4
for(k=0;k<n;k++)
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			 e[i][j] = min(e[i][j],e[i][k]+e[k][j]);

状态转移方程的含义就是从i点出发，借助k点到达j点是否会更快。

Bellman-Ford

Bellman-Ford优化了Floyd固定\(O(N^3)\)的时间复杂度，继承了Floyd的动态规划的思想，得到了这么一个可以控制最大中转次数的路径算法：

1
2
3
4
5
6
7
 for(int k=1;k<=n-1;k++){				//控制中转的次数
	for(int i=1;i<m;i++){				//遍历所有的边
		  if(dis[v[i]]>disp[u[i]]+w[i]){
			dis[v[i]]=dis[u[i]]+w[i];
		}
	}
}

其中v[i]是第i条边的入点，u[i]是第i条边的出点，w[i]是第i条边的权值。它的单点时间复杂度是\(O(N\times M )\)，可能会比Dijkstra还慢，但是优点是算法简单，而且可以解负权边。

如何解负权边呢？

Dijkstra不能解负权边的原因是遇到负环的时候会出现死循环，而Bellman-Ford可以控制中转次数。如果一条路的中转次数超过了n-1次那说明它一定经过了负权环。也就是说，只要我们再进行一次内部循环，如果dis数组发生了改变，说明存在负权环。下面是判断负权环的代码其实很简单

1
2
3
4
5
6
7
8
9
10
11
12
for(int k=1;k<=n-1;k++){
    for(int i=1;i<m;i++){
        if(dis[v[i]]>disp[u[i]]+w[i]){
            dis[v[i]]=dis[u[i]]+w[i];
        }
    }
}
for(int i=1;i<m;i++){
    if(dis[v[i]]>disp[u[i]]+w[i]){
        printf("存在负权环\n");
    }
}

SPFA

SPFA全称shortest path fast algorithm，也就是快速最短路算法，非常正宗的中国式英语。
SPFA尝试解决Bellman-Ford每次都需要遍历n-1条边的问题，用队列维护需要修改的点。具体做法如下：
首先把起点放入队列，然后松弛相邻的所有点，如果松弛成功而且点不在队列中，则把点加入队列。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
//spfa参考代码
void spfa(int u){
    queue<int>q;
    q.push(u); vis[u]=1;
    while(!q.empty()){
        int x = q.front();
        q.pop(); vis[x]=0;
        for(int i=first[x];i!=-1;i=next[i]){
            int y = v[i];
            if(dis[x]+w[i]<dis[y]){
                dis[y]=dis[x]+w[i];
                if(!vis[y])vis[y]=1,q.push(y);
            }
        }
    }
}

同样的，spfa也能够判断负环的存在，但是需要我们记录每次使用的点，只要某一个节点我们使用的次数超过了n-1次，那么也能说明存在负环。这里就不给代码了。

差分约束

差分约束是最短路的一种问题表达形式，本质上是同一个问题。
2024-04-05 15:13:11 星期五后续再更新差分约束

第二课——线段树

上一节课讲了树状数组，也介绍了树状数组的优点与不足，这里简单回顾一下。
**优点：**树状数组的代码非常简短，易于实现，被刘老师亲切的称为IO选手的”HelloWorld!“，就是因为代码短。
**缺点：**树状数组的缺点也非常的明显，只能处理单点修改区间查询或者区间修改单点查询的问题（以较高的效率）。而区间修改区间查询的问题没有办法很优雅的解决，于是引出了线段树。

线段树

先来看一个问题：

7-1 张煊的金箍棒（2）
张煊的金箍棒升级了！

升级后的金箍棒是由几根相同长度的金属棒连接而成（最开始都是铜棒，从1到N编号）；

张煊作为金箍棒的主人，可以对金箍棒施以任意的变换，每次变换操作就是将一段连续的金属棒（从X到Y编号）改为铜棒，银棒或金棒。

金箍棒的总价值计算为N个金属棒的价值总和。其中，每个铜棒价值为1；每个银棒价值为2；每个金棒价值为3。

现在，张煊想知道多次执行操作后的金箍棒总价值。
输入格式:

输入的第一行是测试数据的组数（不超过10个）。

对于每组测试数据，第一行包含一个整数N（1 <= N <= 100000），表示金箍棒有N节金属组成，第二行包含一个整数Q（0 <= Q <= 100,000），表示执行变换的操作次数。

接下来的Q行，每行包含三个整数X，Y，Z（1 <= X <= Y <= N，1 <= Z <= 3），它定义了一个操作：将从X到Y编号的金属棒变换为金属种类Z，其中Z = 1代表铜棒，Z = 2代表银棒，Z = 3代表金棒。
输出格式:

对于每组测试数据，请输出一个数字，表示操作后金箍棒的总价值。

每组数据输出一行。
输入样例:

1
2
3
4
5
1
10
2
1 5 2
5 9 3

输出样例：

1
24

可以看到题目中非常明显的区间修改+区间查询的意图，这也是线段树的一道入门题目。接下来我来介绍这个神奇的数据结构。

构成

线段树由一个四倍原数组长的数组组成，对于数组中的元素也有着特殊的含义，但是比起树状数组来说要好理解多了。
首先我们不从数组层面来看这个数据结构，而是从一个二叉树，一棵完全二叉树。假设我们的原数组有8个数。那么线段数和原数组的关系就像这样：
image
在线段树上的每一个节点表示对应区间的某个属性值，只要这个属性值满足区间的加法即可。
举个例子，这个属性值可以是区间和，可以是区间最值等等，这些具体的属性由题目来决定了，由于属性值的自由度极高，导致线段树在非常多的场合可以用于加速。

见过了线段树的二叉树形状，接下里给线段树一个数组的表示方式。这个也非常的简单：
image
和《数据结构》中一致，从根节点开始为 1 ，宽度优先搜索的顺序升序标号。有了标号，我们就能用数组来存储这棵二叉树了。可是为什么我们需要四倍的原数组空间呢？

这里我们从长度为5的数组开始，来探讨一下这个问题。
image
原数组长度为5，那么理论上黑色的节点已经够用了，但是我们使用的静态数组，一般会选择直接把完全二叉树所需的空间开出来，所以会用到最多四倍的空间。

左右子结点的访问

学过《数据结构》的读者可以跳过这块内容。
这部分比较简单，假设根的标志是1，那么左右子结点分别可以用以下两个函数访问：

1
2
3
4
5
6
int left(int d){
	return data[d<<1];
}
int right(int d){
	return data[d<<1|1];	//等价于 d * 2 + 1
}

稍微解释一下访问右节点的操作，一个二进制数在左移动后最低位一定是0，那么这时候可以用1与该数位或，就能得到乘2加一的效果。

树的初始化

首先定义一下线段树的结构（代码层面）

1
2
3
4
5
6
7
8
9
class SeqTree{
public:
	//方法定义
private:
	struct Data{
		int val;
	}data[N<<2];
}seqTree;			//线段树类
int arr[N];			//原数组

树的初始化是从原数组构造我们的线段树，以前面提到的题目为例子。节点的属性是区间和。

1
2
3
4
5
6
7
8
9
10
11
12
13
//arr为原数组、l为区间左边、r为区间右边、rt为线段树上的位置
void build(int*arr,int l,int r,int rt){
	if(l==r){							//到了叶子节点，直接赋值
		data[rt].val = arr[l-1];
	}
	int m = (l+r)>>1;					//寻找左右子结点的区间边界
	build(arr,l,m,rt<<1);				//递归构造两边的线段树
	build(arr,m+1,r,rt<<1|1);
	pushUp(rt);							//利用两边的子节点更新当前节点
}
inline void pushUp(int rt){
	data[rt].val = data[rt<<1].val + data[rt<<1|1].val;
}

可以发现线段树的构造是非常容易理解的。由于二分的存在它的复杂度也只是O(NlogN)。

单点修改

线段树的修改，相当于修改最下层的某个节点，它会影响到上层的非常多节点，依照树的初始化的想法，我们可以很容易的写出修改代码，这里不提供。

区间查询

首先有一个理论保障：线段树的每次查询不会超过O(logN)的复杂度。为什么呢？

  • 任一连续区间至多由\(2log_2^N\)个子区间组成
    • 原因：任一区间不在线段树同一层出现两个子区间，并且树高不超过\(logN\)
      • 原因的原因：因为区间连续，所以如果在同一层出现了两个子区间，那么这两个子区间一定可以合成上一层的一个区间。

所以查询的复杂度有了保障。于是我们来讲查询的思路。
对于一个区间查询\([L,R]\)，我们从根节点[0,4N]出发，进行二分查找，并把符合要求的区间上的节点都进行修改。Idea is pool show me the code!

1
2
3
4
5
6
7
//区间查询[R,L]
int query(int R,int L,int r,int l,int rt){
	if(R>l||L<r)return 0;
	if(R<=r&&L>=l)return data[rt].val;
	int m = (l+r)>>1;
	return query(R,L,r,m,rt<<1)+query(R,L,m+1,l,rt<<1|1);
}

是不是超级简单？哈哈，刘老师说：“当年我们没有人教，没有题目刷的时候，学会了线段树就开始大杀四方，当时觉得是很稀奇的东西。你们今天倒好，随便就能学到如此有意思的算法。”

区间修改

这个就厉害了，不仅实现了区间修改，还引入了最高效的偷懒方式——lazy
思路是这样的：我们修改一个区间的时候，如果要把值给到每个受影响的节点，会非常的麻烦，并且涉及到多次修改时，程序的复杂度会较高。但是仔细想想，我们线段树上的节点不是能代表属性么，那是不是也可以记录修改的属性呢？于是lazy诞生了。
修改一个区间的时候，我们不修改对应的叶子节点，而是在最上层的区间节点上记录本次修改，并在查询的时候应用。
我们首先修改一下数据结构体：

1
2
3
4
5
6
7
8
9
10
class SeqTree{
public:
	//方法定义
private:
	struct Data{
		int val;
		int lazy;	//lazy标志
	}data[N<<2];
}seqTree;			//线段树类
int arr[N];			//原数组

然后重写之前的各个方法：

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
void build(int*arr,int l,int r,int rt){
	if(l==r){
		data[rt].val = arr[l-1];
	}
	int m = (l+r)>>1;
	data[rt].lazy = 0;										//给lazy初始化值
	build(arr,l,m,rt<<1);
	build(arr,m+1,r,rt<<1|1);
	pushUp(rt);
}
//区间查询[R,L]
int query(int R,int L,int r,int l,int rt){
	if(R>l||L<r)return 0;
	if(R<=r&&L>=l)return data[rt].val;
	int m = (l+r)>>1;
	pushDown(rt,m-r+1,l-m);									//新加了一个应用lazy的函数
	return query(R,L,r,m,rt<<1)+query(R,L,m+1,l,rt<<1|1);
}
//区间修改把[R,L]修改为C
void update(int R,int L,int C,int r,int l,int rt){
	if(R>l||L<r){
		return;
	}
	if(R<=r&&L>=l){
		data[rt].val=C*(l-r+1);								//更新节点值
		if(r<l)
			data[rt].lazy=C;								//查询到此，继承lazy值
		return;
	}
	int m = (l+r)>>1;
	pushDown(rt,m-r+1,l-m);									//应用lazy
	update(R,L,C,r,m,rt<<1);
	update(R,L,C,m+1,l,rt<<1|1);
	pushUp(rt);												//这里有一个细节，应用lazy要在向上计算value之前
}

下面是应用lazy的函数的实现：

1
2
3
4
5
6
7
8
9
inline void pushDown(int rt,int rn,int ln){
	if(data[rt].lazy){
		data[rt<<1].val=data[rt].lazy*rn;
		data[rt<<1].lazy=data[rt].lazy;
		data[rt<<1|1].val=data[rt].lazy*ln;
		data[rt<<1|1].lazy=data[rt].lazy;
		data[rt].lazy=0;
	}
}

到这里就讲完了，线段树我似乎没有进行多少理论的分析，大部分都是show you the 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
#include <iostream>
#include <algorithm>
#include <vector>
#include <string>
#include <string.h>
using namespace std;

#define N 100000

class SeqTree{
public:
	inline void clear(int size){
		memset(data,0,sizeof(Data)*(size<<2));
		for(int i=1;i<=(size<<2);i++){
			data[i].val = 1;
		}
	}
	inline void pushUp(int rt){
		data[rt].val = data[rt<<1].val + data[rt<<1|1].val;
	}
	inline void pushDown(int rt,int rn,int ln){
		if(data[rt].lazy){
			data[rt<<1].val=data[rt].lazy*rn;
			data[rt<<1].lazy=data[rt].lazy;
			data[rt<<1|1].val=data[rt].lazy*ln;
			data[rt<<1|1].lazy=data[rt].lazy;
			data[rt].lazy=0;
		}
	}
	void build(int*arr,int l,int r,int rt){
		if(l==r){
			data[rt].val = arr[l-1];
		}
		int m = (l+r)>>1;
		data[rt].lazy = 0;
		build(arr,l,m,rt<<1);
		build(arr,m+1,r,rt<<1|1);
		pushUp(rt);
	}
	//区间查询[R,L]
	int query(int R,int L,int r,int l,int rt){
		if(R>l||L<r)return 0;
		if(R<=r&&L>=l)return data[rt].val;
		int m = (l+r)>>1;
		pushDown(rt,m-r+1,l-m);
		return query(R,L,r,m,rt<<1)+query(R,L,m+1,l,rt<<1|1);
	}
	void update(int R,int L,int C,int r,int l,int rt){
		if(R>l||L<r){
			return;
		}
		if(R<=r&&L>=l){
			data[rt].val=C*(l-r+1);
			if(r<l)
				data[rt].lazy=C;
			return;
		}
		int m = (l+r)>>1;
		pushDown(rt,m-r+1,l-m);
		update(R,L,C,r,m,rt<<1);
		update(R,L,C,m+1,l,rt<<1|1);
		pushUp(rt);
	}
	void debug(int size){
		cout<<"############## debug ##############\n";
		for(int i=1;i<=(size<<2);i++){
			cout<<data[i].val<<" "<<data[i].lazy<<"\n";
		}
		cout<<"############## debug ##############\n";

	}
private:
	struct Data{
		int val;
		int lazy;
	}data[N<<2];
}seqTree;

int main(){
	int b,n,q;
	cin>>b;
	while(b--){
		cin>>n>>q;
		seqTree.clear(n);
		int x,y,z;
		for(int i=0;i<q;i++){
			cin>>x>>y>>z;
			seqTree.update(x,y,z,1,n,1);
		}
		cout<<seqTree.query(1,n,1,n,1)<<"\n";
		// seqTree.debug(n);
	}
}

好好领悟线段树的节点属性吧。

第一课——树状数组

前缀和算法可以计算某一个区间的累记和，但是出现修改的时候，前缀和的效率便得不到保障。于是数状数组出现了。出现原因总结——需求从单纯的区间查询变为了单点修改 + 区间查询。

树状数组

本文不探讨树状数组的开发过程，这里先给出树状数组的结构：
image
树状数组的设计非常巧妙，它让下标为\(i\)（从1开始）的位置存储原数组的部分和。比如下标为2的位置，存储了前两个数据的和，而下标为4的位置存储了前四个数据的和。
但是也有些特殊的位置，比如6。你会发现，虽然它是偶数下标，但是它并没有存储前6个数据，而是只存储了5、6两个数据。下面给出树状数组的核心机制\(lowbit\)。

lowbit()

\(lowbit\)的定义是：一个二进制数的低位零的个数。
比如，2 的二进制是 10 ，那么 \(lowbit(2)\) = 1 。
而 6 的二进制是110，所以和2一样 \(lowbit(6)\) = 1。
于是我们可以给出树状树组下标为 \(i\) 的位置的定义：

1
2
TreeArr[i] = 0;
for(int j = 0; j <= lowbit(i) ; j++)TreeArr[i] +=data[i-j];

这里我直接写了C++的代码，但是阅读应该没有困难。其中data是原数组的数据，TreeArr是我们构造的树状数组。

lowbit()函数的实现

\(lowbit\)函数有两种实现方式，其中第一种是比较容易理解的：

1
2
3
int lowbit(int i){
	return x - (x & (x - 1));
}

第二种是比较抽象的，但是我个人比较喜欢，因为它更加的简洁优美。

1
2
3
int lowbit(int i){
	return i & -i;
}

这个有兴趣的朋友可以自己推导一下，不过也不会很复杂的。

update(int i,int data)函数

然后我来看一下单点更新如何在树状数组这样奇怪的数据结构上实现。
首先这个操作传入两个参数，一个是在原数组的位置，另一个是修改后的数值。再来看一下树状数组的结构：
image
假设我给data[3]加上一个1，那么树状数组中受到影响的节点有3、4、8，不难发现，我们可以从底部的3出发，自下而上的找出所有被影响的点，而4 = 3 + lowbit(3)、8 = 4 + lowbit(4)，是不是非常的妙？前推和后推都回到了\(lowbit\)上，不然怎么数\(lowbit\)是树状数组的核心机制呢？
有了这个理论基础，我们就能轻松的写出\(update\)函数的代码了。

1
2
3
4
5
6
int update(int i ,int data){
	while(i < n){				// n 是数组的长度
	TreeArr[i] += data;
	i += lowbit(i);
	}
}

不难看出这是一个**O(lgN)**的更新操作。

Sum(int j)函数

Sum(int j)函数实现了原数组中前j个数据的求和。
前面提到过，TreeArr[i]包括了从i开始的往前数lowbit(i)个数据的和，那么在求前 j个数据的和时，我们可以利用和更新中类似的方法，每次计算当前lowbit(j)个数据的和，然后去到j \- lowbit(j)的位置继续求前面的值。代码如下，也是非常的简洁

1
2
3
4
5
6
7
8
int Sum(int j){
	int rest = 0;
	while(j > 0){
		rest += TreeArr[j];
		j -= lowbit(j);
	}
	return rest;
}

不难看出这是一个**O(lgN)**的求和。

应用：单点更新区间查询

P3374 【模板】树状数组 1
https://www.luogu.com.cn/problem/P3374

区间更新 + 单点查询

这时，我们的需求改变了，我们不再需要区间查询了，而只要单点查询，但是需要实现区间修改。这时我们可以使用到一个数学概念——差分。使用差分作为树状数组存储的内容，可以让树状树组从单点修改 + 区间查询变为区间修改 + 单点查询。

差分的定义

假设d[n]是一个差分数组，那么：

\[d[n] = data[n] - data[n-1] \]

非常好理解，如果我要修改全部的数据，比如把所有数据加 1，那么我们只需要在第一个位置加上一就好了，因为差分数组的性质，其他的位置值并不需要修改。
那么如果我们要进行单点的查询，比如 data[n]（原数据），就需要计算前d数组的前n项和。这一就完美地完成了从单点修改 + 区间查询变为区间修改 + 单点查询的过程。

树状数组的不足

当我们的问题变成区间修改 + 区间查询时，树状树组便不能完成这个工作了（一维），我们需要更好的数据结构，下节课——线段树完美解决树状数组的问题。
PS：其实并不是，树状树组的内存是(N)，而线段树需要(N<<2)也就是原数组四倍的内存空间。

B-Spline

B-Spline

注：
本文只是作者自己对B样条的理解和实现，如有错误欢迎指出，谢谢。
作者邮箱：zhywyt@hdu.edu.cn
B-Spline 是图形学中非常基础的一个曲线，我的导师徐岗徐老师，让我们从 B样条的绘制开始，无疑是一个很好的开始。

同 Bezier 曲线，B样条也是一个由伯恩斯坦基函数加权得到的曲线，因此它和 Bezier 曲线有这很多的相似点，这里我们把 Bezier 曲线和 B 样条放在一起来研究。

贝塞尔曲线：

\[\vec{P(t)}=\sum_{i=0}^{n-1}\vec{p_i}B_{i,n}(t),t \in [0, 1] \]

B 样条曲线：

\[\vec{P(t)}=\sum_{i=0}^{n-1}\vec{p_i}B_{i,d}(t)，其中 t_{min}\leq t \leq t_{max}，2\leq d \leq n \]

可以看到两个曲线的区别就在于基函数的定义：Bezier 曲线的基函数只与曲线的控制点数有关，曲线的控制点数就确定了曲线的阶数；而 B 样条的基函数中，n 换成了 d ，这也意味着曲线的控制点数不再决定曲线的阶数，而是独立的一个自由度。

下面是我的一些符号规定：

\[\begin{array} {cll} \hline 参数名 & 符号 \\ \hline 参数 & t\\ 顶点数 & n+1 \\ 曲线阶数 & p \ \&\ d \\ 控制顶点 & P_i \\ 下标为\ i \ 的节点向量 & u_i \\ 基函数的次数 & p-1\\ 第\ i+1\ 个控制顶点的\ j\ 阶基函数 & B(i,j)\\ 头节点重数 & h \\ 节点向量个数 & n+p \\ 自由度 & n,p,h\\ \hline \end{array} \]

解释：曲线阶数在数学中我一般使用\(d\)，而在代码中我喜欢使用\(p\)

为什么要修改基函数？

众所周知，Bezier 曲线的一大缺陷就是曲线的定义是全局的，移动一个控制点就会影响整条曲线，但是设计上往往希望能够局部修改。于是产生了分段 Bezier 曲线，但是新的问题出现了：分段 Bezier 曲线的连续性很难保证，不能达到很好的效果。B 样条因此诞生，以修改基函数的方式达到目的。

为什么修改基函数可以达到这个效果？

B 样条曲线：

\[\vec{P(t)}=\sum_{i=0}^{n-1}\vec{p_i}B_{i,d}(t)，其中 t_{min}\leq t \leq t_{max}，2\leq d \leq n \]

再来理解一下曲线的定义，每个控制点，在参数 t 处的对曲线影响的权值由基函数控制。要控制每个点的影响范围，就需要定义基函数的非零区间，在非零区间里，该控制点就对曲线有影响，而在零区间，就没有任何影响。这样的话，就能控制每个点对曲线影响的范围了。
在此之前，先介绍一下节点向量的概念。

B样条的节点向量

节点向量描述了B样条曲线上某一参数 t对应基函数的哪一部分。你可以把它想象成一部电影的时间线，给定一个时间就能得到相应的画面，而B样条就是给定一个参数 t就能得到对应的坐标。
下面是一个例子：

  • \(P_i:\) 第 i 个控制点
  • \(u_i:\) 下标为 i 的节点向量的值
    下图是一个有五个控制顶点的(\(n+1=5\))，四阶B样条(\(p=4\))。\([u_0 , u_8]\)就是节点向量了，而上面画的曲线则是某个控制点对下面的区间的参数的权值。给定一个节点向量，就能唯一确定一个基函数的对应关系，比如最简单的，\(u_0\)对应参数的点，只被\(P_1\)这个点所影响，而\(u_3\)却被\(P_1，P_2，P_3，P_4\)四个点所影响。当然，你会发现区间\([u_0 , u_3]\)其实是被小于曲线阶数的顶点所影响的，它的权值和不等于一，这不满足B样条曲线的性质，实际上B样条的定义不从\(u_0\)开始，\(u_0\)的存在仅仅是为了其他区间的计算权值。这也就反映了B样条参数的定义域了。不是\([0,1]\)而是根据节点向量的\([u_{p-1},u_{n+1}]\)在下面这个例子中是区间\([u_3 , u_5]\)。
    image
    可以很清楚的看到，\(P_1\)控制点不会影响到\(u_5\)及以后的曲线。每一段曲线比如\([u_3, u_4]\)只由\(P_1 - P_4\)四个控制点所影响，实现了曲线的局部控制。

这样的基函数如何定义？

从一个简单的一阶基函数定义看起：

\[B_{k,1}(u) \begin{cases} 1 , u[k]\leq u\leq u[k+1]\\ 0 , else \end{cases} \]

它只有在自己的定义区间内才有值 1，其余部分都为 0。
画出来也非常的清晰：
image
那么如果要进行升阶，我们就需要使用两个定义在相邻区间的一阶基函数：
image
通过一个简单的线性组合就能得到二阶的基函数：
image
可以看到这个基函数的定义域成为了两个低阶基函数定义域的并集，从定义上看这个线性组合就是：
image
先不要害怕两个奇怪的权值，观察上面的一阶到二阶的加权不难发现，这两个权值其实是曲线在参数为 u 处的简单加权。也就是上面出现过的一次加权直线的函数取值。但是你可能会觉得奇怪，当\(u\leq u_{k+1}\)的时候，第二项的权值好像大于一了？！这不对吧？

实际上在\(u\leq u_{k+1}\)的时候，右边这一项\(B_{k+1,d-1}(u)\)的值等于零了。

看下面这个图，计算三阶基函数会显得更清晰：
image
于是我们得到完整的 B 样条数学定义：

\[\begin{cases} \begin{cases} 1 , u_k \leq u\leq u_{k+1}\\ 0 , else \end{cases} \quad ,k = 1 \\ B_{k,d}(u)=\frac{u-u_k}{u_{k+d-1}-u_k}B_{k,d-1}(u)+\frac{u_{k+d}-u}{u_{k+d}-u_{k+1}}B_{k+1,d-1}(u) \quad ,k \ge 2 \end{cases} \]

B样条的自由度

说到 B 样条，就不得不说说它的自由度了。相比于 Bezier 曲线仅有的一个自由度（顶点数），均匀B 样条的自由度可多达三个（个人理解）：顶点数、曲线阶数、头节点重数设置（非均匀B样条修改为节点向量的设置即可）。姑且允许我这么随便的称呼自由度，可能不是那么标准，但是我是想让你明白 B 样条需要设置哪些参数。这里，我们讨论的是（准）均匀节点的 B 样条，也就是节点向量只会在头和尾相等而在中间以相同速率增加。比如：

1
2
3
4
5
6
7
8
9
10
11
12
13
n+1	: 5			number of control points.
p	: 4			rank of BSpline.
h	: 1			repeat of head and tail.
nv	: 9			number of node vector.
head vector:
0 0.111111 0.222222 0.333333 0.444444 0.555556 0.666667 0.777778 1 
//change para of B-Spline
n+1	: 5			number of control points.
p	: 3			rank of BSpline.
h	: 3			repeat of head and tail.
nv	: 8			number of node vector.
head vector:
0 0 0 0.25 0.5 1 1 1

相信你们顶点数和曲线阶数一定都能理解意思，那么我来解释一下头节点重数的定义：

头节点重数的作用

头节点重数： 头节点重数是用来描述节点向量头和尾的重合个数的。首先，一个节点向量是单调不减的，而我们讨论的又是均匀节点的 B 样条曲线，所以头节点重数就可以完整的控制节点向量的生成。

而节点向量就能定义控制点对曲线的影响。

下面举两个直观的例子来带你了解头节点重数的作用。
首先我们绘制一个正常的 B 样条曲线，设置它的头节点重数为一：
image
然后我们修改它的头节点重数，使得头节点重数等于基函数阶数。
image
会发现，头节点重数可以使得曲线插值端点。下面我用两个图来从基函数的层面带你理解头节点重数为什么能够影响曲线。首先是一张有五个控制顶点的4阶B样条的基函数关系图。
image
如果现在让头和尾重合一部分。比如 \(h = p\) 也就是头节点重数等于曲线的阶数。节点向量就会变成下面这样：

\[\begin{cases} u_0 = u_1 = u_2 = u_3 = 0 \\ u_8 = u_ 7 = u_6 = u_5 = 1 \\ u_4 = 0.5 \\ \end{cases} \]

让我们把节点向量的值写到基函数的关系图上去。红色部分的基函数由于区间为零，值恒为零。所以可以看到不为零的基函数构成了伯恩斯坦基函数，也就是Bezier曲线同款基函数。
image

B样条退化为Bezier曲线

再想象一下，如果控制点的个数减少一个，也就是 \(n+1 = p\) ，控制点的个数等于曲线阶数，并且满足头节点重数等于曲线阶数的时候，中间非零的节点区间只剩下了一个，这个时候，B样条真正的退化为了Bezier曲线。总结一下，B样条退化为Bezier曲线的条件：

  • 1.控制顶点个数等于曲线阶数 \(\quad\) \(n+1 = p\)
  • 2.头节点重数等于曲线阶数 \(\quad\) \(h = p\)

B样条的实现

这部分其实不难，如果完全理解了B样条真的是非常简单的实现，但是如果不清楚自己是否完全理解了B样条的思想，可以动手实现一下，这里提供实现思路和源代码。这个项目中我用到了几个绘制算法，你无需理解他们是如何工作的，只需要像调用函数一样调用它即可。如果想要使用我的框架进行B样条曲线的学习，请确保自己有看懂C++简单代码的能力，以及稍微阅读一下函数接口的使用。这篇文章会介绍B样条曲线实现用到的所有接口。如果对算法感兴趣的同学也可以去看一下这几个绘制算法的具体实现。
直线的绘制\(\quad\)圆形的绘制\(\quad\)橡皮筋技术

准备工作

请先下载项目源代码，并保证根据README的指导可以让程序跑起来。这个项目提供了linux和windows两个操作系统的编译方案。windows使用VS进行编译，linux使用cmake进行编译。
B-Spline
其中CMakeLists.txt文件的内容如下：

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
#最低版本
cmake_minimum_required(VERSION 2.8)

#项目信息
project(BSpline)

#查找当前目录下所有源文件
#并将文件名保存在DIR_SOURECE
aux_source_directory(. DIR_SOURCE)
add_executable(${PROJECT_NAME} ${DIR_SOURCE})

find_package(GLUT REQUIRED)
find_package(OpenGL REQUIRED)

if(GLUT_FOUND)
    include_directories(${GLUT_INCLUDE_DIR})
    target_link_libraries(${PROJECT_NAME} ${GLUT_LIBRARIES})
else(GLUT_FOUND)
    message(FATAL_ERROR ”GLUT library not found”)
endif(GLUT_FOUND)
if(OPENGL_FOUND)
    include_directories(${OPENGL_INCLUDE_DIR})
    target_link_libraries(${PROJECT_NAME} ${OPENGL_LIBRARIES})
else(OPENGL_FOUND)
    message(FATAL_ERROR ”OpenGL library not found”)
endif(OPENGL_FOUND)

#指定生成目标

b_spline_basis

这将是我们要实现的第一个函数，b_spline_basis，函数原型如下

1
GLdouble b_spline_basis(int i, int p, GLdouble u, const vector<GLdouble> &nodeVec)

接受 （i，p，u），返回下标为i的控制点的p阶基函数在参数u处的函数值。
首先编写一阶的基函数。

1
2
3
4
5
6
if (p == 1){
	if (nodeVec[i] <= u && u <= nodeVec[i + 1])
		return 1;
	else
		return 0;
}

在一阶之外的基函数，首先计算两个权值的分母：

1
2
GLdouble len1 = nodeVec[i + p - 1] - nodeVec[i];
GLdouble len2 = nodeVec[i + p] - nodeVec[i + 1];

然后利用递归计算两个区间的子基函数

1
2
3
4
GLdouble alpha, beta;
alpha = (u - nodeVec[i]) / len1;
beta = (nodeVec[i + p] - u) / len2;
return (alpha)*b_spline_basis(i, p - 1, u, nodeVec) + (beta)*b_spline_basis(i + 1, p - 1, u, nodeVec);

由于出现了对len1和len2的除法，所以需要检测区间的长度，避免浮点误差导致的离群点。而且u的区间也不能超过基函数的定义区间，需要把上面代码修改为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
if (u < nodeVec[i] || u > nodeVec[i + p])			//限制参数区间
	return 0;
GLdouble len1 = nodeVec[i + p - 1] - nodeVec[i];
GLdouble len2 = nodeVec[i + p] - nodeVec[i + 1];
GLdouble alpha, beta;
if (len1 <= 1e-10){									//避免浮点误差导致的错误
	alpha = 0;
}
else{
	alpha = (u - nodeVec[i]) / len1;
}
if (len2 <= 1e-10){
	beta = 0;
}
else{
	beta = (nodeVec[i + p] - u) / len2;
}
return (alpha)*b_spline_basis(i, p - 1, u, nodeVec) + (beta)*b_spline_basis(i + 1, p - 1, u, nodeVec);

然后就能得到完整的基函数计算方法了。下面是完整代码的展示。

b_spline_basis

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
GLdouble b_spline_basis(int i, int p, GLdouble u, const vector<GLdouble> &nodeVec){
	if (p == 1){
		if (nodeVec[i] <= u && u <= nodeVec[i + 1])
			return 1;
		else
			return 0;
	}
	if (u < nodeVec[i] || u > nodeVec[i + p])
		return 0;
	GLdouble len1 = nodeVec[i + p - 1] - nodeVec[i];
	GLdouble len2 = nodeVec[i + p] - nodeVec[i + 1];
	GLdouble alpha, beta;
	if (len1 <= 1e-10){
		alpha = 0;
	}
	 else{
		  alpha = (u - nodeVec[i]) / len1;
	}
	if (len2 <= 1e-10){
		beta = 0;
	}
	else{
		beta = (nodeVec[i + p] - u) / len2;
	}
	return (alpha)*b_spline_basis(i, p - 1, u, nodeVec) + (beta)*b_spline_basis(i + 1, p - 1, u, nodeVec);
}

draw

然后可以开始实现绘制算法了。
绘制算法首先要计算参数对应点的坐标，这个坐标等于每一个控制点在参数u处的贡献。下面代码很好理解，我就只提一点，GLUT的使用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
GLdouble basis;
vector<GLdouble> rx(node, 0), ry(node, 0);
GLdouble u_begin = m_nodeVec[p-1];					//注意开始参数与结尾参数
GLdouble u_end = m_nodeVec[m_n+1];
GLdouble dis = (u_end - u_begin) / node;
GLdouble u = u_begin;
for (int i = 0; i < node; i++) {
	 for (int j = 0; j <= m_n; j++){
		basis = b_spline_basis(j, p, u, m_nodeVec);
		rx[i] += m_points[j].x * basis;
		ry[i] += m_points[j].y * basis;
	}
u += dis;
}
glBegin(GL_POINTS);
for (int i = 0; i < node; i++){
	glVertex2d(rx[i], ry[i]);
}
glEnd();
glFlush();

本项目中，只使用了GLUT的点绘制功能，虽然它内置的算法大多比我写的高效，但是毕竟是学习目的，了解绘制算法的底层对自己没有坏处。glVertex2d(int x,int y)可以在屏幕上绘制出一个在(x, y)的点，非常的简单。
好了就讲这么多了，其余的还请大家自己看源代码了，核心代码就这么多了，下次见。白白

后言

如果有小伙伴发现了本文的错误，请您第一时间联系我，我会感激不尽。
作者邮箱：zhywyt@hdu.edu.cn
或者：zhywyt_6.10@qq.com
又或者您可以直接到我的gitee项目中提出issue，我有时间一定会继续更新代码的。

CG

花期内的花的数目

2251.花期内的花的数目
看到题目的第一想法是桶排序，但是想想肯定会超时，在题解区看到了这么一种解法，感觉很有意思，就记录一下。
要统计某一时间内多少花开放，也就是统计某一时间有多少花开放在它之前，结束在它之后。因为一朵花的开始一定是比结束早的，所以并不需要关心匹配问题。
使用两个数组分别存储花开与花谢的日期，排序后使用二分法即可快速得到某一日期的花开数量。

  • 优点 ：不用存储花期的对应关系。
  • 缺点 ： 不能确定是哪些花在这天会开放。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
class Solution {
public:
    vector<int> fullBloomFlowers(vector<vector<int>>& flowers, vector<int>& persons) {
        int n = flowers.size();
        vector<int> starts(n), ends(n);
        for (int i = 0; i < n; ++i) {
            starts[i] = flowers[i][0];
            ends[i] = flowers[i][1];
        }
        sort(starts.begin(), starts.end());
        sort(ends.begin(), ends.end());
        int m = persons.size();
        vector<int> ans(m);
        for (int i = 0; i < m; ++i) {
            int x = upper_bound(starts.begin(), starts.end(), persons[i]) - starts.begin();
            int y = lower_bound(ends.begin(), ends.end(), persons[i]) - ends.begin();
            ans[i] = x - y;
        }
        return ans;
    }
};

简单的拓扑排序

[OI WiKi]什么是拓扑排序?
简单来说，拓扑排序要解决的问题是给一个有向无环图的所有节点排序。
使用一个队列维护入度为零的节点，取出队列中的节点，存入答案，并把该节点的后续节点入度减一，得到新的有向图。

例题一 ： 标准拓扑

课程表II

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
class Solution {
public:
    vector<int> findOrder(int numCourses, vector<vector<int>>& prerequisites) {
        vector<vector<int>>g(numCourses);
        vector<int>indgree(numCourses,0);
        for(auto&i:prerequisites){
            g[i[1]].push_back(i[0]);
            indgree[i[0]]++;
        }
        queue<int>q;
        for(int i=0;i<numCourses;i++){
            if(indgree[i]==0){
                q.push(i);
            }
        }
        vector<int>ans;
        while(!q.empty()){
            int i=q.front();
            q.pop();
            ans.push_back(i);
            for(int j:g[i]){
                if(--indgree[j]==0){
                    q.push(j);
                }
            }
        }
        return ans.size()==numCourses?ans:vector<int>();
    }
};

例题二 ： 拓扑+贪心

课程表III

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
class Solution {
public:
    int scheduleCourse(vector<vector<int>>& courses) {
        //按ddl排序
        sort(courses.begin(),courses.end(),[](const auto&a1,const auto&a2){
            return a1[1]<a2[1];
        });
        priority_queue<int>q;
        int total=0;
        for(const auto& cour:courses){
            int ti=cour[0],ddl=cour[1];
            if(total+ti<=ddl){
                //可以加入,直接加
                total+=ti;
                q.push(ti);
            }
            else if(!q.empty()&&q.top()>ti){
                //有安排的情况下，最后安排的事情耗时大于当前事件
                total-=q.top()-ti;
                q.pop();
                q.push(ti);
            }
        }
        return q.size();
    }
};

例题三 ： 拓扑 + 搜索

课程表IV

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

class Solution {
public:
    vector<bool> checkIfPrerequisite(int numCourses, vector<vector<int>>& prerequisites, vector<vector<int>>& queries) {
        vector<vector<int> >g(numCourses);
        vector<int>indgree(numCourses,0);
        vector<vector<bool> >isPre(numCourses,vector<bool>(numCourses,false));
        for(auto&p:prerequisites){
            indgree[p[1]]++;
            g[p[0]].push_back(p[1]);
        }
        queue<int>q;
        for(int i=0;i<numCourses;i++){
            if(indgree[i]==0){
                q.push(i);
            }
        }
        while(!q.empty()){
            auto cur=q.front();
            q.pop();
            //cur的出口
            for(auto&out :g[cur]){
                isPre[cur][out]=true;
                //父节点传递
                for(int i=0;i<numCourses;i++){
                    isPre[i][out]=isPre[i][out] | isPre[i][cur];
                }
                --indgree[out];
                if(indgree[out]==0)q.push(out);
            }
        }
        vector<bool>ans;
        for(auto&que:queries){
            ans.push_back(isPre[que[0]][que[1]]);
        }
        return ans;
    }
};
|