同步于 博客园洛谷博客

由于各网站之间 LaTeX\LaTeX 及 Markdown 渲染方式的差异,故不保证除博客园外的排版完好,对于排版有明显不正确的地方请以博客园为主。

质数 约数 同余 组合

初等数论及组合数学入门复习概览

矩阵乘法

前置知识

  1. 矩阵:在 OI 中,一个 n×mn\times m 的矩阵可视为一个 n×mn\times m 的二维数组,其数学本质为一个向量,对于一个 n×mn\times m 的矩阵 AA,表示为:

A=[a1,1a1,2a1,ma2,1a2,2a2,man,1an,2an,m]A = \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,m} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,m} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n,1} & a_{n,2} & \cdots & a_{n,m} \end{bmatrix}

矩阵加法 矩阵减法:设 A,B,CA,B,C 均为 n×mn\times m 的矩阵,则有

C=A±Bi[1,n],j[1,m],Ci,j=Ai,j±Bi,jC=A\pm B\Leftrightarrow \forall i\in [1,n],\forall j\in [1,m],C_{i,j}=A_{i,j}\pm B_{i,j}

矩阵数乘:设 A,BA,B 均为 n×mn\times m 的矩阵,kRk\in \mathbb{R},则有

B=kAi[1,n]j[1,m]Bi,j=k×Ai,jB=k\cdot A\Leftrightarrow \forall i \in [1,n] \forall j \in [1,m] B_{i,j}=k\times A_{i,j}

矩阵乘法:由向量的数量积演变而来,若 AAn×mn\times m 的矩阵且 BBm×pm\times p 的矩阵,则对于 C=A×BC=A\times BCCn×pn\times p 的矩阵,且:

i[1,n],j[1,p],Ci,j=k=1mAi,k×Bk,j\forall i\in [1,n], \forall j\in [1,p], C_{i,j}=\sum_{k=1}^{m} A_{i,k}\times B_{k,j}

矩阵乘法模拟过程

一个模拟矩阵乘法的网站

注意若对于 n×mn\times m 的矩阵 AAp×qp\times q 的矩阵 BBmpm\neq p,则 AABB 无法进行矩阵乘法

单位矩阵:一个 n×mn\times m 的矩阵 II,且 i[1,min(n,m)],Ii,i=1\forall i\in [1,\min(n,m)], I_{i,i}=1,形如:

In=[100010001]I_n = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}

对于任意一个 n×mn\times m 的矩阵 AA,均存在 A×I=AA\times I=A(在 AA 的大小范围内相等)。

运算性质

  1. 结合律:(A×B)×C=A×(B×C)(A\times B)\times C=A\times (B\times C)

  2. 不存在交换律

矩阵乘法应用

加速线性递推

满足优化所需条件

  1. 在每个单位时间内每个阶段中只发生一次变化。(即递推)

  2. 变化的形式是一个递推式。

  3. 递推式在每个阶段都可以作用于不同的数据上,但本身保持不变。

  4. 未知的阶段数很大,但递推式长度很小。

例题Fibonacci 数列(斐波那契数列)

Fn={1 (n2)Fn1+Fn2 (n3)F_n = \left\{\begin{aligned} 1 \space (n \leqslant 2) \\ F_{n-1}+F_{n-2} \space (n\geqslant 3) \end{aligned}\right.

Fnmod109+7F_n\bmod {10^9+7},其中 1n<2631\leqslant n < 2^{63}

题解:设转移矩阵为 basebase,通项矩阵为:

[FiFi1]\begin{bmatrix}F_i & F_{i-1}\end{bmatrix}

对于每一项,有如下式:

Fi=Fi1×1+Fi2×1F_i=F_{i-1}\times 1+F_{i-2}\times 1

Fi1=Fi1×1+Fi2×0F_{i-1}=F_{i-1} \times 1+F_{i-2}\times 0

由系数得:

base=[1110]base=\begin{bmatrix}1 & 1\\ 1 & 0\end{bmatrix}

则有:

[FiFi1]=base×[Fi1Fi2]=[1110]×[Fi1Fi2]=[Fi1×1+Fi2×1Fi1×1+Fi2×0]=[FiFi1]\begin{aligned}\begin{bmatrix}F_i & F_{i-1}\end{bmatrix} &=base\times \begin{bmatrix}F_{i-1} & F_{i-2}\end{bmatrix}\\ &=\begin{bmatrix}1 & 1\\ 1 & 0\end{bmatrix}\times \begin{bmatrix}F_{i-1} & F_{i-2}\end{bmatrix}\\ &=\begin{bmatrix}F_{i-1}\times 1+F_{i-2}\times 1 & F_{i-1}\times 1+F_{i-2}\times 0\end{bmatrix}\\ &=\begin{bmatrix}F_i & F_{i-1}\end{bmatrix} \end{aligned}

故对于第 nn 项 Fibonacci 数列,只需求 basen1base^{n-1}(在 base0=Ibase^0=I 时可用此)或 [11]×basen2\begin{bmatrix}1 & 1\end{bmatrix}\times base^{n-2},最终矩阵 Ans1,1Ans_{1,1} 即为答案。

时间复杂度 O(23logn)\mathcal{O}(2^3\log n)

扩展 Fibonacci

对于形似

Fi={p1,p2,,1imFi1+Fim,im+1F_i=\begin{cases}p_1,p_2,\dots,1\leqslant i \leqslant m\\ F_{i-1}+F_{i-m},i\geqslant m+1\end{cases}

的递推式,都可以通过矩阵乘法加速线性递推在 O(m3logn)\mathcal{O}(m^3\log n) 的时间复杂度内求出 FnF_n 的值。

广义矩阵乘法

AAn×mn\times m 的矩阵,BBm×pm\times p 的矩阵,则对于 C=A×BC=A\times BCCn×pn\times p 的矩阵,有:

i[1,n],j[1,p],Ci,j=k=1m(Ai,kBk,j)\forall i\in [1,n],\forall j\in [1,p],C_{i,j}=\bigoplus_{k=1}^{m}(A_{i,k}\otimes B_{k,j})

且下列运算律成立:

交换律ab=baa\otimes b=b\otimes a

结合律(ab)c=a(bc)(a\otimes b)\otimes c=a\otimes (b\otimes c)

分配律a(bc)=(ab)(ac)a\otimes (b\oplus c)=(a\otimes b)\oplus (a\otimes c)

则称其为广义矩阵乘法。

其中 \oplus 常为 min\minmax\max\otimes++

代码

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
struct Matrix{
int n,m;
int a[N][N];
inline Matrix(int x,int y){memset(a,0,sizeof(a));n=x,m=y;}
inline void setI(){
int k=min(n,m);
memset(a,0,sizeof(a));
F(i,1,k) a[i][i]=1;
return ;
}
inline void print(){
F(i,1,n){
F(j,1,m) printf("%lld ",a[i][j]);
puts("");
}
return ;
}
inline Matrix operator+(const Matrix& x) const{
Matrix ret(n,m);
F(i,1,n) F(j,1,m)
ret.a[i][j]=a[i][j]+x.a[i][j],ret.a[i][j]%=p;
return ret;
}
inline Matrix operator-(const Matrix& x) const{
Matrix ret(n,m);
F(i,1,n) F(j,1,m)
ret.a[i][j]=a[i][j]-x.a[i][j];
return ret;
}
inline Matrix operator&(const int& k) const{
Matrix ret(n,m);
F(i,1,n) F(j,1,m)
ret.a[i][j]=a[i][j]*k,ret.a[i][j]%=p;
return ret;
} // 数乘
inline Matrix operator*(const Matrix& x) const {
Matrix ret(n,x.m);
F(i,1,n) F(j,1,x.m) F(k,1,m)
ret.a[i][j]+=a[i][k]*x.a[k][j],ret.a[i][j]%=p;
return ret;
}
inline Matrix operator^(const int &x) const {
int b=x;
Matrix a=*this,I(a.n,a.m);
I.setI();
while(b){
if (b&1) I=I*a;
a=a*a,b>>=1ll;
}
return I;
}
};
// Matrix a(n,m); 以新建一个 n*m 的全 0 矩阵 a。

0/1 分数规划

定义:类似于给定 22 个含 nn 个整数的序列 a,ba,b,求一组解 xi (i[1,n],xi[0,1])x_i\space (\forall i\in [1,n],x_i\in[0,1])(即选或不选)使得下式最大化,给出其最大值:

i=1nai×xii=1nbi×xi\frac{\sum_{i=1}^{n}a_i\times x_i}{\sum_{i=1}^{n}b_i\times x_i}

求解

二分

二分答案 midmid

i=1nai×xii=1nbi×xi>mid\frac{\sum_{i=1}^{n}{a_i\times x_i}}{\sum_{i=1}^{n}{b_i\times x_i}} > mid

i=1nai×ximid×i=1nbi×xi>0\sum_{i=1}^{n}{a_i\times x_i}-mid\times \sum_{i=1}^{n}{b_i\times x_i} > 0

i=1nxi×(aimid×bi)>0\sum_{i=1}^{n}{x_i\times (a_i-mid\times b_i)} > 0

二分左式最大值即可,时间复杂度 O(nlog2n)\mathcal{O}(n\log^2 n)

Dinkelbach 算法

取一个值 xx 并将其多次迭代使其逼近于答案。

此方法不常用,时间复杂度玄学。

高斯消元

定义:通过初等行变换把增广矩阵变为简化阶梯形矩阵的线性方程组求解算法。

增广矩阵:在线性代数中系数矩阵右边添上常数列所得到的矩阵。

线性方程组转化为增广矩阵:对于一个 nn 元一次方程组:

{a1,1x1+a1,2x2++a1,nxn=p1a2,1x1+a2,2x2++a2,nxn=p2an,1x1+an,2x2++an,nxn=pn\begin{cases} a_{1,1}x_1+a_{1,2}x_2+\dots+a_{1,n}x_{n}=p_1\\ a_{2,1}x_1+a_{2,2}x_2+\dots+a_{2,n}x_{n}=p_2\\ \vdots\\ a_{n,1}x_1+a_{n,2}x_2+\dots+a_{n,n}x_{n}=p_n \end{cases}

可转化为一个 nnn+1n+1 列的增广矩阵:

[a1,1a1,2a1,np1a2,1a2,2a2,np2an,1an,2an,npn]\left[\begin{array}{cccc|c} a_{1,1} & a_{1,2} & \dots & a_{1,n} & p_1\\ a_{2,1} & a_{2,2} & \dots & a_{2,n} & p_2\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ a_{n,1} & a_{n,2} & \dots & a_{n,n} & p_n \end{array}\right]

初等行变换:平时解方程组的方法,在增广矩阵上可概括为以下三种操作:

  1. 用一个非零的数乘某一行。

  2. 把其中一行的若干倍加到另一行上。

  3. 交换两行的位置。

阶梯形矩阵:形如

[a1,1a1,2a1,np10a2,2a2,np200an,npn]\left[\begin{array}{cccc|c} a_{1,1} & a_{1,2} & \dots & a_{1,n} & p_1\\ 0 & a_{2,2} & \dots & a_{2,n} & p_2\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & \dots & a_{n,n} & p_n \end{array}\right]

的矩阵,其系数矩阵部分被称为上三角矩阵

简化阶梯矩阵:形如

[a1,100p10a2,20p200an,npn]\left[\begin{array}{cccc|c} a_{1,1} & 0 & \dots & 0 & p_1\\ 0 & a_{2,2} & \dots & 0 & p_2\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & \dots & a_{n,n} & p_n \end{array}\right]

的矩阵,由阶梯形矩阵简化而来,该矩阵直接给出了方程的解。

自由元:若 i[1,n],xi\exists i\in [1,n],x_i 可以取任意值,即方程组有无数组解,则称 xix_i 为自由元。

主元:不为自由元的元即为主元。

代码核心逻辑

  1. 对于每一个 xix_i,找出最大的 ak,i,(k[1,n])a_{k,i},(k\in [1,n])ak,i0a_{k,i}\neq 0 所属的 kk,将该行移动到第 ii 行(若找不出 kk 则说明方程组无解或存在自由元)。

  2. 使当前行的首项变为 11,当前行其它项除原首项系数(即将当前行约分至首项系数为 11)。

  3. 将当前行下面所有行的当前列都消成 00,同时更新正在处理的行后面的系数。

  4. 若找不出 kk,处理完剩余行后遍历所有行,若 i[1,n],ai,n+10\exists i \in [1,n],a_{i,n+1}\neq 0,则说明方程组无解(未消元的行系数必为 00),反之则存在自由元,方程有无数组解。

  5. 反之,则将阶梯形矩阵转化为简化阶梯矩阵,ai,n+1a_{i,n+1} 即为 xix_i 的对应解。

代码

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
const double eps=1e-6;
int n;
double a[N][N];
inline int Gauss(){
int i=1; // 当前行
for (int j=1;j<=n;++j){ //当前列
int maxh=i;
F(x,i,n)
if (abs(a[x][j])>abs(a[maxh][j])) maxh=x; // 找 k
if (abs(a[maxh][j])<eps) continue; // 若找不出 k,则说明方程组无解或存在自由元
F(x,j,n+1) swap(a[maxh][x],a[i][x]); // 将该行移动到第 i 行
FJ(x,n+1,j) a[i][x]/=a[i][j]; // 使首项为 1,其它项需除原首项系数
F(x,i+1,n) if (abs(a[x][j])>eps)
FJ(y,n+1,j) a[x][y]-=a[i][y]*a[x][j]; // 处理后面行的系数
++i;
}
if (i<=n){ // 若存在 i 找不出 k
F(x,i,n)
if (abs(a[x][n+1])>eps) return -1; // 系数项不为 0,则方程组无解
return -2; // 存在自由元
}
FJ(x,n,1) F(y,x+1,n)
a[x][n+1]-=a[x][y]*a[y][n+1]; // 转化为简化阶梯矩阵
return 0;
}
signed main(){
n=read();
F(i,1,n) F(j,1,n+1) a[i][j]=read();
int ret=Gauss();
if (ret==-1) puts("No Solution"),exit(0);
else if (ret==-2) puts("Multiple Solutions"),exit(0);
F(i,1,n)
printf("x_%lld=%.2f\n",i,abs(a[i][n+1])<eps?0:a[i][n+1]); // 防止输出 -0.00
// 存在唯一解,则 a_{i,n+1} 为 x_i 的对应解
}

高斯-约旦消元

相比高斯消元的优势

  1. 更高的精度。

  2. 代码更简洁。

缺点

  1. 判无数组解的时候比较麻烦。

代码核心逻辑

  1. 找出一个此前未被找过的未知数作为主元。

  2. 将当前主元移到当前行。

  3. 使当前行的首项变为 11,后面的项约分。

  4. 因为这样处理后 ai,ia_{i,i} 即为该主元的系数,所以若 i[1,n],ai,i=0\exists i\in [1,n],a_{i,i}=0,则说明要么无解或有无数组解。

  5. 如果无解或有无数组解,则遍历每个未知数,若存在常数项不为 00,则说明方程组无解,若不存在则说明方程组有无数组解。

  6. 若有唯一解,则 ai,n+1ai,i\frac{a_{i,n+1}}{a_{i,i}} 即为 xix_i 的对应解。

代码

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
inline int GaussJordan(){
int i=0;
F(j,1,n){
++i; // 实际上 i=j,这里仅为了便于区分行列
int maxh=i;
F(x,1,n)
if(x>j||abs(a[x][x])<eps) // 防止 a[x][x]=0 而被用上
if (abs(a[x][j])>abs(a[maxh][j])) maxh=x; // 找主元
F(x,1,n+1) swap(a[i][x],a[maxh][x]); // 移到第 i 行,使后面找到的主元均为之前未选过的
if (abs(a[i][i])<eps) continue;
F(x,1,n){
if (x!=i){
double cur=a[x][j]/a[i][j];
F(y,j+1,n+1) a[x][y]-=a[i][y]*cur; // 约分
}
}
}
bool flg=0;
F(i,1,n) if (abs(a[i][i])<eps){
if (abs(a[i][n+1]/a[i][i])>eps) return -1;
flg=1; // 无解比无数组解优先级高
}
return flg?-2:0;
}
signed main(){
n=read();
F(i,1,n) F(j,1,n+1) a[i][j]=read();
int ret=GaussJordan();
if (ret==-1) puts("-1"),exit(0);
else if (ret==-2) puts("0"),exit(0);
F(i,1,n)
printf("x_%lld=%.2f\n",i,abs(a[i][n+1]/a[i][i])<eps?0:a[i][n+1]/a[i][i]);
return 0;
}

求解异或方程组

对于一个 nn 元一次异或方程组:

{a1,1x1xora1,2x2xorxora1,nxn=p1a2,1x1xora2,2x2xorxora2,nxn=p2an,1x1xoran,2x2xorxoran,nxn=pn\begin{cases} a_{1,1}x_1\operatorname{xor}a_{1,2}x_2\operatorname{xor}\dots\operatorname{xor}a_{1,n}x_{n}=p_1\\ a_{2,1}x_1\operatorname{xor}a_{2,2}x_2\operatorname{xor}\dots\operatorname{xor}a_{2,n}x_{n}=p_2\\ \vdots\\ a_{n,1}x_1\operatorname{xor}a_{n,2}x_2\operatorname{xor}\dots\operatorname{xor}a_{n,n}x_{n}=p_n \end{cases}

依旧构造增广矩阵,因为异或本质就为不进位加法,所以仅需要在原高斯消元代码的基础上稍加修改即可,详见代码。

需要注意的是因为自由元只能取 0011,故出现自由元时解的数量并不是无数组,而是 2自由元个数2^{\text{自由元个数}}

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
inline int XorGauss(){
int i=1,ret=1;
F(j,1,n){
int maxh=i;
F(x,i+1,n)
if (a[x][j]) maxh=x;
if (!a[maxh][j]) continue;
F(x,j,n+1) swap(a[maxh][x],a[i][x]);
// 因为都为 0 或 1 也就不需要处理后面行系数了
F(x,i+1,n) FJ(y,n+1,j)
a[x][y]^=a[x][j]&a[i][y]; // 用异或替换减
++i;
}
if (i<=n){
F(x,i,n){
if (a[x][n+1]) return -1; // 出现 0=1,无解
ret<<=1ll; // 统计解的个数
}
}
return ret; // 因为自由元可为 0 或 1,故解的个数为 2^{自由元个数}
}

线性空间

前置知识

  1. 向量加法a+ba+b,其中 aabb 均为向量。

  2. 标量乘法(数乘)kaka,其中 aa 为向量,kk 为标量。

表出:若向量 bb 能够由向量 a1,a2,,aka_1,a_2,\dots,a_k 经过向量加法和标量乘法得出,则称 bb 能被向量 a1,a2,,aka_1,a_2,\dots,a_k 表出。

线性空间:若干个向量 a1,a2,,aka_1,a_2,\dots,a_k 能表出的所有向量所构成的集合。

生成子集:对于一个线性空间来说,表出其所有向量的若干个向量 a1,a2,,aka_1,a_2,\dots,a_k 为此线性空间的生成子集。

线性相关:任意选出线性空间中的若干个向量,若其中存在一个向量能够被其它向量表出,则称这些向量线性相关,反之则称这些向量线性无关

基底:对于一个线性空间来说,线性无关的生成子集被称为线性空间的基底,简称为

维数:对于一个线性空间来说,基所包含的向量个数被称为线性空间的维数。(一个线性空间所有基包含的向量个数都相等)

与矩阵的联系

对于一个 n×mn\times m 的矩阵,其每一行都可看作一个长度为 mm 的向量,称其为行向量nn 个行向量所能表出的所有向量构成一个线性空间,其维数被称为矩阵的行秩

类似地,有列向量列秩,由于行秩一定等于列秩,故合称为

将该矩阵变为系数列均为 00 的增广矩阵进行高斯消元,得到简化阶梯矩阵。该简化阶梯矩阵所有非零行向量都线性无关,其构成的生成子集为该线性空间的一个基,这个基的个数为该矩阵的秩。

异或空间

表出:若整数 bb 能够由整数 a1,a2,,aka_1,a_2,\dots,a_k 经过异或运算得出,则称 bb 能被 a1,a2,,aka_1,a_2,\dots,a_k 表出。

异或空间:若干个整数 a1,a2,,aka_1,a_2,\dots,a_k 能表出的所有数字所构成的集合。

其它的定义和对于线性空间的定义类似,其中异或空间的基又被定义为异或空间的极大线性无关子集

对于 nn 个整数 aia_iai[0,2k1]a_i\in [0,2^k-1],则可以构造一个 nnmm 列的 0101 矩阵,对于第 ii 行,从右到左依次为 (ai)2(a_i)_2 从低位到高位,不够的在左边补 00

比如,对于 n=5,k=4,a=[5,12,2,7,9]n=5,k=4,a=[5,12,2,7,9],得到的矩阵为:

[01011100001001111001]\begin{bmatrix} 0 & 1 & 0 & 1\\ 1 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 1 & 1 & 1\\ 1 & 0 & 0 & 1 \end{bmatrix}

不够的在右边补 00,然后通过异或高斯消元求出它的简化阶梯矩阵:

[10010101001000000000]\begin{bmatrix} 1 & 0 & 0 & 1\\ 0 & 1 & 0 & 1\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix}

即该异或空间的基为 [9,5,2][9,5,2],其数学含义为:

i[1,n],S,S,xorxS=ai\forall i \in [1,n],\exists S\subseteq \text{基},S\neq \varnothing,\operatorname{xor}_{x\in S} = a_i

数据结构求解线性基

以异或线性基为例,用 Ai(i[1,log2(max{ai})])A_i(i\in[1,\left \lceil \log_2(\max\{ a_i \})\right \rceil]) 保存基底。

对于插入线性基操作 Insert(x)\text{Insert(x)},有:

  1. xx 最高位 11 在第 ii 位。

  2. Ai=0A_i=0,则使 Ai=xA_i=x

  3. 反之,则使 x=xxorAix=x\operatorname{xor} A_i,然后继续 Insert(x)\text{Insert(x)},直到 x=0x=0

AiA_i 的非 00 项即为基底。

代码

1
2
3
4
5
6
7
8
9
inline void Insert(int x){
FJ(i,51,0){ // 从高到低找,找到的第一个 1 必定为最高位的 1
if (x>>i){
if (bs[i]) x^=bs[i];
else {bs[i]=x;return;}
}
}
return ;
}

1
2
3
4
5
6
inline void Insert(int x){
if (!x) return ;
int pos=__lg(x); // 找最高位的 1
if (bs[pos]) Insert(x^bs[pos]);
else {bs[pos]=x;return;}
}