山东游记 - 五一数学专题

CHANGE LOG

• 2023.5.4 更新扩展中国剩余定理代码与解释。
• 2023.5.5 添加线性预处理逆元解释说明。
• 2023.5.7 重构初等矩阵部分。
• 2023.5.11 添加快读模板。
• 2023.5.24 更新高斯消元法代码，添加高斯 - 约旦消元法。
• 2023.6.3 补充逆元部分概念，修改时间复杂度事实性错误。
• 2023.6.9 补充逆元一般性求法，增加扩展卢卡斯定理（exLucas）。
• 2023.7.29 优化求逆元部分文章结构，修改扩展欧几里得事实性错误，对其代码进行优化。

写在前面

本博客中所有 ll 代指 long long，ksm 代指 快速幂 代码，如下：

template <typename _Tp> _Tp ksm(_Tp base,int ind,int mod=mod) {
	_Tp ans=1;
	while(ind) {
		if(ind&1) ans=ans*base%mod;
		base=base*base%mod,ind>>=1;
	}
	return ans%mod;
}

请注意，当前写法是每次 乘法后取模，未拆解大数乘法为加法，这意味着 int 型快速幂计算 不能涉及过程量大于 46341（即 \(\lfloor\sqrt{2147483647}\rfloor\)）的数字，大数快速幂请自觉送参 long long。

本快速幂模板运用了运算符 * 与 %，并使用了转换构造函数，呼吁广大自编结构体完善构造函数和重载 operator * 与 operator %。

同时，推荐使用如下快读模板：

template <typename _Tp> inline void read(_Tp &x) {
    int f=1; char ch=getchar(); x=0;
    while(!isdigit(ch)) { if(ch=='-') f*=-1; ch=getchar(); }
    while(isdigit(ch)) x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    x*=f;
}

---

从接触竞赛到现在，不知不觉 8 个月了；如今，春回大地，万物复苏——是时候出去看看了！

于是——

TSOI2022 进军山东！

4.28 启程

上午八点，TSOI2022 师生 5 人从唐一出发，开始了征程。

上午十点半左右正式沿秦滨高速跨入山东，久违地看见了浪漫的浓积云~

下午五点左右安置完毕，晚上还去了海边。在辛安河特大桥边到沙滩上跑了跑。

归程下起了细雨，意外的有趣。

这两天还了解了一些东西，写一下：

• bool 其实是 char，所以理论上是 true 和 false，事实上是 0 或 非 0，是可能出现 0 和 1 之外的数字的。

• 数组越界后操作会在 内存中已有的位置 进行，所以在你数组越界后，其他变量是有可能遭殃的。

• 结构体内的函数可以 随意调用，即使是 顺序在当前函数之后。这意味着你 无需提前声明。

• 结构体和类可以用 *this 来代表当前所在的变量。

---

基础前置

4.29 - Day 1

上午

见到了钟皓曦学长，NOI2013 金牌，清华姚班大佬。
讲了一些比较基础的东西，大部分还都有了解，阿弥陀佛。

那就说几个点吧，一些零碎的知识：

构造函数

在结构体或类中，每次声明变量时执行的函数。声明的 等号 = 和 括号 其实是在调用构造函数。如：

struct high {int len,x[2023];};
high a (0,{0});//显式调用初始化构造函数
high b={0,{0}};//隐式调用初始化构造函数
high c=a;//隐式调用构造函数
high d;//调用默认构造函数
high e=1;//调用转换构造函数
a=1;//赋值
b=a;//赋值

其实都是调用了构造函数。

由上面代码可知，构造函数分为四类：

1. 默认构造函数：声明的时候啥也不说。

2. 初始化构造函数：使用传参构造。

3. 复制构造函数：用同类构造同类。

4. 转换构造函数：使用其他类型的变量构造。

在我们不显式编写构造函数时，系统会 自动生成隐式构造函数。但一旦我们自己编写了构造函数，所有 隐式构造函数全部 失效。

所以上面代码中，就应在结构体中编写：

struct high {
	int len,bit[2023];
	high() {len=1,memset(bit,0,sizeof bit);}//默认构造函数 
	high(int given_len,std::vector <int> given_bit) {
		len=given_len;
		std::copy(given_bit.begin(),given_bit.begin()+given_bit.size(),bit);
	}//初始化构造函数 
	high(int x) {
		int now=0;
		while(x) bit[++now]=x%10,x/=10;
		len=now;
	}//转换构造函数 
};

隐式的复制构造函数始终存在，可以不写。

转换函数

顾名思义，就是用于将当前变量转换为其他类型变量的函数。如：

operator int() const {
		int now=0;
		for(int i=len;i>0;--i)
			now=now*10+bit[i];
		return now;
	}

在刚才的结构体 high 中加入这段代码，就可以实现 high 转 int 了。

取模运算的修正与优化

负数取模后还是负数。正确的做法是，对于 int 型变量 a 和正整数 b，使用 (a%b+b)%b 得出正确结果。

值得一提的是，乘法的速度比除法和取模快。所以，当多次取模时，因为你两个加数 a 和 b 本身是取过模的，所以对于结果 ans 和取模数 mod，可以直接 ans=a+b>mod?a+b-mod:a+b; 总之，面对几乎极限的时间要求时，取模是有可能可以榨出一些油水的。

5.2 补：等比数列快速幂

一种使用类似快速幂求 \(\sum \limits_{i=0}^{n-1} p^i\) 的方法，我愿称此为 阶乘快速幂（fac_pow）。

出现在 杰杰的女性朋友 一题中，后来好不容易研究明白了。

我们设 \(S_n\) 为前 \(n\) 项的和，即 \(\sum \limits_{i=0}^n p_i\)，就有

\[S_{\lfloor n/2 \rfloor -1}=1+p+p^2+\cdots+p^{\lfloor n/2 \rfloor -1}, \]

\[S_{n-1}\begin{cases}n \equiv 0 \pmod 2 & = 1+p+p^2+\cdots+p^{n-1}\\& = 1+p+\cdots+p^{n/2-1}+p^{n/2}+p^{n/2+1}+\cdots+p^{n-1}\\& = (1+p^{n/2})S_{n/2-1}\\n \equiv 1 \pmod 2 &= 1+p+p^2+\cdots+p^{n-1}\\& =1+p+\cdots+p^{\lfloor n/2 \rfloor -1} + p^{\lfloor n/2 \rfloor} + p^{\lfloor n/2 \rfloor +1} +\cdots+ p^{n-1}\\& =(1+p^{\lfloor n/2 \rfloor })S_{\lfloor n/2 \rfloor -1} + p^{ n -1}\end{cases}. \]

至此，我们找到了 \(n\) 和 \(n/2\) 的关系。然后，我们就可以考虑采用倍增的方法，采用类似快速幂的方式计算。

如果采用右推，即完全类似于快速幂的方法，以对指数 \(ind\) 的右移 >>=1 为基本操作，每次操作之间并不是严格的二倍关系，需要记录大量的过程量，这对复杂度并不友好，所以我们考虑采用正推。

左推，即以对 \(1\) 对左移操作 <<=1 为基本操作，直至左移到与指数相同为止。由于左移就是严格的二倍关系，只在当前位为 \(1\) 时加上 \(p^{ n-1 }\) 即可，而 \(p^{n-1}\) 正是 \(p^{2{\lfloor n/2 \rfloor}}\)，同样满足倍增关系，所以一切就变得非常好办了。

template <typename _Tp> _Tp fac_pow(_Tp base,int ind,int mod=INT_MAX) {
	_Tp ans=0,now=1; int st[100] {0};
	while(ind) {st[++st[0]]=ind&1,ind>>=1;}
	while(st[0]) {
		ans=(ans+(now*ans%mod))%mod,now=now*now%mod;
		if(st[st[0]]) ans=(ans+now)%mod,now=(now*base)%mod;
		st[0]--;
	}
	return ans;
}

其实快速幂也可以同理左推实现，同时，我们也可以通过寻找倍增关系快速求出 \(\sum \limits_{i=0}^{n} p^i\) 和 \(\sum \limits_{i=1}^n p^i\)， 他们都是异曲同工的。

---

下午

其实第一天下午就讲了很多线性代数啊……

讲了素数筛法，矩阵等等。大部分还都算是见过的，比较好理解。

一起上课的很多同学年龄都比我们小很多，和我们听同样的课程，山东的竞赛确实比河北要强啊。

矩阵快速幂优化

矩阵运算都比较简单，需要注意的是，矩阵满足结合律但 不满足交换律。（abel 群笑话 hhhh）

矩阵乘法在 OI 中有很多用处：

1. 加速递推：数列的前几项用向量表示，使用矩阵快速幂加速计算过程，\(O(n) \to O(\log n)\)。

2. 加速一维相同的二维数组计算：比如下面这个问题：

   给定一个有向图，从 \(A\) 点走 \(k\) 步恰好走到 \(B\) 点的方案数？

   这道题用 \(E(x,y) \gets 1\) 表示 \(x\) 和 \(y\) 之间有一条有向边，可得对于 \(C=A \times A\) 有 \(C(i,j)=\sum A(i,k) \times A(k,j)\)，把 \(k\) 作为中继点，不难看出其实这表达的就是从 \(i\) 走两步恰好到 \(j\) 的方案数。所以，矩阵快速幂，直接解决。

   可以看一把 P4159 [SCOI2009] 迷路 ，这道题用前缀优化建边和矩阵快速幂来解决。

数组遍历顺序优化

在计算机中，所有临时数据将优先存储在缓存中。遇到数组时，缓存中会将当前访问的数组下标的后几个数同时写入缓存，所以在遍历时，按顺序逐个遍历确实是最快的。

比如经典的矩阵乘法：

for(int i=1;i<=ans.n;++i)
	for(int k=1;k<=m;++k)
		for(int j=1;j<=ans.m;++j)
			ans.x[i][j]+=x[i][k]*a.x[k][j];

i,k,j 的遍历顺序就是速度最快的。

---

初等数论

4.30 - Day 2

上午

一觉醒来啊…7：05，爽。

逆元

众所周知，加减乘法满足 过程中多次取模 和 对结果取模 等效。我愿称此为 模运算的盖斯定律。然而，除法并不满足模运算的盖斯定律。我们可以计算 \(12 \div 4 \bmod 5\)，却对 \(2 \div 4 \bmod 5\) 无计可施。所以，我们提出了 逆元 解决这一问题：

如果 \(\gcd(a,m)=1\) 且存在唯一的 \(b\) 使得 \(a \times b \equiv 1 \pmod m\) 且 \(1 \leq b < m\)，则 \(b\) 为 \(a\) 在 \(\bmod m\) 意义下的 逆元。

事实上，逆元是一个极广阔的概念，他是指一个可以取消另一给定元素运算的元素。在数学里，逆元素广义化了加法中的加法逆元和乘法中的倒数。模意义下的逆元又称为 数论倒数，这也很容易理解。

求解逆元，我们有如下几种方法：

1. 费马小定理

   对于任意一个 质数 \(p\)，若 \(a,p\) 互质，有 \(a^{p-1} \equiv 1 \pmod p.\)

   由费马小定理，我们可以推出 \(a^{p-2} \equiv \dfrac{1}{a} \pmod p\)，快速幂求逆元即可。

2. 欧拉定理

   对于任意一个 自然数 \(m\)，若 \(a,m\) 互质，有 \(a^{\varphi(m)} \equiv 1 \pmod m.\) 其中，\(\varphi(m)\) 为欧拉函数。

   由于对于任意质数 \(p\) 有 \(\varphi(p)=p-1\)，费马小定理其实是欧拉定理在 \(m\) 为质数下的 特殊情况。

   接下来，我们就可以得出 \(\dfrac{1}{a} \equiv a^{\varphi(m)-1} \pmod m\)，于是将 \(\div a\) 转换为 \(\times \dfrac{1}{a}\)，最后转换成 \(a^{\varphi(m)-1}\)，就得到了 \(a\) 的逆元。

   Q: 怎么求欧拉函数？
   A: 你往下看。

3. 线性预处理序列逆元

   如果我们需要且仅需要求一组数的逆元，那么这还算是个不错的方法。

   我们可以线性处理前缀积 \(fac\)，（用这个名称是因为，当这一组数是从 \(1\) 开始的公差为 \(1\) 的等差数列时，前缀积就是 阶乘 了。）费马小 / 欧拉 / 扩欧求前缀积末项逆元 \(ifac_n\) ，向前逆推求出所有前缀积逆元 \(ifac\) ，再用 \(inv_i = fac_{i-1} \cdot ifac_i\) 的方法求出从 \(1\) 到 \(n\) 的所有逆元。本质是不断 约分。

   以从 \(1\) 开始的公差为 \(1\) 的等差数列为例，\(Code\) 如下：

   void cal_inv() {
   	fac[0]=1;
   	for(int i=1;i<=n;++i) fac[i]=1ll*fac[i-1]*i%p
   		ifac[n]=ksm(fac[n],p-2,p);
   	for(int i=n;i>=0;--i) ifac[i]=1ll*ifac[i+1]*(i+1)%p;
   	for(int i=1;i<=n;++i) inv[i]=1ll*fac[i-1]*ifac[i]%p;
   }
   

   这个方法非常巧妙，我们可以具体的板演一下。

   对于需要求逆元的序列 \(a_1,a_2,a_3,\cdots, a_n\)，设

   \[f_x= \prod \limits_{i=1}^{x} a_i, \]

   再设

   \[g_x=\dfrac{1}{f_x}, \]

   则有

   \[g_i=\dfrac{1}{a_1 a_2 a_3 \cdots a_i} = \dfrac{1}{a_1 a_2 a_3 \cdots a_i a_{i+1}} \times a_{i+1} = g_{i+1} \times a_{i+1}, \]

   故求 \(g_n\) 逆元，后可一直逆推至 \(g_1\) 逆元。后可知

   \[\dfrac{1}{a_i} = \dfrac{1}{a_1 a_2 \cdots a_{i-1} a_i} \times a_1 a_2 \cdots a_{i-1} = f_{i-1} \cdot g_i \]

   就求出逆元了。

   建议练习 洛谷P5431 【模板】乘法逆元 2。

4. 另一种线性预处理逆元方法（不常用）

   我们将 \(p \bmod i\) 做如下变换：

   设 \(k=\lfloor{p \div i}\rfloor,r=p \bmod i,\)

   得：

   \[p=ki+r \]
   \[0 \equiv ki+r \pmod p \]
   \[-r \equiv ki \pmod p \]
   \[\dfrac{r}{i} \equiv -k \pmod p \]
   \[\dfrac{1}{i} \equiv \dfrac{p-k}{r} \pmod p \]

   逆元可求。

   \(Code\)：

   inv[1]=1;
   void cal_inv() {
   	for(int i=2;i<=n;++i) {
   		int k=p/i,r=p%i;
   		inv[i]=1ll*(p-k)*inv[r]%p;
   	}
   }
   

   UPD. 2023.5.5 关于预处理逆元，应用场景有限，特殊时候特殊方法对待。在 \(p\) 为质数且不连续的时候，首选 费马小 + 快速幂，\(p\) 不为质数……祝你好运。

5. 扩展欧几里得

   UPD. 2023.6.9 隆重介绍！应用场景更广泛且复杂度比肩费马小的 扩欧求逆元！

   根据逆元的定义，易得

   \[a \cdot x \equiv 1 \pmod p, \]

   其中 \(x\) 为 a 逆元。
   那么对于任意整数 \(k,\) 就有

   \[k \cdot p +1 \equiv 1 \pmod p \]
   \[k \cdot p + a \cdot x \equiv 1 \pmod p, \]

   由于逆元定义写到 \(\gcd (a,p) =1,\) 那么我们就可以用下文中提到的 exgcd 求解方程

   \[p \cdot k + a \cdot x = 1. \]

   \(x\) 即为所求。

   这可真是太强了，即使 \(p\) 不是素数仍然可以这么求解，所以对于 \(p\) 不是素数的情况 认定他就对了。（稍微会比费马小慢，所以可以是质数的时候可以考虑费马小。）

欧拉函数

众所周知对于质数 \(p\) 有 \(\varphi(p) = p-1\)。

后可证：

\[\varphi(p^2) = p(p-1) \]

\[\varphi(p^n)=\dfrac{p-1}{p}\times p^n \]

\[\varphi(n)=n - \dfrac{n}{p_1} - \dfrac{n}{p_2} + \dfrac{n}{p_1 p_2} \]

\[\varphi(n) = n \times \dfrac{p_1-1}{p_1} \times \cdots \times \dfrac{p_k-1}{p_k} \]

其中第三式 \(n\) 为质数 \(p_1\) 与 \(p_2\) 的积，利用第二式结合容斥原理推出；第四式在第三式基础上因式分解，然后合理外推，\(n\) 表示任意自然数，\(p_i\) 表示 \(n\) 的质因数。

至此，我们可以 \(O(\sqrt{n})\) 计算欧拉函数了。

int getphi(int x) {
	int ret=x,tmp=x;
	for(int i=2;i*i<=x;++i) {
		if(tmp%i==0) {
			ret=ret/i;
			ret=ret*(i-1);
		}
		while(!tmp%i) tmp=tmp/i;
	}
	if(tmp>1) ret=ret/tmp*(tmp-1);
	return ret;
}

Miller-Rabin 素性测试

如果 \(n\) 是质数，则以下两条至少一条成立：

取 \(a<n\)，设 \(n-1=d\times 2^r\)，

\[a^d \equiv 1 \pmod n \]

\[\text{or} \]

\[\exists 0 \leq i < r, \text{s.t. } a^{d\times2^i} \equiv 1 \pmod n \]

但是，有一些合数也是有记录通过 Miller-Rabin 测试的，所以我们可以多次测试，提高概率。

利用 Miller-Rabin，我们可以 \(O(\log n)\) 判断质数。

\(Code\)：

bool miller_rabin(int n,int a) {
	int d=n-1,r=0;
	while(!d%2) d=d/2,r++;
	int x=ksm(a,d,n);
	if(x==1) return 1;
	for(int i=0;i<r;++i) {
		if(x==n-1) return 1;
		x=1ll*x*x%n;
	}
	return 0;
}
bool is_prime(int n) {
	if(n<2) return 0;
	for(int i=1;i<=23;++i) {
		int a=rand()%(n-1)+1;
		if(!miller_rabin(n,a)) return 0;
	}
	return 1;
}

或者，为提高成功几率，我们可以制定一个 \(prime\_list\)，不用太大，几个数即可，然后用这些数进行检验，只要能保证他在 int 范围内正确就行。

int prime_list[]={2,3,5,7,13,23,37,73};
bool miller_rabin(int n,int a) {
	int d=n-1,r=0;
	while(!d%2) d=d/2,r++;
	int x=ksm(a,d,n);
	if(x==1) return 1;
	for(int i=0;i<r;++i) {
		if(x==n-1) return 1;
		x=1ll*x*x%n;
	}
	return 0;
}
bool is_prime(int n) {
	if(n<2) return 0;
	for(int i=0;i<8;++i) {
		if(n==prime_list[i]) return 1;
		if(!n%prime_list[i]) return 0;
		if(!miller_rabin(n,prime_list[i]%n)) return 0;
	}
	return 1;
}

扩展欧几里得 - exgcd

给定 \(a,b\)，已知 \(\gcd(a,b)=g\)，求 $x,y,\text{ s.t. } $

\[xa+yb=g. \]

如何解决？

之前学到的辗转相除：

\[gcd(a,b)=gcd(b,a \bmod b)=g \]

\[xa+yb=x'b+y'(a \bmod b) \]

\[xa+yb=x'b+y'(a-b \cdot \lfloor \dfrac{a}{b} \rfloor) \]

\[xa+yb=x'b+y'(a-b \cdot \lfloor \dfrac{a}{b} \rfloor) \]

\[xa+yb=y'a+(x'-y' \cdot \lfloor \dfrac{a}{b}\rfloor)b \]

\[x=y',y=(x'-y' \cdot \lfloor \dfrac{a}{b}\rfloor) \]

问题解决。

\(Code\)：

int exgcd(int a,int b,int &x,int &y) {
	if(!b) {x=1,y=0; return a;}
	int g=exgcd(b,a%b,x,y),tmp=x;
	x=y,y=tmp-(a/b)*y;
	return g;
}

注意，这样递归的边界条件是 \(x \cdot (a,b) + y \cdot 0 = (a,b)\)，所以在判断 !b 后 \(x \gets 1,y \gets 0\)，否则，\(y\) 不影响结果（废话，因为是 \(y \cdot 0\)），\(x\) 赋的值是几 \(ax+by\) 就是几倍 \(\gcd (a,b)\) 。

UPD. 2023.7.29 谢邀，人在 ZROI。以前对 \(x,y\) 我们都是手动 swap，其实可以直接动个小手脚：

int exgcd(int a,int b,int &x,int &y) {
	if(!b) { x=1,y=0; return a; 	}
	int g=exgcd(b,a%b,y,x);
	return y-=a/b*x,g;
}

一个意思。

裴蜀定理

害，就是你可以一步走 \(x\) 米，也可以一步走 \(y\) 米，问你离原点最近几米，就是 \(\gcd(x,y)\)。

我要提出 碘酊定理 ！

\[\forall x \in N^*, 1 \equiv 1 \pmod x \]

伟大！伟大啊！！！

成功解决中国剩余定理。

下午

跟大家说一声，这个 \(\LaTeX\) 不念雷泰克斯，

补：关于 exgcd

我们使用扩欧解决 P1082 [NOIP2012 提高组] 同余方程，找到的解不一定是最小的解。我们用找到的解除以他们的 \(\gcd\)，寻找到一组 互质 的解 \(x,y\)，就可以找到当前方程的所有解了：

\[\begin{cases}x'=x+kb\\y'=y+ka\end{cases},k \in R \]

（扩展）中国剩余定理 - CRT & exCRT

解决

\[\begin{cases}x \equiv a_1 \pmod {p_1}\\x \equiv a_2 \pmod {p_2}\\\cdots\\x \equiv a_k \pmod {p_k}\end{cases} \]

的问题。

不难得出，事实上同余方程组的解是一个同余方程。所以，其实中国剩余定理解决的是 同余方程的合并问题。

我们先以解

\[\begin{cases} x \equiv a_1 \pmod {p_1}\\ x \equiv a_2 \pmod {p_2}\\ \end{cases} \]

为例，介绍两种解法：

1. 大数翻倍法（暴力）

   《小学奥数题常用骗分方法》

   枚举第一个同余方程的解 \(a_1+k\cdot p_1\)，直至找到一个解满足第二个同余方程，或枚举至 大于 \(\operatorname{lcm}(p_1,p_2)\)，枚举至此意味着该方程组 无解。

   就是优化点的暴力呗。时间复杂度 \(O(p_2)\)。

   之所以叫 大数 翻倍法，就是枚举 \(p\) 大的那个方程的解。这样到达无解的极端情况的时间更快。即 \(O(\min (p_1,p_2))\)。

   有人说“卡常小妙招”？ZHX：事实上，他还没被成功卡过。整体来说，他的复杂度是 \(\Theta((\sum_{i=1}^{k} p) -\max\limits_{i=1}^{k}\{p_i\})\)。之所以没被卡，是因为大部分时候都会出 \(p\) 是质数的情况，而如果要保证 \(\prod p\) 在 long long 范围内，那么其实最大的 \(p\) 对于这个积来说还是蛮小的，不是很好卡。

2. 扩展中国剩余定理（正解）

   正解好闪，拜谢正解。

   将方程写成：

   \[x=k_1\cdot p_1 + a_1=k_2\cdot p_2 + a_2 \]
   \[k_1\cdot p_1 - k_2\cdot p_2=a_2 - a_1 \]

   扩欧！

   没错！扩展欧几里得求 \(k_1,k_2\)。

   ZHX 调整后代码：

   #include<bits/stdc++.h>
   #define int long long
   const int N=1e5+5;
   int n;
   int p[N],a[N];
   int exgcd(int a,int b,int &x,int &y) {
   	if(!b) {x=1,y=0; return a;}
   	int g=exgcd(b,a%b,x,y),tmp=x;
   	x=y,y=tmp-(a/b)*y;
   	return g;
   }
   void combine(int x,int y) {
   	int k1,k2;
   	int g=exgcd(p[x],p[y],k1,k2);
   	int k=(a[y]-a[x])/g; k1*=k;
   	int z=k1*p[x]+a[x];
   	p[y]=p[x]/g*p[y]; a[y]=(z%p[y]+p[y])%p[y];
   }
   int solve() {
   	for(int i=2;i<=n;++i)
   		combine(i-1,i);
   	return a[n];
   }
   signed main() {
   	std::ios::sync_with_stdio(0);
   	std::cin.tie(0);std::cout.tie(0);
   	std::cin>>n;
   	for(int i=1;i<=n;++i)
   		std::cin>>p[i]>>a[i];
   	int ans=solve();
   	std::cout<<ans;
   	return 0;
   }
   

   这种做法被称为 扩展中国剩余定理（exCRT）。扩展中国剩余定理不要求 \(p\) 互质，nice。

   UPD. 2023.5.4 关于 ZHX 的代码，他 死 了。原因是过程量爆 long long 了。所以我们需要一个很新的乘法，就是类似于快速幂的方法，然后处处取模。

   #include<bits/stdc++.h>
   #define ll long long
   const int N=1e5+5;
   int n;
   ll p[N],a[N];
   template <typename T> void read(T &x)
   {
   	int f=1;x=0;char ch=getchar();
   	while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();}
   	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
   	x*=f;
   }
   ll mul (ll a,ll b,ll mod) {
   	ll ans=0;
   	while(b) {
   		if(b&1) ans=(ans+a)%mod;
   		a=(a+a)%mod; b>>=1;
   	}
   	return ans;
   }
   ll exgcd(ll a,ll b,ll &x,ll &y) {
   	if(!b) {x=1,y=0; return a;}
   	ll g=exgcd(b,a%b,x,y),tmp=x;
   	x=y,y=tmp-(a/b)*y;
   	return g;
   }
   ll excrt() {
   	ll x,y,g,k,tmp,ans=a[1],P=p[1];
   	for(int i=2;i<=n;++i) {
   		ll B=p[i],C=(a[i]-ans%B+B)%B;
   		g=exgcd(P,B,x,y),tmp=B/g;
   		if(C%g) return -1;
   		x=mul(x,C/g,tmp),ans+=x*P,P*=tmp,ans=(ans%P+P)%P;
   	}
   	return (ans%P+P)%P;
   }
   int main() {
   	std::ios::sync_with_stdio(0);
   	std::cin.tie(0);std::cout.tie(0);
   	read(n);
   	for(int i=1;i<=n;++i)
   		read(p[i]),read(a[i]);
   	std::cout<<excrt();
   	return 0;
   }
   

质数筛法

嗯，TSOI 暑假入门内容。
先讲的埃氏筛：

for(int i=2;i<=n;++i)
	if(!not_prime[i])
		for(int b=a+a;b<=n;b+=a)
			 not_prime[b]=1;
for(int i=2;i<=n;++i)
	if(!not_prime[i]) prime[++cnt]=i;

根据调和级数，这个玩意复杂度 \(O(n \log n)\)，加第二行优化是 \(O(\log (\log n))\)，Vergil 讲了但我不会证。

然后是线性筛。

for(int i=2;i<=n;++i) {
	if(!not_prime[i]) cnt++,prime[++cnt]=i;
	for(int j=1;j<=cnt;++j) {
		int x=prime[j]*i;
			if(x>n) break;
		not_prime[x]=1;
		if(!i%prime[j]) break;//保证每个数只被最小的质因数筛掉
	}
}

就是把埃氏筛的内层拉出来放外层，然后加上关键性的加注释的那一行优化。

积性函数

\(\forall\) 互质的 \(a,b\)，若有 \(f(a) \cdot f(b) = f(ab)\)，则 \(f(x)\) 是 积性函数。

\(\forall a,b\)，若有 \(f(a) \cdot f(b) = f(ab)\)，则 \(f(x)\) 是 完全积性函数。

易证欧拉函数是 积性函数，所以可以在线性筛的同时计算欧拉函数。

for(int i=2;i<=n;++i) {
	if(!not_prime[i]) cnt++,prime[++cnt]=i,phi[i]=i-1;
	for(int j=1;j<=cnt;++j) {
		int x=prime[j]*i;
		if(x>n) break;
		not_prime[x]=1;
		phi[x]=phi[prime[j]]*phi[i];//如果x与prime[j]互素
		if(!i%prime[j]) {
			phi[x]=prime[j]*phi[i];//如果x与prime[j]不互素
			break;
		}
	}
}

今天，ZHX说他 \({\color{RoyalBlue}\text{#临时}}\) 出了一道 \({\color{RoyalBlue}\text{#凑数题}}\)，是\({\color{RoyalBlue}\text{#T3}}\)，结果 \({\color{RoyalBlue}\text{#得分率最低}}\)，是 \({\color{RoyalBlue}\text{#斐波那契}}\)，靠。

5.1 - Day 3

上午

大步小步算法 - Baby Step Giant Step

AKA 北上广深 (BSGS)。

给定质数 \(a,b,p\)，求解方程 \(a^x \equiv b \pmod p\)，保证 \(a,b,p\leq 10^9\)。

考虑使用暴力，枚举至已经出现过的模数，即枚举一个循环内是否有解。由于 \(p\) 是素数，\(a^{p-1} \equiv 1 \pmod p\)，即一个循环的开始必然为 \(1\)，所以枚举至取模后为 \(1\) 即可。

int solve(int a,int b,int p) {
	int v=1;
	for(int x=0;x<p-1;++x) {
		if(x==b) return x;
		v=1ll*v*a%p;
	}
	return -1;
}

srds，他复杂度 \(O(p)\)，算了。

所以，介绍 Baby Step Giant Step 算法。

将一个循环（\(0 \sim p-2\)）分为一组 \(s\) 个，如果

\[\exists x \in \left[ (j-1) s, j s \right) \text{, s.t. } a^x \equiv b \pmod p , \]

那么令 $x=(j-1)s+k (k \in [0,s) \bigcap \mathbf{Z} ),\ $ 就有

\[\begin{aligned} a^{(j-1)s+k} & \equiv b \\ a^k & \equiv b \cdot a^{-(j-1)s} . \end{aligned} \]

即如果第 \(j\) 组有解可以满足该同余方程，则 \(b \cdot a^{-(j-1)s}\) 必然在第 \(1\) 组出现。最小非负整数解即为当前组的基数加上在第 \(1\) 组出现的位置，即 \((j-1)s+k.\)

int BSGS(int a,int b,ll p) {
	int s=sqrt(p),v=1;
	std::map<int,int> mp;
	for(int i=0;i<s;++i) mp[v]=i,v=(ll)v*a%p;
	for(int i=0;(ll)i*s<=p;++i) {
		int chk=(ll)b*ksm(ksm((ll)a,(ll)i*s,p),p-2,p)%p;
		if(mp.count(chk)) return mp[chk]+i*s;
	}
	return -1;
}

复杂度 \(O(\max(\dfrac{p}{s} ,s))\)，所以 \(s=\sqrt{p}\) 时均摊复杂度最低。

\[\color{Red}\textsf{初等数论结课！} \]

---

组合数学

好在主播有所了解。

排列数与组合数

排列数：在一个含 \(n\) 个元素的集合中选取 \(m\) 的元素，考虑顺序（同样的元素不同顺序视作不同集合），组成的子集的数量，记作 \(P(n,m)\) 或 \(A(n,m)\)。

\[P(n,m)=\dfrac{n!}{(n-m)!} \]

组合数：在一个含 \(n\) 个元素的集合中选取 \(m\) 的元素，不考虑顺序（同样的元素不同顺序视作相同集合），组成的子集的数量，记作 \(C(n,m)\)。

\[C(n,m)= \dfrac{n!}{m!(n-m)!} \]

下文中对排列数和组合数用下标表示 \(n,m\)，即 \(P_n^m\) 和 \(C_n^m\)。

一些有意思的结论：

\[C_n^0=C_n^n=1 \]

\[C_n^m=C_n^{n-m} \]

当然，还有就是组合数的 递推式：

\[C_n^m=C_{n-1}^{m-1}+C_{n-1}^m \]

应该都不难理解吧。

根据递推式，我们不难发现 杨辉三角 就是 组合数下三角阵，也易得可以通过 动态规划 计算组合数。也可以用 预处理阶乘 和 预处理组合数 计算排列数，完美解决 除法不满足取模盖斯定律 的问题。

还有一些性质：

\[C_n^0+C_n^1+C_n^2+\cdots+C_n^n=2^n \]

对于 \(n \geq 1\)，有

\[C_n^1 - C_n^2 + C_n^3 - \cdots \ C_n^n=0 \]

证明如下：

\[C_n^1=C_{n-1}^0+C_{n-1}^1 \]

\[C_n^2=C_{n-1}^1+C_{n-1}^2 \]

\[\cdots \]

\[C_n^n=C_{n-1}^{n-1} \]

一直枚举，刚好可以抵消。

二项式定理

将 \((x+y)^n\) 展开可以发现，其各项系数刚好就是杨辉三角的第 \(n+1\) 行。所以其实组合数也正好是 二项式系数。方便书写，我们也将 \(C_n^m\) 写作 \(\dbinom{n}{m}\)，代表二项式系数。

于是有：

\[(x+y)^n=\sum \limits_{i=0}^n \dbinom{n}{i}x^{n-i}y^i \]

同时，如果将 \(\dbinom{n}{m}\) 展开 \(k\) 次，也有：

\[\dbinom{n}{m}=\sum_{i=0}^k\dbinom{k}{i}\cdot \dbinom{n-k}{m-k+i} \]

令 \(j=k-i\)，得：

\[\begin{aligned} \sum_{i=0}^k\dbinom{k}{i}\cdot \dbinom{n-k}{m-k+i} & =\sum_{i=0}^k\dbinom{k}{k-j}\cdot \dbinom{n-k}{m-j} \\ & =\sum_{i=0}^k\dbinom{k}{j}\cdot \dbinom{n-k}{m-j} \end{aligned} \]

使用换元法，将 \(j\) 换成 \(i\)，得到最终结果：

\[\dbinom{n}{m}=\sum_{i=0}^k \dbinom{k}{i} \cdot \dbinom{n-k}{m-k+i} =\sum_{i=0}^k\dbinom{k}{i}\cdot \dbinom{n-k}{m-i} \]

\(6.\)

组合数的 6 种求法

\(ZHX\) 显出极高兴的样子，将两个指头的长指甲敲着讲台，点头说，“对呀对呀！……组合数有六样求法，你知道么？”

接下来我们一个一个说：

\(\textbf{Question. } \dbinom{n}{m} \bmod k\)

1. 特殊的模数

\(\textbf{Sol. } \alpha - n,m \leq 10^{18},k=1\)

\(k=1\) 啊兄弟！就是 \(0\)！

2. 递推式

\(\textbf{Sol. } \beta - n,m \leq 10^{3},k \leq 10^9\)

使用杨辉三角，\(O (n^2)\) 的复杂度。

for(int i=0;i<=n;++i) {
	C[i][0]=1;
	for(int j=1;j<=i;++j)
		C[i][j]=(C[i-1][j-1]+C[i-1][j])%p;
}

3. 预处理阶乘逆元

\(\textbf{Sol. } \gamma - n,m \leq 10^{6}, k \leq 10^9\text{ and } k \text{ is prime.}\)

预处理阶乘，然后除法使用逆元（\(k\) 是质数，费马小），定义式求解组合数。

fac[0]=1;
for(int i=1;i<=1000000;++i)
	fac[i]=(ll)fac[i-1]*i%p;
C(n,m) =(ll)fac[n]*ksm(fac[m],p-2,p)%p*ksm(fac[n-m],p-2,p)%p;

4. 定义式约分

\(\textbf{Sol. } \delta - n \leq 10^{9},m\leq 10^3,k \leq10^9\)

ZHX OI 社会学定律：没有思路看数据范围，一切尽在数据范围最诡异的数中。

本范围的 1e3 就够诡异。于是，我们尝试利用一下这个 1e3。

我们可以将组合数约分，消去 \((n-m)!\) ，就有

\[\begin{aligned}\dbinom{n}{m} & = \dfrac{n!}{m!(n-m)!}\\& = \dfrac{n \cdot (n-1) \cdot (n-2) \cdot ... \cdot (n-m+1)}{1 \cdot 2 \cdot ... \cdot m},\end{aligned} \]

这串式子的上下都只有 \(m\) 项，复杂度是可以压住的；考虑组合数必定是整数，我们就可以每次寻找分子与分母的 \(\gcd\) 进行约分，暴力但有效地解决这一问题。

for(int i=1;i<=m;i++)
	deno[i]=i,nume[i]=n-i+1;
for(int i=1;i<=m;++i)
	for (int j=1;j<=m;++j) {
		int g=gcd(nume[i],deno[j]);
		nume[i]/=g;
		deno[j]/=g;
	}
int ans=1;
for(int i=1;i<=m;++i)
ans=(ll)ans*nume[i]%p;

5. Lucas 定理

\(\textbf{Sol. } \epsilon - n,m \leq 10^9 ,k \leq 100 \text{ and k is prime.}\)

隆重介绍 卢卡斯定理：

对于质数 \(p\)，有：

\[\dbinom{n}{m} \bmod p=\dbinom{\lfloor n/p \rfloor}{\lfloor m/p \rfloor} \cdot \dbinom{n \bmod p}{m \bmod p} \bmod p \]

所以，我们可以将 \(n,m\) 转换为 \(p\) 进制，然后对于 \(n_{(p)}\) 的每位 \(a_i\) 和 \(m_{(p)}\) 的对应为 \(b_i\)，计算

\[\prod \dbinom{a_i}{b_i} \bmod p \]

即可。

int solve() {
	while(n) x[++x[0]]=n%p,n/=p;
	while(m) y[++y[0]]=m%p,m/=p;
	int ans=1;
	for(int i=1;i<=std::max(x[0],y[0]);++i) {
		if(i>x[0]) x[i]=0;
		if(i>y[0]) y[i]=0;
		if(x[i]<y[i]) { ans=0; break; }
		ans=(ll)ans*fac[x[i]]%p*ksm((ll)fac[y[i]]*fac[x[i]-y[i]]%p,p-2,p)%p;
	}
	return ans;
}

6. 扩展卢卡斯定理 - exLucas

\(\textbf{Sol. } \zeta - n,m \leq 10^9 ,k \leq 100\)

\(k\) 不是质数，怎么办呢？

把 \(k\) 拆分成质因数 \(k=p_1^{a_1} \cdot p_2^{a_2} \cdot \dots \cdot p_k^{a_k}\)，然后就会有同余方程组：

\[\begin{cases}\dbinom{n}{m} \bmod p_1^{a_1} = x_1\\ \dbinom{n}{m} \bmod p_2^{a_2} = x_2\\ \cdots\\ \dbinom{n}{m} \bmod p_k^{a_k} = x_k\\ \end{cases} \]

解 中国剩余定理 即可。

具体过程详见 \(dalao\) 题解。

\(\textbf{AC Code:}\)

#include<bits/stdc++.h>
#define ll long long
#define ld long double
#define fsp(x) std::fixed<<std::setprecision(x)
#define forE(u) for(ll p=head[u],v=E[p].to;p;p=E[p].next,v=E[p].to)
const int N=1e5+5;
ll cnt,P[N],A[N];
template <typename _Tp> _Tp ksm(_Tp base,ll ind,ll mod) {
	_Tp ans=1;
	while(ind) {
		if(ind&1) ans=ans*base%mod;
		base=base*base%mod; ind>>=1;
	}
	return ans%mod;
}
ll exgcd(ll a,ll b,ll &x,ll &y) {
	if(!b) return x=1,y=0,a;
	ll g=exgcd(b,a%b,x,y),tmp=x;
	x=y,y=tmp-(a/b)*y;
	return g;
}
inline ll inv(ll a,ll p) { ll x,y; exgcd(a,p,x,y); return (x+p)%p; }
ll f_pk(ll n,ll p,ll pk) {
	if(!n) return 1;
	ll ans=1;
	for(ll i=1;i<=pk;++i)
		if(i%p) ans=ans*i%pk;
	ans=ksm(ans,n/pk,pk);
	for(ll i=pk*(n/pk);i<=n;++i)
		if(i%p) ans=ans*(i%pk)%pk;
	return f_pk(n/p,p,pk)*ans%pk;
}
ll g_p(ll n,ll p) { return n<p?0:n/p+g_p(n/p,p); }
inline ll C_pk(ll n,ll m,ll p,ll pk) { return f_pk(n,p,pk)*inv(f_pk(m,p,pk)*f_pk(n-m,p,pk),pk)%pk*ksm(p,g_p(n,p)-g_p(m,p)-g_p(n-m,p),pk)%pk; }
inline ll excrt() {
	ll x,y,g,tmp,ans=A[1],p=P[1];
	for(ll i=2;i<=cnt;++i) {
		ll B=P[i],C=(A[i]-ans%B+B)%B;
		g=exgcd(p,B,x,y),tmp=B/g;
		x=C/g*x%tmp,ans+=x*p,p*=tmp,ans=(ans%p+p)%p;
	}
	return (ans%p+p)%p;
}
inline ll exLucas(ll n,ll m,ll p) {
	for(ll i=2;i*i<=p;++i) {
		ll tmp=1;
		while(!(p%i)) tmp*=i,p/=i;
		if(tmp==1) continue;
		P[++cnt]=tmp,A[cnt]=C_pk(n,m,i,tmp);
	}
	if(p!=1) P[++cnt]=p,A[cnt]=C_pk(n,m,p,p);
	return excrt();
}
int main() {
	std::ios::sync_with_stdio(false);
	std::cin.tie(nullptr); std::cout.tie(nullptr);
	ll n,m,p;
	std::cin>>n>>m>>p;
	ll ans=exLucas(n,m,p);
	std::cout<<ans<<'\n';
	return 0;
}

tips：关于非常大的 \(\dbinom{n_1}{m_1},\dbinom{n_2}{m_2}\) 比大小，可以转换成 \(\log \dbinom{n_1}{m_1},\log \dbinom{n_2}{m_2}\) 进行比较，开 double 就行。

求和变形

遇到多层求和 \(\sum \sum\) 的时候，复杂度压不住啊。

于是有了求和变形。

1. 增加枚举量

2. 交换枚举顺序

3. 分离无关变量

4. 换元法

比如 P3746 [六省联考 2017] 组合数问题，他想求：

\[\left( \sum_{i = 0}^\infty C_{nk}^{ik + r} \right) \bmod p, \]

我就可以

\[\begin{aligned} & \quad \ \sum_{i = 0}^\infty C_{nk}^{ik + r} \\ & = \sum_{i=0}^\infty \dbinom{nk}{ik+r} \\ &= \sum_{i=0}^{\infty} \sum_{j=0}^k \dbinom{k}{j}\cdot \dbinom{nk-k}{ik+r-j} & \mathit{1}\\ & = \sum_{j=0}^k \sum_{i=0}^{\infty} \dbinom{k}{j}\cdot \dbinom{nk-k}{ik+r-j} & \mathit{2}\\ & =\sum_{j=0}^k \dbinom{k}{j} \sum_{i=0}^{\infty} \dbinom{(n-1)k}{ik+r-j} & \mathit{3} & ,\\ \end{aligned} \]

其中 \(\mathit{1}\) 式使用了增加枚举量，\(\mathit{2}\) 式使用了交换枚举顺序，\(\mathit{3}\) 式使用了分离无关变量。那么，到了这里如何解呢？

我们使用第 4 招——换元法，设

\[f(n,r)=\sum \limits_{i=0}^{\infty} \dbinom{nk}{ik+r}, \]

就有

\[f(n-1,r-j)=\sum \limits_{i=0}^{\infty} \dbinom{(n-1)k}{ik+r-j}, \]

接下来，添加辅助纬度，将原式表达为

\[f_n(1,r)=\sum \limits_{i=0}^{\infty} \dbinom{nk}{ik+r}, \]

再使用多项式定理，就有

\[\begin{aligned} f_n(1,r) & = \sum \limits_{i=0}^{\infty} \dbinom{nk}{ik+r}\\ & = \sum \limits_{i=0}^{\infty} \sum \limits_{j=0}^k \dbinom{k}{j} \dbinom{(n-1)k}{ik+r-j}\\ & = \sum \limits_{j=0}^k \sum \limits_{i=0}^{\infty}\dbinom{(n-1)k}{ik+r-j} \dbinom{k}{j}\\ & = \sum \limits_{j=0}^k f_{n-1}(1,r-j)\dbinom{k}{j}, \end{aligned} \]

怎么样，是不是有点意思？

到了这儿，其实复杂度并没有降，所以是时候想个方法了——矩阵快速幂。但是现在还不满足条件，所以我们开一个新的矩阵 \(D\)，并使

\[D(r-j,r)=\dbinom{k}{j}, \]

于是就有了

\[f_n(1,r)=\sum_{j=0}^{k}f_{n-1}(1,r-j)\cdot D(r-j,r), \]

完全满足矩阵乘法。

所以，我们就可以使用

\[f_n=f_{n-1} \times D \]

结合矩阵快速幂完美解决这道题，也就是说我们使用矩阵快速幂，然后

\[\begin{aligned} \quad \ \sum_{i = 0}^\infty C_{nk}^{ik + r} & =\sum_{j=0}^k \dbinom{k}{j}\sum_{i=0}^{\infty} \dbinom{(n-1)k}{ik+r-j}\\ & =f_n(1,r). \end{aligned} \]

泰裤辣！

抽屉原理

有 \(n+1\) 个物品和 \(n\) 个抽屉，把这些物品放进抽屉，至少有 一个抽屉 里面有 多于一个 物品。

没了。

还在？

真没了。

就把其他的题抽象成这个模型就行了，没了。

容斥原理

最基本的，

\[\left| A_1 \cup A_2 \right|=| A_1 | + |A_2| - |A_1 \cap A_2|, \]

\[\begin{aligned}|A_1\cup A_2 \cup A_3| & =|A_1|+|A_2|+|A_3|\\ & -|A_1 \cap A_2|-|A_1 \cap A_3|-|A_2 \cap A_3|\\ & + |A_1 \cap A_2 \cap A_3|, \end{aligned} \]

反正就是，一层加，一层减，然后就有一个你看不懂但是可以意会其精神的式子：

\[\left| \bigcup_{i=1}^n S_i \right|= \sum_{m=1}^n (-1)^{m-1} \sum_{a_i<a_{i+1}} \left| \bigcap_{i=1}^m S_{a_i} \right|, \]

理论就结束了。

全靠思考，抽象模型。

在容斥原理的运用中，我们可以将较强的条件（苛刻的）转换为较弱的条件（如不满足其中一个条件），然后使用容斥原理求解。

5.2 - Day 4

上午

容斥定理解决春季赛 T2

没错！其实春季赛 T2 也是一道 容斥原理 的题。

[春季测试 2023] 幂次

题目描述

小 Ω 在小学数学课上学到了“幂次”的概念：\(\forall a, b \in \mathbb{N}^+\)，定义 \(a^b\) 为 \(b\) 个 \(a\) 相乘。

她很好奇有多少正整数可以被表示为上述 \(a^b\) 的形式？由于所有正整数 \(m \in N^+\) 总是可以被表示为 \(m^1\) 的形式，因此她要求上述的表示中，必须有 \(b \geq k\)，其中 \(k\) 是她事先选取好的一个正整数。

因此她想知道在 \(1\) 到 \(n\) 中，有多少正整数 \(x\) 可以被表示为 \(x = a^b\) 的形式，其中 \(a, b\) 都是正整数，且 \(b \geq k\)？

输入格式

第一行包含两个正整数 \(n, k\)，意义如上所述。

输出格式

输出一行包含一个非负整数表示对应的答案。

考虑使用枚举法，依次枚举 \(b=k,k+1,k+2 \cdots\) 时所有的数，可以推出在 \(n\) 以内这样的数有 \(\lfloor \sqrt[b] n \rfloor\) 个。当然，其中会有重复计数。比如，\(a^6\) 就会被 \(a^2\) 和 \(a^3\) 同时计数。所以，我们需要扣去一次。

那么，如何确定当前幂次是应该加上还是减去呢？理论上，这需要用到 莫比乌斯反演。但是显然这对目前阶段的我们来说可能还是有点超纲了。所以我们采用一种更加暴力的方法。

使用 \(num^b\) 表示 \(a^b\) 被计数了几次，根据 \(num\) 决定当前幂次是加是减。核心部分如下：

for(int a=2;a<=64;a++)
	num[a] = 0;
for(int a=2;a<=64;a++) {
//	long long v=pow(n,1.0/a)-1; 可能导致精度问题
	long long v=(long long)floor(exp(log(n)/a))-1; //精度优化
	int d=1-num[a];
	ans+=v*d;
	for(int b=a;b<=64;b+=a)
		num[b]+=d;
}

冷知识:

\[2^{\frac{\log_2 n}{a}}=\sqrt[a]{n} \]

可以解决精度问题。

---

线性代数

高斯消元

现在我们要解决

\[\begin{cases} a^1_1x_1+a_1^2x_2+\cdots+a_1^nx_n=b_1\\ a^1_2x_1+a_2^2x_2+\cdots+a_2^nx_n=b_2\\ \cdots\\ a^1_nx_1+a_n^2x_2+\cdots+a_n^nx_n=b_n \end{cases} \]

这么一个 \(n\) 元一次方程组，how？

高斯消元。

其实就是把咱们手解方程的过程转换成代码了。

void gauss() {
	for(int i=1;i<=n;++i) {
		for(int j=i;j<=n;++j)
			if(a[j][i]!=0) {
				for(int k=1;k<=n+1;k++) swap(a[i][k],a[j][k]);
				break;
			}
		for(int j=i+1;j<=n;++j) {
			if(a[j][i]==0) continue;
			int l=a[i][i]/gcd(abs(a[i][i]),abs(a[j][i]))*a[j][i];
			int ratioi=l/a[i][i];
			int ratioj=l/a[j][i];
			for (int k=1;k<=n+1;++k)
				a[j][k] = a[j][k] * ratioj - a[i][k] * ratioi;
		}	
	} 
	for(int i=n;i>=1;--i) {
		for(int j=i+1;j<=n;++j)
			a[i][n+1]-=a[i][j]*x[j];
		x[i]=a[i][n+1]/a[i][i];
	}
}

还有一个优化一点：

bool gauss() {
	for(int i=1;i<=n;++i) {
		int r=i;
		for(int j=i;j<=n;++j)
			if(fabs(a[j][i])>fabs(a[r][i]))
				r=j;
		if(fabs(a[r][i])<eps) return 0;
		if(i!=r) std::swap(a[i],a[r]);
		for(int j=i+1;j<=n;++j) {
			double ratio=a[j][i]/a[i][i];
			for(int k=1;k<=n+1;++k)
				a[j][k]-=a[i][k]*ratio;
		}
	}
	for(int i=n;i>=1;--i) {
		for(int j=i+1;j<=n;++j)
			a[i][n+1]-=a[i][j]*x[j];
		x[i]=a[i][n+1]/a[i][i];
	}
	return 1;
}

最后是升级版的 高斯 - 约旦消元法：

bool jordan() {
	for(int i=1;i<=n;++i) {
		int r=i;
		for(int j=i;j<=n;++j)
			if(fabs(a[j][i])>fabs(a[r][i])) r=j;
		if(fabs(a[r][i])<eps) return 0;
		if(i!=r) std::swap(a[i],a[r]);
		for(int j=1;j<=n;++j) {
			if(i!=j) {
				double ratio=a[j][i]/a[i][i];
				for(int k=i+1;k<=n+1;++k)
					a[j][k]-=a[i][k]*ratio;
			}
		}
	}
	for(int i=1;i<=n;++i) x[i]=a[i][n+1]/a[i][i];
	return 1;
}

年轻人！珍惜约旦优美的"\(\setminus\)"型吧。

单位矩阵

就是这样一个矩阵：

\[\begin{bmatrix} 1 & 0 & \cdots \\ 0 & 1 & \ddots \\ \vdots & \ddots & \ddots \end{bmatrix} \]

有这么一个性质：

\[\begin{bmatrix} x & 0 & \cdots \\ 0 & x & \ddots \\ \vdots & \ddots & \ddots \end{bmatrix} \times A = xA \]

然后就可以把数乘转化成矩阵乘法了，可以运用这么一个构造函数对矩阵赋 int 型的初值。

matrix(int z=0,int a=0,int b=0) {
	n=a,m=b,memset(x,0,sizeof x);
	for(int i=0;i<=24;++i)
	x[i][i]=z;
}

那个 24 是我开的矩阵的数组的长度，不同情况不同对待。

单位矩阵只能是 \(n \times n\) 的，即 行数与列数相同。

初等矩阵

一个奇妙的矩阵：

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

你会发现他可以这样：

\[\begin{bmatrix} 0 & 1\\ 1 & 0 \end{bmatrix} \times \begin{bmatrix} 1 & 2\\ 3 &4 \end{bmatrix} = \begin{bmatrix} 3 & 4\\ 1 & 2 \end{bmatrix} \]

\[\begin{bmatrix} 1 & 2\\ 3 & 4 \end{bmatrix} \times \begin{bmatrix} 0 & 1\\ 1 & 0 \end{bmatrix} = \begin{bmatrix} 2 & 1\\ 4 & 3 \end{bmatrix} \]

这个矩阵是将单位矩阵的第 \(1\) 行和第 \(2\) 行交换，比如还有这个矩阵：

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

他相比单位矩阵，第 \(2\) 行和第 \(3\) 行做了交换。他乘矩阵 \(X\) 后，结果就是 \(X\) 第 \(2\) 行和第 \(3\) 行 交换。

而如果单位矩阵的第 \(i\) 行第 \(j\) 列增加 \(x\)，将该矩阵乘矩阵 \(X\)，就会将 \(X\) 的第 \(j\) 行的 \(x\) 倍加到第 \(i\) 行。

这些矩阵都是 初等矩阵，他们是由单位矩阵经过 初等变换 得来的。

初等变换共有三种：

1. 交换两行或两列；

2. 用一个数 \(k\) 乘某一行；

3. 用一个数 \(k\) 乘某一行后加到另一行去（乘行不变）。

哎，不行就看看 这篇 博吧，反正对于 OIer 来说这些应该已经够了。

逆矩阵

对一个 \(n\) 阶矩阵 \(A\) ，如果存在另一个 \(n\) 阶矩阵 \(B\)，它们满足：\(AB = BA = E\)（\(E\) 为单位矩阵），那么 \(A\) 和 \(B\) 互为逆矩阵。

如果你把

\[\begin{cases} a^1_1x_1+a_1^2x_2+...+a_1^nx_n=b_1\\ a^1_2x_1+a_2^2x_2+...+a_2^nx_n=b_2\\ ...\\ a^1_nx_1+a_n^2x_2+...+a_n^nx_n=b_n \end{cases} \]

写成

\[\begin{bmatrix}a^1_1 & a_1^2 & \cdots & a_1^n\\ a^1_2 & a_2^2 & \cdots & a_2^n\\ \vdots & \vdots & \vdots & \vdots\\ a^1_n & a_n^2 & \cdots & a_n^n \end{bmatrix} \times \begin{bmatrix} x_1\\ x_2\\ \vdots \\ x_n \end{bmatrix} = \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n \end{bmatrix}, \]

然后就不难发现，解方程组其实可以看作求逆矩阵的过程。

同时，也有一些题会问到你，让你去求逆矩阵，可以把原矩阵 \(A\) 和单位矩阵 \(E\) 放在一起进行高斯消元，就会有

\[A \times B = E \]

\[E \times B = Y \]

然后这个矩阵 \(Y\) 就是你要求的逆矩阵了。

---

初等概率

下午

嗨害，直接就是概率啊。

\[P(A)\cdot P(B)=P(AB) \]

• 随机试验：

  1. 不能预先确知结果；

  2. 试验之前可以预测所有可能结果或范围；

  3. 可以在相同条件下重复实验。

• 样本空间：随机试验所有可能结果组成的集合。

• 事件：样本空间的任意一个子集。

• 事件发生：在一次试验中，事件的一个样本点发生。

• 必然事件：样本空间全集。

• 不可能事件：空集。

• 条件概率：\(B\) 发生时，\(A\) 发生的概率。记作 \(P(A|B).\)

  \[P(A|B)=\dfrac{P(AB)}{P(B)} \]
  \[P(A|B)P(B)=P(AB) \]

• 独立事件：\(A\) 是否发生与 \(B\) 无关。

• 期望：我懂了，但我不好说。大概是是每种结果贡献与概率的乘积的和。

  期望有一个非常重要的性质，就是

  \[E[x_1+x_2]=E[x_1]+E[x_2], \]

  即 期望的和等于和的期望。

最后的最后，我们来一道题，实战一下：

小胡站在原点，手里拿着两枚硬币。抛第一枚硬币正面向上的概率为 \(p\)，第二枚正面向上的概率为 \(q\)。
小胡开始抛第一枚硬币，每次抛到反面小胡就向 \(x\) 轴正方向走一步，直到抛到正面。
接下来小胡继续抛第一枚硬币，每次抛到反面小胡就向 \(y\) 轴正方向走一步，直到抛到正面。
现在小胡想回来了，于是他开始抛第二枚硬币，如果小胡抛到正面小胡就向 \(x\) 轴的负方向走一步，否则小胡就向 \(y\) 轴的负方向走一步。
现在小胡想知道他在往回走的时候经过原点的概率是多少呢？

OK，我们可以把计算概率的式子列出来：

\[\sum \limits_{x=0}^{\infty} \sum \limits_{y=0}^{\infty} (1-p)^x p \cdot(1-p)^y p \cdot q^x(1-q)^y \dbinom{x+y}{x} \]

因为有可能走到无限远，所以求和从 \(0\) 到 \(\infty.\) \((1-p)^x p\) 是走到 \((x,0)\) 的概率，\((1-p)^y p\) 是在此基础上走到 \((x,y)\) 的概率，\(q^x(1-q)^y \dbinom{x+y}{x}\) 是再走回 \((0,0)\) 的概率。接下来，进行史诗级变形：

\[\begin{aligned} & \ \sum \limits_{x=0}^{\infty} \sum \limits_{y=0}^{\infty} (1-p)^x p \cdot(1-p)^y p \cdot q^x(1-q)^y \dbinom{x+y}{x}\\ = & \ \sum \limits_{x=0}^{\infty} \sum \limits_{y=0}^{\infty} p^2 (1-p)^{x+y} \cdot q^x(1-q)^y \dbinom{x+y}{x}\\ = & \ \sum \limits_{t=0}^{\infty} \sum \limits_{x=0}^{t} p^2 (1-p)^{t} \cdot q^x(1-q)^{t-x} \dbinom{t}{x}\\ = & \ p^2 \sum \limits_{t=0}^{\infty} (1-p)^t\sum \limits_{x=0}^{t} \cdot q^x(1-q)^{t-x} \dbinom{t}{x}\\ = & \ p^2 \sum \limits_{t=0}^{\infty} (1-p)^t (q+1-q)^t\\ = & \ p^2 \sum \limits_{t=0}^{\infty} (1-p)^t\\ = & \ p^2 \cdot \dfrac{1-(1-p)^{\infty+1}}{1-(1-p)}\\ = & \ p^2 \cdot \dfrac{1}{p}\\ = & \ {\color{Red}p} . \end{aligned} \]

我知道你现在的心情，一方面是像疯了一样，震惊于这个变形绝妙的美感；一方面是像傻了一样，尝试用毕生所学去理解这个式子。接下来，我们一步一步来说。

第一步，合并同类项，即 \(p^2\) 和 \((1-p)^{x+y}\)；

第二步，换元法，设 \(t=x+y\)，更换枚举变量；

第三步，分离无关变量，改变了 \(p^2\) 和 \((1-p)^{x+y}\) 的位置；

第四步，运用 二项式定理，换 \(\sum \limits_{x=0}^{t} \cdot q^x(1-q)^{t-x} \dbinom{t}{x}\) 为 \((q+1-q)^t\)；

\[\textbf{二项式定理：} (a+b)^n = \sum \limits_{x=0}^n a^x \cdot b^{n-x} \cdot \dbinom{n}{x}. \]

第五步，化简该式；

第六步，运用 等差数列求和公式，化 \(\sum \limits_{t=0}^{\infty} (1-p)^t\) 为 \(\dfrac{1-(1-p)^{\infty+1}}{1-(1-p)}\);

\[\begin{aligned} \textbf{等差数列求和公} & \textbf{式推导：}\\ x&=1+a+a^2+\cdots +a^n\\ ax&=a^1+a^2+a^3+\cdots+a^{n+1}\\ (1-a)x&=1-a^{n+1}\\ x&=\dfrac{1-a^{n+1}}{1-a}. \end{aligned} \]

第七步，运用极限思想，认为 \((1-p)^{\infty+1}\) 无限趋近于 \(0\)，将 \(\dfrac{1-(1-p)^{\infty+1}}{1-(1-p)}\) 化简为 \(\dfrac{1}{p}\);

最后一步，化简得出结果。

没想到吧！如此繁杂的式子，变形化简之后竟只剩一个 \(p\)！

---

结语

五一数学专题集训就到此为止了，但是探索的旅途永远不会结束。不论是在数学上，还是在信竞中，人类正以昂扬的姿态谱写着新的历史——而 \(TSOI\) 的新的历史，将由我们书写。

王重阳，2023.5.4 署

\[\Large\color{Red}\textbf{- The End. -} \]